Spatio-Temporal Dynamics of Benthic Macrofaunal Communities in Relation to the Recovery of Coastal Aquaculture Operations Following the 2011 Great East Japan Earthquake and Tsunami

The 2011 Great East Japan Earthquake and tsunami wiped out over 1200 shellfish and ascidian culture long-lines and ~120 salmon farm cages that comprised the entire aquaculture installations in Onagawa Bay, Japan, and severely altered the associated ecosystem. A year later, we launched a coordinated monitoring programme to measure the extent of the damage caused by the disaster and monitor the change in the state of the marine ecosystem. As part of this effort, we conducted multi-seasonal sampling to characterise spatio-temporal variation in benthic macrofaunal community and a range of environmental parameters across Onagawa Bay between March 2012 and January 2018. The 492 total macrofaunal species recorded included Polychaeta (38.8 %), Bivalvia (13.2 %), Amphipoda (10.8 %), Decapoda (9.6 %), Gastropoda (9.3 %) and Echinodermata (4.3 %). At the outermost reference site, macrofaunal abundance, biomass, and species diversity were all consistently high throughout the study period. Inside Onagawa Bay, macrofauna metrics increased steadily from the lowest values at the beginning of the study to the highest over time. During the same period, the spatial extent of aquaculture facilities for long-lines and fish cages recovered steadily to within 60.8 % and 74.8 % of the original state, respectively. The significant variables identified by multivariate analysis to explain spatio-temporal variability in benthic macrofaunal communities were: (1) proximity to the nearest aquaculture facilities; (2) wind fetch length (exposure); (3) sediment grain size; (4) the total area of aquaculture facilities. This study suggests that coastal aquaculture operations may strongly influence the occurrence and distribution of benthic macrofaunal communities and thereby influence the recovery of seafloor biota at ecosystem scales following a catastrophic natural disaster.

On 11 March 2011, the Great East Japan Earthquake and subsequent tsunami devastated extensive areas of the Pacific coast of the Tohoku region in northern Japan. The major disturbances caused by the disaster included the immediate physical force of the tsunami and associated negative impacts on intertidal biota (Urabe et al., 2013), large amounts of terrestrial debris washed out into the ocean, a severe accident at the Fukushima Daiichi nuclear power plant that caused leakage of radioactive contaminants (Sohtome et al., 2014), spills of toxic chemicals (e.g., heavy oil), and associated negative impacts on subtidal benthic fauna (Abe et al., 2015). Coastal regions also suffered from loss of seagrass/seaweed beds with their associated fauna (Takami et al., 2013), abrupt accumulation of tsunami deposits and resuspended mud on the seafloor (Seike et al., 2013;Toyofuku et al., 2014;Kanaya et al., 2015;Seike et al., 2016), and widespread land subsidence. In order to measure the extent of the damage caused by the disaster and to monitor the subsequent change in the state of the marine environment, we launched the "Tohoku Ecosystem-Associated Marine Sciences" (TEAMS) project in March 2012. As part of this study, we initiated a coordinated monitoring program to obtain a holistic view on the spatio-temporal dynamics of the marine ecosystem of Onagawa Bay (Figure 1). To this end, we monitored a range of environmental parameters including biogeochemistry, physical measurements, sediment properties, and community dynamics of benthic and pelagic components potentially impacted by the disaster. This study presents the benthic component and associated environmental factors.
Prior to the earthquake and tsunami, fishing and aquaculture were the primary industries for many towns and villages along the Tohoku Pacific coast. For example, Onagawa Bay supported over 1200 long-lines for culturing scallops, oysters, and edible ascidians as well as around 120 fish cages for farming salmonids. However, the 2011 earthquake and tsunami wiped out entire aquaculture installations. Rehabilitation of coastal aquaculture therefore attracted considerable interest from both marine scientists and policy makers. Recent studies demonstrate varying, but often negative, effects of aquaculture on the environment through changes in the physical, chemical, and biological attributes of sediments below and the water column around aquaculture installations (e.g., Tovar et al., 2000;Borja et al., 2009;Forchino et al., 2011;Sarà et al., 2011;Farmaki et al., 2014;Tomassetti et al., 2016). Most studies demonstrate reduced impacts of shellfish aquaculture than finfish aquaculture because shellfish are grown at comparatively lower intensity and require no external feed (Shumway et al., 2003;McKindsey et al., 2006;Weise et al., 2009). However, shellfish farms typically cover a much greater area than finfish farms and may affect the surrounding environment differently (Weise et al., 2009).
Benthic macrofauna have long been used as indicators of environmental health and trends in coastal marine ecosystems (Borja et al., 2009;Forchino et al., 2011;Tomassetti et al., 2016). We therefore investigated spatio-temporal dynamics of benthic macrofaunal community structure in relation to changes in multiple environmental factors including the rehabilitation of coastal aquaculture facilities in Onagawa Bay following the 2011 Tohoku earthquake and tsunami. Our study addresses four lines of investigation: (1) characterize spatio-temporal dynamics of benthic macrofaunal community structure; (2) describe spatio-temporal variation in major environmental variables; (3) map spatio-temporal distribution of coastal aquaculture facilities to estimate changes in the location, number, and total area of aquaculture facilities; (4) examine the relationship between macrofaunal community structure and the explanatory environmental variables described above.
In addition, we extracted a subset of an external survey data collected near our sampling locations before the 2011 disaster as a baseline for comparison. Our overall aim is to provide insight into the environmental and anthropogenic factors most strongly related to the recovery of seafloor biota and coastal environmental health in this region.

