Drivers of the Distribution of Soil Organic Carbon Stocks in the Topsoil Under Pinus Hartwegii Lindl. Forest, Along an Elevation Gradient in Central Mexico

Background: Mountain forest soils ( ≥ 2,500 m a.s.l. [1] ), where elevation is crucial in the ecosystem dynamic, have a great capacity to capture and preserve carbon for a long time. The aim of this research was to determine the role of elevation combined with soil, climate, and vegetation variables on soil organic carbon (SOC) stocks distribution under P. hartwegii in the Nevado de Toluca Volcano, a protected area in Mexico. Topsoil samples (0-15 cm depth) collected every 100 m in elevation (3,400-4,000 m) were chemically and physically examined in statistical analysis together with vegetation structure and climate variables. Derived from eld forest conditions, elevation plots were additional analyzed as logged (3,400-3,800 m) vs no-logged (3,900-4,000 m). Results: SOC stocks followed a signicant linear trend (r 2 = 0.70; p= 0.02) along the elevation gradient, being highest at 4,000 m (173.1 ± 5.2 Mg C ha -1 ) and lowest at 3,700 m (146.8 Mg C ha -1 ). Multiple regression analyses showed that SOM, BD, and mean annual temperature (MAT) were the main abiotic drivers of SOC stocks variability (94.5 %) along the elevation gradient. Meanwhile, from the logistic multiple regression, higher tree shrubs and herbs density, in addition to lower tree height and grass cover at lower elevations, indicate a signicant effect of logging on soil traits and vegetation structure, depending on the elevation plot. Conclusions: This research evaluated the SOC stocks and the potential effect of current warming over mountain soils, using the elevation as a proxy for both environmental and human drivers of SOC stocks in mountain forest ecosystems. Higher SOM and grass cover, larger-diameter trees together with low temperatures and logging restrictions in the high elevation range suggest a slower decomposition rate and SOC long-term stability. Despite the fact that we do not know about the intensity and cutting cycle, nor the size of the clearings derived from the trees removed, results from this research show how logging could exacerbate these effects and diminish SOC stocks, as well as the capacity of mountain soils to mitigate the effects of climate change. m (Perry 1991), and a decrease of tree height and diameter as elevation increases towards the treeline ecotone (from 3,980 to 4,130 m) has been shown at the northwest and east-southeast side of the volcanic cone of the Nevado de Toluca, as per the ndings of Alfaro-Ramírez et al. (2017). Although this last study was performed on a different side of the Nevado de Toluca (northeastern side) where the ecotone occurs at a much higher elevation, what was reported by Alfaro-Ramírez and collaborators, is an indicator of the limiting conditions to the P. hartwegii growth as elevation increases. Thus, the fact that we found the largest trees and highest SOC stocks at the high elevation plots, and not as expected at 3,700 m, indicates that other external factors have modied the P. hartwegii forest structure as well as the SOM input and the SOC residence time. Change in forest


