The emergence of longevous populations

Data and Analyses for Human Populations.

The data for the hunter-gatherer populations were obtained from published data on two populations, the Hadza and Ache (15, 16) (see Table 3 for summary information). Life tables for preindustrial 18th-century Sweden (1751–1759), modern Sweden (2000–2009), and modern Japan (2012) were drawn from the “Human Mortality Database” (HMD) (19). For these three datasets we extracted death counts by age (0–109 y), calendar year, and sex and corresponding exposure (Table 3). An additional 15 life tables for females and 15 life tables for males were obtained for Fig. 4A from the HMD for Sweden in different periods, Ukraine (1933), and Iceland (1882); from the “World Health Organization” (WHO) (54) for Nigeria, India, Russia, China, and the United States all for the year 2013; and from published sources for Liberian migrants (1820–1843) (17, 18), Trinidad slaves (1813–1816) (55), and England (1600–1725) (56). Additionally, we obtained Siler mortality parameters estimated by ref. 57 for acculturated hunter-gatherers for both sexes combined. For Fig. 4B, we obtained 8,198 yearly period life tables from the HMD for 44 countries.

For the six main human populations we carried out a bootstrap analysis to estimate all our pace–shape measures. For most life tables except Liberia, we assume that the observed number of deaths was the outcome of a random Poisson process in a population with given exposures N(x) and the given age-specific death rate, μ(x), as the “true” signal. For each jth bootstrap step, we randomly drew a new set of number of deaths for age x, D_j(x), from a Poisson distribution with parameter λ_x = μ(x) N(x). From the D_j(x)s we estimated the resulting life table, using standard methods (43) and the subsequent pace and shape measures. We ran 20,000 bootstrap steps, which allowed us to calculate mean, SEs, and confidence intervals for each measure. In the case of Liberia, we assumed that the observed number of deaths was the outcome of a random binomial process with parameters q(x) for the probability of death and the number of individuals entering the interval [x, x + Δx], namely P(x). We then drew D_j(x) values from a binomial distribution with parameters q(x) and P(x). The next steps were equivalent to those for the Poisson bootstrap.

The values of life expectancy and lifespan equality plotted for other populations in Fig. 4 A and B were calculated from data from various sources, as indicated in Table S4. In some cases data were available for 5-y age categories: We estimated values for a single year of age by linear interpolation. When life tables were available for single years of time but we needed tables for multiyear periods of time, we then took simple averages over the period of the single-year values of life expectancy and lifespan equality. We applied these same methods to all other life tables in Fig. 4, except for the acculturated hunter-gatherers for which we used the estimated Siler mortality parameters provided by Gurven and Kaplan (57).

Smoothing of Life-Table Data for Fig. 1.

To display the densities in Fig. 1 we smoothed the death rates for the hunter-gatherer populations using P-splines, because of small sample sizes in these populations, following the proposal of Eilers and Marx (58). We used the implementation by Camarda (59, 60) in R, which has been tailored for smoothing mortality in a Poisson framework. We also smoothed the data for historic Sweden (1751–1759). The reason was not insufficient data, as was the case for the hunter-gatherer population: Our dataset for historic Sweden included ∼8 million person years for each sex when pooling data across years (Table 3). Instead, the data for Sweden from that period suffer from strong age heaping (19, 61). We suppressed those artificial fluctuations over age by having a strong penalty term λ for the P-spline smoothing. Death rates at ages below 80 y for contemporary Sweden (2000–2009) and Japan (2012) were not modified or smoothed in any way. Mortality at ages 80 y and above for historic (1751–1799) and contemporary (2000–2009) Sweden as well as for the population of contemporary Japan (2012) was estimated using a logistic mortality model given byμ¯(x)=aebx1+aγb(ebx−1),[1]

which expresses the population hazard µ(x) as a mixture of Gompertz-distributed hazards with Gamma-distributed “frailty” among individuals (62). We used the same parametric model to estimate mortality for historic Sweden, contemporary Sweden, and Japan at ages 80 y and higher.

We estimated the parameters in a standard Poisson maximum-likelihood framework, given byln ℒ(α,β,γ|D(x)N(x))=∑x=80ωD(x)lnμ¯(x)−N(x)μ¯(x),[2]

