# Forest carbon sink neutralized by pervasive growth-lifespan trade-offs

Sep 8, 2020

### Tree-ring data

We used tree-ring records from over 210,000 trees of 110 species, distributed globally in habitats ranging from the tropics to the Arctic region over more than 70,000 sites (Supplementary Fig. 1, Supplementary Table 1). The largest publicly available data source from which we used data is the International Tree-Ring Data Bank (ITRDB, https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). These were complemented with other datasets to maximize the number of records for each species and to fill in spatial gaps. A particularly large tree-ring dataset used in our analyses is the National Forestry Inventory data from the Ministère des Forêts de la Faune et des Parcs from Quebec, Canada52,53 (hereafter NFI-Quebec). This data consists of a complete set of ring-width data from 156,711 trees from 79.381 sites across the province of Quebec. Field tree-ring data were collected according to specific standard protocols52,53, which consisted of selection of up to nine trees in each plot, with 3–5 trees (>91 mm diameter at breast height, DBH) randomly selected, 1–2 selected from the largest trees, and 1–2 from trees closest to the mean tree diameter of the plot52,53. From selected trees one core per tree was collected. Tropical tree-ring data were compiled from the ITRDB and from unpublished and published records26,27,54,55. For species with larger sample sizes, we distinguished between tree-ring data from trees that died naturally before the moment of sampling and trees that were alive at the moment of sampling, allowing us to test the assumption that living tree ages can be used to estimate trees natural lifespans (see section “Trade-off estimates and assessment of possible artefacts”). Part of these dead tree data were obtained from the ITRDB by selecting tree-ring records of which the last measured ring width was dated to before AD1900. We assumed that most of these trees must have been dead at the time of sampling as no records were collected before 1900. In addition, we compiled published dead tree data from refs. 21,56, and used subfossil tree-ring data from refs. 57,58. Supplementary Table 1 provides an overview of the datasets, and full details of each dataset are available online as supplemental info.

Various data controls and selection procedures were used to assure high confidence in our dataset. Where possible we tried to identify duplicate records, i.e., multiple cores taken from the same individual tree. This is only a problem for ITRDB, but not for the NFI-Quebec dataset where only one core per tree existed, or for datasets from co-authors. We thus merged ITRDB records that had identical ID′s except for the last character of their ID (e.g. 01a, 01b, or ID1-1, ID1-2, etc). From the ITRDB, we only used species for which we could obtain data from a minimum of 3 different sites with at least 20 records each, and only selected species that had a minimum total of 100 separate ring width series. We excluded those sites from the ITRDB that showed relative even age structures, and are thus unlikely to represent old-growth populations that provide robust estimates of trees’ maximum lifespans. To this end, we calculated for each ITRDB site the coefficient of variation in tree ages (CVAge = StandDevAge/ MeanAge × 100) and excluded sites with a CVAge lower than 10%. A large subset of ITRDB-data from 46 species has previously been inspected for data quality by co-author S. Voelker59. In this subset of data, each ring width series was manually re-aligned by cambial age (i.e., ring number from pith), providing more reliable estimates of tree ages. For the datasets that were not acquired from the ITRDB, we used slightly different criteria. From NFI-Quebec, we used all available sites, excluding those that were classified with evidence of recent management (commercial thinning or clear-cutting) and where fire or insect disturbances destroyed more than 25% of the forest cover.

For estimates of species-level early growth rates and lifespans (cf. Fig. 1a), we included only species with a minimum of 30 records, as lower sample size is unlikely to provide good approximations of tree life spans. This resulted in the inclusion of 110 species, with a median sample size of 305 trees and 12 sites per species. To assess within-species relationships between early growth and tree lifespan, we included 82 species. As a general rule, we included only species with more than a total of 150 trees and from at least 3 sites. About half of our species had more than 300 tree records (see Supplementary Information).