Study Site
Onagawa Bay has an area of 27.15 km 2 , located at the southernmost part of a ria coastline, known locally as the "Sanriku Coast, " on the Pacific Ocean of northern Japan (Figure 1). Geologically, this stretch of coast was formed by a combination of sea level rise and the submergence of mountains and river valleys that resulted in an extremely irregular and indented coastline. Onagawa Bay is also a natural sheltered deepwater harbor with a mean depth of 33 m, which makes it suitable for coastal aquaculture.

Onagawa Bay Monthly Survey on RV Suikoh
Most of our data were collected during the Onagawa Bay Monthly Survey on board RV Suikoh between March 2012 and January 2018. A total of 12 sampling stations (i.e., St.1-8, 11, 12, 16, and 17) were established to monitor multiple environmental parameters (Figure 1 and Table 1). Benthic and physical measurements were made at all 12 stations only in March 2012 (Supplementary Figure 1). Thereafter, benthic sampling was conducted mainly at St.1,3,4,6,8,11,12,16, and 17 (nine stations) on a monthly basis between May 2012 and March 2014, which we reduced to every 2 months from March 2014 onward.
We sampled a total of four sediment samples at each station using a Smith-McIntyre grab (sampling area: 0.05 m 2 ). For benthic macrofauna, we sieved and merged three of the grab samples on a 1.0 mm mesh with filtered sea water on board. All organisms retained were preserved in 10% buffered formalin for  subsequent sorting, species identification, counting, and biomass measurements. Species were identified to the lowest possible taxonomic level using a stereomicroscope and a compound microscope where necessary, expressing abundance and biomass (wet weight) as ind./m 2 and g/m 2 , respectively. We then extracted subsamples of the top 5 and 1 cm of the sediment to determine particle size composition and chemical oxygen demand (COD), respectively. Particle size composition was determined using both dry sieving and sedimentation analysis and categorized into three size ranges defined here as "gravel" (GRAV: >2 mm), "sand" (SAND: <2 mm, >75 µm), and "mud" (MUD: <75 µm) ( Table 1). COD (mg/g) was determined using a standardized solution of potassium permanganate as the oxidant (Table 1).
At each sampling station, we recorded profiles of water temperature (TEMP: • C), salinity (SALI), chlorophyll-a concentration (CHLa: µg/l), dissolved oxygen (DO: mg/l), and turbidity (TURB: FTU) using a CTD RINKO-Profiler (JFE Advantech Co., Ltd., Japan) including measurements immediately above the seabed (Table 1). Distance to the nearest land (km) (D.LAND) and wind fetch distance (km) (FETCH) were calculated as indicators of the exposure of coastal locations ( Table 2). FETCH was determined as an average length (distance to the nearest land) of fetch vectors in 36 equiangular directions with a maximum distance set to 200 km if unobstructed by land, following Burrows et al. (2008) and Seers (2018).