where D(x) is the number of people dying at age x and corresponding exposure times N(x); and µ(x) refers again to the population-level hazard of Eq. 1.

We then estimated life tables and life expectancy, using standard methods (43). We obtained the mean survival time at age x of those who die at age x, typically denoted as a(x), from the HMD for contemporary Japan as well as for historic and contemporary Sweden. Based on life table information for the number of survivors at age x, l(x), and the number of person years lived L(x) available for the Hadza, we used values of a(0) = 0.3, a(1) = 0.4, and a(x) = 0.5 for x > 1 for all hunter-gatherer populations.

Data for Nonhuman Primates.

We obtained data for six species of nonhuman primates from the Primate Life History Database (PLHD) (11), which includes longitudinal life-history data for known individuals in six study populations (Table 1). All populations in the PLHD are living in the wild and, with few exceptions, no provisioning or interventions have occurred in any of these populations (11). For all of these primate species one or both sexes undergo natal dispersal (first dispersal after birth), whereas four (baboons, gorillas, sifakas, and capuchin monkeys) have secondary or higher-order dispersals (Table S5).

For each of the six primate populations, the study population was defined as the set of social groups in which individuals were continuously monitored for life-history events. With one exception (muriquis), all study populations were embedded in larger continuous populations that were not isolated from the social groups being studied; consequently, individuals could immigrate into and emigrate from the study population, resulting in a situation in which some individuals were not observed throughout their entire lives. In the case of the muriquis, research expanded to encompass all four social groups in the entire isolated study population in the early 2000s, about 20 y into the research project. Thus, for muriquis subsequent to 2003 individuals could leave the study population but there was no immigration into the study population. All told, the database included nine different types of individual records in two broad categories (Fig. S2): (i) natal individuals, i.e., individuals that were live born in the study population, and (ii) immigrants, i.e., individuals born in social groups that were not part of the study population before their appearance in a study group, who emigrated from their natal social group and immigrated into a study group. For all six studies there were a considerable number of individuals with unknown fate, for both natals and immigrants, because an individual that belongs to the dispersing sex, has reached dispersal age, and disappears could have either died or attempted to disperse into a social group outside the study population.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S2.

Circles represent times of entry; solid circles are known times of birth and open circles are entries after birth. Squares are departure times where solid squares are known times of death and open squares are out-migration. Solid triangles indicate individuals known to be alive until the end of the study and vertical bars indicate unknown type of departure.

Two further aspects of the data required estimation procedures. First, most immigrants had estimated birthdates (i.e., their ages were estimated at immigration). Individuals also had estimated birthdates if they were present in the study population at the time monitoring of the study population first began and they were first individually identified. Second, in most populations, some infants died before their sex was ascertained, resulting in deaths of infants with unknown sex.

Because of these sources of uncertainty, we constructed a Bayesian model to estimate mortality in these populations. Our model included estimation procedures for age and sex in cases where these were not known exactly and incorporated the probability of emigration for both natals and immigrants.

Bayesian Model for Nonhuman Primate Data with Sex-Specific Dispersal.

We used an extension of the Bayesian approach proposed by Barthold et al. (13) to model age- and sex-specific mortality for species where one or both sexes undergo natal dispersal. The model builds upon the Bayesian survival trajectory analysis framework (63, 64) to model sex-specific mortality. We extended the model to include species for which one or both sexes undergo higher-order dispersal as well as the estimation of unknown times of birth (for a full model description see refs. 13 and 14).

This approach allowed us to estimate the essential latent states of sex, age, and dispersal state and hence the parameters for mortality and out-migration from each study population (hereafter termed out-migration). We defined the random variables X for ages at death, Y for ages at natal out-migration, and Z for ages at immigrant out-migration. We defined a dispersal state D that assigns 1 if an individual i out-migrates in its last detection age, xiL=tiL−bi, and 0 if otherwise. We treated D as a latent variable for all individuals with unknown fate. We also defined a variable S that assigns 1 if an individual is female and 0 otherwise.

The mortality function or hazard rate isμ(x)=limΔx→0Pr(x<X<x+Δx|X>x)Δx.

We assumed that the mortality function in all of the primate species analyzed was well described by the Siler hazards rate (65) given byμ(x|θ)=exp[α0−α1x]+κ+exp[β0+β1x], x≥0,[3]

