Global reconstruction of life‐history strategies: A case study using tunas

Abstract Measuring the demographic parameters of exploited populations is central to predicting their vulnerability and extinction risk. However, current rates of population decline and species loss greatly outpace our ability to empirically monitor all populations that are potentially threatened. The scale of this problem cannot be addressed through additional data collection alone, and therefore it is a common practice to conduct population assessments based on surrogate data collected from similar species. However, this approach introduces biases and imprecisions that are difficult to quantify. Recent developments in hierarchical modelling have enabled missing values to be reconstructed based on the correlations between available life‐history data, linking similar species based on phylogeny and environmental conditions. However, these methods cannot resolve life‐history variability among populations or species that are closely placed spatially or taxonomically. Here, theoretically motivated constraints that align with life‐history theory offer a new avenue for addressing this problem. We describe a Bayesian hierarchical approach that combines fragmented, multispecies and multi‐population data with established life‐history theory, in order to objectively determine similarity between populations based on trait correlations (life‐history trade‐offs) obtained from model fitting. We reconstruct 59 unobserved life‐history parameters for 23 populations of tuna that sustain some of the world's most valuable fisheries. Testing by cross‐validation across different scenarios indicated that life‐histories were accurately reconstructed when information was available for other populations of the same species. The reconstruction of several traits was also accurate for species represented by a single population, although credible intervals increased dramatically. Synthesis and applications. The described Bayesian hierarchical method provides access to life‐history traits that are difficult to measure directly and reconstructs missing life‐history information useful for assessing populations and species that are directly or indirectly affected by human exploitation of natural resources. The method is particularly useful for examining populations that are spatially or taxonomically similar, and the reconstructed life‐history strategies described for the principal market tunas have immediate application to the world‐wide management of these fisheries.


| INTRODUC TI ON
Monitoring a population's status and risk of extinction requires detailed knowledge of life-history characteristics, including survival, growth and reproduction (Pearson et al., 2014;Salguero-Gómez et al., 2016;Winemiller, 2005). However, our understanding of these traits is fragmented, error-prone and incomplete for many species (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000). For example, more than 10% of the 87,967 animals, plants and fungi assessed by the International Union for Conservation of Nature (IUCN) are classified as Data Deficient, such that their status and vulnerability cannot be determined (IUCN, 2017). Consequently, identifying reliable ways to detect declines and assess the extinction risk of data-limited populations has been highlighted as a key ecological problem facing policy makers (Kindsvater et al., 2018;Sutherland et al., 2006).
One approach of assessing data-limited populations involves judiciously filling-in missing values of life-history parameters with surrogate data from similar species or other populations of the same species. However, using surrogate parameters concurrently with observed life-history traits neglects the trade-offs, that is, budgetary compromises, that connect the different aspects of a population's demography, constraining the range of possible life-history strategies that can evolve (Stearns, 1992). The biases and imprecisions generated by employing this approach are difficult to quantify, creating a considerable margin for parameter misuse and population mismanagement. The process of using surrogate life-history values to overcome missing data has become increasingly formalized with the widespread availability of meta-analyses for key parameters (e.g. Denney, Jennings, & Reynolds, 2002;Hutchings, Myers, Garcia, Lucifora, & Kuparinen, 2012;Thorson, Taylor, Stewart, & Punt, 2014). Predictive approaches have also been developed to estimate missing values across large groups of species based on correlations between available life-history data (Freckleton & Jetz, 2009;Jetz & Freckleton, 2015;Thorson, Munch, Cope, & Gao, 2017). These studies predict life-history parameters by grouping species based on phylogeny and fine-scale environmental determinants of life-history traits, such as temperature.
Commercial fisheries exploit populations at rates that greatly outpace our ability to accurately assess them. Populations with missing life-history data are particularly problematic in this context because, by definition, the data necessary to develop evidence-based management and policy action are incomplete. At present, less than 20% of the global fish catch originates from populations that are managed by formal population assessment (Costello, 2012), and therefore, threats from over-exploitation are rarely identified before species have suffered large population declines (Burgess, Polasky, & Tilman, 2013). Correlative approaches have been used to estimate specific life-history traits for data-limited fish species (Beverton, 1992;Hordyk, Ono, Sainsbury, Loneragan, & Prince, 2014;Nadon & Ault, 2016), incorporating phylogenetic and environmental determinants of traits to access large numbers of species (Thorson et al., 2017).
However, phylogeny becomes less informative when considering population-level strategies, as opposed to species-level strategies, and when examining species within the same genus. Furthermore, populations that inhabit large geographic areas and migrate over broad latitudinal gradients preclude a straightforward assignment of fine-scale environmental determinants for life-history strategy.
Predictive approaches informed by theoretically motivated constraints that align with life-history theory (Kindsvater et al., 2018) offer a new avenue for providing informative distinctions between populations that are closely placed taxonomically and spatially.  (Costello, 2012;FAO, 2014). However, these populations have also experienced long-term declines in adult biomass since the 1950s (Juan-Jordá et al., 2011), such that 43% are currently overfished with biomass levels below management missing life-history information useful for assessing populations and species that are directly or indirectly affected by human exploitation of natural resources. The method is particularly useful for examining populations that are spatially or taxonomically similar, and the reconstructed life-history strategies described for the principal market tunas have immediate application to the world-wide management of these fisheries.