To assess what minimum sample size is needed to get a representative estimate of the true maximum age of a species or a site, and to evaluate how sample size affected estimates of trade-offs between early growth and longevity, we randomly resampled 500 times varying sample sizes—from 25 to 600 trees—from a subset of 11,752 Picea mariana trees from NFI-Quebec sites located north of 50.7°N. Comparison of the maximum ages of these random subsets of trees with the true observed maximum ages shows that a sample size of 100 trees results in 99.4% of the cases in maximum age estimates larger than the 95th percentile of the original dataset, and in 67.2% of the cases in ages larger than the 99th percentile original age (Supplementary Fig. 2a, b). As more than 70% of the species had at least 100 trees we thus assumed that for most species, the estimates of their lifespan were close to true lifespans. We used this same approach to assess how sample size affected the estimate of the trade-offs (i.e., estimation of the negative exponential decay constant; see next section for details). This analysis showed that sample size of 300 trees (corresponding to median sample sizes for trade-off analysis), leads to mean errors in the estimated slope of 12% (Supplementary Fig. 2c). Thus, for most species we achieve relatively accurate estimates of the trade-off strength. Low sample sizes for some species will nevertheless result in small errors of the mean slope, but we expect that positive and negative errors will cancel out against each other. Indeed, we do not observe a specific bias towards over- or under-estimation for low sample sizes, as the mean exponential decay constant for a simulated sample size of 150 trees is very similar to that observed (i.e., −0.409 versus −0.399).

### Trade-off estimates and assessment of possible artefacts

The strength of the trade-offs between growth and tree lifespan was assessed for each species using a 95th quantile regression between mean early growth rate and the natural logarithm of age using the QUANTREG package in R60, as

$$begin{array}{*{20}{c}} {{{log}}(A_{(95{mathrm{th}},{mathrm{quantile}})}) = a + b cdot overline {{mathrm{RW}}} } \ {{mathrm{or}}} \ {{mathrm{Lifespan}} approx A_{left( {95{mathrm{th}},{mathrm{quantile}}} right)} = {mathrm{exp}},left( {a + b cdot overline {{mathrm{RW}}} } right)} end{array}$$

(1)

where A is age of the tree, (overline {{mathrm{RW}}}) is the mean ring width over the first 10 years. The constant b describes the negative exponential decay constant (i.e., exponential rate of decrease of tree lifespan with increasing early growth rate). This quantile regression fit results in similar estimates of the maximum ages of trees as the 95th percentile ages in binned early growth rate categories (see Supplementary Fig. 3a–c). Note that in contrast the maximum diameter does not vary strongly between slow and fast-growing trees (Supplementary Fig. 3d). We chose a relative short period, the first ten years, for estimating early growth as our study included some relative short-lived species. Previous studies have found similar results when using longer periods (50 years)56, and we expect no substantial difference using different early growth periods as tree growth is usually strongly auto-correlated in time61.

To assess trade-off strengths within species, we calculated the mean decay constant (b, Eq. 2) for each species using relative age, A/max(A), and relative mean early ring width, (overline {{mathrm{RW}}})/max((overline {{mathrm{RW}}})). Maximum of A and (overline {{mathrm{RW}}}) are species level maxima. The mean slope calculation across all species was weighted by the cube root of the sample size to account for the large differences between species in sample size, and confidence of the trade-off estimates.

While these relationships suggest true trade-offs, they may also be affected (or even driven by) the approaches or analytical methods used here. In particular, we here evaluate the effect of the following four possible artefacts on our results; (1) the use of living trees to estimate tree lifespans, (2) effect of recent growth increases on early growth-age relationship, (3) effects of pith offsets and wood decay on early growth-age relationship, (4) sampling artefacts, such as disproportionate selection of large trees.

1. (1)

