Mature Andean forests as globally important carbon sinks and future carbon refuges

Tropical and subtropical forests sequester roughly 70% of the carbon taken up by the world’s forests, but most estimates are based on lowland ecosystems. Extrapolating lowland-based findings to montane forests is problematic due to greater environmental heterogeneity, steep climatic gradients, and distinct disturbance histories in mountains. The Andes, the world’s longest mountain range and a biodiversity hotspot, have sparse and poorly quantified estimates of carbon stocks at regional scales. Extensive human occupation and ongoing disturbances, coupled with global warming, are driving compositional changes (thermophilization) via upslope range shifts, potentially increasing mortality of large trees at the warm edges of species’ ranges and causing biomass losses that may not be immediately offset by recruitment and growth of thermophilic species. Biotic factors such as symbiotic root associations and phylogenetic diversity may also modulate aboveground carbon (AGC) stocks and productivity, but their roles are underexplored in Andean forests. The study synthesizes data from 119 long-term plots across five Andean countries to: (1) determine whether Andean forests are current AGC sinks or sources; (2) identify dominant abiotic and biotic drivers of AGC productivity, mortality, and net change; and (3) estimate changes in total AGC stocks across subtropical and tropical Andes over the past two decades.

Prior work has emphasized the large contribution of tropical forests to global carbon sinks, but focused primarily on lowland Amazonia and Africa. Montane forests have been considered to show declining AGC with elevation due to colder, harsher climates, yet emerging studies suggest complex patterns driven by both abiotic and biotic factors. Belowground symbiotic root associations (arbuscular vs ectomycorrhizal) are increasingly linked to forest dynamics and soil carbon at large scales, but Andean assessments have relied on coarse global models. Evolutionary history and phylogenetic diversity can influence AGC via niche complementarity or selection effects, with evidence of large-statured clades dominating biomass in certain elevations. Thermophilization of tree communities due to warming has been documented across the Andes, raising concerns about increased mortality at mid-elevations. IPCC default rates for montane carbon dynamics are uncertain, being derived from smaller datasets. Overall, there is a need to integrate climate, disturbance, biotic interactions, and evolutionary context to understand Andean carbon dynamics.

Study area and plots: 119 mature forest inventory plots (73 tropical, 46 subtropical) spanning 7.1°N to 27.8°S, 79.5°W to 63.8°W, and 500–3511 m asl. Mean annual temperature 7.3–23.8 °C and precipitation 608–4313 mm y−1. Plots ranged 0.32–1.28 ha (79% ≥1 ha), total sampled area 104.4 ha, >63,000 trees (DBH ≥10 cm). Censused at least twice between 1991 and 2017 (intervals 2–9 years). Secondary forests and plantations were excluded. Measurements and taxonomy: All trees and palms (DBH ≥10 cm) tagged, mapped, measured; DBH growth, recruitment, mortality recorded at re-census. Species names standardized via TNRS v3.0, with manual curation. Aboveground biomass and carbon (AGB/AGC): Tree AGB estimated using Chave et al. (2014) allometry: AGB = 0.0673 × (WD × DBH² × H)^0.97, where WD from literature at species/genus/family or plot average for unidentified, H estimated from local H:DBH models. AGB converted to AGC using 1 kg AGB = 0.456 kg C. AGC per ha computed by summing trees and scaling to 1 ha. Tree height modeling: Height measured on 44,442 trees; plot-level H:DBH models selected among four candidates (two log-log polynomials, Weibull, Michaelis–Menten) using lowest RSE and bias (BIOMASS R package). Where height data absent (32 Argentine plots), a country-level model was used. Pantropical/country models can bias heights along elevation; plot-level models preferred when possible. AGC dynamics components: For each inter-census interval, computed: (i) AGC mortality = sum of AGC of individuals that died divided by elapsed time; (ii) AGC recruitment = sum of AGC of recruits into DBH ≥10 cm minus AGC of a 9.99 cm stem, per year; (iii) AGC growth = sum of AGC increment of surviving stems per year; (iv) AGC net change = (AGC_final − AGC_initial)/years. Analyses focus on aboveground carbon (exclude soils and belowground tissues). Climate data: Extracted CHELSA bioclim variables (30 arc-sec) for temperature and precipitation. Performed PCA separately on temperature (11 vars; PCA_temp1, PCA_temp2 retained) and precipitation (8 vars; PCA_prec1, PCA_prec2 retained). Due to high collinearity, subsequent models used only PCA_temp1 (elevation-related) and PCA_temp2 (latitude-related). Biotic predictors: Thermophilization Rate (TR; °C y−1) estimated from GBIF occurrence-derived species thermal optima (CHELSA Bio1 at collection sites), computing Community Temperature Index (CTI) weighted by species basal area for each census; TR is annualized CTI change. Phylogenetic diversity (PD) computed from a dated phylomatic tree; standardized effect size PDz via null models (independent swap, 999 randomizations). Symbiotic root associations (SRA) assigned as arbuscular (AM) or ectomycorrhizal (EcM) at genus/family, restricted to American taxa, with manual checks; SRA variable is ln(AM/EcM) weighted by stem number. Disturbance/size-dependent mortality parameter β from logistic regression of death probability vs DBH (logit(P) = α + β×DBH); β<0 indicates higher small-tree mortality (post-disturbance competitive thinning), β>0 indicates higher large-tree mortality (active disturbances). Mean quadratic diameter (Dq) and stem density changes used to characterize thinning vs disturbance. Statistical analyses: Generalized additive models (GAMs) examined latitudinal and elevational patterns of AGC stocks and dynamics; ANOVA tested country differences. Drivers of AGC dynamics analyzed via Structural Equation Modeling (SEM) to quantify direct/indirect effects of climate (PCA_temp1,2), AGC_initial (AGC1), TR, SRA, PDz, and β on AGC net change, productivity, and mortality; variables standardized; Satorra–Bentler scaled chi-square used for model fit; R² reported for endogenous variables. Complementary Information-Theoretic (IT) model selection used multimodel inference (AIC, ΔAIC ≤4) with natural model averaging; variables standardized using partial SDs to adjust for multicollinearity. Two-phase approach: (i) abiotic only (PCA_temp1,2 + AGC1) and (ii) abiotic + biotic (add TR, SRA, PDz, β). Forest cover and AGC stock change: Forest cover and loss estimated for Andes (11°N–27.3°S; 82°W–56°W) from Hansen et al. v1.6 via Google Earth Engine, restricting to forest cover ≥70%, MAP ≥700 mm (WorldClim 2) and ecoregion masks. Time frame aligned to plot census medians: initial 2003, final 2014. Summarized by elevational bands: 500–1200, 1200–2000, 2000–2800, 2800–3500 m. Mean AGC stocks per band derived from plot initial (≤2009) and final (>2009) censuses; bootstrapped 95% CIs. Total AGC stocks computed as mean AGC × forest area per year; overall weighted mean AGC across bands; net AGC balance = (Total 2014 − Total 2003)/11 years. Data and code: plot data available via Dryad; code upon request.

