Soil represents the most important terrestrial carbon sink on Earth. Using data from 9 forest sites located across an elevational gradient (range: 167–2928 m a.s.l.) in northwestern Colombia, which spans a wide range of climate, forest, and soil types, we aim to answer the following research questions: (i) How do SOC stocks change along the elevational gradient or between lowlands and highlands? (ii) What are the main drivers of SOC stocks along elevational gradients? Using structural equation modeling (SEM), we found that climate and soil fertility explained 67% of the variation in SOC stocks across the entire dataset, with SOC stocks declining under warmer and wetter conditions, but increasing with higher soil N:P ratios. Because soil mineralogy was closely correlated with elevation (all lowland were kaolinites while all highland sites were allophanes), SEM models for lowland and highland forests were also run separately. Lowland forests were dominated by trees associated with arbuscular mycorrhizas and nitrogen-fixing symbiont root associations, which could also increase decomposition rates, and thus, reduce the SOC stocks. This finding suggests that greater soil P availability stimulated microbial activity and decreased SOC. In highland forests, which had a wider range of climatic conditions and a greater proportion of trees with ectomycorrhizal associations, decreased temperatures as well as ectomycorrhizal modification of soil N:P ratios slowed soil C cycling, resulting in a greater accumulation of SOC. In conclusion, the increase in SOC stocks across either lowlands or highlands in the Northern Andes was driven by different combinations of abiotic and biotic factors. Since increases in temperature are expected to modify forest functioning and composition along elevational gradients, which in turn depends on soil conditions, improving our understanding on the likely fate of the large amount of C stored in soils should be seen as a priority in tropical montane ecosystems.

In each of the 9 study sites, a 100 m × 100 m (1-ha) plot was established in areas covered by natural forests. Each plot was divided into 25 contiguous subplots of 20 m × 20 m (0.04 ha). To assess soil organic carbon (SOC) in each 1-ha plot, we used a systematic sampling design that selected 15 subplots distributed in three transects of 20 m × 100 m each, separated by 20 m between them (Figure S1). In each subplot, we took a sample centrally located with a soil auger. For each soil core, we removed the organic litter layer and extracted a soil sample centered around 10 cm depth, which represents the first 20 cm of soil. Soil samples were taken using a 4 cm wide and 5 cm long (62.83 cm³) cylinder. The samples were preserved in paper bags and returned to the laboratory. One sample from one plot was missed, meaning that only 134 samples were used for the analyses.

Soil cores were oven-dried at 105 °C to a constant weight. The samples were then crushed and sieved to 2 mm to remove all roots and stones and measured again. A 5-10 g subsample from every core was ground and both C and N concentrations determined after dry combustion with a FLASH 2000 Elemental Analyzer. P concentration was determined by employing spectrophotometry based on a Bray II /L-Ascorbic acid dilution. All soil samples were determined to be carbonate free. Analyses were made at the National University of Colombia at Medellin. Bulk soil densities (soil particles < 2 mm) were quantified using the stone-free dry weights and the sampling cylindric volumes corrected for the stone volumes and compaction (g cm^-3^). Soil organic carbon (SOC in Mg ha^-1^; see Supplementary Information), N and P stocks were then calculated according to the soil densities. Then, we obtained SOC, N and P stock values for each subplot and the N:P ratio calculated as an indicator of soil nutrient availability.

Mean annual temperature (MAT °C) and Annual Precipitation (AP mm) were extracted from the CHELSA database employing the geographic coordinates of each plot. The CHELSA database is a topographically weighted interpolation of weather stations. The climatic variable values are the same for all subplots because of the coarse resolution of CHELSA database. As expected, MAT was highly correlated with elevation (*r *= 0.99; *P *< 0.001). In contrast, AP was unrelated to elevation (r= 0.02; P= 0.76).

All stems tallied in each subplot, when possible, were assigned to the following symbiont root associations (SRA): ectomycorrhizas (EM), arbuscular mycorrhizas (AM) and nitrogen-fixation (NF), based on the genus -or family- level designations provided in Steidinger et al., (2019). Overall, the number of genera assigned was >90% (227/249).