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.
Due to data availability, Monthly
SST data in our study can only obtained from two sources as follows:
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
where 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 length-frequency data with the Electronic Length
Frequency Analysis (ELEFAN) approach using the Tropfishr package
(Mildenberger et al. 2017). Size-frequency of the commercial fishing
catch during 1970 – 1982 are not
available but we gathered 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) coinciding
in time and space with the available data. We evaluated L. crocealife-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 (n=2074)
and weight (n=853) of the collected L. crocea individuals were
measured, length-weight relationship was estimated based on the equationw =aLb (n=853). 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-history 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 a
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 at different life 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 to
predict overwintering habitat suitability of L. crocea in the
ECS, which is a type of specie
distribution models (SDMs) used for evaluating organisms-habitat
relationships based on limited data or expert knowledge. 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 andL. 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 value adjusted for
small sample size (AICc) (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 AICcvalue (Supporting information), thus was selected as the final HSI
model. Correspondingly, fitting based SI model was selected as the best
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, average, 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.