• Andean forests are strong aboveground carbon (AGC) sinks. Mean AGC net change across 119 plots: 0.67 ± 0.08 Mg C ha−1 y−1 (equivalent AGB net change 1.44 ± 0.18 Mg ha−1 y−1), a proportional increase of 1.01 ± 0.13% y−1.
• Drivers of AGC net change: With abiotic variables only, PCA_temp2 (latitude-related axis; higher mean temperatures, lower annual temperature range, higher precipitation and lower rainfall seasonality) significantly explained AGC net change (R²_SEM = 0.18). With abiotic + biotic variables, the size-dependent mortality parameter β was consistently selected (SEM and IT), and AGC1 was significant in SEM (R²_SEM = 0.32). Elevation per se did not significantly predict AGC net change.
• AGC productivity: Abiotic-only models identified AGC1 and PCA_temp2 as significant predictors (R²_SEM = 0.46). With biotic variables, productivity was negatively associated with phylogenetic diversity (PDz) and positively with symbiotic root associations (higher ln(AM/EcM)), yielding R²_SEM = 0.50.
• AGC mortality: Not significantly correlated with latitude/elevation when assessed across all plots. Abiotic-only models selected AGC1 (R²_SEM = 0.19). Adding biotic variables, β was the strongest and most consistent driver (R²_SEM = 0.47). IT models also selected TR and PDz; SEM indicated SRA as a direct significant driver (TR indirectly via β). When β was excluded, TR became a key predictor of mortality (IT).
• Post-disturbance recovery contribution: Focusing on plots with minimal signal of competitive thinning (β between 0.25 and 1.0 quartiles), AGC net change remained positive at 0.48 ± 0.09 Mg C ha−1 y−1, suggesting ≈30% of the overall Andean AGC sink relates to recovery from past disturbance; the remainder likely reflects CO₂ fertilization and warming effects.
• Spatial patterns: AGC stocks decrease with elevation at both plot and landscape scales; nearly half of total Andean AGC (1.71 Pg C) resides in foothills (500–1200 m asl).
• Continental carbon balance 2003–2014: Weighted mean rate of AGC change 0.66 (0.34–0.96) Mg C ha−1 y−1, consistent with plot-based estimate. Total AGC stock increased from 3.83 Pg C (3.31–4.34) in 2003 to 4.12 Pg C (3.43–4.81) in 2014, implying an AGC sink of 0.027 Pg C y−1 (0.011–0.042) despite a 4.2% forest cover loss (12,687 km²) releasing 0.33 Pg CO₂-eq (0.28–0.37).
• Comparative strength: The Andean AGC sink (0.67 Mg C ha−1 y−1; 1.01% annually) is stronger than estimates for lowland Amazonian forests over a similar period (~0.42 Mg C ha−1 y−1).
• Biotic mechanisms: Productivity is associated with selection effects of large-statured clades differing by elevation (e.g., Fabaceae, Lauraceae, Moraceae at <2000 m; Cunoniaceae, Melastomataceae, Clusiaceae >2800 m), negative relationship between PDz and productivity, and positive association with AM-dominated SRA. Thermophilization (positive TR) is linked to increased mortality of large, lower wood density species at mid-elevations (1000–1800 m), potentially transient as thermophilic species migrate upslope.
• Policy implication: Reducing deforestation, especially at 500–1800 m, and enhancing restoration and connectivity will increase AGC stocks and facilitate upslope species migrations.

