Life histories determine divergent population trends for fishes under climate warming
Compilation of population life-history data
We compiled population-specific life-history data for Indo-Pacific fishes from primary literature using a combination of systematic and opportunistic searches. We obtained references primarily from Google Scholar and the web page of Fisheries and Oceans Canada (DFO, http://www.isdm-gdsi.gc.ca/csas-sccs/applications/publications/index-eng.asp), using a combination of species name and the following keywords: Pacific, temperature effects, von Bertalanffy growth, age at maturation, or natural mortality. Species names for 107 common teleosts in the Indo-Pacific Ocean were obtained from the global-capture production database of the Food and Agricultural Organization (FAO, ISSCAAP groups 31-37; available at: http://www.fao.org/fishery/statistics/global-capture-production/query/en). Although similar data are available in FishBase, a global and freely available database on fish, we opted not to use data from FishBase because temperature data associated with population records are often missing and difficult to restore due to lacking spatial coordinates in FishBase.
In total, we screened around 8000 references published in 1958–2017. We extracted data from 440 references representing peer-reviewed journal papers, stock assessment reports, governmental documents, graduate theses, and conference proceedings. Our data include 332 species, with conspecific life-history data of ≥2 records in 206 species.
We extracted available population data of seven life-history traits related to growth and maturation for each population: i.e., the von Bertalanffy growth coefficient K (yr-1) and asymptotic length L∞ (cm), age and length at 50% maturity (A50 (yrs) and L50 (cm)), allometric exponent (b) of length–weight relationship (W = aL^b), lifespan (Amax (yrs)), and natural mortality (yr-1). To make the length traits (L∞ and L50) comparable among populations, we standardized different length measures (e.g., standard length or fork length) into total length using available regression models (Supplementary Table 7). To account for between-sex differences in traits values, we coded life-history data of each population as sex-specific or combined-sexes. When data of single- and combined-sexes were both available, we primarily used combined-sexes data for our analysis. When only sex-specific data were available, we used female data.
In our selected references, the von Bertalanffy growth coefficients are primarily estimated via fitting length-at-age data (Lt, where L is length and t is age) with one of these methods: the Ford–Walford method, length-based method, and nonlinear regression fitting. As all these fitting methods involve the same model form (L_t = L_infty left( {1 – e^{ – Kleft( {t – t_0} right)}} right)) and similar optimization processes (e.g., minimizing residual errors), we consider that the growth coefficients estimated with these methods are generally comparable. However, we excluded the growth coefficient data derived from other model forms, as transforming those coefficients could introduce biases (e.g., Gompertz model or weight-based form of von Bertalanffy model). The population natural mortality data were derived from various models with distinct assumptions (Supplementary Table 8); thus, these natural mortality estimates may not be comparable. To investigate the temperature effects on natural mortality, we re-estimated natural mortality estimates for individual populations using a life-history invariant method (see Methods: Derivation of M).
We recorded additional descriptive variables for each population: e.g., ecological group (coded based on habitat types following definitions in FishBase, Supplementary Table 5), family, species, latitude, and longitude of sampling site (in decimal degrees; for the studies with multiple sampling sites, the minimum and maximum of latitudes and longitudes of sampling sites are recorded), minimum and maximum sampling depths, and sample size. Latitudes and longitudes of sampling sites were derived from available spatial coordinates, names of sampling ports, or the maps in the references.
Some references were excluded to avoid ambiguity. For example, we excluded references of fish in estuary and artificial habitats (e.g., aquaculture species). To alleviate noise in the datasets, we rejected data lacking clear descriptions of study area (e.g., studies that did not define sampling sites or those that aggregated samples across a very broad range) and those lacking explicit information on fitting procedures for the growth or maturation models (e.g., lacking information on sample size, suitable range of length-at-age data, independent samples, description of the fitting methods, a plot of data with the fitted line, or other means for assessing the fit). Similarly, we excluded studies without descriptions of fitting methods for estimating natural mortality (as suggested by 14). For stock assessment reports and review papers, we scrutinized the summarized data based on original publications and rejected data lacking clear sampling or fitting information. Lastly, we removed duplicated data cited in multiple reports.
Compilation of population temperature data
We used available in situ data of mean decadal sea temperature as habitat data. These temperature data were obtained from the World Ocean Atlas 2013 (NOAA Atlas NESDIS 73; 43; hereafter referred to as WOA13; available at: https://www.nodc.noaa.gov/cgi-bin/OC5/woa13). Sources of the WOA13 temperature data include temperature profiles measured by various instruments43. The mean decadal temperature profile was calculated by averaging six decadal datasets spanning 1955–201243. For our analysis, we used the mean decadal temperature profile data in the Indo-Pacific region, with a horizontal resolution of 0.25° and depth segments of 5 m from surface to 100 m and depth segments of 25−100 m for >100 m (maximum depth = 5500 m)43. As rates of temporal changes in the ocean temperature were slow (e.g., 1.48 °C per century; 25), we did not account for the small temperature changes due to differences in time of measurements between the temperature and life-history data.
We matched the mean decadal temperature profile data with spatial coordinates for each of these populations. For populations with single sampling sites, we assumed that their potential habitats centered at the sampling sites and extended for 0.5° in the north, south, east, and west directions. Further, for populations with multiple sampling sites, we assumed that their habitats were bounded by the ranges of latitudes and longitudes of sampling sites. To account for habitat depths, we estimated the minimum and maximum depths of populations based on ranges of sampling depths or description of depths of species habitats in the references. When depth information was unavailable, we used the maximum depth of the species in FishBase as the maximum depth for a population. For each population, we derived two habitat variables: SST and BT. SST is the temperature at 0 m (WOA13 data; 43). BT is the temperature at the maximum depths of populations. Finally, we calculated the minimum, mean, maximum, and coefficients of variation of each of the habitat indices for each population.
Derivation of natural mortality (M)
Population-specific natural mortality (hereafter referred as M) provides insight into population sustainability under fishing or environmental changes (i.e., high natural mortality indicates high degree of sustainability13,15). However, because estimates of M from different models are not necessarily comparable, we re-estimated M for each population using the model II of Gislason et al.14, which builds on well-established empirical relationships44,45 and has some theoretical backing46,47. The equation is:
$$ln left( M right) = 0.55 – 1.61ln left( L right) + 1.44ln left( {L_infty } right) + ln left( K right),$$
(1)
where M is natural mortality and L is the midpoint of length range (cm) of a population, respectively14. As the length range data were not available to us, we constrained L to be the length at age-at-50% maturity (A50) to estimate M, following the invariant relationship among mortality-at-age at maturation, L50, L∞ and K13,48. Further, because A50 data are available for fewer number of populations (n = 119, about 8.5% of total populations), we restored the missing values of A50 based on a theoretical linear relationship between population A50 and (frac{{L_{50}}}{{K times L_infty }})49. With the available data of A50, L50, L∞, and K for 70 populations, we fitted this linear relationship (F = 204.6, df = 1,68, R2 = 0.75, P 25 yr-1). Range of M’s for the remaining 1387 populations is 0.05–21.8 yr-1 (Supplementary Fig. 1). These Gislason M estimates were used in the subsequent regression and life-table modeling analysis.
Evaluation of the relationships between temperature and life-history traits
Because life-history variation is nested within phylogenetic levels, we used the linear mixed-effect model (LME) to explore the relationships between each life-history trait and each temperature index, simultaneously accounting for species- and family-related variance as random effects. We used the lme4 and lmetest packages in R (www.r-project.org50;) to construct these LME models. Given the eight different temperature indices (e.g., four descriptive metrics nested within two temperature variables), we constructed eight LME models for each life-history trait. For each model, a natural-log transformed life-history trait is the response variable, a species mean-centered habitat index is a single fixed-effect variable, and species and family are two random-effect variables. We considered four alternative model structures for random effects: i.e., with either family or species as random intercepts, and with either family or species as random intercepts and slopes.
We evaluated strength of each of fixed- and random-effect variables using the likelihood ratio test51. Specifically, we estimated the log likelihood between a pair of full and alternative models (with one less fixed or random variables), deriving the test statistics (i.e., 2 times the difference between log likelihood of two models) and P value. If the two models were not different significantly (e.g., P > 0.05), we selected the more parsimonious model. Otherwise, significant difference between these two models (e.g., P ≤ 0.05) indicates pronounced variation in the response variable due to the additional fixed or random effect in the full model. Furthermore, when neither family nor species accounts for significant variance in the response variable, we reduced the model to a simple linear regression with a habitat index as the sole fixed-effect predictor. We selected the best model structure to depict the relationship between each of the eight habitat indices and a life-history trait.
We observed pronounced intra-specific variability in the empirical K and L∞ for some of our study species; e.g., ≥3-fold differences in these traits within some species. Because large variability in K and L∞ within a species may be partially due to that the estimation of these traits depends on one another, we conducted PCA with these two traits and evaluated the ratios of maximum-to-minimum scores of the first PC1 among species. We found 34 species with the absolute values of the ratios of PC1 scores ≥3, whereas 285 species had |ratios| 0 indicate positive warming effects on population growth rates). To derive (frac{{R_0}}{G}), we built an age-structured model (i.e., a life-table33;), incorporating functions for growth increment (in length and weight), maturity states, fecundity, and survivorship probability. The length-at-age data (Lt) were calculated using the von Bertalanffy growth model (Eq. (3)):
$$L_t = L_infty left( {1 – e^{ – Kleft( {t – t_0} right)}} right),$$
(3)
where t denotes age, and L∞, K, and t0 are the von Bertalanffy growth coefficients. Weight-at-age data (Wt) were estimated using an allometric function of Lt (Eq. (4)):
$$W_t = alpha L_t^beta.$$
(4)
Also, fecundity (mt) is an allometric function of Wt (Eq. (5)23:
$$m_t = delta W_t^gamma.$$
(5)
Because of lacking data on the length-weight and fecundity-weight relationships for many populations, we assumed constant intercepts and slopes for Eqs. (4) and (5) (i.e., α = 0.02, β = 3, γ = 1.18, δ = 2930). However, L∞ and K data are available for most populations in our data (n = 1,387). As a result, we used these data to derive length-at-age for each population (assuming t0 = 0). Also, we used the model-derived A50 estimates (Supplementary Fig. 7) to account for maturity state for mt for each population (i.e., mt = 0 for t More
