Constructing age-structured matrix population models for all fishes

Data

To construct an age-structured matrix population model, data from three sources were combined. Here, the term “data” refers to the statistics available from databases and publications rather than the original raw data. The first source was FishBase (Froese & Pauly, 2024), which is an online database of fishes. It includes life history and other data for 35,600 species. The amount of data varies substantially among species and depends on the number of researchers who contribute to the database. The data from FishBase were directly downloaded using the rfishbase package (v 2.0; Boettiger, Lang & Wainwright, 2012). The second source was the FishLife R package (Thorson, 2024). The package takes the life history data from FishBase and predicts missing life history parameters for other species based on their taxonomic relatedness (Thorson et al., 2017). A newer version of the FishLife package also uses data from the database of stock-recruitment relationships developed by Myers, Bridson & Barrowman (1995) to predict density-dependent parameters (Thorson, 2019). Finally, the third source of data was the general mass-dependent mortality of fish derived by Lorenzen (1996). This data was used instead of the mortality estimates from FishBase, which often only provides mortality estimates independent of size or age (typically a point estimate for adults). Table 1 summarizes the life history data used in constructing the population matrix.

Construction of age-structured matrix population models

The data were assembled into age-structured matrix population models as follows. More details for each step are described later in this section.

1. The size of the population matrix was determined based on the maximum age.

2. The age starting with 1 was converted into length using the von Bertalanffy growth equation.

3. The length was converted into mass using length-mass allometric relationships.

4. The mass was converted into the instantaneous mortality rate.

5. The instantaneous mortality rate was converted into the age-specific finite survival rate, which was entered into the subdiagonal elements of the age-strcutured matrix population model (see further details below).

6. For ages greater than the age of maturation, the fertility rate was calculated assuming it is proportional to the mass of individuals for the age and scaled such that the asymptotic growth rate is 1. (Density-Independent Matrix Population Model)

7. The density-dependent parameters were calculated from the maximum annual spawner biomass per spawner biomass in excess of replacement, age of maturity, age-specific individual mass, and age-specific survival rate. (see Density-Dependent Matrix Population Models)

Determining the size of the population matrix

The size of a population matrix was determined by the maximum age (tmax), obtained from the FishLife package, rounded up to the nearest integer. Hereafter, this integer value is denoted as t_max. This assumes that the time unit of the model is one year. Therefore, only cases in which the maximum age is greater than 1 were considered. In other words, annual populations are not addressed in this paper because they lead to unstructured population models.

Age to length conversion

Age t was converted into length L (t) using the von-Bertalanffy equation: (1)${\displaystyle L\left(t\right)=L_{\infty }\left(1-e^{\left(-\kappa \left(t-t_0\right)\right)}\right)}$

where L_∞ is the asymptotic length, κ is the parameter determining how fast individuals grow in size, and t₀ is an extrapolated age at which the length would have been 0. L_∞ and κ were obtained from the FishLife package. t₀ was calculated by solving Eq. (1) for t₀ and substituting the age of maturity and the length of maturity (both also obtained from the FishLife Package) along with L_∞ and κ. When the obtained value for t₀ was greater than −0.1, −0.1 was used instead; it is the default value used in Thorson (2019). The unit of length in the package was in cm, and it assumes the total length (rather than standard length or fork length).

Length to mass conversion

The age-specific mass W (t) was used in determining the instantaneous mortality rate and was converted from the length L (t) of individuals using the following allometric relationship: (2)${\displaystyle W=a{L}^b}$

where a and b are the allometric coefficients obtained from FishBase. Note that Eq. (2) omits the age to reduce clutter, but both length and mass are a function of age. Here, the unit of mass is in grams.

The length in the length-mass relationship may be measured in terms of total length, fork length, or standard length. If the total length was used, the allometric parameters were directly used in Eq. (2). When fork length or standard length was used, the total length L (t) was converted into the corresponding length measure using the length-length relationship obtained from FishBase and substituted into Eq. (2). The length-length relationship is expressed as a linear equation. For example, the total length L_T is converted into standard length L_S as: (3)${\displaystyle L_S=\alpha +\beta L_T}$

where α and β are the intercept and slope coefficients, respectively, for the length-length relationships. Note that the original database uses symbols a and b, but they were changed into α and β to avoid using the same symbols for different quantities in this paper. The total length can be converted into the fork length L_F in the same manner. Then, L_S or L_F was substituted into Eq. (2) instead of L_T when the allometric relationship was for standard length or fork length, respectively.