K E Y W O R D S
Bayesian imputation, data limited, demography, fecundity, life-history theory, missing data, principal market tuna, Scombridae recommendations (ISSF, 2018). Furthermore, population-specific life-history data are missing for more than a third (36%) of maturity and fecundity parameters (Figure 1; see Supporting Information Table S1). For example, age of maturity and annual fecundity are poorly resolved, and are only reported for 10 and two of the 23 populations respectively. Consequently, population assessments often rely on age of maturity estimates from other populations, and spawning stock biomass is generally used as an index of fecundity when fitting stock recruitment models (e.g. ISC, 2016, Rice, Harley, Davies, & Hampton, 2014.
Here, we develop a Bayesian hierarchical approach that considers a multi-population, multispecies dataset and reconstructs a full life-history strategy for the 23 populations of principal market tuna that are recognized globally. We employ a multivariate approach that determines the direction and strength of life-history trade-offs among traits in order to resolve the similarity between species. This framework provides access to traits that are difficult to measure directly, such as age of maturity and annual fecundity. Furthermore, by considering multiple sources of information and including error in the observed data, it also generates a measure of uncertainty associated with each parameter. The model outlined here is applied to tunas; however, the approach could be applied to a wide range of taxa and we explore the predictive capacity of the method through a series of crossvalidation experiments.

| Life-history traits for the principal market tunas
We considered a seven trait life-history strategy for the principal market tunas. Data on six of these traits were collated from Juan-Jordá, Mosqueira, Freire, Ferrer-Jordá, and Dulvy (2016). As the principal market tunas are sexually dimorphic, we removed trait estimates based on males only. Traits included adult survival (Φ) derived from maximum age (T max ; Equation 1); age of maturity (A, age that 50% of the sampled individuals are mature); spawning frequency (S, number of days between each batch of eggs, oocytes, being produced); spawning duration (D, proportion of the year that the population is actively spawning); absolute batch fecundity (F, mean number of oocytes produced per individual per batch); and somatic growth rate (K) from the von Bertalanffy growth function.
Three additional studies on batch fecundity published after the original dataset was compiled were also included in order to bolster populations that were missing data across the four fecundity traits (Supporting Information Table S1). To resolve a full set of fe-

Southern Bluefin
Pacific Bluefin Atlantic Bluefin S1). Finally, as broad similarities in life-history strategy can be drawn between species of tuna that inhabit the same biome (Juan-Jordá, Mosqueira, Freire, & Dulvy, 2013, each population was assigned a binary habitat variable (H) delineating tropical or temperate preference (Figure 1). To promote convergence across life-history parameters with different scales (Congdon, 2003), each trait was standardized to a mean of 0 and a standard deviation of 1.
Rates of adult survival (Ф) were based on the following estimation of natural mortality (M) from maximum age (T max ), where 4.3 is an estimated constant (Kenchington, 2014): Mortality estimates that are based on maximum age reportedly perform better than those based on other parameters, such as body size . Using a maximum age estimator also provided the largest sample size for rates of adult survival from the starting dataset, and allowed the estimation of this trait to be independent from metrics of body growth that would be used simultaneously to impute mortality in later analysis. Natural mortality is more commonly included as a demographic trait in fisheries assessment, compared to survival. However, we include the con-  Table S1), and few studies provided information on more than one or two life-history traits.
Thus, examination of trade-offs at the individual-study level was not possible.