Use of living trees: our analysis includes mostly trees that were sampled when still alive, and may thus not be representative for the true lifespan trees may achieve. To assess to which degree use of living trees may affect our results, we analyse and compare the trade-off strengths of trees that died before 1900 to living trees for 12 species with sufficient data availability (minimum of 150 dead and 150 living trees). As the slopes of dead and living trees do not differ significantly (Supplementary Fig. 5f, g, paired t-test exponential decay coefficient, t = −0.1095, p = 0.915, n = 12), we conclude that 95th quantile regressions on living trees can be used to approximate tree lifespan.

2. (2)

Effect of recent growth increases: recent growth stimulation of trees due e.g. to CO2 fertilisation, warming in higher latitudes, and/or nitrogen deposition, may result in observation of a trade-off. This is because recent increases in growth will lead to higher early growth rates for young trees compared to old trees, resulting in a negative relationship between early growth and tree age. The comparison of trade-off strength of dead versus living trees provides strong evidence that this effect does not drive the trade-off. In addition to this, we use a data driven forest simulation (see section “Examining the effects of growth stimulation on forest dynamics”) to assess how growth increases affect estimations of the trade-off strength. In this simulation, we used the actual tree-ring data to simulate realistic growth increases of Picea mariana tree-ring trajectories in response to high latitude warming. By sampling from these trajectories at the end of the growth increase period (i.e., year 350), and in a period without any recent growth increases (i.e., year 600, see Fig. 3e), we establish that growth increases result in only a small over-estimation of the trade-off, decreasing the exponential decay coefficient from −0.37 to −0.44 (see Supplementary Fig. 9c). Thus, it is unlikely that recent growth stimulation is the cause for the negative relation between early growth and tree lifespan.

3. (3)

Pith offset: tree-ring data, especially those acquired from ITRDB, may miss the innermost sections due to incomplete cores, decayed centres, or imperfect increment borer alignment. Missing rings will result in underestimation of tree ages and inaccurate early growth rates estimation and could thus affect the estimated relationship between early growth and lifespan. However, ring widths in most species decrease with tree age and size17, and even trees showing constant wood production with age, will show decreasing ring width because of geometry. Thus early growth in these samples will underestimate true growth rates and would most likely weaken the observed trade-off, rather than strengthening it. A comparison of species present in both the NFI-Quebec and the ITRDB datasets confirms this. The NFI-Quebec dataset was less affected by pith offset problems, as the trees were carefully screened and trees with substantial differences between cumulative ring widths and field diameters were excluded. Yet, we find that slopes were more negative for NFI-Quebec compared to ITRDB (mean b of −0.25 for Quebec vs. −0.10 for ITRDB, two-sided paired t.test, t = 2.49, p = 0.047, n = 7 species) and pith offsets thus do not explain the relationship. This comparison also shows that estimates of the strength of the trade-offs between early growth and longevity inferred from ITRDB data are probably conservative, as the Quebec data can be considered to be of higher quality, and were collected according to standard protocols. In contrast, data from the ITRDB may contain incomplete series and were collected for unknown purposes, and these issues probably weaken trade-offs in the ITRDB.

4. (4)