where θ⊤=[α0,α1,κ,β0,β1] is a vector of mortality parameters to be estimated, with α₀, β₀ ∈(−∞, ∞) and α₁, κ, β₁ > 0. From the hazard rate in Eq. 3 we define the survival functionl(x|θ)=Pr(X>x) =  exp[−∫0xμ(t|θ)dt]=exp[eα0α1(e−α1x−1)−κx+eβ0β1(1−eβ1x)],[4a]

the probability density function (PDF) of ages at deathd(x|θ)=μ(x|θ)l(x|θ), x>0,[4b]

and the cumulative distribution function (CDF) for ages at death given by F(x|θ) = 1 − l(x|θ).

For out-migration, we assumed that the age at natal out-migration was Y∼GY(y|γ) for ages y > 0, where G_Y(y) is the Gamma distribution function with parameter vector γ = [γ₁, γ₂]. This distribution yields the PDF of ages at natal out-migration given bygY(y|γ)={γ1γ2Γ(γ2)(y−yd)γ2−1e−γ1 (y−yd)if y≥yd0if y<yd,[5]

where y_d is the minimum age at natal dispersal and γ₁, γ₂ > 0.

For immigrants in species with higher-order dispersal we assumed that the age at immigrant out-migration was Z∼GZ(z|λ) for ages z > 0, where λ = [λ₁, λ₂] is the vector of parameters to be estimated. The PDF of ages at immigrant out-migration is thengZ(z|λ)={λ1λ2Γ(λ2)(z−zd)λ2−1e−λ1 (z−zd)if z≥zd0if z<zd,[6]

with λ₁, λ₂ > 0 and where z_d is the minimum age at higher-order dispersal. Here we assumed that z_d = y_d.

The full Bayesian model is given byp(du,su,bu,θ,γ,λ|dk,sk,bk,tF,tL)∝p(dk,sk,bk,tF,tL|du,su,bu,θ,γ,λ)︸likelihood×p(d)p(s)p(x)︸priors for states×  p(θ)p(γ)p(λ)︸priors for parameters,[7]

where d is the vector of dispersal states (i.e., d_i = 1 if the individual out-migrated and 0 otherwise), s is the indicator vector for females (s_i = 1 if female and 0 if male), b is a vector of ages at birth, and t^F and t^L are the vectors of times at first and last detection, respectively, where x^F = t^F – b, and x^L = t^L – b. Each of these vectors has two subsets represented by the subscripts u for unknown and k for known.

To estimate the parameters and latent states we used a Markov chain Monte Carlo (MCMC) algorithm to fit the model in Eq. 7 that uses a Metropolis-within-Gibbs sampling framework (66, 67) (for details of the estimation, see refs. 13 and 14). In addition to the implementation in refs. 13 and 14, we extended the model to estimate unknown times of birth. The conditional posterior for unknown times of birth, b_i, isp(bi|di,xiL,xiF, θ, γ, λ)∝p(di,xiL,xiF, θ, γ, λ|bi) p(bi|biL,biU),[8]

where the first term in Eq. 8 corresponds to the likelihood function and the second term is a prior for times of birth. These priors are either uniform or normally distributed, both with upper and lower bounds provided with the dataset (i.e., b^L and b^U).

Estimation of Prior Parameters.

To estimate the prior parameters for the mortality and out-migration parameters, we used a combination of published data, expert information provided by the PLHD coauthors, and an agent-based model designed to simulate natal and immigrant out-migration. The model required the parameters for adult mortality provided in Bronikowski et al. (12). Because these were calculated only for adults, we fixed the prior parameters that control juvenile mortality in the first exponential term of the Siler model in Eq. 1 at α0∗=−2 and α1∗=0.3. We then used nonlinear least-squares estimation to find the remaining three parameters of the Siler model, using as reference the adult age-specific mortality constructed with the parameters found in ref. 12. For each species we calculated the sum of squares asQ(θp∗|θb)=∑i∈v[μ(xi|θp∗)−μg(xi−xm|θb)]2,[9]

where θp∗ was the vector of prior mortality parameters to be estimated, µ(x|θ^*) was the Siler hazards rate constructed with the prior parameters, µ_g(x – x_m | θ_b) = a exp[b (x − x_m)] with a, b > 0 was the Gompertz mortality function where θ_b = [a, b] was the vector of parameters estimated in ref. 12, and v⊤=[x1,x2, …,xT] was a sequence of T = 10 equally distanced ages starting at the age at maturity x₁ = x_m and with last age x_T, such that the Gompertz survival with parameters θ_b was S_g(x_T |θ_b) = 0.05. We used the R built-in function nlm to find the prior parameters θ_p that minimized the sum of squares in Eq. 9 (see Table S6 for mean values for prior parameters).

To estimate the priors for out-migration parameters γ and λ, we implemented an agent-based model. We used published information on dispersal behavior for each species as well as expert information provided by the PLHD researchers. Then, for each species we simulated a hypothetical population of dispersing individuals that could disperse between study groups, represented by set G, or to nonstudy groups, defined by set A, with universal set S = {G ∪ A} that represented all of the possible groups for each population. The number of groups in G and A varied among the six studies, and any movement between groups (either study or external groups) is considered to represent dispersal in the population.

To parameterize the density functions of ages at natal and higher-order dispersal we used information on minimum, maximum, and average ages at natal and higher-order dispersal provided in the literature and by the PLHD researchers. We calculated mean, y¯1 for natal, and y¯2 for higher-order dispersal and variances σ12 for natal and σ22 for higher-order dispersal and estimated the corresponding Gamma parameters asαj=y¯j2σj2 and βj=y¯jσj2,

for j = 1,2, where α_j and β_j are the shape and rate parameters for the natal and higher-order dispersal (see Table S7 for resulting parameters).

We used an individual-based model to simulate mortality and dispersal, to find the prior parameters for the distribution function of ages at natal and immigrant out-migration. The step-by-step algorithm for first dispersal is as follows: (i) For every individual i simulate age at death (x_i) from a Siler mortality function with parameters θ_p; (ii) simulate ages at natal dispersal (yi1) by randomly sampling from a Gamma distribution with parameters {α₁, β₁}; (iii) individuals where xi>yi1 successfully dispersed, and otherwise they died before dispersing for the first time; (iv) all dispersers move with the same probability to any of the groups within the study population (G) or to one of the areas outside the study population (A); (v) those individuals that moved to set A have out-migrated and their age at out-migration y_i is stored, and those that moved to set G have remained in the population; (vi) in species with higher-order dispersal, those in set G disperse again at a further age (yi2) by randomly sampling from a Gamma distribution with parameters {α₂, β₂}; (vii) individuals where xi>yi2 can disperse to any group in G or A, and otherwise they died before dispersing again; and (viii) repeat steps vi and vii until all individuals are either dead at age x_i or have out-migrated at age y_i.

All individuals that had out-migrated were assigned the indicator d_i = 1 and those that died before migrating were assigned d_i = 0. We then took the results to construct the log-likelihoodln ℒ(γp|y,x)=∑di=1ln[g(yi|γp)]+∑di=0ln[1−G(xi|γp)],[10]

where γ_p is a vector of natal out-migration prior parameters to be estimated. These parameters were then used in the main model as uninformative priors (with large SEs) for natal out-migration in Eq. 7.

We followed a similar procedure to estimate the prior parameters for secondary out-migration applied to immigrants in the main model. To do this, we used all individuals that had out-migrated from the previous results and assumed that they were moving into the study population. Here is the algorithm for secondary out-migration (which applies only in species that show secondary out-migration): (i) All immigrants are potential secondary dispersers; (ii) simulate age at dispersal (yit) by randomly sampling from a gamma distribution with parameters {α₂, β₂}; (iii) if xi>yit, allow the individual to move to any group in G or in A with same probability; (iv) all dispersers that moved to A have out-migrated and are assigned an age at dispersal y_i, and those that moved to set G have remained in the study population; and (v) repeat ii and iii with all individuals that moved to G until all either are dead at ages x_i or have out-migrated at ages y_i.

All individuals that had out-migrated were assigned the indicator d_i = 1 and those that died before migrating were assigned d_i = 0. We then took the results to construct the log-likelihood for ages at secondary out-migrationln ℒ(λp|y,x)=∑di=1ln[g(yi|λp)]+∑di=0ln[1−G(xi|λp)],[11]

where λ_p is the vector of immigrant out-migration prior parameters to be estimated. These parameters were then used in the main model as informative priors (with low SEs) for immigrant out-migration in Eq. 7.

Regressions Between Life Expectancy and Lifespan Equality.

In Fig. 4A we depict two regression lines between lifespan equality and life expectancy. The gray line corresponds to the regression over the 22 human populations. Because the data included eight time spans for Sweden, we performed the Durbin–Watson test (68) to determine the level of serial autocorrelation in the residuals. The test was not significant (D > D_u with D = 1.48, upper critical value D_u = 1.174 for alpha level 0.01).

For the primate regression line (i.e., yellow line in Fig. 4A) we used all six nonhuman primate estimates and averaged the estimates for Sweden 1751–1759 and 1800–1809 and England 1600–1725, as well as the Ache, Hadza, and acculturated hunter-gatherer estimates. Thus, we had a single pair of values for each of the seven species, with the values for humans being averaged over six populations. We then implemented a weighted phylogenetic generalized least-squares regression model on the female points, with weights given bywi=1Ni∑j=1Ninj, for i=1,…,7,

where N_i is the total number of datasets for the data point i, and n_j is the sample size for the jth population. The phylogenetic component of this regression model was the variance–covariance (VCV) matrix, which we based on the phylogenetic relationships (and therefore statistical nonindependence) of the seven species. To account for phylogenetic signal we transformed this VCV matrix using Pagel’s lambda (69), which we optimized using maximum likelihood. We obtained a consensus tree for the seven primates from the 10kTrees website (70), but substituted Brachyteles arachnoides for Brachyteles hypoxanthus because phylogenetic data for the latter were not available.

The resulting coefficients of this phylogenetic generalized least-squares regression were ε^0i=− 0.18+0.014 e0i, where ε^0i denotes the estimated lifespan equality for the ith population and e0i life expectancy (slope: t = 3.34, P = 0.02, df = 7), with an optimized value of Pagel’s lambda of λ = 0.861 (i.e., strong phylogenetic signal).

Upper Limit for Density Curves in Fig. 1.

As a general rule the ages to plot the probability density functions in Fig. 1 range from age 1 y to the age when l_x = 0.01; this is when 1% of the cohort is still alive. The only exceptions are the curves for the muriqui due to the large amount of censoring in the data and the resulting uncertainty in the curves. In addition, for the nonhuman primates, we shaded the polygons and used dashed lines after the maximum estimated age in our populations to highlight that after those ages the curves correspond to predicted values.

Data Preparation for Additional Mammal Species and Analysis.

We obtained life tables for females of nine additional mammal species from seven taxonomic families and ranging in body size from the ∼50-g tundra vole (Microtus oeconomus) to the 800-kg gaur (Bos gaurus). These data were obtained from the DatLife database (www.demogr.mpg.de/en/laboratories/evolutionary_biodemography_1171/projects/datlife_the_demography_of_aging_across_the_tree_of_life_database_744.htm), and the ultimate sources are given in Table S3. For each life table we calculated life expectancy at birth (e₀) and lifespan equality (ε₀), using the equations given in Box 1, as for humans and other primates.

We then implemented a phylogenetic generalized least-squares regression model on the female points for the six primate species and the nine additional mammals. The phylogenetic component of this regression model was the VCV matrix, which we based on the phylogenetic relationships (and therefore statistical nonindependence) of the 15 species. To account for phylogenetic signal we transformed this VCV matrix using Pagel’s lambda (69), which we optimized using maximum likelihood. We obtained a consensus tree for the 15 mammals from ref. 71 hosted at the Evo10 website (https://www.evoio.org/wiki/File:Bininda-emonds_2007_mammals.nex). Because the phylogenetic data for B. arachnoides and Papio cynocephalus were not available, we substituted them for B. hypoxanthus and Papio hamadryas, respectively.

The resulting coefficients of this phylogenetic generalized least-squares regression were ε^0i=− 0.2 + 0.018 e0i,, where ε^0i denotes the estimated lifespan equality for the ith population and e0i life expectancy (slope: t = 1.99, P = 0.07, df = 13), with an optimized value of Pagel’s lambda of λ = 0 (i.e., no phylogenetic signal).

All statistical analyses were performed in the open-source free package R (72).