| Model of life-history traits
The Bayesian hierarchical model comprised a demographic component with parameters describing species-specific relationships between life-history traits, and an observation component that integrated observed data for each population with trait-specific uncertainty associated with measurement error and process variability in the data (e.g. spatial and temporal variation). The code and data file for the model are provided in Appendices S1 and S2 of the Supporting Information, and key model parameters are listed in Supporting Information Table S2. All models were implemented in JAGS (v. 4.3.0) via the "jagsUI" library (v 1.4.9) for program r (v. 3.4.1). Models were fitted by running three Monte Carlo Markov chains (MCMC) for 1.5 × 10 7 iterations and retaining every 200th step in order to increase the effective MCMC sample size for the same amount of computer memory. The first 5,000 MCMC draws were removed as burn-in, and each chain was initialized at different points in the parameter space. Convergence of the chains was confirmed using the Brooks-Gelman-Rubin diagnostic tool (all values r ≤ 1.01) and the effective sample size for each parameter (we required this to be >500).
The demographic component of the model (Equation 2) contained a function for each life-history trait. Each trait was modelled in relation to a binary habitat term (H) delineating tropical and temperate species, and two random-effect terms (γ and ω) that represent the correlated residuals (i.e. trade-offs, Charnov, 1993) between traits at the species-and population-level, respectively: Here, subscripts i and j denote population and species respectively. Lower case symbols (ϕ, a, s, d, f, v, k) represent simulated life-history traits that are modelled on the linear scale; that is, logit for survival and spawning duration, log for all other variables. The trait-specific coefficients for the effect of habitat Consequently, all methods and results relate to the model that included growth rate.
To capture the trade-offs between the seven life-history traits at the between-and within-species level, we modelled two random-effect terms (γ and ω respectively) using multivariate normal distributions. The covariance structures (Σ i and Σ j : Equation 3) parameterizing these distributions allow the correlations, or lifehistory trade-offs, that connect the different aspects of a population's demography to emerge during model fitting: The mean of the multivariate normal distributions (μ, Equation 3) had normal priors. Standardizing the input data to a mean of zero and a standard deviation of one meant that all traits were able to receive priors centred on zero. For the species-level term (γ), being centred on zero reflects the standardized trait mean. For the population-level term (ω), being centred on zero allows the population mean to operate as a classic random effect. The standard deviation in the species-level priors was set to 0.25 (i.e. a precision of 15) to allow flexibility, and this was reduced to a standard deviation of 0.14 for the population-level (i.e. a precision of 50). The variance-covariance matrices in the multivariate normal distributions (Σ i and Σ j , Equation 3) were assigned inverse Wishart prior distributions with an identity matrix for the scale matrix (Ω) and eight degrees of freedom (df): This equation is repeated in the model, such that c is either species (i) or population (j). Setting the df parameter for the Wishart prior to one more than the dimensions of the variance-covariance matrix (n = number of traits) achieves a uniform prior distribution on the individual correlation parameters; that is, an equally likely probability between −1 and 1 (Gelman & Hill, 2006). A scaling parameter (ξ) assigned from a uniform prior: U(0,5), was also incorporated to overcome constraints on the scale parameters associated with an inverse-Wishart model (Gelman & Hill, 2006).
Observed life-history data for each trait (Φ, A, S, D, F, V, K) were incorporated through an observation component that reflected a combination of measurement error and process variability in the data (e.g. spatial and temporal variation). Upper case notation is used to denote observed, population-specific data. Data were included at the study level (p) for each species (i) and population (j) : Variance was assumed to be Gaussian, where the recon- Equation 2, and the precision (τ ϕ , τ a , τ s , τ d , τ f , τ v , τ k ) was set from a normal distribution centred on 20: N(20,0.1). This assumes a reasonably high precision on the reconstructed value for data-limited traits; however, for data-rich traits, this prior distribution will allow the precision to travel into the tails, that is, become high or low in line with the observed variance in the data. Thus, in addition to imputing missing values, posterior distributions for the observed data points were also obtained. This acknowledges that even observed life-history data are a sample from a theoretical population of possible trait-specific values whose mean and credible intervals can be obtained. Although this exercise is similar to the traditional approach of sampling distributions, the joined inference across all data prevents the resulting parameter posteriors from being restricted by small sample sizes.

| Model validation
We  Table S1). For the second scenario, we removed the life-history data for a species represented by a single population, retaining somatic growth rate only; here, we chose Southern bluefin tuna (Figure 1; Supporting Information Table   S1). We also ran an additional six model  Figure 2B-j).
The reconstructed life-history strategies (Supporting Information   Table S3) of the principal market tunas appeared along a thermalgrowth gradient, whereby tropical species had faster rates of somatic growth compared to species that spend a large portion of their annual cycle in temperate waters (Figure 2). In addition, there was a gradual increase in survival rates from the fast-growing tropical species to the slower growing temperate species ( Figure 3A). The (4) tropical species also matured earlier compared to the temperate species and returned predominantly larger credible intervals for annual fecundity ( Figure 3B and F). The relationship with maturation was further demonstrated by the coefficient of the relationship with biome (97.5 credible interval for ψ a , Equation 2: −2.000 to −0.003: Figure 3B). Based on the median values, spawning duration was longer in tropical species ( Figure 3D). In contrast, spawning frequency and batch fecundity did not clearly differ between the two habitat groupings ( Figure 3C and E), although batch fecundity increased and was more variable for the three slow-growing bluefin species ( Figure 3E).
Large intra-population variation in the raw data for spawning duration was not reduced through model fitting (see Supporting Information Figure S1). Furthermore, the large disparity between the mean published values of maturation for the two populations of Atlantic bluefin tuna was maintained ( Figure 3B, Supporting Information Figure S1). Intra-species variability appeared as different median estimates of traits within species, for example, in survival ( Figure 3A) and somatic growth ( Figure 3G), as well as differences in the credible interval for each trait, for example, in age of maturity ( Figure 3B), spawning duration ( Figure 3D), spawning frequency ( Figure 3C) and annual fecundity ( Figure 3F; Supporting Information Table S3). Within the three bluefin tuna species, the breeding strategy of Pacific bluefin tuna appeared more distinct, compared to the other two species (Figure 3).