Sampling biases: one potential bias in our dataset may arise due to the tendency of tree-ring studies to sample predominantly large trees in the field (i.e., big tree selection bias62,63,64). This may result in a negative relationship between early growth and tree age, as young slow-growing trees tend to be underrepresented in the tree-ring sample (i.e., have not reached the field minimum size threshold yet), compared to fast-growing young trees that are much larger, and therefore more likely to be sampled. This effect would reduce the number of trees with slow early growth and young ages in the tree-ring sample (i.e., trees in the lower left-hand corner of the early growth-lifespan graphs, cf. Supplementary Fig. 6a), and results in overestimation of the 95th percentile age estimates for slow-growing trees. Our approach to estimate to which degree this bias affected our estimates of growth-lifespan trade-offs was as follows. We first used the tree-ring NFI-Quebec data of Picea mariana, combined with plot data from Quebec to reconstruct a new artificial tree-ring dataset with a size frequency distribution identical to the population size frequency distribution for this species in Quebec (Supplementary Fig. 6b). For each tree of Picea mariana sampled for their tree-rings we know the early growth rate and age, and also the complete diameter- and age-trajectory up to the year of sampling. From these data, we resampled for each size class (in bin widths of 2 cm) the same number of trees as that observed in the field. By doing this we filled in the lacking data of trees smaller than 91 mm, and created a new artificial tree-ring dataset that had an identical size structure to that observed in the field. We know the mean growth rate over the first ten years and the age at which each individual tree reached the diameter of their respective size class, and could thus reconstruct the early growth rate versus tree age graphs for the full population, including the smaller size classes which were missing from our original tree-ring sample. We then compared the early growth -lifespan relationship for the complete population to that of the trees larger than 91 mm, mimicking the NFI-Quebec field collection protocol. This shows that the exponential decrease is marginally larger (b = −0.505 compared to −0.470 for trees >91 mm) and that the intercept is lower (159 years compared to 220 for trees >91 mm) when sampling all trees compared to only trees with diameters >91 mm (Supplementary Fig. 6). Hence, the use of a minimum size threshold (91 mm) in the NFI from Quebec results in a slight underestimation of the trade-off (by ~7%) for the Picea mariana dataset. We also resampled from this artificial dataset the 10% largest trees, to mimic a hypothetical standard tree-ring sampling scenario that only samples the largest trees. Such a sampling scenario resulted in a decay constant of −0.432, thus again causing a small underestimation of the true trade-off. This simulation proves that the trade-off is not a result of a sampling bias.

### Possible environmental drivers of the trade-offs

We evaluated whether the observed trade-off between early growth and tree lifespans could be caused by covariance of growth and lifespan with climate, soil or competition. Temperature variation for example reduces tree growth and lifespan in various species (cf. Fig. 2a, b). To this end, we calculated site-level mean early growth rates and the maximum tree age for a set of species covering different geographic regions (North America, Europe and Quebec). For Quebec, we combined multiple nearby locations to obtain a minimum of 30 trees per site, as sample sizes were low for each location. Site-level mean annual temperature and precipitation was obtained from WorldClim65. We then assessed for nine different species the effect of temperature and precipitation on site-level mean early life growth and maximum tree age using major axis regression from the package smart-366. These analyses confirm that early life growth is positively related to temperature for all nine species studied, and that lifespan decreases significantly with temperature for seven out of nine species (see Fig. 2a, b). Using linear mixed effect models with species as random factor (nlme-package-R67), we find that across all nine species, mean early life growth increases on average by 0.11 mm for each degree temperature increase, while lifespan decreases by 13 years for each degree temperature increase. Precipitation has no significant effect on early life growth or tree lifespan.

To disentangle whether lifespan decreases are a direct effect of temperature increases, or due to increases in early life growth, mixed effect models were run for all nine species that simultaneously included temperature and mean early life growth rate as explanatory variables for variation in tree lifespan. To account for species differences in growth and age, we used relative mean early ring width ((overline {RW})/max((overline {RW}))), and relative maximum age (A/max(A)), and used species as a random factor with random intercepts for both early life growth and temperature. This analysis shows that mean early life growth is a stronger predictor of tree lifespans than temperature (t-value early growth = −7.2, p < 0.001, t-value temperature = −2.5, p = 0.012). A similar analysis for Picea mariana alone confirms that the primary driver of lifespan is the mean early life growth rate and not temperature, or other environmental variables. For example, analysis of early life growth vs. lifespan in different temperature classes (of 2 °C) remains strong, while the relationship between temperature and lifespan breaks down when this is analysed in growth rate classes (Fig. 2d). Similarly, trade-offs remain strong for Picea mariana even when analysing the data in crown cover classes (i.e., an indication of the stand level competition), or when analysing data in different soil classes (Supplementary Fig. 8).

### Examining the effects of growth stimulation on forest dynamics