The study demonstrates that mature Andean forests are globally significant carbon sinks with a strong net AGC gain that not only matches but exceeds lowland tropical sink rates in some regions. By integrating climate gradients with biotic factors and disturbance signals, the analyses show that size-dependent mortality (β) and associated post-disturbance recovery are central to explaining current AGC gains, while climate variability along latitude (PCA_temp2) modulates both net change and productivity. The mortality signal linked to thermophilization suggests that warming drives increased losses where species approach the hot edge of their thermal niches, particularly between 1000–1800 m asl; however, these effects may be transient as thermophilic, large-statured species colonize higher elevations given sufficient connectivity. The negative association between productivity and phylogenetic diversity, alongside the dominance of a few productive clades at different elevations, supports selection effects over complementarity in these montane systems. The positive role of AM-type symbioses in productivity underscores the importance of nutrient acquisition strategies in cold, nutrient-limited montane soils, though such advantages may diminish under future CO₂-enriched, N-limited conditions. At the regional scale, increasing AGC stocks in remaining forests have more than offset emissions from deforestation over 2003–2014, highlighting the Andes as emerging carbon refuges as lowland sinks weaken. The findings emphasize that conserving mid-elevation forests, curbing deforestation, and maintaining elevational connectivity are essential not only for biodiversity but also for sustaining and enhancing carbon storage under climate change.

Mature subtropical and tropical Andean forests currently function as strong aboveground carbon sinks (~0.67 Mg C ha−1 y−1; ~1% per year), with total AGC stocks increasing from 3.83 to 4.12 Pg C between 2003 and 2014 despite measurable forest loss. Size-dependent mortality dynamics (β) linked to disturbance and post-disturbance recovery, climate gradients (PCA_temp2), symbiotic root associations, and evolutionary composition (PDz and dominant clades) jointly govern AGC productivity, mortality, and net change. Approximately one-third of the sink strength is attributable to recovery from past disturbances, with the remainder likely reflecting CO₂ fertilization and warming. Thermophilization currently elevates mortality at mid-elevations but may be transient if connectivity enables upslope migration of thermophilic, productive species and their symbionts. To secure and enhance the Andes’ role as future carbon refuges, policies should prioritize halting deforestation—especially at 500–1800 m—and scaling restoration to improve connectivity. Future research should expand sampling intensity and spatial coverage, incorporate belowground and soil carbon, explicitly quantify regrowth contributions, refine mycorrhizal mapping and nutrient constraints, and track species and symbiont migrations to improve forecasts of montane carbon dynamics.

• Carbon accounting focuses on aboveground tree carbon (DBH ≥10 cm), excluding soils and belowground biomass; thus total ecosystem carbon dynamics are underestimated.
• Tree height data were unavailable in 32 Argentine plots, necessitating use of country-level H:DBH models that can homogenize AGC estimates along elevation; pantropical/country allometries may bias AGB at high/low elevations.
• Sampling intensity and plot size can bias mortality detection and AGC dynamics estimates; many plots are within forests <100 km² and may be influenced by edge effects.
• Asynchrony in plot census dates required using initial (≤2009) and final (>2009) stocks to estimate band-level changes, potentially introducing temporal mismatch.
• High collinearity between temperature- and precipitation-derived climate axes limited models to temperature PCA axes, potentially omitting independent precipitation effects.
• SRA assignments at genus/family level and regional restriction may introduce classification uncertainty; AM/EcM proportions use stem counts rather than basal area/biomass.
• Thermophilization (TR) estimates rely on GBIF occurrence data and CHELSA-derived thermal optima, which may be biased by sampling effort and spatial errors.
• Forest regrowth areas (~500,000 ha, 2001–2014) were not included in total AGC stock estimates, likely underestimating net regional AGC gains.
• β and AGC mortality are both size-dependent, leading to inherent correlation that complicates causal inference; distinguishing competitive thinning from exogenous disturbance remains challenging.