Natural Mortality

The natural mortality was obtained by using the Lorenzen mass-dependent mortality rate: (4)${\displaystyle m=-\mu W^{\upsilon }}$

where m is the instantaneous mortality rate. The two coefficients µand υ are those obtained from Lorenzen (1996). The estimated values are µ= 3.00 and υ =−0.288. In the original paper, these two coefficients were denoted by M_u and b, respectively, but to reduce confusion with other parameters, the symbols were changed in this paper. The unit of the instantaneous mortality rate is per year (year⁻¹).

Age-specific finite survival rate

The age-specific natural mortality rate m was converted to the age-specific finite survival rate s_t as follows: (5)${\displaystyle s_t=\frac{exp\left({\int }_1^{t+1}m\left(\tau \right)d\tau \right)}{exp\left({\int }_1^{t}m\left(\tau \right)d\tau \right)}.}$

The finite survival rate is the proportion of individuals that survive from age t to t +1, and it is calculated for t =1, …, t_max-1, where t_max is the size of the population matrix. The finite survival rate was entered into the <t +1, t>elements of the population matrix.

Fertility rate

The fertility rate for individuals of age t (F_t) was assigned to a population matrix based on the following assumptions:

1. Individuals start contributing to reproduction once they reach the age of maturity.

2. The first age class is for age 1.

3. The fertility rate for age t is proportional to the mass at age t.

4. The finite population growth rate (λ) is 1.

Therefore, the <1,t>element of the population matrix is given by F_t =cW (t) for t greater than or equal to the age of maturity (t_mat). Here, c is a proportionality constant, and the same value is applied to all age-classes of the same population. To obtain the age of maturity (t_mat), the value denoted by tm in the FishLife package was rounded up to the nearest integer. Finally, the value for c was numerically searched such that the dominant eigenvalue (λ) of the population matrix was 1.

Finally, 0 was entered into empty elements of the population matrix to complete the density-independent population matrix.

Density-dependent population matrix

The density-dependent population matrix was developed by modifying the density-independent population matrix based on the following assumptions:

1. All density-dependent processes take place between ages 0 and 1.

2. Density dependence operates through the total biomass.

3. The age of maturity is the age of recruitment.

4. Density dependence takes the form of compensatory density dependence,

5. The total number of mature fish at the equilibrium point is 10,000 individuals per unit area (an arbitrary value, as the habitat area for a population is often unknown).

Therefore, the density-dependent fertility rate F_t,x for age t at time x is given as: (6)${\displaystyle F_{t,x}=\frac{D_a{W}_t}{1+D_b{B}_x}}$

where D_a and D_b are density-dependent constants and (7)${\displaystyle B_x=\sum _{t=t_{mat}}^{t_{max}}{n}_{t,x}{W}_t}$

where n_t,x is the density of individuals of age t at time x.

The coefficient D_a was determined from the maximum annual spawner biomass per spawner biomass in excess of replacement (MASPS) in the FishLife package as follows: (8)${\displaystyle D_a=\frac{p}{M_{spawner}{l}_{mat}}}$

where p is the maximum lifetime spawner biomass per spawner biomass (Myers, 2001) and given as, (9)${\displaystyle p=1+\frac{\gamma }{1-\overline{S}}}$

where $\overline{S}$ is the mean survival rate after maturation and γ is the MASPS (see Eq. (7) in Thorson, 2019). The survivorship from age 1 to maturation is denoted by l_mat and is calculated as ${\prod }_{t=1}^{t_{mat}-1}{s}_t$. The survivorship is included in the denominator of Eq. (8) because the MASPS calculation assumes that density dependence occurs somewhere between reproduction and recruitment, whereas the model in this paper assumes that the density dependence occurs somewhere between reproduction and age 1. Therefore, the age-structured population matrix in the current study includes survival from age 1 to the age of maturity separately. The spawner biomass per recruit is denoted by M_spawner, and it is calculated as ${\sum }_{t=t_{mat}}^{t_{max}}\left(\frac{l_t}{l_{mat}}{W}_t\right)$ where l_t is the survivorship from age 1 to age t. Finally, the value of D_b is determined numerically, by setting ${\sum }_{t=t_{mat}}^{t_{max}}{n}_t=10,000$, which is arbitrarily determined as it depends on the habitat in which a population lives (i.e., rarely known parameter).

Matrix population model analysis

Once the density-independent and density-dependent population matrices were obtained, various analyses were performed. Table 2 shows the seven quantities calculated in this study.

