Lei Xing

Yong Chen
email: ychen@maine.edu

Robert Boenish
email: robert.boenish@gmail.com

Kisei R Tanaka
email: kisei.tanaka@noaa.gov

Nicolas Barrier
email: nicolas.barrier@ird.fr

Yiping Ren

Evaluating the impacts of fishing and migratory species in a temperate bay of China using the ecosystem model OSMOSE-JZB

Small-scale fisheries (SSFs) play a vital role in the sustainability of local economies. Migratory species moving into and out of an ecosystem may influence the dynamics of local fish communities and SSFs. We used the end-to-end model, OSMOSE-JZB (Object-oriented Simulator of Marine ecOSystEms), to evaluate the impacts of fishing and a migratory shrimp (Trachypenaeus curvirostris) on the ecosystem of Jiaozhou Bay, China. Increased fishing intensity (i.e., annual fishing effort) resulted in the decline of four ecological indicators, including the total biomass of the community, mean trophic level of the community, inverse fishing pressure, and large fish index. Compared to managing fish stocks under uniform fishing mortality over the fishing season, landings and com-munity biomass were higher when a "race to fish" (i.e., large catches in a short period) occurred. The results suggested that managing fishing seasonality (i.e., temporal allocation of fishing effort) could mitigate the negative impact of fishing intensity. Two resident high-trophic-level fishes were sensitive to changes in fishing intensity and fishing seasonality. The changes in trophic interactions had larger impacts on species at low trophic levels than fishing. Pearson's correlation analysis showed that T. curvirostris biomass was negatively correlated with the biomass of resident species and positively correlated with the biomass of other migratory species. We also found that fishing changed the impact of varying T. curvirostris migratory biomass on other species. Resident species were more sensitive to changes in fishing and T. curvirostris migratory biomass than other migratory species. We argue that SSFs management can benefit from the development of temporal fishing strategies and consideration of trophic interactions stemming from migratory species dynamics.

Introduction

Small-scale fisheries (SSFs) constitute an important source of global fisheries landings, employ millions of people, and support food and nutrient security for more than a billion people (Food and Agriculture Organisation (FAO, 2018). SSFs have been increasingly recognized for their significant socio-economic importance to the local economy and development [START_REF] Teh | Contribution of marine fisheries to worldwide employment[END_REF][START_REF] Said | Aligning the sustainable development goals to the small-scale fisheries guidelines: a case for EU fisheries governance[END_REF]. Concerns about the viability and sustainability of SSFs have grown in recent years [START_REF] Salas | Viability and Sustainability of Small-Scale Fisheries in Latin America and The Caribbean[END_REF]. The multispecies and multigear nature of SSFs requires managing fish stocks at an ecosystem level [START_REF] Herrón | Toward ecosystem-based assessment and management of small-scale and multi-gear fisheries: insights from the tropical eastern Pacific[END_REF]. SSFs management remains less well studied compared to large-scale industrial fisheries. Data limitations impede sound stock assessment and decision-making processes for SSFs management [START_REF] Ramírez | Improving stock assessment and management advice for data-poor small-scale fisheries through participatory monitoring[END_REF].

Overfishing threatens the sustainable development of SSFs [START_REF] Gough | Evidence of overfishing in small-scale fisheries in Madagascar[END_REF]). Fishermen's behaviors are inextricably tied to stock status [START_REF] Kiyama | The impact of stock collapse on small-scale fishers' behavior: evidence from Japan[END_REF]. The importance of understanding the behavioral response of fishermen to varying resource conditions is becoming widely recognized for achieving sustainable management [START_REF] Fulton | Human behaviour: the key source of uncertainty in fisheries management[END_REF]. Large-scale environmental changes not only influence spatio-temporal variations in natural resources and fishery productivity (Moullec et al., 2019) but also systematically increase uncertainty for SSFs management [START_REF] Finkbeiner | Local response to global uncertainty: insights from experimental economics in small-scale fisheries[END_REF]. Some work has suggested that changes in migratory species biomass may alter the stability of local ecosystems [START_REF] Mormul | Migratory benthic fishes may induce regime shifts in a tropical floodplain pond[END_REF][START_REF] Stockwell | Habitat coupling in a large lake system: delivery of an energy subsidy by an offshore planktivore to the nearshore zone of Lake Superior[END_REF]. Shifts in ecosystem function further challenge the competing goals of balancing the short-term needs and long-term conservation of natural resources in SSFs management.

Marine ecosystems are becoming fragile as a result of decadal overfishing [START_REF] Pauly | The fisheries of the South China Sea: major trends since 1950[END_REF]. China's central government and provincial governments have implemented management tools to mitigate harmful fishing impacts in the last decade, including marine protected areas (MPAs), a summer fishing moratorium, and minimum mesh size regulations [START_REF] Su | Evolution of marine fisheries management in China from 1949 to 2019: How did China get here and where does China go next?[END_REF]. Nonetheless, reports of depleted fish stocks have increased in recent years, and species at lower trophic levels are becoming more dominant in Chinese coastal waters [START_REF] Szuwalski | High fishery catches through trophic cascades in China[END_REF]. Fishermen often exhibit a "race to fish" behavior where substantial fishing effort is concentrated immediately following a seasonal fishing closure or during a spawning season [START_REF] Wang | The Analysis on Trajectory Simulation of Single Otter Trawl Vessels and Its Production Estimates As Well as Spatio-Temporal Characteristics in Zhoushan Fishing Ground in Autumn and Winter[END_REF][START_REF] Sys | A temporal race-for-fish: the interplay between local hotspots of flatfish and exploitation competition between beam trawlers after a seasonal spawning closure[END_REF]. The warming climate influences migratory fish stocks, such as Trichiurus japonicas and Larimichthys polyactis [START_REF] Yuan | Impact of climate change on nekton community structure and some commercial species in the offshore area of the northern East China Sea in winter[END_REF]. These factors further complicate the management of SSFs, which greatly depend on local fish stocks.

China's government launched the "Blue Granary" concept aiming to safeguard food security in China [START_REF] Yang | Development ideas and implementation approaches of blue granary scientific and technological innovation in China[END_REF]. Ecosystem-based fisheries management (EBFM) is emphasized as a national direction and pathway for fisheries management. EBFM seeks a broadening consideration of trade-offs between managing fisheries for socio-economic profit and ecosystem stability [START_REF] Pikitch | Ecosystem-based fishery management[END_REF]. Ecosystem modeling approaches are recognized as a necessary complement to single-species modeling approaches in moving towards EBFM [START_REF] Hilborn | Future directions in ecosystem based fisheries management: a personal perspective[END_REF]. Ecosystem models, such as Ecopath and size-spectrum model, have been used to advise Chinese fisheries management [START_REF] Zhang | An evaluation of implementing long-term MSY in ecosystem-based Fisheries management: incorporating trophic interaction, bycatch and uncertainty[END_REF][START_REF] Zhang | Evaluating fishing effects on the stability of fish communities using a size-spectrum model[END_REF][START_REF] Zhang | Dynamic monitoring and analysis of number of fishing vessel voyages in Hainan Province based on Beidou position data[END_REF][START_REF] Dai | Current researching status of marine fishery ecosystem in China outlined basing on Ecopath model[END_REF]. To date, few studies have examined how variations in the temporal allocation of fishing effort and migratory species biomass influence the conservation of natural resources at an ecosystem level [START_REF] Chen | Bibliometric-based analysis of research progress in science of fishery resources[END_REF].

In this study, we simulated a small-scale mixed species fishery in Jiaozhou Bay, China. Jiaozhou Bay is a semi-enclosed bay located in the southeastern Shandong Peninsula, China. The fish community structure in Jiaozhou Bay has markedly shifted from large demersal species to small pelagic species [START_REF] Ma | Comparative study of Jiaozhou Bay ecosystem based on an Ecopath model[END_REF]. As a remediation strategy, a seasonal fishing closure has been implemented for many years [START_REF] Wang | A Strategic Study on Development of Marine Ecology Fishery in Blue Economic Area of Shandong Peninsular[END_REF], but it has been ineffective in avoiding further depletion [START_REF] Zhu | Research on the Effects of China's Summer Fishing Moratorium -A Perspective of Institutional Analysis[END_REF]. Apart from seasonal fishing closures, there are few regulations governing the temporal allocation of fishing effort. The biomass of species that migrate into and out of Jiaozhou Bay has also shown changes in ecosystem productivity [START_REF] Xu | Tempo-spatial distribution characteristics of fish resources in Jiaozhou Bay[END_REF]. We used the end-to-end model, OSMOSE-JZB (Object-oriented Simulator of Marine ecOSystEms, OSMOSE), which was developed to simulate the full ecosystem dynamics of Jiaozhou Bay [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF], as the operating model in our simulations. We address the following three issues regarding SSFs management in China: (1) impacts of variability in intensity and seasonality of fishing effort on fish stocks;

