Study area
This study was conducted using tree census data collected from 119 forest inventory plots (73 tropical, 46 subtropical) situated across a latitudinal range of 7.1°N (Colombia) to 27.8°S (Argentina), a longitudinal range of 79.5° to −63.8° W, and an elevation range of 500–3511 m asl (Fig. 1). The mean annual temperature (MAT) of plots ranged from 7.3 to 23.8 °C (mean = 16.7 ± 4.1 °C; mean ± SD) and mean annual precipitation (MAP) of the plots ranged from 608 to 4313 mm y−1 (mean = 1405.0 ± 623.9 mm y−1) (External Databases 1). The number of plots sampled in each country was: Argentina = 46, Bolivia = 26, Peru = 16, Ecuador = 21, and Colombia = 10 (Fig. 1). The 119 forest plots ranged in size from 0.32 to 1.28 ha and represent a cumulative sample area of 104.4 ha (horizontal areas corrected for slope) that containe more than 63,000 trees with a diameter at breast height (DBH, 1.3 m) ≥10 cm (External Database 1). Ninety-four of the plots (79.0%) were ≥1 ha in size. Neither secondary forests nor plantations were included. However, only seven of the plots (five in Argentina and two in Bolivia) were located in forests >100 km2 in extent41, which suggests that at least the edges and borders of some plots could have experienced some degree of disturbance or degradation. All plots were censused at least twice between 1991 and 2017 (census intervals ranged between 2 and 9 years).
In each plot, we tagged, mapped, measured, and collected vouchers of all trees and palms (DBH ≥ 10 cm). DBH was measured 50 cm above buttresses or aerial roots when present (where the stem was cylindrical). During the second or subsequent set of censuses, DBH growth, recruitment, and mortality were recorded. In cases where the recorded DBH growth of the second census was less than −0.1 cm y−1 or greater than 7.5 cm y−1, the DBH of the second census was augmented/reduced in order to match these minimum/maximum values42. To homogenize and validate species names of palms and trees recorded in each country and plot, we submitted the combined list from all plots to the Taxonomic Name Resolution Service (TNRS; http://tnrs.iplantcollaborative.org/) version 3.0. Any species with an unassigned TNRS accepted name or with a taxonomic status of ‘no opinion’, ‘illegitimate’, or ‘invalid’ was manually reviewed. Families and genera were changed in accordance with the new species names. If a full species name was not provided or could not be found, the genus and/or family name from the original file was retained.
Aboveground carbon stocks
The aboveground biomass (AGB) of each tree was estimated using the allometric equation proposed by Chave et al43., defined as: AGB = 0.0673 × (WD × DBH2 × H)0.976 where AGB (kg) is the estimated aboveground biomass, DBH (cm) is the diameter of the tree at breast height, H (m) is the estimated total height, and WD (g cm−3) is the stem wood density. To estimate WD, we assigned the WD values available in the literature44 to each species found in each plot. In cases where we could not assign a WD value at the species level, we used the average value at the genus- or family level. For unidentified individuals, we used the average WD value of all other species in the plot. Tree height (H) was estimated (see below) based on the heights measured on a subset of the individual stems in each plot using digital hypsometers or clinometers. The estimated AGB of each tree was then converted to units of aboveground carbon (AGC) by applying a conversion factor of 1 kg AGB = 0.456 kg C45. The AGC per ha was then determined by converting kg to Mg, summing the values for all trees in a plot, and extrapolating or interpolating to a sample area of 1 ha.
Estimates of AGB and AGC are highly dependent on tree height. Unfortunately, tree height was difficult or impossible to measure on all stems due to physical and logistical constraints. Therefore, we estimated the height of each stem based on allometric relationships between DBH and tree height that we developed for each plot based on height and DBH measurements taken on a subset of individuals. Although the AGB/AGC estimates are only for trees with DBH ≥ 10, we used trees with DBH ≥ 5 cm to construct the H:DBH models when possible in order to be as comparable as possible with the existing pantropical H:DBH models46. In total, 44,442 trees had their heights measured in the field and were employed to construct the H:DBH models. The percentage of trees with direct field measurements of H (DBH ≥ 5 cm) in each country was: Argentina = 19%, Bolivia = 98%, Peru = 96%, Ecuador = 97%, and Colombia = 46%. In Argentina, 32 of 46 plots did not have any field measurements of H, while all plots in all other countries had field measurements of H for at least a subset of trees.
We tested and compared the expected effects of using H:DBH models constructed using the local (plot), country, or pantropical (regional) level data. To select the best model to estimate H from DBH at the plot and country level, we used the function modelHD available in the BIOMASS package for R47. We chose the best allometric model from four candidate models (two log-log polynomial models, the three-parameter Weibull model, and a two-parameter Michaelis-Menten model (Supplementary Table 7)) by selecting the model with the lowest RSE and bias (Supplementary Table 8). At the regional level, we used a pantropical model46. The use of country and pantropical H:DBH allometries underestimates tree heights in the lowlands and overestimates tree heights in highlands, thereby homogenizing AGB estimates along elevational gradients10,48 (Supplementary Figs. 11, 12, 13). Using plot level allometries eliminates this problem. However, in the 32 plots in Argentina where we had no information about tree height, we used the country-level H:DBH model developed with the data available in the remaining 14 plots to estimate the height of each tree, which could have homogenized the AGC estimates along the Argentinian elevational gradient (Supplementary Figs. 11, 12, 13).
Aboveground carbon dynamics
The AGC dynamics of each plot was estimated from the annualized values of AGC mortality, AGC productivity (AGC change due to recruitment + growth), and AGC net change3. The calculations of the separate AGC dynamic components was performed as follows: (i) AGC mortality (Mg ha−1 y−1) = the sum of the AGC of all individuals that died between censuses divided by the time between measurements. (ii) AGC recruitment (Mg C ha−1 y−1) = the sum of the AGC of individuals that recruited into DBH ≥ 10 cm between censuses divided by the time between measurements. However, for each tree recruited (DBH ≥ 10 cm), we subtracted the corresponding AGC associated with a tree of 9.99 cm (i.e. just below the detection limit) in order to avoid overestimations of the overall increase in AGC due to recruitment49. (iii) AGC growth (Mg ha−1 y−1) = the sum of the increase in AGC of all individuals with DBH ≥ 10 cm that survived between censuses divided by the time between censuses. (iv) AGC net change (Mg ha−1 y−1) = the difference between AGC stock in the last census (AGCfinal) and AGC stock in the first census (AGC1) divided by the elapsed time (t; in years) between measurements [(AGC net change = AGCfinal − AGC1)/t]. We recognize that these methods exclude C stored in soils or in belowground tissues9,48; however, quantifying just aboveground C stocks and fluxes provides valuable information about the overall status of these forests as net C sinks or sources.
Climate
Climate variables at each plot location were extracted from the CHELSA28 bioclimatic rasters at a resolution of 30-arcsec (~1 km2 at the equator). The climate variables extracted were: Mean Annual Temperature (MAT), Mean Diurnal Range (MDR), Isothermality (Isoth), Temperature Seasonality (TS), Maximum Temperature of Warmest Month (MaxTWarmM), Minimum Temperature of Coldest Month (MinTCM), Temperature Annual Range (TAR), Mean Temperature of Wettest Quarter (MeanTWarmQ), Mean Temperature of Driest Quarter (MeanTDQ), Mean Temperature of Warmest Quarter (MeanTWetQ), Mean Temperature of Coldest Quarter (MeanTCQ), Mean Annual Precipitation (MAP), Precipitation of Wettest Month (PWetM), Precipitation of Driest Month (PDM), Precipitation Seasonality (PS), Precipitation of Wettest Quarter (PWetQ), Precipitation of Driest Quarter (PDQ), Precipitation of Warmest Quarter (PWarmQ), Precipitation of Coldest Quarter (PCQ). We separated all variables associated with temperature (°C) from those associated with precipitation (mm y−1) and applied a Principal Component Analysis (PCA) to the 11 variables associated with temperature (PCAtemp) and a separate PCA to the eight variables associated with precipitation (PCAprec). The first two principal components of both PCAtemp and PCAprec (four PCA axes in total) were selected for use in subsequent analyses. Plot elevations were estimated based on their coordinates and the SRTM 1 ArcSec Global V3 (https://lta.cr.usgs.gov) 30 m resolution digital elevation model (DEM).
PCAtemp1 (Supplementary Fig. 1a) explained 53.0% of the total variance of the temperature variables and had high loading from Isothermality and Maximum Temperature of Warmest Month, which was primarily associated with changes in elevation (r = −0.97, p < 0.001) (Supplementary Fig. 3b). PCAtemp2, explained 45.2% of the total variance of the temperature and had high loading of Temperature Annual Range and Minimum Temperature of Coldest Month. The PCAtemp2 was primarily associated with changes in latitude (r = 0.78, p < 0.001) (Supplementary Fig. 3c). PCAprec1 (Supplementary Fig. 1b) explained 68.9% of the total variance of the precipitation variables was highly loaded by Mean Annual Precipitation and Precipitation Seasonality, and was primarily associated with changes in latitude (r = 0.68, p < 0.001) (Supplementary Fig. 3e). PCAprec2 (Supplementary Fig. 1b) explained 26.3% of the total variance of the precipitation variables was highly loaded by Precipitation Seasonality and Precipitation of Driest Month, and was also primarily associated with changes in latitude (r = 0.64, p < 0.001) (Supplementary Fig. 3g).
PCAprec1 was negatively correlated with PCAtemp1 (−0.55, p < 0.001), indicative of the fact that precipitation generally decreased from the lowlands to highlands. Likewise, PCAprec2 was negatively correlated with both PCAtemp1 (−0.46, p < 0.001) and positively with PCAtemp2 (0.58, p < 0.001). These associations suggest that precipitation seasonality decreases towards the equator, but increases with elevation (Supplementary Fig. 2). Therefore, we only used the PCAtemp1 and PCAtemp2 axes as explanatory variables in the subsequent models to represent an overall climate gradient of elevation and latitude, respectively (Supplementary Fig. 3). The use of PCA axes instead of raw climate variables minimizes multicollinearity among variables and better represents the multidimensional gradient of climate that occurs across our study plots.
Thermophilization rate
We used the Thermophilization Rate (TR; °C y−1)13,14 as a metric of the direction and rate of compositional changes in the tree communities in relation to species’ thermal optima through time. A positive TR indicates an increase in the relative abundances of lowland, thermophilic species (as expected due to the upward shift of species ranges under increasing temperatures) and a negative TR indicates an increase in relative abundances of less-thermophilic, highland species.
To calculate the TR of each plot, we downloaded all available herbarium collection records from the forested regions of tropical Andean countries through the Global Biodiversity Information Facility (GBIF) data portal (www.gbif.org; accessed September 2019). We estimated the mean annual temperatures (MAT) at the collection locations of all specimens by extracting the Bio1 values from the CHELSA climate rasters28. For each of the tree species found in the plots and that were represented by 10 or more GBIF records, we estimated their thermal optimum as the mean temperature across all available collection locations within the forested Andean region50. For species with <10 available records, we estimated the thermal optimum as the mean temperature across the collection locations of all congeneric individuals within the forested Andean region. Then, for each plot census, we calculated the Community Temperature Index (CTI; °C) as the mean thermal optima of all constituent species weighted by their relative basal areas. Finally, we calculated the annual TR (°C y−1) of each plot as the annualized change in the CTI over the entire census period of each plot13,14.
Taxonomic and phylogenetic diversity
All species and genus names were checked and standardized using the Taxonomic Name Resolution Service51. In the dataset, 91.3% of stems were identified to species level, 7.3% to genus, 0.8% to family, and 0.6% remained unidentified. To estimate species diversity while accounting for differences in plot size and stem numbers, we used the rarefied Species Richness (SR)52 at the minimum stem number found per plot (86 individuals). To calculate the phylogenetic diversity sensu stricto (PD)53, we first generated a phylogenetic tree of hypothesized relationships using phylomatic for the complete species list as recorded across all 119 plots (excluding unidentified taxa). We used the bladj algorithm to date the phylogenetic tree by adjusting phylogenetic branch lengths to respective fossil ages54. The PD of each plot community was then calculated as the total sum of the phylogenetic branch lengths connecting the co-occurring species in each plot along the minimum spanning path to the root of the tree.
The observed PD was compared to a null distribution to control for the sampling effects and differences in regional diversity. The null model used an independent swap algorithm that maintained the frequency and richness of species in each plot while randomizing community composition55. The standardized effect size of the PD (PDz) was then calculated by subtracting the expected mean PD derived from the null distribution of 999 random draws to the observed PD value in each plot, divided by the standard deviation of the null distribution. This metric was estimated using the Picante package in R56. SR and PDz were negatively correlated (r = −0.59, p < 0.001). This negative correlation is mainly due to the increased mixture of temperate-affiliated and tropical-affiliated species, both at higher southern latitudes and at higher elevations in the tropics57,58. However, SR was never a significant explanatory variable in the models due to the greater relative importance of climate variables, and we only used PDz in the models as the surrogate of complementarity effects.
Symbiotic root associations
Symbiotic root associations (SRA), particularly those involving mycorrhizal fungi, are increasingly recognized as important factors influencing plant community productivity and dynamics16,17. To incorporate the potential contributions of SRA to AGC stocks and dynamics in subtropical and tropical Andean forests, individuals were assigned an SRA status either as arbuscular mycorrhizal (AM) or ectomycorrhizal fungi (EcM) based on genus- or family-level designations18. Overall, the number of genera assigned an SRA status was >90%. We chose these two taxonomic levels for three reasons: (1) assigning species-level taxonomy was difficult for some individuals in our plots; (2) using these higher taxonomic levels greatly increases the ability to provide SRA assignments; and (3) SRA is largely conserved at the genus and family level59. Here, we restricted matches for our genera and families to only those present in North and South America in the compiled reference list18. Any genera or families lacking symbiotic root assignment were manually checked and, when possible, assigned SRA on the basis of primary literature searches.
We next calculated the relative representation of the two SRA types (AM vs EcM) within each plot. In calculating the relative abundances, we weighted SRA types by the proportional stem number. Since the association between tree species is generally specialized to a single mycorrhizal fungal type (either AM or EcM), the relative representation of the two groups is a good indicator of different modes of nutrient acquisition and soil C stocks60. Overall, a high AM to EcM ratio is expected in forest soils that are characterized by faster and inorganic-dominated nutrient cycling, whereas low AM-to-EcM ratios are expected in forest soils with slower and organic-dominated nutrient cycling. Additionally, as the AM-to-EcM ratio decreases, there is generally more C per unit N present in soil, particularly in the upper layers61. We used the log-transformed ratio of AM to EcM (e.g. ln (AM/EcM)) as an explanatory variable in our statistical analyses.
Forest cover and forest loss in the Andes
We estimated forest cover and forest loss for the Andes mountains between 11°N and 27.3°S and between 82°W and 56°W between the years 2003 and 2014 from Hansen et al. v1.630 using Google Earth Engine31. We excluded pixels with forest cover <70% from the analysis. We further limited our study area by masking out forests with less than 700 mm of annual rainfall to represent tropical mountain evergreen conditions62 using Worldclim 2.0 dataset63, and forests inside terrestrial ecoregions64. The time frame of forest change assessed was defined by the median of the year of plot censuses carried out between 1991 and 2009 (2003) and the median of the year plot censuses carried out between 2010 and 2019 (2014). The patterns of forest cover and forest loss were then summarized within four elevation bands that represent recognized habitat zones:25 (i) 500–1200 m, (ii) 1200–2000 m, (iii) 2000–2800 m, and (iv) 2800–3500 m.
To define the initial (2003) mean AGC carbon stock in each elevational band we used the initial AGC stocks (first plot census) before 2009 (including it). To define the final (2014) mean AGC carbon stock in each elevational band we used the final AGC stocks (last plot census) after 2009 (excluding it). We did not use the plot AGC net change rates due to the high asynchrony among census’ dates. That said, in some countries the last (final) census in a particular plot could almost always be before 2010 (e.g. Peru), while in others the initial census was always after 2010 (e.g. Ecuador).
We used bootstrapping to assess the mean and 95% confidence intervals (CI) of AGC stocks in each elevation band for both the initial (2003) and final (2014) years defined to assess changes in forest cover. In each elevational band, total AGC stocks were calculated as the product of the mean AGC stock in the initial/final census and the forest cover at 2003/2014, respectively. The overall mean AGC stock for Andean forests along the whole elevational gradient was quantified as the forest cover weighted mean of AGC stocks across elevational bands. Finally, the net AGC (Mg C y−1) balance was calculated as the difference between the total AGC stocks in the final census (2014) minus the total AGC stocks in the initial census (2003) divided by the elapsed time (i.e. 11 years) (Table 1).
Statistical analysis
Disturbance and competitive thinning
We used the size-dependent parameter of mortality (β) derived from the logistic regression to differentiate plots along a gradient of disturbance that ranges from sites strongly influenced by competitive thinning post internal disturbance (low β) to sites more influenced by active disturbances (high β). The size-dependent mortality parameter (β) was calculated as the probability of tree death (P) as a function of DBH according to logistic regression (logit(P) = a + β × DBH; if β < 0, smaller trees have a higher probability of dying due to competitive thinning during post internal disturbance, and if β > 0, larger trees have a greater probability of dying as expected under active disturbances29. The β parameter was also compared with the mean square diameter (Dq) of each plot. Dq was calculated by (Dq = sqrt {frac{{{sum} {{mathrm{DBH}}_i^2} }}{n}}); where DBH = diameter at breast height, and n = the total number of stems in each plot. As a separate means of characterizing competitive thinning, we graphically evaluated the temporal changes in mean square diameter (Dq) and stem density. An increase in the number of individuals along with a decrease in Dq is expected in recently disturbed plots where the large trees have been lost and the recruitment of small stems is increasing; an increase in Dq but a decrease in the number of individuals is expected under competitive thinning29.
Latitudinal and elevational patterns of aboveground carbon stocks and dynamics
We used generalized additive models (GAMs) to evaluate the pattern of change of AGC stocks and AGC dynamics (AGC net change, AGC productivity, and AGC mortality) along the latitudinal and elevational gradients. A One Way Analysis of Variance (ANOVA) was employed to evaluate differences in both AGC stocks and dynamics among countries.
Drivers of aboveground carbon dynamics
To select the most important explanatory variables of AGC dynamics along both latitudinal and elevational gradients in the 119 South American forest plots, we used Structural Equation Modelling (SEM)26 to assess and visualize the likely multiple correlations between explanatory variables and their exogenous or latent influence on AGC dynamics. Overall, the SEM approach was employed to analyze the direct and indirect effect of climate (PCAs), symbiotic root association (SRA), phylogenetic diversity (PDz), directional compositional change (TR), initial aboveground carbon (AGC), and the size-discriminant β parameter associated to disturbance on AGC net change, AGC productivity, and AGC mortality. All variables were standardized before their inclusion in the SEM. The SEM included climate and SRA as exogenous variables. Climatic variables affected all variables in the model. In contrast, SRA only affected AGC1 and AGC dynamics. PDz was included as an explanatory variable of both AGC1 and dynamics. We controlled for the effect of climate and TR on the variation of PDz. AGC1 was affected by all explanatory variables, and affected both β and AGC dynamics. The β parameter directly influenced AGC dynamics, and was also affected by TR. We used a Satorra–Bentler scaled chi-square test statistic to determine whether the covariance matrix observed in our data significantly deviated from that predicted by the SEM26. Finally, we used estimates of standardized coefficients in each path of the model and R2 for each endogenous variable to assess the importance of each set of explanatory variables for determining changes in the AGC dynamics along the latitudinal and elevation gradients in the Andean mountains of South America.
As a complementary analysis, we also used an information-theoretic (IT) approach27,65. The IT approach employs model-based inferences to generate a set of candidate models that represent competing hypotheses including different sets of explanatory variables. The IT model selection uses the relative Kullback-Leibler information to compute the difference in information loss between each model and reality. The relative Kullback-Leibler information is in turn assessed by the Akaike Information Criterion (AIC)65, which identifies the best model as the one that minimizes the overall difference in information between the model and the ecological reality66.
We used a natural model-averaging technique to select the most important explanatory variables of AGC net change, AGC productivity, and AGC mortality. The explanatory variables were associated with our hypotheses (Supplementary Table 1) on the relative importance of climate (PCAs), symbiotic root association (SRA), phylogenetic diversity (PDz), compositional change (TR), initial aboveground carbon (AGC1), and the size-discriminant β parameter as determinants of AGC dynamics. We assessed the nested effect of plots within countries as a means of differentiating geographic subregions, but it was not significant (p > 0.05). The set of models representing competing hypotheses were those models that had a net difference in AIC ≤ 465. The natural model-averaging technique allows us to infer and predict AGC net change, AGC productivity, and AGC mortality based on a set of alternative models. The natural model-averaging technique calculates the average of each variable´s parameter estimates over the models where the variable was selected66,67. The explanatory variables were previously standardized to have mean = 0 and standard deviation = 1, and then each parameter was standardized by the partial standard deviations68. The use of the partial standard deviations68 aims to correct for the likely effect of multiple correlation among variables. The partial standard deviation (s*xi) is calculated as follows:
$$S_{xi}^ ast = S_{xi}sqrt {left( {{mathrm{VIF}}_{xi}^{ – 1}} right)} times sqrt {left( {frac{{left( {n – 1} right)}}{{left( {n – p} right)}}} right)}$$
where Sxi is the standard deviation of the Xi variable; VIF is the variance inflation factor; n is the number of observations; p the number of predictors in the model. The standardized coefficients (βi*) are then transformed as βi* = βi S*xi, where βi is the unstandardized coefficient. All the information-theoretic (IT) analyses were performed with the MuMIn package in R.
We used both SEM and IT independently in a two-phase analysis. In the first phase, we assessed the importance of only the abiotic explanatory variables (PCAtemp1,2 + AGC1) in determining AGC dynamics. By doing so, we allow regional projections based on widely used and available abiotic information. In the second phase, we assessed the importance of both the abiotic and the biotic variables (TR, SRA, PDz, and β) in determining AGC dynamics. Note that we only used the temperature-derived PCA axes (PCAtemp1,2) as climatic explanatory variables in these analyses due to high collinearity with the PCAprec1,2 axes (Supplementary Fig. 6, Supplementary Table 9).
Source: Ecology - nature.com