We examine the effect of the observed growth longevity trade-off on forest dynamics (growth, mortality, and standing stocks) for a realistic growth stimulation as expected from changes in temperature. We do this by using a data driven forest simulation approach in which we use observed tree-ring data and realistic estimates of size-related mortality rates that result in a trade-off between early growth and lifespan very similar to the observed trade-off. The growth stimulation of trees was estimated using the temperature sensitivity of tree-ring data from the large temperature gradient for Quebec (cf. Supplementary Fig. 7). The full approach is described below.

As demonstrated in our main manuscript, tree-ring data from Quebec reveal a tree growth-longevity trade-off, which may be mediated by rapid increase of mortality with increasing diameter. Not only tree cores, but also detailed forest census data exist for this region, allowing reconstruction of size-dependent mortality relationships. This analysis was done using only data from Quebec forest inventories north of 50.7°N. First, we estimate mortality, μ, as a function of tree diameter using both forest census data inventories and individual tree-ring records from Picea mariana from NFI-Quebec. We assume a stationary state of the tree diameter (D) number distribution, N(D). Thus, for a diameter class with D in the interval [D,D+δD]

$$0 = frac{{dN(D)}}{{dt}} = I(D) – L(D) – mu cdot (D).$$

(2)

Here I(D) is the number of individuals growing into diameter class [D,D+δD] per unit time, and likewise L(D) the number of individuals leaving the same diameter class per unit time. Thus, mortality can be estimated from the steady diameter distribution, obtained from plot data (Supplementary Fig. 7b), and number of in-growers and leavers, estimated from the diameter distribution and ensembles of randomly selected growth trajectories, as

$$mu left( D right) = frac{{I(D) – L(D)}}{{N(D)}}.$$

(3)

The result is shown in Fig. 3a, together with a fourth order polynomial model of the form