Aquaculture Facilities in Onagawa Bay
Aquaculture operations in Onagawa Bay commonly use two main methods. Hanging culture uses a long-line system to suspend cultured organisms such as oysters, scallops, and ascidians on vertical ropes, typically down to around 15-20 m in depth below surface floats, in contrast to floating fish cages used for farming salmonids. We digitized and mapped these floating aquaculture structures using time-series of high resolution satellite imagery to create distribution maps for the long-lines (polylines) and fish cages (polygons) for the whole of Onagawa Bay from 2011 March to 2017 March (Table 3). Each polyline represents a single unit of the long-line system which, in fact, comprises two parallel lines of the same length (typically ∼100 m), approximately 1.5 m apart and connected at regular intervals by floats and/or ropes. We estimated the sea surface area occupied by each long-line by creating a buffer of 0.75 m around each polyline and then computing the area of the buffer zone (polygon).

Data Analysis
We identified structural and spatio-temporal trends in the macrofauna data in Onagawa Bay using multivariate community analysis in PRIMER v6 (Clarke and Warwick, 2001). Following fourth-root transformation of macrofauna counts to downweight the influence of highly abundant species, we performed cluster analysis (group-average linkage) on a resemblance matrix of the transformed data based on the Bray-Curtis similarity index. A similarity profile (SIMPROF) permutation test on the resemblance matrix identified statistically significant clusters (significance level, p < 0.01). SIMPER (similarity percentage) analysis then identified those species most responsible for the similarity within the groups determined by the cluster analysis.
We assessed general temporal trends by constructing time series for the 16 environmental variables (Table 1) and normalized them to generate a resemblance matrix based on Euclidean distance for conducting RELATE routines and BIO-ENV stepwise (BEST) analysis. RELATE routines examined potential differences in macrofaunal community structure between the clusters in relation to each of the 16 environmental variables. The BEST procedure examined rank correlations between the multivariate patterns of macrofauna and environmental factors to identify the subsets of variables that best explained the overall pattern. We then examined the biological metrics of each cluster group (i.e., abundance, biomass, species richness, and Shannon's diversity index H ) using box plots. Further, based on the BEST procedures, we assessed those environmental variables that significantly explained the macrofaunal variation using the same approach.
In addition, we extracted a subset of external data collected before the 2011 disaster from annual measurements of warm drainage water from the Onagawa nuclear power plant (Miyagi Prefecture and Tohoku Electric Power Company, 2010, 2012. This program has monitored a wide range of physical and biological data in Onagawa Bay since 1981. The program conducted benthic surveys twice a year in February and August at 18 sampling stations, of which 5 sites, namely TE1, 4, 6, 8, and 18, were close to sampling stations of St.1, 17, 11, 6, and 8, respectively, in our study (Figure 1), reporting the following information: (1) total abundance of macrofauna (ind./m 2 ); (2) species richness of macrofauna; (3) SAND (%); (4) MUD (%); (5) GRAV (%); and (6) COD (mg/g) (coding of environmental variables as in Table 1). The survey reports include details on sampling methods and protocols (Miyagi Prefecture and Tohoku Electric Power Company, 2010, 2012. We extracted data for August and February in each year from 2008 to 2011 (pretsunami: n = 6 per variable at TE1, 4, 6, 8, and 18), for comparison of mean values with our corresponding data for July and January from 2015 to 2018 in the present study (post-tsunami: n = 6 per variable at St.1, 17, 11, 6, and 8).
We constructed distribution maps for the aquaculture facilities using ESRI ArcGIS 10.5 and calculated values for FETCH using "fetchR" package (Seers, 2018). All the graphics relating to the multivariate analysis were constructed using R version 3.4.1 package (R Development Core Team, 2017).

Benthic Macrofaunal Community Structure
We collected a total of 382 benthos samples and identified 492 species comprising 191 from the class Polychaeta (38.8%), 65 from the class Bivalvia (13.2%), 53 from the order Amphipoda (10.8%), 47 from the order Decapoda (9.6%), 46 from the class Gastropoda (9.3%), 21 from the phylum Echinodermata (4.3%), and 69 others (14.0%). Macrofaunal abundance, species richness, and Shannon index H were lowest soon after the 2011 disaster, but increased rapidly over subsequent years to peak values toward the end of the study period (Figures 2A,C,D). Macrofaunal biomass decreased slightly shortly after sampling began, but then increased toward the end of the study period ( Figure 2B).
Cluster analysis and a SIMPROF test revealed 11 distinctive macrofaunal groups (Groups A-K) and outliers (OUT) that further subdivided into 40 statistically significant groups (Figure 3). While Groups A, D, B, and G comprised only 5 (1.3%), 5 (1.3%), 7 (1.8%), and 8 (2.1%) of the samples, respectively, the majority of samples were categorized into K (100 samples: 26.2%), I (77 samples: 20.2%), H (50 samples: 13.1%), and C (43 samples: 11.3%) groupings (Figure 3). Spatio-temporal patterns in benthic groups appeared largely random during the first year (Groups A-H except for Group C) (Figure 4). SIMPER indicated lowest average similarities in Groups A,B,and D (32.4,28.4,and 31.7%,respectively). Relatively few species, e.g., Paraprionospio coora (Polychaeta), Capitella sp. (Polychaeta), and Theora fragilis (Bivalvia), were predominantly responsible for similarity within these groups (Figure 5). Groups H-K then began to occur consistently over space and time from the second year onward (Figure 4). The highest average percent similarity occurred in Group K (45.3%) and, except for Group C, the number of species responsible for similarities generally increased  from Group A toward Group K (Figure 5). Group C occurred almost exclusively in the outermost part of Onagawa Bay (St.8) and remained essentially unchanged from the beginning to the end of our study (Figure 4). Despite relatively low percent similarity (35.3%), Group C contained high numbers of species, including several crustaceans, e.g., Byblis japonicus (Amphipoda), Ampelisca naikaiensis (Amphipoda), and Iphinoe sagamiensis (Cumacea), that contributed to average similarity within this group (Figure 5).

Environmental Parameters
Temperatures (TEMP), chlorophyll-a concentration (CHLa), and DO clearly demonstrated an annual (seasonal) cycle at all nine stations (Figures 6A,C,D). TEMP varied little among stations (Figure 6A), whereas CHLa was disproportionately higher at times at several locations (e.g., St.1, 3, 12, and 17) ( Figure 6C). Further, at St.1 only, summer minimum DO became much lower in recent years ( Figure 6D). Values for both salinity (SALI) and turbidity (TURB) fluctuated between 33-34 and 0-10, respectively, although these ranges were exceeded in several instances (Figures 6B,E). COD and sediment fraction of both sand (SAND) and mud (MUD) varied distinctively among stations, with St.1 and St.8 often at opposite ends of the spectrum (Figures 6F,H,I). COD at St.1 was almost always the highest and St.8 the lowest, with some distinct peaks for most stations at the beginning and toward the end of the study period ( Figure 6F). Sediment grain size at St.8 consistently contained a very high proportion of sand and little mud (Figures 6H,I). Nonetheless, a higher fraction of mud and lower fraction of sand characterized all stations at the beginning of our study, but stations gradually became increasingly sandy toward the end (Figures 6H,I). Relatively high fractions of gravel recurred at St.17 and, less frequently, at St.1 and 11 ( Figure 6G).

Aquaculture Facilities in Onagawa Bay
In March 2011 (just before the 2011 disaster), the 1217 longlines and 119 fish cages occurred mainly along low-energy coastal waters and inlets in Onagawa Bay (Figure 7A). After the complete removal of these facilities by the 2011 tsunami, 135 long-lines and 56 fish cages were re-constructed by February 2012 (Figures 7B,C), corresponding to 11.1 and 47.1% recovery, respectively (Figures 8A,B). By March 2017, the total surface area occupied by the long-lines and fish cages became 0.115 and 0.016 km 2 , i.e., 60.8 and 74.8% recovery, respectively (Figures 7D-H, 8C,D).

Macrofaunal Community Structure and Environmental Drivers
The RELATE tests showed that the macrofauna community resemblance matrix correlated most strongly with distance to the nearest long-line culture (D.LINE, ρ = 0.485, p < 0.001) followed by fetch length (FETCH, ρ = 0.476, p < 0.001) and distance to the nearest fish cage culture (D.CAGE, ρ = 0.464, p < 0.001) ( Table 4). The other aquaculture-related variables, namely A.LINE and A.CAGE, also correlated significantly with macrofauna community structure ( Table 4). D.LAND, sediment fraction of both sand (SAND) and mud (MUD), and water depth (DEPTH) also correlated significantly with the macrofauna resemblance matrix ( Table 4). The BEST analysis showed that the combination of FETCH, SAND, and A.LINE produced the highest correlation, followed by FETCH, MUD, and A.LINE (Table 4). Although the variables of D.LINE, FETCH, D.CAGE, SAND, MUD, A.LINE, DEPTH, and A.CAGE, appeared repeatedly in the five best subsets selected in the BEST procedure, no other environmental parameters in Table 4 were selected in a single combination.
We plotted values for macrofaunal abundance, biomass, species richness, and Shannon index H against cluster groups A-K identified by the SIMPROF test (Figures 9A-D). Groups C and G both showed markedly higher values for these biological metrics. Apart from these two groups, all biological metrics generally increased from Group A to Group K which also appeared to correspond to sampling date from the oldest to the newest (Figures 9A-D). In terms of proximity to the nearest aquaculture facilities, Group C mostly contained the farthest stations, whereas Group G stations occurred in closer proximity (Figures 9E,F). Apart from these groups, values for D.LINE and D.CAGE generally decreased from Group A to Group K with Group K associated with sites closest to the aquaculture facilities (Figures 9E,F). Groups C, I, J, and K tended to coincide with sites with increasing or highest cover of aquaculture facilities, whereas Groups A, B, and G characterized low aquaculture sites (Figures 9G,H). High fractions of sand characterized Group C and, to a lesser extent, Group G (Figures 9I,J). Apart from these two groups, values for SAND generally increased but MUD decreased from Group A to Group K, respectively, with Group K as somewhat intermediate (Figures 9I,J). Groups C and J were both significantly associated with exposed and deeper sites, whereas Group I was constrained to more sheltered and shallow locations (Figures 9K,L). Both Groups G and K were largely associated with stations with intermediate exposure (Figures 9K,L).
Finally, mean values for macrofaunal abundance and species richness between 2015 and 2018 did not differ significantly from before the 2011 disaster (2008)(2009)(2010)(2011) at most of the sampling stations except for St.6/TE8 (Figures 10A,B). The mud fraction (MUD) became significantly higher for some stations after the 2011 tsunami (Figures 10C,D). The gravel fraction was substantially higher around St.17/TE4 before the 2011 disaster (Figure 10E), and values for COD were significantly higher around St.1/TE1 in more recent years compared to pre-2011 values.

DISCUSSION
Our multivariate assessment of benthic macrofaunal community structure following the 2011 Tohoku earthquake and tsunami showed general increases in macrofaunal abundance, biomass, and species diversity during the study period, indicating substantial recovery of seafloor biota and coastal environmental health within 7 years. These community changes correlated with a combination of both natural and anthropogenic factors, adding insights into the mechanisms of community recovery from a catastrophic event.
Of the 11 macrofaunal cluster groups identified, exceptionally high abundances, biomass, and species diversity characterized Groups C and G. Group C occurred almost exclusively at the outermost site of St.8, some distance offshore from the outer boundary of Onagawa Bay. Physical conditions at St.8, particularly sediment grain size and exposure, differed greatly from the sampling stations inside Onagawa Bay. High species diversity and higher proportions of crustaceans (including amphipods and cumaceans), polychaetes, and bivalves characterized Group C community structure; these taxa strongly resembled those typical of high-energy sandy bottom environments. Group C persisted at St.8 throughout our study, suggesting that St.8 could act as a reference site with minimal or negligible physical impacts of the 2011 tsunami and the re-establishing process of aquaculture facilities.
In contrast, Group G occurred only in the earliest stage of our study up until September 2012 and only in the presumably low-energy waters at St.11, 16, and 17. Despite their inshore locations, however, sediment granulometry in Group G varied, characterized by high fraction of sand, which was remarkably consistent with Group C at St.8. This resemblance suggests that Group G stations may also be hydrodynamically high-energy sites (Takahashi et al., 2017) that could favor establishment of rich macrofaunal communities. In addition, Group G occurred near sites where aquaculture facilities were first re-installed in Onagawa Bay. Long-lines and fish cages were both re-constructed at and expanded from these locations (i.e., around St.11,16,and 17). Coastal aquaculture operations can alter pelagic-benthic energy fluxes, enhancing flux of organic matter to the bottom and altering local sediment characteristics and benthic community composition (Crawford et al., 2003). Environmental effects of shellfish and ascidian aquaculture on seafloor environments may differ from those associated with finfish culture in cages (Weise et al., 2009). Cultured shellfish and ascidians feed entirely on naturally occurring phytoplankton in the water column, and captured food and nutrients return to the environment as undigested waste or feces, which falls to the seafloor and may become important food source for benthic deposit feeders (Shumway et al., 2003;McKindsey et al., 2006). Chemicals and excess nutrients from food and feces associated with finfish farming may have similar effects but intense operations can also disturb benthic communities (Borja et al., 2009;Forchino et al., 2011;Tomassetti et al., 2016). For Group G, such organic enrichment generated by both types of aquaculture may have fuelled benthic recovery, particularly if benthic communities were impoverished from the immediate effects of the 2011 disaster via, e.g., vigorous influx of tsunami deposits/resuspended mud (Seike et al., 2013(Seike et al., , 2016Toyofuku et al., 2014;Kanaya et al., 2015) or spills of heavy oil/toxic chemicals (Abe et al., 2015). Group G occurred near the aquaculture facilities (D.LINE and A.LINE), likely explaining the markedly high abundance and biomass in this group.
One key finding of our study was that all the biological metrics for the cluster groups (i.e., abundance, biomass, species richness, and Shannon index H ) increased from Group A to Group K (except for Groups C and G). Comparison with the pre-tsunami data (2008)(2009)(2010)(2011) confirmed that the benthic macrofaunal community largely recovered over the duration of this study. A combination of distance to or total area of the aquaculture facilities, wind/wave fetch length, and sediment granulometry best explained variability in macrofaunal community structure. D.LINE and D.CAGE generally decreased from Group A to Group K suggesting significant, positive effects of the physical presence of aquaculture facilities on benthic macrofaunal populations, acknowledging that aquaculture prior to the disaster may have altered the benthos relative to a pristine condition. A.LINE and A.CAGE generally increased from Group A to Group K suggesting that increasing extent of aquaculture facilities can also positively influence benthic macrofaunal populations. As of March 2017, the total surface area occupied by the long-lines and fish cages were 0.115 and 0.016 km 2 , representing only 0.42 and 0.06% of the total area of Onagawa Bay, respectively. However, recent studies suggest that bio-deposits produced from aquaculture operations can disperse and spread through the water column and thereby affect benthic community at a much broader spatial scale than the immediate vicinity of aquaculture facilities (up to several kilometers) (Borja et al., 2009;Weise et al., 2009;Yokoyama, 2010;Sarà et al., 2011). If sphere of influence extends distances of 50, 100, and 150 m from each of the long-lines as suggested, we estimate the potential footprint of aquaculture operations on the seafloor to be 4. 70, 6.93, and 8.96 km 2 , or 17.3, 25.5, and 33.0% of the total area of Onagawa Bay, respectively. This estimate points to a need to assess accurately the spatial extent of benthic impact resulting from the re-establishment of aquaculture facilities and associated fishery operations in this region.
With respect to changes in macrofauna and sediment grain size following the earthquake, the sand fraction (SAND) and the mud fraction (MUD) clearly increased and decreased, respectively, from Group A to Group K. In estuarine and coastal soft-bottom habitats, sediment granulometry critically influences benthic macrofaunal community (e.g., Constable, 1999;Thrush et al., 2004;Anderson, 2008). The trend of gradually For the BEST analysis, we present the five best combinations of environmental variables that generated the highest Spearman's rank correlations (ρ).  Table 1).
increasing sand and decreasing mud across Onagawa Bay over time suggests these granulometric changes also affected the dynamics of macrofaunal community structure. Before the 2011 disaster, the sand fractions were significantly higher at some stations prior to anomalous deposition of mud across the whole bay following vigorous erosion and deposition of fine-grained sediments associated with a strong tsunami-induced current (Seike et al., 2016). Abrupt increase in fine-grained sediment could adversely affect benthic community through smothering and clogging feeding structures of suspension-feeders (Thrush et al., 2004), which may have wiped out benthic communities immediately after the 2011 tsunami. Despite gradual decreases in mud fractions, the sand fractions have not yet returned to pre-tsunami levels at some locations. Different rates of removal of fine-grained sediments could reflect different hydrodynamic processes operating at different locations. Our study used wind fetch length (FETCH) as an indicator of exposure and our analysis indicated that Groups I and J were constrained to more sheltered (low-energy) and more exposed (high-energy) locations, respectively, whereas Group K occurred at locations of intermediate exposure (i.e., FETCH: Group J > Group K > Group I). However, the dominance of sand fractions associated with Group K > Group J > Group I suggests that hydrodynamics alone did not control changes in the sediment particle size composition. Links between sediment granulometry and benthic macrofaunal communities have yet to be fully critically evaluated because coastal aquaculture operations can enhance flux of organic matter to bottom sediments which may also alter local sediment characteristics (Crawford et al., 2003). Overall, the relative influence of hydrodynamic processes and the flux of organic matter generated from the aquaculture operations on local variation in sediments is unclear, pointing to a need for future research on how distribution of aquaculture operations relates to sediment characteristics of seafloor habitat and how other environmental parameters such as coastal topography and/or hydrodynamic processes influence benthicpelagic dependencies across Onagawa Bay. Our study demonstrated the primary importance of coastal aquaculture operations, along with sediment granulometry and coastal topography (exposure), in determining the occurrence and distribution of benthic macrofauna and thereby influencing recovery in abundance, biomass, and diversity of benthic community at an ecosystem-scale following even a catastrophic natural disaster such as the 2011 Tohoku earthquake and tsunami in Japan. The potential links between the extent of coastal aquaculture operations and spatio-temporal distributional responses of benthic macrofauna as demonstrated in our study illustrate the need to improve understanding of the influence of multiple coastal anthropogenic operations on the structure and dynamics of the marine bentho-pelagic interactions that, in turn, affect the functioning of the whole marine ecosystem. Application of such research to other locations can inform marine policy issues relating to post-tsunami and/or aquaculture operations.

AUTHOR CONTRIBUTIONS
TF led the overall data analysis, writing of the paper, and the statistical analyses. AKi and KK coordinated the project and the survey. KK, TF, MK, YN, AKa, DT, YG, and HA performed the sampling at sea. KK, MK, and YN led the processing and management of the survey data. TF and KK led the analysis of the macrofauna and environmental data. HM, CY, and TF performed the analysis of the satellite imagery.

FUNDING
This study was performed as part of the "TEAMS" project and supported by grants-in-aid for scientific research by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT/JSPS).