| Validation testing
Life-history traits were recreated convincingly for a popula-  Figure 4B). However, maturity was imputed to be lower than the original data and the median value reconstructed by a model that incorporated all available data. Furthermore, the credible intervals of reconstructed traits increased considerably, compared to the model including all available data (Supporting Information Table   S4). Inference was only improved by refitting the model with the respective life-history data incorporated (Supporting Information Figure S2).
Examination of the coefficient parameters for the effect of habitat preference demonstrated that age of maturity was the only trait with a step difference between biomes, as opposed to a more continual change from the fastest growing tropical species to the slowest growing temperate species. Removing the remaining habitat terms from the demographic functions generated marginal changes to the posterior credible intervals for annual fecundity and the median values for spawning duration. Furthermore, the remaining traits were largely unchanged (Supporting Information Figure S3).
However, the overall fit of a model without a complete set of habitat terms was reduced compared to a model including these terms (minimum effective sample for model without habitat terms = 2,604, all r ≤ 1.13; minimum effective sample size for model including habitat for the majority of populations, max. 7%; Supporting Information Table S5). The imputed median life-history traits were also highly correlated with those estimated by the model that included these additional four metrics of fecundity (Supporting Information Figure   S4).

| D ISCUSS I ON
In this study, we applied a formal hierarchical Bayesian approach For this species, the median posterior values were relatively accurate for traits where the original data appeared close to the interspecies linear regression with growth rate (e.g. survival and annual fecundity). However, the observed values of maturation and batch fecundity are both higher than expected for its growth rate, that is, these values have large residuals from the linear regression, thus reducing predictive capacity ( Figure 3B and E). Furthermore, Southern bluefin tuna exists at the edge of the parameter space where the likelihood surface may be less well defined. This is likely to have con- Southern bluefin contributing factor in their overexploitation (Juan-Jordá et al., 2015). where λ and β are parameters. The former is a measure of the maximum per capita productivity, that is, related to annual fecundity, and the latter a measure of the strength of density dependence. Steepness (h) is defined to be the reproduction when spawning biomass is 20% of its unfished level (B 0 ) relative to reproduction at the unfished level; h = R(0.2B 0 ) R(B 0 ) (see Mangel et al., 2010Mangel et al., , 2013. For the BH-SRR, steepness is (Mangel et al., 2010(Mangel et al., , 2013: where W is the average spawning biomass of a female that has survived to age a (Φ(a)), and the probability of being mature at age a (p m (a)) and mass at age a (W(a)) according to W = ∑ a (a)p m (a)W(a). Our methodology provides values for survival and maturation and can easily be adapted to include mass as an additional life-history trait.
Application of Equation 6 to determine steepness requires knowing λ, which is the product of mass-specific fecundity δ (eggs, or larvae for those species with live birth, per unit spawning biomass) and survival from the egg or larval stage until recruitment to the population ϕ EL , thus, λ = δϕ EL . Survival to recruitment into the population can be modelled using well-established allometries between size and mortality (e.g. Brodziak, Mangel, & Sun, 2015;Kai & Fujinami, 2018;Mangel, Brodziak, & DiNardo, 2010;Simon et al., 2012), and the described hierarchical model of life-history traits allows the computation of total fecundity. In the Supporting Information Appendix S4, we show how to obtain mass-specific fecundity from total fecundity under appropriate assumptions. Thus, once mass-specific fecundity is known, the parameter λ can be determined for a data-poor stock.
Our analyses reveal that by considering multiple life-history traits, populations and species simultaneously in a Bayesian framework, it is possible to statistically maximize the utility of sparse datasets, quantify the trade-offs that connect different aspects of an organism's life-history and access traits that are difficult to measure empirically. The percentage of missing lifehistory data will inevitably be higher for other, commercially less important groups of fishes, such that a Bayesian hierarchical approach informed by the principles of life-history evolution, including metabolic theory, could further assist in population-specific trait reconstruction.