2.2 Environmental data
Three environmental factors were collected for L. croceaoverwintering habitat suitability modelling: depth (m), sea surface temperature (oC) and sea surface salinity. The Bathymetry data (30 arc-second spatial resolution) was obtained from Gridded Bathymetry Chart of the Oceans (GEBCO) (https://www.gebco.net/data_and_products/gridded_bathymetry_data/, accessed on September 2020) to represent depth. Historical monthly SST and SSS data (0.5° spatial resolution) between 1971 and 2001 was acquired from the Simple Ocean Data Assimilation (SODA) reanalysis dataset (Carton and Giese 2008) (from the Climate Data Library:http://iridl.ldeo.columbia.edu/, accessed on September 2020). Monthly SST data (4 km spatial resolution) between 2002 and 2019 was downloaded from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the satellite Aqua platform provided by the Ocean Biology Processing Group at NASA Goddard Space Flight Center (https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/, accessed on September 2020). To investigate the environmental data from when fishing occurred, environmental data were resampled by the mean value of each month to 0.5° spatial resolution, then matched with the fisheries data.
2.3 Life-history parameters ofL. crocea in the ECS
To understand the life-history parameters of L. crocea in the study area, we analyzed the length-based data with electronic length frequency analysis (ELEFAN) using the TropfishR package (Mildenberger et al. 2017). Size-frequency of the commercial fishing catch during 1970 – 1982 are not available but we gather the life-history information from published literature (Yu and Lin 1980, Xu et al. 1984a, b, Liu and de Mitcheson 2008, Ye et al. 2012) consisted with the study period and study area. We evaluated L. crocea life-history parameters during 2018 – 2019 using body length data collected by otter trawl cruises, seasonally from November 2018 to November 2019 (e.g. November 18th, 2018; January 20th, 2019; April 18th, 2019; July 17th, 2019; September 28th, 2019; November 18th, 2019). we selected a total of 2074L. crocea individuals between 2018 to 2019: the length and weight of the collected L. crocea individuals were measured (n=853), length-weight relationship was estimated based on the equationw =aLb . The length frequency was calculated for each season for Electronic Length Frequency Analysis (ELEFAN) in TropfishR package.
We fit the von Bertalanffy growth function (VBGF) through the length frequency and life-history data to estimate life parameters (e.g. maximum length, weight at length, length at maturity (Lmat ) and VGBF parameters including von Bertalanffy growth constant (K ) and asymptotic length (Linf ), age at zero length(t0 ) and et al.)(Pauly and David 1981, Brey and Pauly 1986, Sparre and Venema 1998) and estimated mortality parameters (e.g. total mortality (Z ), natural mortality (M ) and fishing mortality (F )) (See details in Supporting Information). Specifically, the VBGF and mortality parameters were estimated following four approaches, including: (i) K-Scan for the estimation of K for a fixed value of Linf(ELEFAN K.S.); (ii) ELEFAN Response Surface Analysis (ELEFAN R.S.A); (iii) ELEFAN with simulated annealing (ELEFAN S.A.); (iv) ELEFAN with a genetic algorithm (ELEFAN G.A.). Among above four workflows, ELEFAN R.S.A, ELEFAN S.A and ELEFAN G.A. allow the simultaneous estimation ofK and Linf (Taylor and Mildenberger 2017), while a advantage of ELEFAN S.A. and ELEFAN G.A. is the possibility to estimate all parameters of the seasonalized VBGF simultaneous (Scrucca 2013, Xiang et al. 2013).Therefore, we estimated life-history parameters using the 60 scenarios with different bin of length (bin=10 and bin=20), move average (MA) value (MA=5, MA=7, MA=9, MA=11) and four different workflows (ELEFAN K.S., ELEFAN R.S.A., ELEFAN S.A., ELEFAN G.A.) following Mildenberger et al. (2017).
For model validation and selection, we used all-subsets model selection based on the fraction of the estimated sum of peaks (Rn) following Gayanilo et al.(1997). Additionally, other ratios of life-history parameters, such as M /K andF /M , were calculated with the estimated parameters. Moreover, Caddy et al. (1998) pointed out that the trophic level of a certain fish would be changing with size. Hence, we used size and trophic levels relationship to estimate the size truncation effect on overall population trophic level between 1980s and recent study (Supporting information).
2.4L. crocea overwintering distribution patterns and overwintering habitat suitability in the ECS
L. crocea is overall habitat specialist during overwintering phase, previous studies have shown that L. crocea has strong depth, temperature, and salinity preferences, while pH, dissolved oxygen, light intensity, sound, water velocity and other factors may affect its distribution pattern, survival and growth in different life-history stages (Liu 2013, Wang et al. 2016, 2017). Unfortunately, detailed data on seasonal environmental data was limited between 1971 and 1982. Hence, only depth, SST and SSS data were available as model inputs to determine suitable habitats for L. crocea . We used HSI model, which is a kind of well-known SDMs (Species Distribution Models) to predict overwintering habitat suitability of L. crocea in the ECS. In our HSI model, we used standardized abundance, HSI, as response variable, and three environmental variables with strongest correlation and the best data availability, depth, SST and SSS, as predictor. Firstly, we constructed both fitting-based (Lee et al. 2019, Hua et al. 2020, Yu et al. 2020) and regression-based (Chang et al. 2019, Jin et al. 2020) suitability index (SI) models to describe the relationship between each environmental variable and L. crocea abundance (Supporting information). Then, we combined two types of SI models into HSI models, respectively. For each type of SI model, we used two empirical HSI models: the arithmetic mean model and the geometric mean model (Supporting information Fig. S1), under different environmental variable combinations (Lee et al. 2019).
For model validation and selection, we used catch data (abundance) from 1971 to 1980 and corresponding environmental data as training data set, catch data from 1981 to 1982 and corresponding environmental data as testing dataset. We assumed a positive linear relationship between predicted HSI values and L. crocea abundance and the evaluated the goodness-of-fit of the above relationship for each HSI model based on \(R^{2}\) and the Akaike Information Criterion (AIC) value (Chang et al. 2013). Fitting-based arithmetic mean model with two variables (e.g. depth and SST) yielded the maximum \(R^{2}\) and the minimum AIC value (Supporting information), thus was selected as the final HSI model. Correspondingly, fitting based SI model was selected as the final SI model. Finally, we retrained the SI model (See parameters and statistical test results in Supporting information) and the HSI model using catch data (abundance) from 1971 to 1982 and corresponding environmental data.
We drew SI curves using the retrained final SI models. We then computed winter (December, January, and February next year) mean SST and used them to predict yearly winter mean HSI distribution between 1971 and 2019 with the retrained final HSI model, then calculated decadal winter mean HSI of the 1970s (1971 – 1980), 1980s (1981 – 1990), 1990s (1991 – 2000), 2000s (2001 – 2010) and 2010s (2011 – 2019). Regions with HSI value > 0.7, 0.7 > HSI value > 0.3, and HSI value < 0.3 were regarded as optimal, common, and poor habitat, respectively. We computed the area of different habitat types in each year and analyzed if there are significant difference in area between decades using Person One-way ANOVA and Scheffe test.