An age-structured population matrix (often denoted by A) is called a Leslie matrix (Leslie, 1977), and it takes the following form: (10)${\displaystyle \mathrm{A}=\left[\begin{array}{ccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill F_t\hfill & \hfill F_{t_{max}}\hfill \\ \hfill s_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill s_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \ddots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill s_t\hfill & \hfill 0\hfill \end{array}\right]}$

where the subdiaconal elements are the age-specific survival rates, and the first row is for the fertility rates (the product of fecundity and survival of offspring until age 1; see Kendall et al., 2019). The fertility rate is assumed to be 0 before individuals mature. In other words, the <1,t>element is 0 if t < t_mat.

Once the density-independent population matrix A is obtained, the asymptotic population growth rate λ₁ is given by the dominant eigenvalue of the matrix. The associated right eigenvector gives the stable age distribution (w), which is often scaled to sum to 1. The associated left eigenvector gives the reproductive value (v), which is scaled in various ways.

The sensitivity matrix (S) is given by the product of the reproductive value (v) and stable age distribution (w) as (11)${\displaystyle \mathrm{S}=\frac{v{w}^T}{v^{T}w}}$

where both v and w are assumed to be column vectors. The <i,j>element (e_ij) of elasticity matrix (E) is given by $e_{ij}=\frac{a_{ij}{\varsigma }_{ij}}{{\lambda }_1}$ where a_ij is the <i,j>element of the population matrix and ς_ij is the <i,j>element of the sensitivity matrix.

The damping ratio measures how quickly transient dynamics dissipate and the population returns to asymptotic growth/decline after a perturbation (Caswell, 2001). It is given by the ratio of the dominant eigenvalue and the second eigenvalue as (12)${\displaystyle \rho =\frac{{\lambda }_1}{\left|{\lambda }_2\right|}.}$

Generation time has several different definitions. Here, it is defined as the expected mean age of the parents of newborn individuals in a cohort. It is calculated as (13)${\displaystyle G=\frac{{\lambda }_1{v}^{T}w}{v^T\mathrm{F}w}}$

where F is a matrix consisting of only the fertility rates (Bienvenu & Legendre, 2015).

Finally, the resilience R is given by the dominant eigenvalue of the linearized population matrix λ_J (Caswell & Neubert, 2005) as (14)${\displaystyle R=-ln\left(\left|{\lambda }_J\right|\right).}$

The linearized population matrix J is given as (15)${\displaystyle \mathbf{J}=\mathbf{A}{|}_{n^{\ast }}+\left[{\left.\frac{\partial \mathbf{A}}{\partial n_1}\right|}_{{\mathit{n}}^{\ast }}{\mathit{n}}^{\ast },{\left.\frac{\partial \mathbf{A}}{\partial n_2}\right|}_{{\mathit{n}}^{\ast }}{\mathit{n}}^{\ast },\dots {\left.\frac{\partial \mathbf{A}}{\partial n_k}\right|}_{{\mathit{n}}^{\ast }}{\mathit{n}}^{\ast }\right]}$

where n^∗ is a population density vector at the equilibrium and n_i is the density of age class i.

Quantification of uncertainties and computation

To quantify the uncertainties associated with the metrics, the data obtained from the FishLife package were simulated using the mean and variance from the package, assuming a multivariate normal distribution. Under each simulation, age-structured matrix population models (both density-independent and -independent models) were constructed, and all metrics were calculated. For the analysis presented in the paper, the data for each species were simulated 1,000 times.

All analyses were performed using R Statistical Software (v4.3.1; R Core Team, 2024) with the tidyverse package (v2.0.0; Wickham et al., 2019) and the MASS package (v 7.3; Venables & Ripley, 2002). The code used for the case study can be found in Fujiwara (2024).

Writing

After the draft of the manuscript was written, it was uploaded to ChatGPT-4o (Open AI, 2024) with memory off. The text was revised using the following prompt:

“I am writing a scientific paper to be published in a peer-reviewed journal in the field of ecology. The following is the draft of the manuscript. Please revise it for grammatical errors. Please also make the tense of verbs consistent with scientific writing practices. Keep the common names of species in singular form and lowercase without any article”.

The manuscript was revised one section/subsection at a time. Then, the abstract was written using the following prompt:

“I am writing a scientific paper to be published in a peer-reviewed journal in the field of ecology. The following is the draft of the manuscript. Please write an abstract of the paper in approximately 300 words”.

After all these processes, the manuscript was further edited manually for accuracy.