(2) detection of species-specific responses to varying migratory species biomass; and (3) evaluation of the combined effect of fishing and migratory species on different species. We believe these issues have broad applicability in Chinese waters and beyond.

Materials and methods

Operating model: OSMOSE-JZB

OSMOSE-JZB is an application of the multispecies individual-based model OSMOSE that simulates the ecosystem dynamics of Jiaozhou Bay from plankton to top predators (Fig. 1). The modeled food web in OSMOSE-JZB contains low-trophic-level (LTL) and high-trophic-level (HTL) groups [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]. The spatio-temporal distribution of LTL biomass is simulated by a biogeochemistry model, NEMURO (North Pacific Ecosystem Model Used for Regional Oceanography; [START_REF] Aita | Interdecadal variation of the lower trophic ecosystem in the northern Pacific between 1948 and 2002, in a 3-D implementation of the NEMURO model[END_REF], coupled to an ocean physics model, FVCOM (Finite Volume Coastal Ocean Model; [START_REF] Chen | An unstructured grid, finite-volume, threedimensional, primitive equations ocean model: application to Coastal Ocean and estuaries[END_REF]. The OSMOSE model [START_REF] Shin | Exploring fish community dynamics through size-dependent trophic interactions using a spatialized individual-based model[END_REF] simulates the full life cycle of HTL groups from egg to adult stages. The outputs of the FVCOM-NEMURO model are used to force the OSMOSE model. Spatio-temporal distributions of LTL biomass are used as prey fields to HTL groups. The LTL groups are classified as small Fig. 1. Impacts of fishing and migratory species in Jiaozhou Bay, China (JZB). A grid map of the study area is shown on the left. phytoplankton (PS), large phytoplankton (PL), small zooplankton (ZS), large zooplankton (ZL), and predatory zooplankton (ZP). The parameter Plank.access is used to calculate the proportion of available LTL biomass to HTL groups. As we focus on the impacts of fishing and migratory species on HTL groups, the configuration of FVCOM-NEMURO is not presented in this study. A detailed description of our model can be found in [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]. Here, we briefly introduce the model structure and parameterization of the OSMOSE model.

The OSMOSE model is a multispecies individual-based model that simulates species life cycles and trophic interactions. The time step of the OSMOSE model is half of a month. The model is built based on data from four stratified random bottom trawl surveys conducted during 2011 in Jiaozhou Bay, China. There are 14 HTL groups modeled in the OSMOSE model, including 7 fishes, 2 crabs, 5 shrimp species, and 2 cephalopods (Table 1). Fish school is the basic modeling unit for each HTL group. There are six modeled ecological processes, including spatial distribution, species migration, predation process, somatic growth, species-specific reproduction, and various mortality sources, such as fishing.

Spatial distribution

Our simulated Jiaozhou Bay consists of 1435 grid cells (0.01 • ×0.01 • ; 35.95 • -36.30 • N; 120.04 • -120.45 • E; Fig. 1). We assumed fish could move throughout Jiaozhou Bay. At each time step, the diffusions of fish schools follow a "random walk" process, and the range of "random walk" is set based on their swimming ability. Fish and other species were set to 2 and 1, respectively (Table A1).

Species migration

The HTL groups can be divided into (1) resident species that spend all their life stages in Jiaozhou Bay and (2) migratory species that move in and out of the simulated area (Table 1). The biomass of migratory species was set to zero at time steps corresponding to migrating out of the simulated area. The timing of migration is set based on survey data and relevant references (Table A1). The flux of species biomass (Bio flux ) migrating into Jiaozhou Bay was estimated from the model calibration. Each age class for Bio flux is assumed to be uniform.

Predation process

The model assumes size-based opportunistic predation processes [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]. Prey availability is dependent on the prey-predator size ratio and spatio-temporal co-occurrence of prey and predators in the same cell at a given time step. The minimum and maximum prey-predator size ratios were set based on stomach content analysis and relevant studies (Table A1).

Somatic growth

The somatic growth of individuals is based on the amount of food obtained in predation processes. Individuals grow after ingesting enough food at a given time step. Alternatively, they stop growing or Note: In a simulated year, there were 24 time steps numbered as 0-23. Migratory species were marked with *. Annual M fishing was the sum of fishing mortality rate at each time step. The annual M fishing (initial value) was the same to that set in [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]. Octopus sp. (SP5*) migrated out of Jiaozhou Bay at a given age. This HTL group occurred in the simulated domain all the year around. Therefore, its time steps for moving out of simulated domain was showed with null.

Fig. 2. Flow chart of the processes for simulation design. The processes of scenario design and impact evaluation are marked with dotted rectangles and solid rectangles, respectively.

starve due to the lack of food. The maximum ingestion rate, set to 3.5 g of food per gram of body weight per year, is used to define the upper limitation of food intake [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]. The critical predation efficiency ξ crit , which is equal to 0.57, determines whether individuals ingest enough food to support somatic growth [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]. If predation efficiency ξ is higher than ξ crit , the growth of individuals is calculated by the von Bertalanffy model. The weight-length relationship formula calculates the individual weight at a given size (Table A1).

Species-specific reproduction

The model considers the species-specific seasonality of reproduction. The number of eggs spawned in the spawning season is related to spawning stock biomass and relative fecundity. No eggs are produced outside of the spawning season. The spawning stock refers to individuals larger than the minimum maturity size in the spawning season. Each species' spawning season and maturity size were derived from relevant studies (Table A1).

Various mortality sources

The five mortality sources included in the model are predation mortality (M predation ), starvation mortality (M starvation ), fishing mortality (M fishing ), larval mortality (M larval ), and additional natural mortality (M natural ). Each mortality source is calculated based on a stochastic mortality algorithm in which five mortalities compete with each other. M predation and M starvation are related to the predation processes and calculated by the model instead of being assigned to invariant values. Additional natural mortality (M natural ) represents the sum of mortality sources that are not explicitly represented in the OSMOSE model. Larval mortality (M larval ) is the sum of the mortality of non-fertilization and export of eggs and first feeding larvae.

The seasonal fishing closure from Jun. 1st to Sep. 1st was considered. No catches occur during the period of closed fishing. The minimum catch size is set to 5 cm so that only individuals larger than the minimum catch size will be harvested during the fishing season. Fishing mortality rates of HTL groups are assumed to be uniform over the fishing season. The partial parameters and equations are summarized in Appendix A.

The unknown parameters, which cannot be directly obtained from survey data and relevant studies, are estimated from the model calibration by fitting the simulated species biomass to observed survey data collected in 2011. An evolutionary algorithm (EA) was employed to find the optimal combination of unknown parameter values (Oliveros-Ramos et al., 2017). There were five types of unknown parameters: the proportion of available biomass of LTL groups to HTL groups (Plank.access), flux of migratory species biomass (Bio flux ), species-specific relative fecundity (RF), larval mortality (M larval ), and annual fishing mortality (M fishing ). The unknown parameters were estimated within a reasonable range at the initiation of model calibration. Model calibration was hierarchical and estimated the parameters sequentially along the three phases: Plank.access and Bio flux were estimated in all phases; RF and M larval were estimated in the last two phases; and M fishing was estimated in the last phase. Note that the RF and M larval of species that did not spawn in Jiaozhou Bay were set to null. Further details are presented in [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

Simulation design

Our study focused on an input control rule that restricted fishing effort. This rule is most commonly applied to Chinese fisheries. We set the fishing mortality rate at each time step to simulate the variability in fishing intensity (i.e., annual total fishing effort) and fishing seasonality (i.e., temporal allocation of fishing effort). Trachypenaeus curvirostris (SP12) is a fast-growing, short-lived migratory shrimp in Jiaozhou Bay. T. curvirostris biomass changed greatly in the water of the south Shandong Peninsula from 2010 to 2017 (see Appendix B). Three scenarios (S1 to S3) were designed to evaluate the possible impacts of fishing and T. curvirostris (Fig. 2): (Table 1). The sum of percentages for each HTL group in a simulated year was equal to 100 %. FSa, FSb, and FSc correspond to managing species under uniform fishing mortality, large catches landing immediately after the end of the seasonal fishing closure, and large catches landing at the start and end of the fishing season, respectively. The fishing mortality rate of one HTL group was set to zero during seasonal fishing closure or when they migrated out of Jiaozhou Bay. Type 1 refers to species that occur in the modeled area year round, including O. oratoria (SP0), small shrimp (SP1), Octopus sp. (SP5), S. schlegelii (SP6), A. hexanema (SP8), J. belangerii (SP11), and C. joyneri (SP13). Types 2 to 5 are species that migrate out of the modeled area at a given time. Type 2 includes C. bimaculata (SP2), Loligo sp. (SP4), and T. kammalensis (SP9). Type 3 includes C. japonica (SP3) and T. curvirostris (SP12). Types 4 and 5 include P. fangi (SP7) and L. tanakae (SP10), respectively. Time steps for moving out of the simulated domain are found in Table 1 and marked with red lines on the X axis.

(1) Scenario S1 quantified the combined effect of fishing intensity and fishing seasonality. Seven fishing intensity levels were evaluated for the HTL groups. The annual fishing mortality rate of each HTL group ranged from 0.85 to 1.15 times their initial values (Table 1) with a step of 0.05. At each level of fishing intensity, three fishing seasonality types (FSa, FSb, and FSc) were considered based on relevant studies [START_REF] Zhu | Research on the Effects of China's Summer Fishing Moratorium -A Perspective of Institutional Analysis[END_REF][START_REF] Wang | The Analysis on Trajectory Simulation of Single Otter Trawl Vessels and Its Production Estimates As Well as Spatio-Temporal Characteristics in Zhoushan Fishing Ground in Autumn and Winter[END_REF][START_REF] Zhang | Evaluating fishing effects on the stability of fish communities using a size-spectrum model[END_REF][START_REF] Zhang | Dynamic monitoring and analysis of number of fishing vessel voyages in Hainan Province based on Beidou position data[END_REF]. The fishing seasonality FSa simulated the situation that the "race to fish" did not occur (i.e., uniform fishing mortality for HTL groups over the fishing season that was the same as [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]). The fishing seasonality FSb simulated the "race to fish" occurring at the beginning of the fishing season (i.e., large catches immediately after the end of the seasonal fishing closure). The fishing seasonality FSc simulated the "race to fish" occurring at the start and end of the fishing season (i.e., large catches in May and September). Fig. 3 shows the percentage distribution of the fishing mortality rate for HTL groups in a simulated year. The fishing mortality rate of each HTL group at each time step was set equal to the product of the annual fishing mortality rate and the corresponding percentage at each time step. Fishing seasonality did not change the value of the annual fishing mortality rate; (2) Scenario S2 investigated the responses of other species to changes in T. curvirostris migratory biomass. We set the flux of migratory biomass (Bio flux ) to simulate the variability in T. curvirostris migratory biomass. The initial Bio flux of T. curvirostris was 196.92 t [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]. Given the varying nature of migratory biomass in the ecosystem, a Monte Carlo simulation approach was employed to randomly generate multipliers for the Bio flux of T. curvirostris from the uniform distribution. The bounds of the uniform distribution and number of multipliers are presented in Fig. 2; (3) Scenario S3 evaluated the impact of T. curvirostris migratory biomass on other species under different fishing scenarios. There were 9 fishing scenarios in S3. Three levels of fishing intensity were considered: low fishing intensity (annual fishing mortality rate of each HTL group set to 0.85 times the initial values, which are shown in Table 1) (Fx0.85), current fishing intensity (annual fishing mortality rate was the same as the initial values) (Fx1), and higher fishing intensity (annual fishing mortality set to 1.15 times the initial values) (Fx1.15). Three fishing seasonality types (FSa, FSb, and FSc) were considered and set the same as those in scenario S1. The fishing mortality rate at each time step was calculated as in S1. In each fishing scenario, the Bio flux of T. curvirostris changed at 15 levels ranging from 0.3 to 1.7 times the initial values set in [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF] with a step of 0.1.

Each scenario ran for a total of 120 years. To stabilize the model, all simulations had an 80-year period of burn-in under the same configuration of the initial OSMOSE-JZB model built by [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]. The final 40 years were used for analysis. Every simulation subscenario in S1 and S3 was repeated 100 times to account for model stochasticity. The model outputs were averaged over the 100 simulation runs. In S2, each of the above simulations was repeated 10 times to account for model stochasticity. There were a total of 2000 (200 × 10) simulations in S2.

The initial OSMOSE-JZB model repeatedly ran 100 times. The outputs simulated by the model configuration of [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF] were used as the reference/base state representing the ecosystem of Jiaozhou Bay in 2011 to compare with other scenarios.

Impact evaluation

In scenario S1, we evaluated the fishing impacts at the community, population, and individual levels. The changes in the community were described by four ecological indicators: total biomass of the community (Bio com ), mean trophic level of the community (mTL com ), inverse fishing pressure (B/L), and large fish index (LFI 20 ). Species biomass (Bio) was used to measure the stock status at the population level. Sebastes schlegelii and Johnius belangerii are slow-growing, long-lived, and hightrophic-level fishes with small populations in Jiaozhou Bay. They are usually exposed to high fishing pressure due to their high commercial value. Amblychaeturichtys hexanema is a small, low-trophic-level fish with abundant biomass. At the individual level, fishing (M fishing ) and predation (M predation ) mortality rates described the responses of individuals of two age groups (age-0 individuals who were below 1 year old and age-1 individuals who were between 1 and 2 years old) for S. schlegelii, A. hexanema, and J. belangerii. M predation described the indirect impact of fishing, which influenced fish individuals by changing predator pressure. The sum of predation and fishing mortality rates (i.e., M sum ) indicated the combined effect of pressure from fishing and predators.

In scenario S2, Pearson's correlation coefficient was used to test the relationship between the biomass of T. curvirostris and other species. In scenario S3, we used a linear regression analysis to quantify the species- specific response to the variability of T. curvirostris biomass under different fishing scenarios [START_REF] Halouani | An end-to-end model to evaluate the sensitivity of ecosystem indicators to track fishing impacts[END_REF]. There were 9 linear regression models for each species in S3. The formula of the linear regression model is:

Relatvie Bio = Bio sp,s Bio sp,ref (1) Relatvie Bio = a × Multiplier + b (2)
where Relative Bio is the ratio of species biomass in scenario s (Bio sp,s ) and that in the reference state (Bio sp,ref ); Multiplier is the ratio of migratory biomass of T. curvirostris in scenario s and the reference state; and a and b are the slope and intercept in the linear regression model, respectively. The slope represented the increasing impact rate of varying T. curvirostris biomass under different fishing scenarios. The P-value determined the goodness-of-fit for the linear regression model. The indicators used in the study are shown in Table 2.

Results

Fishing impacts

In general, fishing intensity had a negative relationship with the four ecological indicators (Fig. 4). The relative changes in Bio com and mTL com were smaller than those in B/L and LFI 20 . B/L was highest when fishing mortality was uniform over the fishing season (FSa) and lowest when large catches landed immediately after the end of the seasonal fishing closure (FSb). LFI 20 was highest with large catches at the beginning and end of fishing seasons (FSc). LFI 20 had a greater increase under FSc after reducing fishing intensity (Fig. 4).

The changes in fishing seasonality altered the impact of fishing intensity on different species (Fig. 5). Except for Pholis fangi, migratory species had smaller changes than resident species when fishing intensity and fishing seasonality changed. The biomass of Oratosquilla oratoria and A. hexanema declined when fishing mortality was below the reference state. A. hexanema biomass was low under FSb. Small shrimp, P. fangi, and Cynoglossus joyneri biomasses were highest under FSb and lowest under FSc. S. schlegelii and J. belangerii were subject to changes in fishing intensity. In contrast to most species, the biomass of S. schlegelii and J. belangerii declined with increasing fishing mortality. S. schlegelii biomass was high under FSc, and J. belangerii biomass was high under FSb. Fishing seasonality had a larger impact on S. schlegelii and J. belangerii under low fishing pressure levels (Fig. 5).

The fishing intensity and fishing seasonality could alter trophic interactions (Fig. 6). The M predation and M fishing of age-0 S. schlegelii were lower under FSc than under FSa and FSb. At a low fishing pressure level (Fx0.85), age-0 J. belangerii underwent high predation mortality under FSa and FSc. In contrast to S. schlegelii and J. belangerii, the M fishing of age-0 A. hexanema showed small changes. The M sum of S. schlegelii and A. hexanema below 1 year old was low under FSc. The M sum of age-0 J. belangerii was low under FSb (Fig. 6a). The M predation of S. schlegelii and J. belangerii between 1 and 2 years old showed small changes in response to changes in fishing intensity and seasonality. The M predation of age-1 A. hexanema declined after increasing fishing intensity. The M fishing of age-1 S. schlegelii was higher under FSc than under the other two fishing seasonality types. The M predation of age-1 A. hexanema declined Fig. 6. Relative changes in predation (M predation , black bars) and fishing mortality (M fishing , red bars) rates for S. schlegelii (SP6), J. belangerii (SP11), and A. hexanema (SP8) under scenario S1. There were two age groups: individuals below 1 year old (Age-0) and individuals between 1 and 2 years old (Age-1). The indicator M sum (black broken line) was the sum of relative changes in M predation and M fishing . The fishing mortality rate was set from 0.85 to 1.15 times the initial values, and there were three fishing seasonality types (FSa, FSb, and FSc). with increasing the fishing mortality of all species. The M predation of age-1 A. hexanema under FSc was higher than that under the other two fishing seasonality types (Fig. 6b).

Impacts of migratory species

Pearson's correlation analysis showed that T. curvirostris biomass had negative correlations with resident species biomass and positive correlations with migratory species biomass (Table 3). Except for Liparis tanakae, the correlation coefficient values for migratory species were higher than 0.75. The correlation coefficient values of O. oratoria and small shrimp (r < -0.9) were much lower than those of other resident species. It indicated that their biomass had strong negative relationships with T. curvirostris biomass. The correlation coefficient values of S. schlegelii, J. belangerii, and C. joyneri were lower than that of A. hexanema. The absolute value of the correlation coefficient for A. hexanema was lowest, indicating a weaker relationship between A. hexanema biomass and T. curvirostris biomass compared to other species (Table 3).

Impacts of migratory species under different fishing scenarios

The absolute values of slopes for O. oratoria, A. hexanema, C. joyneri, and 8 migratory species were lower than 0.1, suggesting that changes in T. curvirostris biomass had small impacts on these species (Fig. 7). The increase of T. curvirostris biomass had negative impacts on small shrimp, S. schlegelii, and J. belangerii (Slope < 0). The slopes of S. schlegelii and J. belangerii were lower than those of small shrimp, suggesting that changes in T. curvirostris biomass had larger impacts on S. schlegelii and J. belangerii. The slopes of small shrimp were high under low fishing pressure. In contrast, the slopes of S. schlegelii and J. belangerii declined after reducing fishing intensity. In terms of fishing seasonality, the slopes of S. schlegelii and J. belangerii were low under FSc and FSb, respectively. The slopes of small shrimp had small changes when the fishing seasonality changed (Fig. 7).

Discussion

Heavy fishing will result in the decline of large predatory fishes and the increase of small pelagic fishes [START_REF] Szuwalski | High fishery catches through trophic cascades in China[END_REF]. The S1 results highlighted the importance of managing the fishing effort to protect large fishes. [START_REF] Han | Evaluation of effects of shellfish aquaculture and capture fishery on a semi-closed bay ecosystem[END_REF] demonstrated that small fishes had an advantage against large fishes under high fishing pressure. Similarly, we found that the proportion of individuals larger than 20 cm (i.e., LFI 20 ) declined after increasing fishing intensity (Fig. 4). Two large resident fishes (S. schlegelii and J. belangerii) were more sensitive to fishing changes than low-trophic-level species (Fig. 5). The low-trophic-level species biomass decreased after reducing the fishing pressure to all species (Fig. 5). This result suggested that changes in trophic interactions had larger impacts on low-trophic-level species than fishing. Fishing can indirectly influence fish stocks by changing trophic interactions [START_REF] Lynam | Interaction between top-down and bottom-up control in marine food webs[END_REF]. In our study, trophic interactions more strongly changed A. hexanema predation mortality (i.e., M predation ) compared to fishing mortality (Fig. 6). This result suggested that A. hexanema was more sensitive to changes in trophic interactions than the direct impact of fishing. Fishing changes had larger impacts on age-1 individuals for S. schlegelii and J. belangerii than age-0 individuals (Fig. 6). In contrast to A. hexanema, the predation mortality of age-1 S. schlegelii and J. belangerii had small changes after fishing changed (Fig. 6b). This result suggested that age-1 S. schlegelii and J. belangerii had an advantage in competing with other species.

Da-Rocha et al. ( 2014) advocated for the management approach of landing large catches in a short time period to let the stock grow over a longer timeframe. Other work suggests that "race to fish" (i.e., harvesting a lot in a short time period) may erode the recruitment of fish stocks [START_REF] Wang | Modeling the effect of the seasonal fishing moratorium on the Pearl River Estuary using ecosystem simulation[END_REF]. Bio com , mTL com , and LFI 20 were highest when the "race to fish" occurred at the start and end of the fishing season (i.e., FSc) (Fig. 4). The "race to fish" seemed to be a potential option to mitigate negative fishing impacts compared to managing fish stocks under uniform fishing pressure (i.e., FSa). At a high fishing intensity level, we noted a convergence of ecological indicators (e.g., Bio com , mTL com , and LFI 20 ) under three fishing seasonality types (Fig. 4). This result suggested that high fishing intensity might undermine the effectiveness of managing fishing seasonality. Past work suggested that the growth of individuals below 1 year old could influence the stocks of S. schlegelii and J. belangerii (Xing et al., 2020a). The "race to fish" that occurred at the beginning of the fishing season (i.e., FSb) favored the growth of age-0 J. belangerii (Fig. 6a). The J. belangerii biomass was highest under the fishing seasonality FSb (Fig. 5). The fishing seasonality FSc was more beneficial to age-0 S. schlegelii, in which the pressure of fishing and predators was low (Fig. 6a). The S. schlegelii biomass was high under the fishing seasonality FSc (Fig. 5). SSFs in China nearly always target multiple species. Seemingly, ideal harvest strategies will be species-specific and fishery-dependent. We note it is a significant challenge to balance optimized harvest strategies for different species simultaneously.

Migratory species can influence local marine ecosystems via trophic interactions [START_REF] Mariani | Trophic impact of Atlantic bluefin tuna migrations in the North Sea[END_REF]. [START_REF] Han | A mass balanced model of trophic structure and energy flows of a semi-closed marine ecosystem[END_REF] reported that a large number of T. curvirostris will add additional trophic competition with O. oratoria and small shrimp. Similarly, the S2 results showed that T. curvirostris had strong negative correlations with O. oratoria and small shrimp (r < -0.9) (Table 3). Moreover, Pearson's correlation analysis showed that increasing T. curvirostris biomass had negative impacts on resident species (r < 0) and positive impacts on migratory species (r > 0) (Table 3). This implied the competition between resident species and migratory species. The changes in T. curvirostris migratory biomass might present additional challenges for SSFs management. The S3 results showed that small shrimp, S. schlegelii, and J. belangerii were more sensitive to changes in T. curvirostris biomass than other species (Fig. 7). Moreover, fishing could influence their responses to changes in the migratory biomass of T. curvirostris. T. curvirostris had a larger impact on small shrimp under high fishing pressure. Fishing seasonality did not seem to influence the impact of T. curvirostris on small shrimp. The impacts of T. curvirostris on S. schlegelii and J. belangerii were larger under low fishing pressure. At the same fishing pressure, S. schlegelii was more sensitive to changes in T. curvirostris biomass under FSc, while J. belangerii was more sensitive to changes in T. curvirostris biomass under FSb. The compound effect of fishing and T. curvirostris migratory biomass further complicated SSFs management in Jiaozhou Bay.

SSFs management in China has become complicated in recent years. Except for the summer fishing moratorium, there are no specific management policies to control the temporal allocation of fishing effort during the fishing season. Such policies result in a "race to fish" each year following the lifting of the summer moratorium. Nationally, "double control" and fishing license systems are designed to restrict fishing effort [START_REF] Shen | An overview of marine fisheries management in China[END_REF]. Emerging from resource instability, "fishing less" and "fishing more" are potentially adaptive strategies by fishermen in response to resource decline [START_REF] Cinner | Responding to change: using scenarios to understand how socioeconomic factors may influence amplifying or dampening exploitation feedbacks among Tanzanian fishers[END_REF]. On the one hand, fisheries resource decline may decrease fishing effort due to increased costs and reduced catch. On the other hand, available fuel subsidies and growing market demand (i.e., price incentives) can incentivize fishermen to take more fishing trips amidst a declining resource [START_REF] Mallory | Fisheries subsidies in China: quantitative and qualitative assessment of policy coherence and effectiveness[END_REF]. Such behaviors certainly introduce uncertainty into the effectiveness of the Chinese seasonal fishing closure, "double control", and fishing license systems.

Our findings highlight the importance of incorporating the temporal distribution of fishing effort into the decision-making process. No temporal allocation of fishing effort was beneficial to all species. Balancing the relationship between focal species and other ecological components is a challenge for SSFs management. There are increasing reports of variations in stocks of migratory shrimp and cephalopods in Chinese waters, which are fast-growing and short-lived species [START_REF] Lv | Assessment of Important Fishery Resources in the South Offshore of Shandong from 2010 to 2017[END_REF][START_REF] Pang | Variability of coastal cephalopods in overexploited China Seas under climate change with implications on fisheries management[END_REF][START_REF] Wu | Inter-annual variation in the community structure of crustaceans in the Bohai Sea during summer[END_REF]. The study demonstrated that the compound effects of fishing and migratory species on the local ecosystem are non-stationary and complex. We recommend that the impact of migratory species should be considered in fisheries management, especially when migratory movement has or is predicted to markedly change.

Additional aspects need to be considered when this approach is used to inform fisheries management. The behaviors of migratory species are subject to environmental changes such as sea temperature anomalies. In addition to focusing on variations in T. curvirostris biomass migrating into Jiaozhou Bay, it is necessary to investigate the impacts caused by changes in migration timing. Migratory species have different responses to environmental variability, and future studies should quantify how varying the biomass of multiple migratory species influences food web stability. We recommend further cooperation between our approach and others. For example, generalized additive models and habitat suitability indices are commonly used to investigate the relationship between biomass and environmental factors [START_REF] Grüss | Producing distribution maps for informing ecosystem-based fisheries management using a comprehensive survey database and spatio-temporal models[END_REF][START_REF] Zeng | Gillnet data enhance performance of rockfishes habitat suitability index model derived from bottom-trawl survey data: a case study with Sebasticus marmoratus[END_REF]. Incorporating species distribution models into our model is a future research direction that will improve the realness of our simulations. We have investigated the impacts of imprecise parameters, such as mortality rates and relative fecundity, on the model performance (Xing et al., 2020b). As fishing and climate change can influence individual growth rates [START_REF] Martino | Using otolith chronologies to understand long-term trends and extrinsic drivers of growth in fisheries[END_REF], there is a need to evaluate the impacts of changing parameter values in growth sub-models.

In this study, we quantified the species-specific response to variations in fishing and migratory species biomass. In addition to controlling fishing intensity, we found that the effectiveness of recovering large predatory species can be improved by managing the "race to fish". Managing fish stocks under uniform fishing pressure may be less effective because the exploitation (high B/L) and community biomass (low Bio com ) are low. Fishing can influence the impacts of migratory species on resident species. The study provides some insights into managing SSFs in China. The modeling approach can help identify potential winwin management strategies that consider both ecosystem conservation and economic profit. Although some of our conclusions may be similar to other studies and the results are based on a particular ecosystem, the approach can be applied to other small-scale fisheries ecosystems to promote ecosystem-based management of fisheries resources.

CRediT authorship contribution statement

Declaration of Competing Interest

The authors report no declarations of interest.

Table A1

Summary of equations and partial parameters modeled in the OSMOSE-JZB model built by [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

Process Description of provisos and equations

Details of partial parameters

Spatial distribution

"Random walk" is employed by the model for mimicking species diffusion during foraging. The range of "Random walk" determines number of adjacent cells which can be used for a given species.

The range of random foraging was set based on their swimming ability. Seven fishes were set to 2, and other species were set to 1.

Migration

The species, which migrate into the simulated area, are characterized by biomass, age, and length. The flux of species biomass (Bio flux ) migrating into Jiaozhou Bay for each age class are uniform. The model allows for setting the timing of species migration. The migratory species will be removed when they migrate out of the simulated domain.

Time steps for moving out of simulated domain were showed in Predation Proviso for biomass of available prey to predator [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF]:

AP i, x, y, t ' = ∑ j ai,j⋅B j, x, y, t ' | Li srmax < Lj < Li srmin ⎧ ⎪ ⎨ ⎪ ⎩ if AP i, x, y, t ' > ri⋅B i, x, y, t ' PBi, j, Δt = ai,j⋅Bj,t * r * Bi,t AP i, x, y, t ' if AP i, x, y, t ' ≤ ri⋅B i, x, y, t ' PBi, j, Δt = ai,j⋅Bj,t
where AP i,x,y,t' represents the total prey biomass available for school i in the cell (x, y) at time step t', B j,x,y,t' is the biomass of prey school j in the cell (x, y) at t', a i,j is the availability coefficient of prey j to predator i, and sr is the suitable predator-prey size ratio. t' corresponds to the sequential time within a time step at which fish schools feed in a random order. PB i,j,Δt represents the biomass of school j preyed upon by school i;

The maximum rate of ingestion r i (r max ) was set for all species to 3.5 g of food per gram of body mass per year [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF].

The suitable predator-prey size ratio of each species (sr min and sr max ) was listed in Table A The predation efficiency ξ i was calculated in the predation process. The critical predation efficiency ξ crit , which was set to 0.57, identified whether ingested food was enough to maintain its basic metabolic rates. If ξ i was lower than 0.57, individual would starve and stop growing.

Parameters of von Bertalanffy model and weight-length relationship were listed in Table 2 of [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

Reproduction

The number of eggs spawned by sexually mature individuals at each time step within a year [START_REF] Shin | Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing[END_REF]:

Neggs,t = Ratiosex * RF * SSBs,t * seasonality
where Ratiosex is the sex ratio, RF is the relative fecundity of species s, SSB s,t is the total biomass of sexually mature fish (larger than minimum maturity size) of species s at time step t, and seasonality is the percentage of eggs produced at each time step within a simulated year.

There were eight HTL groups spawning in the Jiaozhou Bay: O. oratoria (SP0), small shrimp (SP1), C. bimaculata (SP2), S. schlegelii (SP6), P. fangi (SP7), A. hexanema (SP8), J. belangerii (SP11), and C. joyneri (SP13).

Ratio sex values of above species were assigned to 0.5.

The percentage of eggs produced at each time step was set based on relevant references which were listed in Table B.1 of [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

RF of each species was estimated from the model calibration.

Mortality

There are five mortality sources included in the model: predation mortality (M predation ), starvation mortality (M starvation ), fishing mortality (M fishing ), larval mortality (M larval ), and additional natural mortality (M natural ). The starvation mortality is calculated as below [START_REF] Travers-Trolet | Emergence of negative trophic level-size relationships from a size-based, individual-based multispecies fish model[END_REF]: The minimum catch size for each HTL group was set to 5 cm. In other words, fishing mortality occurred to individuals which were larger than 5 cm.

Mξ i = -
M ξmax for each HTL group is set to 0.3.

M fishing and M larval were estimated from the model calibration.

The values of annual M fishing , M larval , and M natural were presented in Table 2 of [START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

Percentage distribution of fishing mortality rate of each HTL group at different time step was shown in Fig. 3.

Fig. 3 .

 3 Fig. 3.Percentage distribution of fishing mortality rate for HTL groups. HTL groups were divided into five types based on the time migrating into/out of Jiaozhou Bay (Table1). The sum of percentages for each HTL group in a simulated year was equal to 100 %. FSa, FSb, and FSc correspond to managing species under uniform fishing mortality, large catches landing immediately after the end of the seasonal fishing closure, and large catches landing at the start and end of the fishing season, respectively. The fishing mortality rate of one HTL group was set to zero during seasonal fishing closure or when they migrated out of Jiaozhou Bay. Type 1 refers to species that occur in the modeled area year round, including O. oratoria (SP0), small shrimp (SP1), Octopus sp. (SP5), S. schlegelii (SP6), A. hexanema (SP8), J. belangerii (SP11), and C. joyneri (SP13). Types 2 to 5 are species that migrate out of the modeled area at a given time. Type 2 includes C. bimaculata (SP2), Loligo sp. (SP4), and T. kammalensis (SP9). Type 3 includes C. japonica (SP3) and T. curvirostris (SP12). Types 4 and 5 include P. fangi (SP7) and L. tanakae (SP10), respectively. Time steps for moving out of the simulated domain are found in Table1and marked with red lines on the X axis.

Fig. 4 .

 4 Fig. 4. Relative changes in ecological indicators under scenario S1. The fishing mortality rate was set from 0.85 to 1.15 times the initial values. There were three fishing seasonality types: uniform fishing mortality over the fishing season (FSa), large catches landing immediately after the end of the seasonal fishing closure (FSb), and large catches landing at the start and at the end of fishing season (FSc).

Fig. 5 .

 5 Fig. 5. Relative changes in species biomass under scenario S1. The fishing mortality rate was set from 0.85 to 1.15 times the initial values, and there were three fishing seasonality types (FSa, FSb, and FSc). The migratory species are marked with *.

Fig. 7 .

 7 Fig. 7. Linear regression model of the relative biomass of species except for T. curvirostris and the corresponding multiplier. The slope values are shown in the figure and indicated by the colors. The P and R 2 are shown in Fig. C.1 to Fig. C.3. The migratory species are marked with *.

Lei Xing :

 Xing Conceptualization, Data curation, Formal analysis, Methodology, Software, Writing -original draft. Yong Chen: Conceptualization, Methodology, Supervision, Writing -review & editing. Robert Boenish: Methodology, Writing -review & editing. Kisei R. Tanaka: Methodology, Writing -review & editing. Nicolas Barrier: Methodology, Software. Yiping Ren: Data curation, Supervision, Funding acquisition.

  and M ξmax are the starvation mortality rate of fish school i and maximum starvation mortality rate, respectively. ξ i and ξ crit are predation efficiency and critical predation efficiency, respectively.

Fig. C2 .

 C2 Fig. C2. Results of linear regression analysis for Loligo sp. (SP4), Octopus sp. (SP5), P. fangi (SP7), T. kammalensis (SP9), and L. tanakae (SP10) in scenario S3. The P values of these five species in the subscenarios were all lower than 0.05.

Fig. C3 .

 C3 Fig. C3. Results of linear regression analysis for small shrimp (SP1), S. schlegelii (SP6), and J. belangerii (SP11) in scenario S3. The P values of the three fishes in the subscenarios were all lower than 0.05.

Table 1

 1 Summary of high-trophic-level (HTL) groups modeled in the initial OSMOSE-JZB built by[START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF].

	HTL group	Category	Species (common name)	Annual M fishing (Initial value)	Time steps for moving out of simulated domain
	SP0 SP1 SP2* SP3* SP4* SP5* SP6 SP7* SP8 SP9* SP10* SP11 SP12* SP13	Mantis shrimp Small shrimp Small crab Large crab Loligo Octopus Large fish Small fish Small fish Small fish Large fish Large fish Medium-sized shrimp Medium-sized fish	Oratosquilla oratoria (Japanese mantis shrimp) Palaemon gravieri (Chinese ditch prawn); Parapenaeopsis tenella (Smoothshell shrimp); Alpheus japonicus (Japanese snapping shrimp) Charybdis bimaculata (Two-spot swimming crab) Charybdis japonica (Japanese swimming crab) Loligo sp. (Squid) Octopus sp. (Octopus) Sebastes schlegelii (Korean rockfish) Pholis fangi (Gunnel) Amblychaeturichthys hexanema (Pinkgray gody) Thryssa kammalensis (Kammal thryssa) Liparis tanakae (Tanaka's snailfish) Johnius belangerii (Belanger's croaker) Trachypenaeus curvirostris (Southern rough shrimp) Cynoglossus joyneri (Red tonguesole)	0.420 0.365 0.314 0.246 0.419 0.350 0.482 0.426 0.229 0.467 0.488 0.411 0.426 0.365	Null Null 0~5 0~11 0~5 Null Null 12~17 Null 0~5 0~5; 12~17 Null 0~11 Null

Table 2

 2 Definitions of indicators used in the study.

	Indicator	Symbol	Description	Source
	Total biomass of the community Mean trophic level of the community Inverse fishing pressure Large fish index Species biomass Predation mortality rate Fishing mortality rate Sum of predation and fishing mortality rates Relative indicators	Bio com mTL com B/L LFI 20 Bio M predation M fishing M sum	Biocom = denotes the biomass of a given ∑ Bioi, where Bio i i species i in the simulated area. mTLcom = ∑ i TLi * Bioi Biocom , where Bio i denotes the biomass of a given species i; Bio com indicates the total biomass of the community calculated as aforementioned computational method; TL i represents the mean trophic level of species i, and its computational method was the same as in Xing et al. (2017). B/L = Biocom Catchcom , where Catch com denotes the total catch of the community; Bio com indicates total biomass of the community calculated as aforementioned computational method. LFI20 = Bio20 Biocom , where Bio20 and Biocom indicate the biomass of individuals larger than 20 cm and total biomass of the community, respectively. The total biomass of a given species in the simulated area. Predation mortality rate provided in outputs of the model. Fishing mortality rate provided in outputs of the model. The sum of predation and fishing mortality rate provided in outputs of the model. Relative EI = EIs , where EI s and EI ref EI ref indicate the ecological indicators predicted in the scenario s and reference state, respectively; Relative Bio = Bios , where Bio s Bio ref and Bio ref indicate the species biomass predicted in the scenario s and reference state, respectively; Relative M = Ms -M ref , where M s and M ref indicate the mortality rates predicted in the scenario s and reference state, respectively.	Travers et al. (2006) Blanchard et al. (2010) Shin et al. (2005)

Table 3

 3 Results of Pearson correlation analysis in S2.

	Type	HTL group	Species	Correlation coefficient (r)	P value (P)
	Resident species Migratory species	SP0 SP1 SP6 SP8 SP11 Sp13 SP2* SP3* SP4* SP5* SP7* SP9* SP10*	O. oratoria Small shrimp (P. gravieri, P. Tenella, and A. japonicus) S. schlegelii A. hexanema J. belangerii C. joyneri C. bimaculata C. japonica Loligo sp. Octopus sp. P. fangi T. kammalensis L. tanakae	-0.901 -0.928 -0.643 -0.362 -0.579 -0.667 0.806 0.953 0.962 0.785 0.841 0.839 0.426	<0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01

Table 1 .

 1 References for speciesspecific migration were listed in Table B.2 of Xing et al. (2017).

	The Bio flux Table A.3 of of each species was estimated from the model calibration. Their values were listed in Xing et al. (2017).

  .1 of[START_REF] Xing | An individualbased model for simulating the ecosystem dynamics of Jiaozhou Bay, China[END_REF] and set based on relevant studies and stomach content analysis of samples collected in 2011. Plank.access represents the proportion of available biomass of LTL groups to HTL groups. Plank. access values were estimated from the model calibration and listed in Table A.2 of Xing et al. (2017). Growth When individuals ingest enough food for, they will grow and incremental length is calculated by von Bertalanffy model ( denote the fish age and mean increase of body length of species s at age a, respectively. L∞ s , Ks, and t0 s are parameters in the von Bertalanffy model. The individual weight at a given size is calculated by weight-length relationship (Shin and Cury, 2004): W

	s )
	) e -Ks (a-t0
	Shin and Cury, 2004): Ls,a = Δ L∞ s ( 1 -where a and Δ L s,a = csL where W s,a,i,t ' s,a,i,t ' bs and L e -Ks Δt s,a,i,t '

s,a,i,t ' represent the body weight and length of school i for species s at age ɑ and at time step t', respectively. The c s and b s are species-specific allometry parameters.

Acknowledgements

This work was funded by Marine S & T Fund of Shandong Province for Pilot National Laboratory for Marine Science and Technology (Qingdao) (No. 2018SDKJ0501-2) and National Key R&D Program of China (2018YFD0900904, 2018YFD0900906). The computational simulation and data analysis were primarily done at Dr. Chen's lab at the University of Maine and supported by the China Scholarship Council, Ocean University of China, and University of Maine. We are extremely grateful to members of Dr. Ren's lab at the Ocean University of China for generously providing data. The authors are grateful to Dr. Punt (Editor) and two anonymous peer reviewers for helpful comments, suggestions, and edits.

As a result of overfishing and environmental change, there are marked variations in natural resources in the water off the south of Shandong Province. A bottom trawl survey was conducted in each summer from 2010 to 2017 to investigate stock status [START_REF] Lv | Assessment of Important Fishery Resources in the South Offshore of Shandong from 2010 to 2017[END_REF]. The survey area is shown in Fig. B1. Trachypenaeus curvirostris is a low-trophic-level shrimp off the coast of Shandong Province. The T. curvirostris population had markedly changed in the last ten years (Fig. B2). The T. curvirostris biomass was highest in 2010 (760.5 kilotons) and lowest in 2011 (28.7 kilotons).


F Forget
email: fabienforget@gmail.com

L Dagorn

B Mérigot

J C Gaertner

J Robinson

PD Cowley

Y Rilwan

M Koonjul

V Mangar

M Taquet

F Ménard

Beta diversity of pelagic assemblages at fish aggregating devices in the open ocean

Keywords: FADs, Indian Ocean, Maldives, Mauritius, pelagic fish diversity, Seychelles, species composition, underwater visual census

Owing to difficulties in accessing the vast open ocean, the beta (β) diversity of pelagic fish assemblages remains poorly studied. We investigated the relationship between assemblage similarity and geographical distance between anchored fish aggregating devices (FADs), sampled by standardised underwater visual censuses in three anchored FAD arrays in the Indian Ocean-at the Maldives, the Seychelles and Mauritius. The use of two complementary indices of β-diversity, based on presence/absence data (Jaccard similarity coefficient) and abundance data (Bray-Curtis index), revealed that geographical distance between sampling sites (from 4 to 257 km) appeared to have no effect on the similarity of fish assemblages associated with FADs within each array. The results of this preliminary study question the generalisation of the paradigm of an increase in β-diversity with geographic distance to the open-ocean fish community. Large-scale studies using a variety of datasets should be conducted to further investigate patterns of β-diversity in the open ocean.

Characterising how diversity changes according to different spatial and temporal scales is a fundamental goal in community ecology (e.g. [START_REF] Legendre | Spatial pattern and ecological analysis[END_REF]. However, the study of patterns in fish diversity in the open ocean, the largest ecosystem on the earth [START_REF] Angel | Biodiversity of the pelagic ocean[END_REF], is still in its infancy. This situation is largely due to the logistical and financial constraints associated with collecting relevant data for the analysis of biodiversity in such a difficult-to-access environment [START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF]). As a consequence, the first attempts to infer worldwide patterns of species richness were restricted to both large predatory fishes and fisheries-dependent datasets [START_REF] Worm | Predator diversity hotspots in the blue ocean[END_REF][START_REF] Trebilco | Mapping species richness and human impact drivers to inform global pelagic conservation prioritisation[END_REF]. Although fisheries data remain an important data source, they have several shortcomings, and fisheries-independent methods are needed to complement the general picture of the diversity in pelagic ecosystems [START_REF] Murphy | Observational methods used in marine spatial monitoring of fishes and associated habitats: a review[END_REF]. In the last decade, several non-destructive and reproducible fisheries-independent methods have been developed to study pelagic fish assemblages, such as midwater baited camera surveys [START_REF] Letessier | Assessing pelagic fish populations: the application of demersal video techniques to the mid-water environment[END_REF][START_REF] Santana-Garcon | Development and validation of a mid-water baited stereo-video technique for investigating pelagic fish assemblages[END_REF] or underwater visual censuses (UVCs) at fish aggregating devices (FADs) [START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF]. Anchored or drifting FADs are artificial floating devices built by fishers to facilitate the capture of targeted pelagic species, usually tunas [START_REF] Hunter | Association of fishes with flotsam in the offshore waters of Central America[END_REF][START_REF] Fonteneau | A worldwide review of purse seine fisheries on FADs/Revue mondiale des pêcheries à la senne sous DCP[END_REF]. Several studies have investigated pelagic fish assemblages at FADs [START_REF] Castro | Fish associated with fish aggregation devices off the Canary Islands (Central-East Atlantic)[END_REF][START_REF] Dempster | Temporal variability of pelagic fish assemblages around fish aggregation devices: biological and physical influences[END_REF][START_REF] Addis | Spatial and temporal changes in the assemblage structure of fishes associated to fish aggregation devices in the Western Mediterranean[END_REF][START_REF] Andaloro | Annual trend of fish assemblages associated with FADs in the southern Tyrrhenian Sea[END_REF][START_REF] Taquet | Characterizing fish communities associated with drifting fish aggregating devices (FADs) in the western Indian Ocean using underwater visual surveys[END_REF]. [START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF] advocated that sampling at FADs allows the study of a well-delimited and reproducible part of the pelagic fish community, which then can help the monitoring of open-ocean fish diversity.

Beta (β) diversity involves investigations of species changes in space and/or time, as well as ecological connectivity [START_REF] Anderson | Navigating the multiple meanings of β-diversity: a roadmap for the practicing ecologist[END_REF]). As such it allows the identification of heterogeneous communities and is thus regarded as a valuable tool for conservation and Introduction ecosystem-based fisheries management [START_REF] Anderson | Beta diversity of demersal fish assemblages in the north-eastern Pacific: interactions of latitude and depth[END_REF]. Generally, β-diversity has been found to increase with respect to geographical distance in both terrestrial and marine environments [START_REF] Soininen | A multivariate analysis of beta diversity across organisms and environments[END_REF][START_REF] Shackell | Spatial scale of similarity as an indicator of metacommunity stability in exploited marine systems[END_REF]. We tested this hypothesis for open-ocean pelagic fish diversity. We assessed changes in β-diversity between pairs of FADs along a spatial gradient within anchored FAD arrays in three areas of the Indian Ocean, considered as three different case studies.

Materials and methods

In each of the anchored FAD arrays (at the Maldives, the Seychelles and Mauritius) (Appendix 1), UVCs were conducted by SCUBA divers on four specific anchored FADs which were sampled on the same day. These UVCs were repeated (conditions permitting) five times in the Maldives, four times in the Seychelles, and three times in Mauritius, between May 2008 and January 2011 (Table 1). All the surveyed FADs were anchored at depths greater than 300 m and were thus considered to be offshore deep-water FADs. The FAD design within each country was identical but could differ between countries. In the Maldives, the floating structure of the FAD consisted of a single large fibreglass buoy (1.5 m in diameter), whereas in Mauritius and the Seychelles a chain of rigid buoys (0.2 m in diameter) was used. The dives were conducted between 08:00 and 15:00, at 15 m depth and for a 30-minute duration, following a standardised UVC protocol [START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF]. The horizontal visibility during the UVCs varied between 15 and 20 m, and thus we considered that the volumes sampled during all the UVCs were comparable. Owing to safety concerns, the UVCs were performed only in calm to moderate sea conditions (i.e. wind speed <30 km h -1 , wave height <1.5 m, and surface current speed ≤1.9 km h -1 ). Tuna species were not included in the censuses as they often occurred outside the visual range of the scientific divers, making estimations of these species particularly inaccurate [START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF].

Within each FAD array, we investigated the effects of the geographical distance between each pair of FADs i and j on the composition of their fish assemblage. Here, low values of similarity indices indicated high β-diversity. We used two complementary similarity indices: Jaccard's coefficient of similarity (S ij ) and the Bray-Curtis similarity index (1-BC ij ); S ij is based on presence/absence data and takes into account the proportion of unshared species between two FADs, whereas 1-BC ij is based on abundance data [START_REF] Anderson | Navigating the multiple meanings of β-diversity: a roadmap for the practicing ecologist[END_REF]. A log transformation, log(n+1), was applied to the abundance data n to minimise the effect of highly abundant species (e.g. schooling fish), prior to calculating the Bray-Curtis index. The pairwise comparisons were computed independently within each FAD array and on the same sampling dates only. This was done to exclude possible seasonal effects as reported in other studies [START_REF] Dempster | Temporal variability of pelagic fish assemblages around fish aggregation devices: biological and physical influences[END_REF]. Furthermore, linear mixed-effects models were used to test the effect of the geographical distance on each index (Y ij ) for each FAD array. Each pair of UVCs was classified according to the date on which they were carried out. The sampling-date effect (µo i ) was thus treated as random variations around a population mean, and the distance between anchored FADs (β 1 dist ij ) was assessed as a fixed continuous covariate, as follows:

Y ij = β 0 + β 1 dist ij × µo i + ε ij

Results

A total of 19 species (12 families), 14 species (9 families) and 21 species (10 families) were observed in the Maldives, the Seychelles and Mauritius, respectively (Figure 1). Carangidae (such as Elagatis bipinnulata, Caranx sexfasciatus and Decapterus macarellus) were the most frequent recorded species in each FAD array (Figure 1a). Figure 1b displays the distribution of fish species richness for each FAD array. The Maldives had a higher mean species richness per UVC (mean 7.3 [SD 2.4]; Kruskal-Wallis test: χ 2 = 20.8, df = 2, p < 0.05) than both Mauritius (mean 4.7 [SD 1.7]) and the Seychelles (mean 3.8 [SD 1.2]). β-diversity expressed through both the S ij and 1-BC ij similarity indices generally exhibited no trend with an increase in geographical distance between FADs in each FAD array (i.e. the geographical distance between FADs was generally not a significant fixed-effect in the mixed models) (Figure 2; Appendix 2). The only significant test result occurred with 1-BC ij (p = 0.04) for the Maldives, and is clearly explained by index values computed for the two distances in this array (39 and 158 km). Finally, both S ij and 1-BC ij similarity values indicate that there is an overall low similarity (mean S ij 0.33 [SD 0.12]; mean 1-BC ij 0.33 [SD 0.14]) of fish assemblages between two pairs of FADs (Figure 2).

Discussion

Standardised UVC observations allowed us to investigate fish assemblage patterns among a delimited component of the epipelagic diversity, and species richness measured here is comparable to that of previous studies measuring pelagic fish diversity using baited remote underwater video systems (BRUVS) [START_REF] Bouchet | Drifting baited stereo-videography: a novel sampling tool for surveying pelagic wildlife in offshore marine reserves[END_REF][START_REF] Rees | Rules of attraction: enticing pelagic fish to mid-water remote underwater video systems (RUVS)[END_REF][START_REF] Santana-Garcon | Pelagic stereo-BRUVs: development and implementation of a fishery-independent technique to study pelagic fish assemblages[END_REF], UVCs [START_REF] Addis | Spatial and temporal changes in the assemblage structure of fishes associated to fish aggregation devices in the Western Mediterranean[END_REF][START_REF] Gaertner | Visual censuses around drifting fish aggregating devices (FADs): a new approach for assessing the diversity of fish in open-ocean waters[END_REF], and fisheries data at FADs [START_REF] Lezama-Ochoa | Biodiversity in the by-catch communities of the pelagic ecosystem in the western Indian Ocean[END_REF][START_REF] Ruiz | Bycatch of the European, and assocated flag, purse-seine tuna fishery in the Indian Ocean for the period 2008-2017[END_REF]). Among the wide variety of available methods to assess β-diversity, [START_REF] Anderson | Navigating the multiple meanings of β-diversity: a roadmap for the practicing ecologist[END_REF] outlined two different conceptual types of β-diversity: (i) a non-directional variation in species' identities or community structure among sample units within a given area or region at a given spatial (or temporal) scale; and (ii) a directional turnover in community structure or species' identities along a spatial, temporal or environmental gradient. We used indices from the second category as we investigated the effect of geographical distance between each pair of FADs on the fish assemblage composition. Surprisingly, in each FAD array, the geographical distance between sampling sites (from 4 to 257 km) generally had no effect on the similarity of assemblages for both indices. Patterns in β-diversity are sensitive to variation in the number of species measured at alpha (α) diversity between two assemblages when sampling effort differs [START_REF] Chase | Using null models to disentangle variation in community dissimilarity from variation in α-diversity[END_REF]. Thus, comparing assemblages when this effort differs automatically induces changes in β-diversity. In such a case, the β-diversity Raup-Crick index, based on presence/absence data, accounts for the effect of differential α-diversity on measuring β-diversity [START_REF] Chase | Using null models to disentangle variation in community dissimilarity from variation in α-diversity[END_REF]). However, sampling effort was standardised in our study, and the Raup-Crick values (not shown) confirmed results from the Jaccard indices.

Our findings contrast with the general pattern in β-diversity usually observed for many communities in marine ecosystems, in which the similarity in species assemblage decreases with distance between sites; this is known as the distance-decay relationship [START_REF] Soininen | A multivariate analysis of beta diversity across organisms and environments[END_REF][START_REF] Shackell | Spatial scale of similarity as an indicator of metacommunity stability in exploited marine systems[END_REF]. Generally, environmental factors play a central role in this relationship that is observed in other ecosystems. However, few studies have investigated the influence of environmental factors at small-scale (~1-20 km) or meso-scale (~20-200 km) on pelagic fish. Sub-mesoscale frontal dynamics (1-10 km) in the Mozambique Channel were observed to shape the distribution of marine top predators (Tew [START_REF] Kai | Top marine predators track Lagrangian coherent structures[END_REF], and in the Humboldt Current System fine-scale ocean dynamics (1-4 km) were found to affect the seascape from zooplankton to predators [START_REF] Bertrand | Broad impacts of fine-scale dynamics on seascape structure from zooplankton to seabirds[END_REF]. However, no studies have investigated the effect of oceanographic conditions on the composition of species assemblages at FADs.

In conclusion, our results question the current perception of fish diversity in pelagic ecosystems, suggesting heterogeneity in assemblages among FADs even at close distances and in a seemingly featureless environment. Pelagic fishes are considered to be highly mobile organisms, with the standard order of magnitude of the horizontal scale in the distribution of communities considered to be large (~10 3 km: [START_REF] Angel | Biodiversity of the pelagic ocean[END_REF]. These features are peculiar to the pelagic realm, and could explain why the distance-decay relationship, which is observed in most marine and terrestrial ecosystems, was observed here in a pelagic ecosystem. Equipping FADs with oceanographic probes would allow the collection of in situ data and the investigation of whether or not environmental changes can occur at small spatial scales.

Additionally, investigation of the behavioural dynamics driving multispecies aggregations could help to interpret our results. One hypothesis is that fish assemblages from FADs that are close to each other (e.g. a few kilometres apart) could have interactions that could shape their composition. For instance, complex behavioural processes within species (e.g. social behaviour) or between species (e.g. agonistic behaviour) could affect species occurrence and abundance at such FADs. Hence, there is a need to collect additional simultaneous in situ data (particularly behavioural data) on fish assemblages to develop a better understanding of the processes that forge pelagic diversity patterns in the open ocean at spatial scales over which such processes take effect. Finally, a meta-analysis using a variety of available datasets should be conducted to further investigate patterns of β-diversity in the open ocean.

Figure 1 :

 1 Figure 1: (a) Overall frequency of occurrence and (b) frequency of species richness observed per underwater visual census at anchored fish aggregating devices in the three study areas (Maldives, Seychelles and Mauritius) between May 2008 and January 2011. S max = maximum species richness

Figure 2 :

 2 Figure 2: Jaccard (S ij ) and Bray-Curtis (1-BC ij ) similarity indices compared at different distances between the anchored fish aggregating devices within each study area (Maldives, Seychelles and Mauritius). The p-values refer to the linear mixed-effects models, and n denotes the number of pairwise comparisons on the same sampling day

Table 1 :

 1 Summary of the sampling by underwater visual census (UVC) conducted at anchored fish aggregating devices (FADs) in three countries in the western Indian Ocean, between May 2008

	and January 2011
	Area	Sampling date (no. UVCs)	Total no. of UVCs	No. of FAD pairs compared	Distance between FADs (km)
	Maldives	Nov 2009 (4)	19	27	39-257
		Feb 2010 (3)
		Apr 2010 (4)
		Oct 2010 (4)
		Jan 2011 (4)
	Seychelles	Jan 2010 (4)	15	21	4-34
		Apr 2010 (4)
		Oct 2010 (4)
		Dec 2010 (3)
	Mauritius	May 2008 (3)	11	15	11-31
		Oct 2010 (4)
		Dec 2010(4)

Acknowledgements -This study was supported by the BIOPS project funded by the 'Fondation pour la Recherche sur la Biodiversité' (No. CD-AOOI-07-013). We would like to thank the Marine Research Centre (MRC) in the Maldives, the Albion Research Institute (ARI) in Mauritius, and the Seychelles Fishing Authority (SFA) and the Desroches Island Resort for their help with collecting data for this study.

ORCID

Shiham Adam: https://orcid.org/0000-0001-7443-1901 Paul Cowley: https://orcid.org/0000-0003-1246-4390 Laurent Dagorn: https://orcid.org/0000-0003-3896-8589 Fabien Forget: https://orcid.org/0000-0002-4845-4277 Jean-Claude Gaertner: https://orcid.org/0000-0001-9248-6260 Frédéric Ménard: https://orcid.org/0000-0003-1162-660X Bastien Mérigot: https://orcid.org/0000-0001-5264-4324