$$mu left( D right) = left{ {begin{array}{*{20}{c}} 0 & {for,D ,le, D_0} \ {b + k cdot left( {D – D_0} right)^4} & {for,D ,> , D_0} end{array}} right..$$

(4)

Here b = 0.025, k = 3 × 10−11 and D0 = 91 mm, which is the minimum sampling diameter used for the NFI-Quebec. The rationale to set the mortality below D0 to zero in this model is that the trees sampled for tree rings did all survive to this diameter (as only trees with diameters >91 mm were cored). Thus, only by setting mortality to zero for trees <91 mm allows for a proper comparison of simulation results with observed data. We ended up using this model as it reproduced the observed trade-off between growth and longevity quite well (see Fig. 3b). The strength of the simulated trade-off is robust with regard to the choice of mortality function as linearly increasing mortality rates equally reproduced the observed trade-off. Increasing mortality rates towards larger size classes are consistent with observations of size-dependent mortality in temperate and tropical trees35,36,40,68, and supports the notion of a species-specific maximum size threshold64. For Picea mariana, the 99th percentile maximum tree size is 353 mm, close to diameters at which we find large increases in tree mortality (cf. Fig. 3a). For comparison purposes, we also calculated an age-dependent mortality rate and performed alternative simulations using this mortality model. The age-dependent mortality relationship was derived in a similar manner to the size mortality curve by repeating the above procedure and calculations (Eqs. 2 and 3) for (yearly) age classes. Age-dependent mortality was parametrized as:

$$mu (A) = left{ {begin{array}{*{20}{c}} 0 & {for,A le A_0} \ {a + b^ast A} & {for,A ,> , A_0} end{array}} right.,$$

(5)

with a = 0.021 and b = 0.0000015. Analogous to the diameter-dependent mortality model, we set mortality to zero for trees younger than 74 years (A0), which is the age at which the average Picea mariana tree reaches 91 mm in diameter (D0).

We next created a 600 year sequence of annually seeded, 1250 member, tree cohorts. Each member of a cohort is a randomly selected tree diameter growth trajectory derived from the tree-ring cores, extended in time to an age of 500 years. Short growth trajectories were extended by using the mean growth of the 10 oldest trees that had a similar ring width in the first ten years of growth. To this end all early mean ring width was grouped into six equal early ring width classes.

All trees of this so-constructed tree cohort set live, by construction, exactly 600 years. Realistic age structures were realised by sequentially (year-on-year and tree-by-tree) assigning death to trees where a random number generator identified those individuals  smaller than μ(D) or μ(A). To test the realism of this procedure, we compared the predicted and observed tree age versus early ring width relationship—or i.e. the growth-longevity trade-off. The slope of the relationship for the diameter-dependent mortality model μ(D) is very close to observed (Supplementary Fig. 9a), and is thus a realistic representation of the observed mortality process and justifies its use to examine the effect of a growth stimulation on standing stocks. In contrast, we find that the age-dependent mortality model μ(A) does not result in a significant trade-off between early growth and tree lifespan (Supplementary Fig. 9b), providing an ideal comparison for models that fail to incorporate the observed trade-offs.

To mimic growth stimulation, we boosted growth of trajectories from year t0 = 300 year onwards of the 600 year sequence of cohorts, while exposing the trajectories over the entire 600 year period to the mortality algorithm just described. We stimulated growth rate, (overline {{mathrm{RW}}}) (mm year−1) from year t0 = 300 year onwards according to

$$overline {{mathrm{RW}}_{stim}} left( t right) = overline {{mathrm{RW}}} left( t right) cdot exp left( {lambda left( A right) cdot delta Tleft( t right)} right),$$

(6)

and

$$delta Tleft( t right) = left{ {begin{array}{*{20}{c}} {frac{{overline {dT} }}{{dt}} cdot left( {t – t_0} right),} & {t – t_0 ,<, tau } \ {{mathrm{const}}}, & {t – t_0 ge tau } end{array}} right.$$

(7)

where δT(t) is a normalized temperature trend (year−1), λ(A) a unitless function of tree age representing that growth sensitivity of young and old trees to temperature may vary with tree age59, and τ is 50 years, the duration of the period of growth stimulus. We calculated the observed temperature trend for Quebec over the past 100 years from CRUTEMP data69, and simulated growth increases from year 300 to year 350 in response to observed warming rate, estimated to be 0.0221 °C year−1. We used a space-for-time substitution approach on the full dataset of Quebec43 to estimate the ring width response of Picea mariana to temperature. We conducted these simulations in age-bands of 10 years (0–10, 10–20, …. 140–150, >150 yrs) as younger trees are more sensitive to temperature increases than older trees (Supplementary Fig. 9d), and use an exponential model of the form

$$overline {{mathrm{RW}}} left( {A,,T} right) = a cdot {mathrm{exp}}left( {lambda left( A right) cdot delta T} right),$$

(8)

where a is a constant, λ(A) is the exponential increase rate for age class [A, A+ δA] and δT = (T0) (°C). We then used the exponential growth increase with temperature for each age band, λ(A), to estimate the relationship between tree age and λ(A) (Supplementary Fig. 9e). In all cases the age modulation of the stimulus is

$$lambda left( A right) = left{ {begin{array}{*{20}{c}} {0.0000132 cdot A^2-0.00291 cdot A + 0.229,} & {A , < ,135} \ {0.07,} & {A , ge ,135} end{array}} right..$$

(9)

Finally, we compared the effect of growth stimulation on forests dynamics against simulation without growth increases (baseline) for the two different mortality models. We also performed one simulation where we multiplied the full growth series by 2, as a representation of the effect of growth stimulation on a faster-growing species. Finally, we evaluated for each simulation the mean ring width growth, the stem mortality rate, the age of the largest (75th percentile) trees that died and the change in the total basal area stock over time for the full population. For the stem mortality and the basal area stocks, we calculate and present the change in dynamics of the growth stimulation scenario relative to the baseline scenario without growth stimulation. All analyses and simulations were performed using R-studio, version 0.99.90370. Maps in SI figures were produced using the ggmap function from the R-package ‘ggplot2’71