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.