Background
Mountain forests (≥ 2,500 m) account for 23 % (9,000,000 km 2 ) of the Earth's forest cover (FAO 2011). These forests have been widely recognized for their role in regulating climate and water cycles, but little is known about the role that low temperatures could have in the carbon (C) cycle. At a global scale, forests are the major component of the C cycle on Earth, storing 1,146 Pg of C, with as much as two-thirds of this C in soils (Dixon et al. 1994). Mountain soils are expected to store more C than soils in low-elevation forests (< 1,000 m) (Spracklen and Righelato 2014; Swetnam et al. 2017), because lower temperatures promote C accumulation through lower decomposition rates (Wiesmeier et al. 2019). However, since mountain forests are considered among the most sensitive terrestrial ecosystems to climatic change (IPCC 2014; Pepin et al. 2015), there is great concern in potential C stock changes in mountain soils as temperature increases (Kirschbaum 2000;Heimann and Reichstein 2008), particularly in the tropics, where warming rates may be more severe at higher elevations (Beniston et al. 1997; IPCC 2014).
The rise of temperature (regional driver) could increase both SOC inputs and outputs in mountain ecosystems (Davidson and Janssens 2006;Wiesmeier et al. 2019) based on the two following opposite theories (Kirschbaum 2000;Berg and Meentemeyer 2002). The rst hypothesis states that an increase in soil temperature will result in higher soil CO 2 emission to the atmosphere (positive feedback) through accelerated SOM decomposition. The second hypothesis assumes an increase in soil nitrogen (N) with faster SOM decomposition, resulting in greater plant productivity with higher litter production and C input to soil, rendering a negative feedback to the atmosphere. However, there is uncertainty on which will be the dominant underlying trend linked to warming, due to the great spatial heterogeneity in mountain forests and the possible in uence of other factors.
Local effects of climatic change on SOC balance might depend on both biotic (plant species) and abiotic (soil type, soil depth, slope, elevation, etc.) conditions, as well as the interaction of both (Salome et al. 2010; Kumar et al. 2012;Bohra et al. 2014;Dymond et al. 2016). Some of these factors might play an important role in mediating SOC stability, i.e., SOC decomposition over the long term (Six et al. 2002;Tian et al. 2016). SOC stability is given through physical protection, as well as biological and chemical reactions. Physical protection of SOC is given by soil structures through aggregate formations, as well as by soil biota such as fungi, bacteria, roots, etc. (Six et al. 2002;Six et al. 2004;Wiesmeier et al. 2019).
Potential impacts of environmental warming on SOC stocks have been determined by using elevation gradients as a surrogate for temperature gradients, since they provide more realistic eld conditions than manipulative experiments in laboratories (Körner 2007;Kumar et al. 2012;Tito et al. 2020). Previous studies using elevation gradients have identi ed different factors as main SOC drivers, including thermal gradients, plant communities, and other local factors For example, Tashi et al. (2016) attributed the direct relationship between SOC and elevation to the reduction in temperature and increase of precipitation, emphasizing temperature as the main driving factor of SOC accumulation along the gradient.
In contrast, Tewksbury and Miegroet (2007) found no interrelation between SOC and elevation, even though C accumulation in the forest oor was more abundant at higher (cooler) elevations, where C turnover was lowest. On the other hand, Shedayi et al. (2016) maintain that elevation controls SOC stocks through vegetation composition and structure, all linked to temperature reduction and increased elevation. In addition, other factors that modi ed vegetation structure such as harvesting and selective logging have showed to have a stronger effect in low forestlands than in higher terrains (Jafari et al. 2013). Neglect effect has also been reported to exert an effect over SOC stocks, despite severely diminishing OC in the topmost centimeters of soil (Nave et  More than 50 % of Mexico's continental territory rises above 1,000 m (Challenger 1998), of which approximately 3.5 % are Andosols (unaged volcanic soils) which contribute signi cantly to the mitigation of regional CO 2 , with high OM content and high SOC sequestration potential (Neall 2006). However, the relationships between SOC stocks and elevation-temperature have not been thoroughly studied in the Mexican mountains where Pinus hartwegii Lindl. is distributed. This mountain pine, with a wide elevation transect (3,500-4,000 m) for a single dominant tree species, provides an unrivaled opportunity for ecological studies along a wide elevation gradient. Under the climate change context, P. hartwegii has been reported as one of the most vulnerable species to the warming projections of climatic change models (Gómez-Mendoza and Arriaga 2007).
A reduction in the environmental suitability of the species has been reported, due to warming coupled with an upward redistribution (Cruz-Cardenas et al. 2016; Manzanilla-Quiñones et al. 2019; Alfaro-Ramírez et al. 2020), in addition to extensive human pressures such as extractive activities (legal or non-legal) and land-use change (Franco et al. 2006). Both events could have an important impact on modifying not only the local, regional, and whole range cover of P. hartwegii forests, but also through many other ecosystem services they provide, such as climate regulation and SOC storage. Therefore, it is important to consider the OC stored in the mountain forest soils in the C national inventories (Santini et al. 2019), and not only the C in the aboveground, as it has been done, as well as the effect of elevation on SOC pools, which would help to know the impact of warming on SOC stocks in these mountain forests.
In this sense, this research aims to evaluate the in uence of elevation on SOC stocks and other soil properties along the elevation distribution gradient of P. hartwegii in the Nevado de Toluca, as well as to determine the main biotic and abiotic variables in uencing SOC stocks in the topsoil of P. hartwegii using elevation as a surrogate for temperature to gain knowledge on the potential effects of warming on SOC stocks in this particularly important model of mountain forest ecosystems.

Study sites
The study was conducted at the Protected Area of Flora and Fauna (PAFF) Nevado de Toluca, an inactive volcano located east of the Transversal Volcanic  García (1990), in the Nevado de Toluca the average annual temperature (MAT) is reported between -2 and 5 °C. In the Tenango del Valle weather station at 2,560 m MAT is 13.23 °C, while at 4,139 m, where the Nevado de Toluca station is located, it is 3.7 °C. Although rainfall can occur throughout the year, the rainy season extends mainly from May to October; July is the wettest month, with rainfall records over 90 mm per day (CONANP 2016). Snowfall is common in winter (December to February) at 3,700 m or higher (Challenger y Soberon 2008).

Field sampling
Randomized strati ed sampling was conducted along an elevation gradient from 3,400 m, the lower limit of P. hartwegii in the area, to 4,000 m. Seven sampling points ("elevation plots") were set every 100 m in elevation along the northeastern face of the Nevado de Toluca Volcano (Fig. 1 C; Table 1).
Elevation plots were selected as homogenous as possible in terms of slope and forest conditions (i.e., no re impacts, pests, or any other type of forest disturbance or tree damage). Accessibility was an important issue because of very steep slopes in some areas along the elevation transect. At each elevation plot, we selected 10 healthy P. hartwegii trees with diameter at breast height (DBH) ranging between 20 and 50 cm, without any sign of physical damage on the trunk (cuts, burn scars, or pests). Two 15 cm-deep soil samples were collected, 50 cm away from the base of the trunk, both above and below the slope of each selected tree, using polyvinyl chloride (PVC) tubes 16 cm in length and 10 cm in diameter. PVC tubes were embedded with the help of two hammers, hitting simultaneously on both sides to insert them homogeneously into the soil. Twenty samples were taken per plot, for a total of 140 soil samples. Additionally, the sampled-tree diameter was identi ed and considered in the analyses as well as the diameter and height of each tree recorded at each plot.

Soil analysis
For each soil sample, the following variables were measured: Bulk Density (BD), Soil Organic Matter (SOM), Soil Organic Carbon (SOC), pH, and texture (reported as clay, silt, and sand content). BD of undisturbed soil samples was determined following the cylinder method. Each soil sample was weighed, after which the soil was taken out from the PVC tubes. Rocks and visible organic material were removed by hand, and each component was weighted. The soil samples were then passed through a 2 mm sieve. SOM content was determined following the oxide-reduction method (Walkley and Black 1934) and subsequently the SOC content was estimated. SOC stocks were obtained using the following equation: SOC stocks = (SOC (%) * BD (g cm -3 ) * Depth (cm) * 100), where Depth was the depth of the extracted sample (15 cm). Results are reported in Mg C ha -1 . Soil pH was measured with 1:2 solutions of 1-N KCl using an OAKTON pH/CON500 potentiometer previously calibrated with pH 4, 7, and 10 buffer solutions. Soil texture was analyzed following the Bouyoucos method (Bouyoucos 1962). The results were analyzed using the United States Department of Agriculture's textural triangle (USDA 2014). In addition, normal climate variables MAT and mean annual precipitation (MAP) were obtained from the ClimateNA model v5. 10 (Wang et al. 2016) for the location of each elevation plot where the soil was sampled, using the site coordinates captured with a GPS device.

Data analysis
In order to identify signi cant differences between the seven elevation plots, a one-way ANOVA was carried out for the soil traits (SOC stocks, SOM, BD, pH, sand, clay, and silt content). Normality tests were previously applied for all continuous variables, con rming adjustment to a normal distribution. Subsequently, linear regression analyses were conducted to identify the in uence of the elevation and its different trends on the SOC stocks and the other soil traits (SOM, BD, pH, sand, clay, and silt) obtained from topsoil samples under P. hartwegii forest along the elevation gradient. The two soil samples taken from each selected tree (above and below the slope) were averaged, using 70 pieces of data for each soil variables in the statistical analysis. The association degree between soil traits, tree height, and DBH of sampled trees was determined by the Pearson correlation coe cient through the data obtained at the elevation plots. Consecutively, in order to describe the variation between trees across elevation plots, a PCA was performed using the 70 pieces of data. To avoid the multicollinearity effect, a variable was removed from the pairs of variables that were highly correlated (r > 0.90). The degree of association between soil, climatic (MAT and MAP), and vegetation structure variables (Table 1) was also determined ( Table 2). In order to identify the main variables in uencing SOC stocks along the elevation gradient, a multiple regression analysis was conducted considering SOC stocks as the dependent variable, while independent variables were arranged in two groups. Group one included mean values of SOM, BD, pH, sand, clay, silt, and DBH of sampled-trees at each sampling point (i.e., the 10-trees each elevation plot). The second group, associated with climatic variables and vegetation structure, included MAT, MAP, trees diameter, tree height, tree density, grass cover, canopy cover, as well as herbs and shrubs density, at each elevation plot.
Finally, due to evidence of timber harvesting on the sampled plots (extraction of trees, pruning, and grassland burning), a dichotomous qualitative variable ("logged" and "not-logged") was used to evaluate logging effects on soil and vegetation throughout a one-way ANOVA and a multiple logistic regression, where soil traits, vegetation variables and climatic variables were associated with this dichotomous qualitative variable. The logged plots were located from 3,400 to 3,800 m in elevation, while the non-logged plots included the 3,900 and 4,000 m elevation plots. All statistical analyses were conducted using the R project version 4.0.3 and the JMP version 8 (SAS Institute, Inc. 2008), with a con dence level of 95 % (α ≤ 0.05), except for Pearson's correlations, which were tested with con dence levels of 99 % (α ≤ 0.01) and 95 %.

SOC stocks and soil properties along the elevation gradient
The elevation plots were signi cantly different for most soil traits. The highest SOC stocks (173.0 ± 5.2 Mg C ha -1 ) were recorded at the highest elevation site (4,000 m), with 16 % less at 3,700 m, where the lowest SOC stocks (146.8 Mg C ha -1 ) were found. In the same way as SOC, the highest amount of SOM (5.13 + 0.06 %) was also found at the highest elevation. The lowest pH was registered at 3,900 m (4.92 + 0.11). Sand content was relatively similar from 3,400 m to 3,700 m, but increased from 3,700 m to 3,800-4,000 m. The clay fraction was lowest at 3,900 m (13.7 + 0.97 %) and highest at 3,400 m (18.7 + 0.74 %); while silt content was lowest to 3,800 m (23.12 + 1.15) and the highest content (26.04 + 0.38 %) at 3,400 m. SOC stocks in topsoil (0-15 cm) under P. hartwegii forest along the elevation gradient were signi cantly related to elevation (r 2 = 0.70, p = 0.02; Fig. 3), i.e., elevation explains 70 % of the spatial variability in SOC stocks with a con dence level of 95 % (Fig. 3). SOM displayed an elevation trend (r 2 = 0.55, p = 0.01) consistent with SOC stocks (Fig. 4 A), except sand content. All other soil variables measured in this study (BD, pH, clay, and silt) did not show any signi cant elevation trend (Fig. 4 B The rst principal component from the PCA explains 38.56 % of the total variability and is made up of tree height and sand, then BD and SOC stocks (Fig. 5 A). The second principal component explains 17.06 % of the total variability and is mostly clay and silt, with pH being less involved (Fig. 5 A). Consecutively, trees located at higher elevations (3,800; 3,900 and 4,000 m) showed greater dispersion, grouping towards PC1, while at the lower elevations (3,400; 3,500 and 3,600 m) trees were less dispersed and clustered towards PC2. Finally, in the center of the gradient (3,700 m) trees were clustered together but separated from the rest (Fig. 5 B).
Relationship between elevation, soil properties, climate, and vegetation structure Considering all variables measured, elevation was signi cantly and positively associated with SOC stocks (r = 0.84, p = 0.04, Table 2 Finally, by associating climate, soil traits and vegetation structure variables with the dichotomous variable of "logging", it can be observed that it presented a statistically signi cant relationship with elevation (p = 0.000, α = 0.05; Fig. 6 A) and MAT (p = 0.000; Fig.6 J). SOC showed a statistically signi cant relationship (Fig. 6 B) with the management of greater stocks at the "no logged" sites, and consistent relationship between logging -SOM (p = 0.000; Fig. 6 C), logging -sand (p = 0.000; Fig. 6 F), and logging -DBH (p = 0.008; Fig. 6 I). The grass cover and the trees height (Fig. 6 N) were higher in not-logging sites, that is, in the highest elevation sites (p = 0.000).

Discussion
Mountain forests and their soils provide multiple bene ts to the inhabitants of both high and lowlands, including the sustainability of the forest itself. However, today these areas are exposed to the impact of climate warming and selective logging. These, in uence the mountain forests capacity to climate warming mitigation by C capture and storage, and these impacts might depend on the elevation at which these forest occur. In this study, SOC stocks in topsoil under P. hartwegii forests increases with elevation, while in the lower elevations the explanation of lowest SOC stocks appears to be confounding, due to environmental and logging factors. Nonetheless, despite we did not measure logging as such and the fact that we do not know the intensity and cutting cycles, results in this study indicate an effect of long-term selective logging not only on the forest structure but also on SOC stocks. Thus, logging has generated not only a greater structural change in P. hatwegii stands between 3,400-3,800 m, which show to be younger tree stands at earlier successional stages (higher density of trees, herbs and shrubs, as well as lower tree diameter and height), but also the SOC stocks decrease in an elevation where the species climax should result in higher accumulation of SOC.
Effect of elevation on soil traits SOC stocks in topsoil under P. hartwegii forests increase with elevation along the 600 m gradient at the Nevado de Toluca; consistent with this, SOM content also raised as elevation increased, indicating, therefore, that this variable as the most signi cant in describing SOC stock trends along the elevation gradient. This association is expected because SOM is the main pool of SOC, which depends on two main drivers. One of them is precisely the input of OM to the soil through litterfall production, which is the dominant way of OM transfer from aboveground to the soil (Six et al. 2004 2016) concluded that SOC residence time (stability of SOC) does not depend on temperature only but on the labile size of SOC. Therefore, the SOC residence time in P. hartwegii mountain forests depends not only on temperature increase by climate change but also on the proportion of labile C in this mountain soils rather than the total SOC. BD, the second variable explaining most of the SOC stocks variability along the elevation gradient in this study, increased signi cantly with elevation. BD is typically low (>0.9 g cm -3 ) in volcanic soils such as Andosols due to the dominance of largely weathered volcanic glass in the sand and silt fractions (Delmelle et al. 2015). However, BD is not only related to the density and arrangement of mineral soil (sand, silt, and clay), but also to the OM particles and SOC (Neall et al. 2006;Delmelle et al. 2015). Thus, the role of soil texture in both SOC content and BD is widely known; for example, sandy soils have shown higher BD and lower SOC, while low BD soils tend to have higher OC (Lukac and Godbold 2011). In this study, sand and silt fraction were positively related to elevation, and negatively related to MAT (Table 2). An increase in sand content as elevation increases is expected in mountain volcanic soils (resulting from volcanic ash) where rock weathering and soil formation processes are less intense, resulting in shallower, coarser, and sandier soils at higher elevations (Simon et al. 2018). . Therefore, it would be important to explore the role of the sand fraction on potential mechanisms of SOC stabilization in these relatively young volcanic ash soils.

Main drivers of SOC stocks on mountain forest soils from the Nevado de Toluca
In this study, the fact that the largest trees were found in the highest elevation plot, consistent with the largest amount of SOM and SOC, suggests litter as the main reason why there were greater SOC stocks recorded at the highest elevation, i.e, there is greater litter production at the highest elevation. However, it is well known that plant growth is reduced at higher elevations due to environmental limitations, such as low partial CO 2 pressure, short growing seasons, soil nutrients limitations, extremely low temperatures, etc., as reported for P. hartwegii (Harsh and Bader 2011; Körner 2015; Alfaro-Ramírez et al. 2017). For this pine species, the climax has been reported at 3,700 m (Perry 1991), and a decrease of tree height and diameter as elevation increases towards the treeline ecotone (from 3,980 to 4,130 m) has been shown at the northwest and east-southeast side of the volcanic cone of the Nevado de Toluca, as per the ndings of Alfaro-Ramírez et al. (2017). Although this last study was performed on a different side of the Nevado de Toluca (northeastern side) where the ecotone occurs at a much higher elevation, what was reported by Alfaro-Ramírez and collaborators, is an indicator of the limiting conditions to the P. hartwegii growth as elevation increases. Thus, the fact that we found the largest trees and highest SOC stocks at the high elevation plots, and not as expected at 3,700 m, indicates that other external factors have modi ed the P. hartwegii forest structure as well as the SOM input and the SOC residence time. Change in forest structure is also supported by the highest herbs and shrub density, and the lowest grass cover, not typical of P. hartwegii forests at low elevations. According to Challenger (1998), typical P. hartwegii forests have a relatively simple structure, with only one tree layer and few or no herbs and shrubs in the understory. Jafari et al. (2013), report similar results with respect to that factors such as harvesting and selective logging, that modi ed vegetation structure, have a stronger effect in lower elevation forestlands than in higher.
In general, changes in the plant community due to natural or anthropogenic causes have a major impact on SOC stocks (de la Cruz-Amo et al. 2020). These mountain forests have been under big anthropogenic pressure for a long time. Timber from P. hartwegii has been in great demand for commercialization of furniture's and construction wood (most of them unauthorized), the largest and highest trees (>35 cm and >20 m, respectively) being the most utilized (Franco et al. 2006). On the other hand, illegal logging occurs in a much lesser degree on the upper area, the buffer zone of the PAFF of the Nevado de Toluca. The logging restrictions allow the permanence of large trees, contrary to the low elevations (up to 3,700 m) where we found the smaller trees and where the forest is more degraded because of the different human activities. Thus, the lack of protection and implementation of the law means that there are no mechanisms to control illegal logging and its impact on the capacity of these forests and their soils to store OC. SOC losses of about 8 % due to logging have been reported in temperate forests (Nave et al. 2010). Logging effects on SOC stocks topsoil can be direct and indirect; tree removal directly reduces OM input to the soil, and indirectly affects decomposition rate, modifying litter quality and microsite conditions ( In summary, SOC stocks along the elevation gradient are controlled through the elevation effect on MAT and soil texture. The results suggest that SOC stocks in high elevation soils occur due to lower temperatures that allow long-term SOC accumulation which is simultaneously promoted through microsite conservation by restrictions on forest exploitation; conversely, at lower altitudes, changes in forest structure strongly reduce SOC stocks. Based on this study, soils under mountain forests not only should be considered in the carbon inventories but also in conservation and differential reforestation programs according to the elevation gradient. Elevation and forest management activities, as important drivers of SOC stocks, could make it possible to keep the greatest SOC reservoirs and long-term mountain forest sustainability.

Conclusions
Elevation had a signi cant positive effect on the distribution of topsoil SOC along P. hartwegii elevation gradient at the Nevado de Toluca. A positive elevation effect was also observed for SOM pools, sand content, and tree diameter, while MAT was found to decrease with elevation. Spatial variability of SOC stocks along the elevation gradient appears to be mainly controlled by SOM, BD, and MAT. However, our results also suggest that this control differs according to the elevation range. In this manner, at low elevations (from 3,500 to 3,800 m) selective timber extraction may have modi ed microsite environment, increasing SOC outputs, while at high elevations (from 3,900 to 4,000 m) lower MAT reduced the decomposition rate and ultimately decreased SOC stocks. Nonetheless, it is important to deepen our knowledge about the control of SOM and SOC dynamics along elevation gradients, as well as to elucidate the role of temperature and other mechanisms in SOC stocks regulation and stabilization in mountain forest topsoil, including tree extraction. In particular, additional research is needed to understand how global warming could impact SOC stocks in alpine forest ecosystems and how unplanned and/or excessive logging could exacerbate its effect.