Soil organic carbon stocks in native forest of Argentina: a useful surrogate for mitigation and conservation planning under climate variability

Background The nationally determined contribution (NDC) presented by Argentina within the framework of the Paris Agreement is aligned with the decisions made in the context of the United Nations Framework Conven‑ tion on Climate Change (UNFCCC) on the reduction of emissions derived from deforestation and forest degradation, as well as forest carbon conservation (REDD+). In addition, climate change constitutes one of the greatest threats to forest biodiversity and ecosystem services. However, the soil organic carbon (SOC) stocks of native forests have not been incorporated into the Forest Reference Emission Levels calculations and for conservation planning under cli‑ mate variability due to a lack of information. The objectives of this study were: (i) to model SOC stocks to 30 cm of native forests at a national scale using climatic, topographic and vegetation as predictor variables, and (ii)


Introduction
Soil organic carbon (SOC) is the main terrestrial C reservoir and contains three times more C than the atmosphere and twice the amount stored in the vegetation (Friedlingsten et al. 2019;Mayer et al. 2020).Increasing C storage in forest ecosystems is widely recognized as having high climate mitigation potential as well as environmental and socio-economic benefits (Lewis et al. 2019;Duarte-Guardia et al. 2020).Also, SOC correlates positively with several soil functions (soil fertility, infiltration, structural stability, porosity) and biodiversity (e.g.plants, insects, birds, and other living things) that support ecosystem services to societies (Mori et al. 2017;Trivedi et al. 2018;Duarte-Guardia et al. 2019;Canedoli et al. 2020;Magnano et al. 2023).
SOC stocks vary as a function of C inputs and outputs resulting from ecosystem factors such as climate, soil, and vegetation, land use change and the mechanisms controlling C residence times (Jobbágy and Jackson 2000; Don et al. 2011;Stockmann et al. 2013;Villarino et al. 2017;Berhongaray and Alvarez 2019;Duarte-Guardia et al. 2020).Different statistical techniques have been applied in the mapping of SOC stocks, including multiple linear regression, linear mixed models, geographically weighted regression, regression-kriging and methods from the machine learning field to map SOC stocks (Vågen and Winowiecki 2013).In Argentina, different simulation models have been used to estimate and project SOC changes.For example, Villarino et al. (2014) used the IPCC (Intergovernmental Panel on Climate Change) empiric Tier 1 and Tier 2 approaches to estimate historical SOC stocks and flows in the Pampa Region and Piñeiro et al. (2006) applied the Century model to simulate historical SOC changes in temperate grasslands.However, there are no antecedents of SOC stocks estimation of native forests in Argentina at the country level due to the lack of field measurements over large forested areas including SOC, bulk density and rock fragment content, and because high spatial variability of SOC requires a very high sampling density to get accurate estimates.
The Working Group I contribution to the IPCC Sixth Assessment Report stated that unless there are deep reductions in global greenhouse gas (GHG) emissions, the goal of limiting warming well below 2 °C and close to 1.5 °C will be out of reach (IPCC 2021).In this context, Argentina has reaffirmed its commitment to the Paris Agreement, extending the country's goal regarding the reduction of GHG by 2030 and announcing a longterm, low-emission development strategy to achieve carbon-neutral development by 2050.Also, the nationally determined contribution (NDC) presented by Argentina within the framework of the Paris Agreement is aligned with the decisions made in the context of the United Nations Framework Convention on Climate Change (UNFCCC) on the reduction of emissions derived from deforestation and forest degradation, as well as forest carbon conservation (REDD+).The Forest Reference Emission Levels (FREL) is one of the pillars of the REDD+ process, which defines a baseline to evaluate the performance of a country in the implementation of REDD+ activities in terms of reducing GHG in the forestry sector.Despite its importance, the SOC stock of native forests has not been incorporated into the FREL calculations.For this, it is necessary to generate accurate estimates of SOC stocks and changes for the major land-use and management conditions.Thus, in the context of international policy agendas on GHG emission mitigation, it becomes relevant to develop reliable tools for SOC stock monitoring at a national scale (Lal 2011).
SOC is the basis for sustain net primary production, which in fact is one of the key functional dimensions of ecosystems to maintain biodiversity.Furthermore, biodiversity is strongly influenced by climatic and environmental variables (Read et al. 2020;Rosas et al. 2023), but also by many other site characteristics, where soils are pointed as a key factor (Khaziev 2011).Recently, a close relationship between SOC stocks and biodiversity was informed (Peri et al. 2019).In other words, changes in abiotic factors as temperature, or radiation, can alter the biogeochemical cycles inducing losses or gains in forest SOC and consequent alteration of NPP, directly affecting biodiversity (Ma et al. 2013;Sardans and Peñuelas 2013;De Frenne et al. 2021;Mina et al. 2021).When the species are not able to adjust to these changes, populations can decline over time due to temporal mismatches (e.g.food and reproduction) (Socolar et al. 2017;Silveira et al. 2021).High spatial variability in phenology (proxy: vegetation greenness) and climate (proxy: soil surface temperature) can reduce the biodiversity loss threat by increasing ecosystem resilience when high inter-annual variability occurs (Nystrom and Folke 2001;Peaucelle et al. 2019).High spatial variability can be related to higher food resources, habitat diversity and climatic refugia, enhancing species survival (Thuiller et al. 2005;Oliver et al. 2010).Inter-annual variability in the vegetation greenness and seasonality of temperature were previously calculated for Argentina (Silveira et al. 2021) using temperature infrared data (Hengl et al. 2012) and time series of satellite-based vegetation indices, as the Enhanced Vegetation Index (EVI) (Hu et al. 2019).
The objectives in this study were: (i) to model SOC stocks to 30 cm of native forests at a national scale using climatic, topographic and vegetation as predictor variables, and (ii) to relate SOC stocks with spatiotemporal satellite-derived indices related to vegetation greenness and climate to determine biodiversity conservation concerns due to threats from high interannual climate variability.We hypothesized that SOC is related to environmental gradients, which together with other factors (e.g.abiotic and biotic) define the potential biodiversity of one particular forest ecosystem.This work aims to assist the Government of Argentina in the quantification of SOC content of the country's native forests to have information in relation to the mitigation measures in the National Forest and Climate Change Action Plan.

Study area and SOC database
The Argentinian territory covers 2.78 million km 2 area from 21° 46ʹ to 55° 03ʹ S latitude and from 53° 38ʹ to 73° 34ʹ W longitude that determines a wide range of ecosystems, where native forest is grouped in six forest regions: Monte forests, Espinal forests, Parque Chaqueño forests, Yungas rainforests, Paranaense rainforests, and Patagonian forests (Tierra del Fuego forests and Andean Patagonian forests) (Appendix 1).Within these regions we masked the native forest lands by using the map of the National Native Forest Inventory (SAyDS 2005) and its update (SAyDS 2019).On this map, we used the current definition of the native forest implemented by the Argentine Government within the framework of Law 26,331/07, which regulates the implementation of institutional policies: "all natural forest ecosystems in different stages of development, of primary or secondary origin, with a tree cover greater than or equal to 20% and a minimum tree height of 3 m, including palm groves" (SAyDS 2019).The study area was the total native forest that covers an area of 0.468 million km 2 .
We used 1040 forest soil samples (0-30 cm depth) taken from 2015 to 2021.The mean soil sampling density was 2.37 samples per 1000 km 2 (Appendix 1).We extracted soil samples (n = 3-9 in randomly selected areas) using a hand soil sampler including the first 30 cm depth below litter layer.The soil sampler has a known volume, making it possible to calculate soil bulk density (SBD).The calculations were conducted with air-dried samples, with any particles > 2 mm removed previously by sieving (roots, stones, coarse woody debris).The SOC stock (kg m −2 ) was calculated according to the equation: where C is the concentration of SOC (kg 100 kg −1 ); BD is bulk density (kg m −3 ), 0.3 is the thickness of the layer (m), CF is the fraction of coarse fragments in the soil sample and 0.01 is a scale factor.

Environmental predictors for SOC estimation
We selected 52 potential predictive environmental covariates, which represent key factors for the spatial distribution of SOC content such as climate, topography, soil and vegetation (McBratney et al. 2003) (Appendix 2).All covariate maps were generated or uploaded to the Google Earth Engine cloud-based computing platform for subsequent modelling.The original covariates' spatial resolution was adjusted to a common resolution of 200 m and then masked the native forest lands of Argentina.

Climatic variables
We included climatic variables (more specifically precipitation and temperature) because climate influences the SOC stock through the control of soil processes such as weathering of parent materials, erosion rate and mineralization/humification of organic matter, but also because it affects the net primary productivity (NPP) that is essential in SOC accumulation (Gaitán et al. 2019).Data on the mean  annual precipitation and monthly mean, minimum and maximum temperature were obtained from the Worldclim global database.Derived from these monthly temperature and precipitation data, we used a group of 17 maps of bioclimatic variables (Hijmans et al. 2005).
Additionally, we included the mean annual  of day and night land surface temperatures, which were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) MOD11A1 V6 product.

Topographic variables
Topography has the potential to control the spatial patterns of soil properties such as SOC and, therefore, models that take topographic variables into account can provide better estimates of SOC stocks (McBratney et al. 2003).This is because topography, through its influence on gravity, solar insolation and micro-climate, determines the water, solute, and sediment fluxes throughout the landscape.A set of 13 morphometric maps (Appendix 2) was generated using the Terrain Analysis in Google Earth Engine (TAGEE, Safanelli et al. 2020).This algorithm uses the digital elevation model (DEM) derived from a radar image from the Shuttle Radar Topography Mission (SRTM, Farr et al. 2007) with a 30 m spatial resolution.

Soil variables
Clay has been shown to actively protect organic matter from decomposition by adsorption and aggregation, slowing turnover and increasing SOC residence times (Paul 1984;Schimel et al. 1985).We also included pH because it controls microbial mechanisms of carbon accumulation (Malik et al. 2018).On the other hand, soil water erosion influences SOC in two ways: redistribution of C within the watershed, and loss of C to the atmosphere (Polyakov and Lal 2004).Therefore, in this study, we used national maps of soil clay content, pH (0-30 cm depth) and an Universal Soil Loss Equation (USLE)-derived map of soil loss rate by water erosion (Gaitán et al. 2017) (Appendix 2).

Vegetation variables
We selected vegetation indices (VIs), based on satellite observations, that are mathematical transformations of reflectance measurements in different spectral bands, especially the visible (usually red) and near-infrared bands, that are widely used to obtain information about vegetation characteristics (Jackson and Huete 1991; Bannari et al. 1995).In this study, we used 10 different VIs and the Blue, Red and NIR bands (Appendix 2) derived from the atmospherically corrected surface reflectance data produced by the collection 2, Tier 1 of the Landsat 8 OLI sensor, available in the GEE platform (GEE collection snippet: ee.ImageCollection of 'LANDSAT/LC08/ C02/T1_L2').To calculate the VIs, we used the CFmask algorithm to mask the cloud and shadow produced, as well as a per-pixel saturation mask (Zhu and Woodcock 2012) and used the mean of the 2015-2021 period of the bands' surface reflectance.Lamichhane et al. (2019) review 120 studies designed to correlate environmental variables with SOC stock by applying different digital soil mapping techniques.Among those methods, Random Forest (RF, Breiman 2001) algorithm performed better than other methods in most comparative studies.Therefore, in our study the RF algorithm was selected to predict and map the SOC stocks in the native forest of Argentina.This algorithm builds a set of regression trees.Each tree predicts the result in each pixel, while the final prediction is obtained by averaging these values (Breiman 2001).Each tree is built from a bootstrap sample of the original data set which allows for robust error estimation with the remaining test set, the so-called Out-Of-Bag (OOB) sample.The excluded OOB samples are predicted from the bootstrap samples and by contrasting predicted versus observed values different metrics can be calculated to assess the performance of the models.For each tree, only a subset of the predictor variables is used, which allows estimation of the variable importance measured by the mean decrease in prediction accuracy before and after permuting a variable.Of our 52 potential covariates, in preliminary simple correlation analyses, several showed a low correlation with stock SOC and a high correlation with each other.Therefore, due to random choices of subsets of variables being made at each split, many noise variables unrelated to the SOC stock may impact model's performance (Bahl et al. 2019).In addition, the presence of correlated predictors impacts RF's ability to identify the strongest predictors by decreasing the estimated importance scores of correlated variables (Gregorutti et al. 2017).

Random forest (RF) modelling for SOC map
To overcome this problem, we applied the Random-Forest-Recursive Feature Elimination (RF-RFE) algorithm (Gregorutti et al. 2017).The RF-RFE selects the subset of covariates that is most predictive for SOC stock by a backward selection algorithm that iteratively eliminates the weakest explanatory variable from the initial 52-covariates model.The performance of the models was evaluated using tenfold cross-validation.For doing this, the SOC dataset was randomly split into ten sets in which 70% of the data were used for calibration of models and 30% for validation.At each step of the RF-RFE algorithm, we calculate the mean of tenfold crossvalidation root-mean-square error (RMSE) and plotted it versus the number of remaining covariates (Appendix 3).Based on this analysis, we finally selected 10 covariates for modelling and mapping SOC stocks (Appendix 2).To evaluate the performance of the final SOC stock model, three indices were calculated and averaged using tenfold cross-validation: root-mean-square deviation (RMSE), mean absolute error (MAE) and the coefficient of determination (R 2 ).Additionally, we computed the lower (Q5) and upper (Q95) percentiles for each pixel, resulting in maps of SOC stocks for lower and upper percentiles.The uncertainty was estimated as the difference between the 5th and 95th percentiles (e.g.90% prediction interval) as proposed by Heuvelink (2014).

SOC and biodiversity conservation concern
To determine the biodiversity threats from high interannual climate variability, we employed the spatio-temporal remotely sensed indices determined by Silveira et al. (2021) for the whole of Argentina based on EVI and land surface temperature (LST) images from Landsat imagery.We used the final products (GRIDS at 250-m resolution): (i) inter-annual variability in phenology based on land surface temperature (IAV-LST); (ii) spatial variability based on land surface temperature (SV-LST); (iii) inter-annual variability in phenology based on vegetation greenness (IAV-PHE); and (iv) spatial variability based on vegetation greenness (SV-PHE).Silveira et al. (2021) assumed that: (i) high inter-annual variability in the seasonality of land surface temperature and phenology of vegetation greenness threatens biodiversity.This high variability can exceed the rate at which organisms are able to adapt, disrupting the synchrony of ecological interactions, and thus, the survival of the species.(ii) High spatial variability in land surface temperature and vegetation greenness enhances resilience to high inter-annual variability.Spatial variability provides a higher variety of resources in close proximity by increasing suitable conditions during times of climate extremes.Biodiversity conservation concerns due to threats from high inter-annual climate variability were considered: (i) low concern when IAV is low and SV is high (L); (ii) medium concern when both IAV and SV are high, or both IAV and SV are low (M); and (iii) great concern when IAV is high and SV is low (H).With these layers, we calculated 3 indices to synthesize the biodiversity conservation concerns due to threats from high inter-annual climate variability: (i) I-LST = integration between IAV-LST and SV-LST, where IAV and LST were classified in high (value = 1 for the area with 50% higher values) and low values (value = − 1 for the area with 50% lower values) for each forest region, and then calculated the index as I-LST = (IAV-LST − SV-LST)/2; (ii) I-PHE = integration between IAV-PHE and SV-PHE, where IAV and LST were classified in high (value = 1 for the area with 50% higher values) and low values (value = − 1 for the area with 50% lower values) for each forest region, and then calculated the index as I-PHE = (IAV-PHE − SV-PHE)/2; and finally (iii) I-GLO = global integration between I-LST and I-PHE, where I-GLO = (I-LST + I-PHE)/2.
We used a mask of native forests and forest regions to generate the final GRIDS: (i) a mask of native forests proposed by Silveira et al. (2022), which included forested areas taller than 5 m in height and 10% in canopy cover, as was defined by Global Forest Change (GFC) data set (Hansen et al. 2013).(ii) Forest cover classified into an adaptation of the proposed forest regions (Peri et al. 2021a) For data extraction, we used the hexagonal binning technique, a spatial methodology that offers the advantage of integrating different pixels (e.g.averaging values for each pixel) within polygonal regions to effectively capture spatial patterns (Battersby et al. 2017).We implemented a hexagonal binning process that involved one spatial matrix dividing the territory of Argentina into hexagonal areas of 5000 ha each (Rosas et al. 2023).Then, the average values of each hexagonal grid were computed.We excluded hexagons that presented less than 10% of native forest cover (e.g.< 500 ha in each hexagon).

Data analyses
Our final database included 21,420 hexagons classified in: (i) forest regions (TDF = 417 hexagons, APF = 1,487 hexagons, MON = 347 hexagons, ESP = 3,347 hexagons, PCH = 13,487 hexagons, YUN = 1,440 hexagons, and PAR = 895 hexagons); (ii) four SOC levels (L = low, M = medium, MH = medium high, H = high) representing quartiles calculated for each forest region; (iii) seven grids related to LST and PHE, and the derived indices.Data were compared through: (i) two-way analysis of variance (ANOVA) with SOC levels and forest regions as fixed factors, analysing IAV-LST, SV-LST, I-LST, IAV-PHE, SV-PHE, I-PHE, I-GLO; and (ii) simple ANOVAs comparing different SOC levels for the same variables, conducted at each forest region individually.We used the Fisher test (F) for the ANOVAs, and the means were compared through the Tukey test at p < 0.05.Finally, we determined one homogeneity index (HI) for each region based on the internal deviation of each SOC level with respect to the average value for each biodiversity conservation concern variable (IAV-LST, SV-LST, IAV-PHE, SV-PHE).The final index was the average of the four variables presented as the percent deviation of the absolutes obtained values.
The RF model explained 69% of the variance in SOC stock (RMSE = 33.99;MAE = 23.07;R 2 = 0.69).The 10 most important variables for predicting SOC stocks across all cases were as follows: day land surface temperature, temperature seasonality (standard deviation × 100), surface reflectance of blue band, mean annual precipitation, maximum temperature of warmest month, monthly maximum temperature, precipitation of warmest quarter, night land surface temperature, soil water erosion rate, and monthly minimum temperature (Fig. 2).Thus, climate variables alone explained most of the variation in SOC stocks.

SOC and biodiversity conservation concern
Biodiversity conservation concern was influenced by the interactions between SOC levels and the forest regions (Table 2).Great SOC contents were related to high values of IAV-LST and SV-LST and resulting in higher I-LST.While high IAV showed a higher biodiversity conservation concern, high SV is associated with higher landscape resilience.However, when both variables were combined, I-LST indicated that high SOC values had higher biodiversity conservation concern.High SOC contents were also related to higher values of IAV-PHE and SV-PHE, and high I-PHE values occurred in the lowest SOC contents.Both IAV and SV presented the same trend (high biodiversity conservation concern and high landscape resilience).However, when both variables were combined, I-PHE indicated that high SOC values showed, in average a low biodiversity conservation concern for PHE.Both indices combined into a single one (I-GLO) did not present significant differences, but had high values (high biodiversity conservation concern and low landscape resilience) associated with high SOC values.
The forest regions presented significant differences in the studied variables (IAV-LST, SV-LST, IAV-PHE, SV-PHE), but not in the indices (Table 2).The lowest values of IAV-LST were found in the driest forest areas (MON and ESP) (lowest threat for biodiversity conservation concern), while the highest values were found in rainforests (YUN and PAR) (highest threat for biodiversity conservation concern).The lowest values of SV-LST were found in the central area of Argentina (ESP and PCH) (lower landscape resilience), while the highest values were found in Patagonian forests (APF) (higher landscape resilience).The lowest values of IAV-PHE were found in the rainforests (YUN and PAR) (lowest threat for biodiversity conservation concern), while the highest values were found in the most austral forests of Argentina (TDF) (highest threat for biodiversity conservation concern).Finally, the lowest values of SV-PHE were found in the driest forest (MON and ESP) and central area of Argentina (PCH) (lower landscape resilience), while the highest values were found in Patagonia and Tierra del Fuego temperate forests (APF and TDF) (higher landscape resilience).These results indicated that the trends for biodiversity conservation concerns varied depending on the studied indicator (LST or PHE).Driest areas (MON, ESP and PCH) showed different trends compared to rainforests (YUN and PAR) when inter-annual variability was considered (LST and PHE).However, the spatial variability of both variables (LST and PHE) had similar trends, where the driest areas (MON, ESP and PCH) presented low landscape resilience compared to southern forests (APF and TDF), which presented the maximum values.
It is important to consider that these analyses showed significant interactions for all the studied factors (Table 2), depending on forest region particularities (see Fig. 3, Appendices 4 and 5).Patagonian forests (TDF and APF) presented high IAV values (high threat for biodiversity conservation concern) at high SOC contents and high SV (enhanced landscape resilience) at middlehigh SOC contents.Dry and central forests of Argentina (MON, ESP and PCH) presented contrasting response in IAV.Generally, higher IAV values (high threat for biodiversity conservation concern) were found at higher SOC contents, except for MON where IAV-LST was associated with low SOC contents and for PCH where IAV-LST was related to middle-high SOC contents.These forests had a similar pattern in SV, where high values (enhanced landscape resilience) were related to high SOC values.Rainforests (YUN and PAR) showed different patterns than other regions: (i) high IAV-LST were related to high SOC contents, but high IAV-PHE was related to low SOC contest, and (ii) high SV was related to high SOC values in YUN and low SOC values in PAR.These relationships can also be observed in the maps of Appendix 4, where the different forest regions presented a significant landscape variation for the studied variables (SOC and

Discussion
The estimation of forest carbon storage at large scales has attracted considerable interests from diverse actors to mitigate climate change and conserve biodiversity (Vågen and Winowiecki 2013;Friedlingsten et al. 2019;Peri et al. 2019;Dante-Guardia et al. 2020).To our knowledge, this work is the first to quantify the SOC stocks in native forests of Argentina at a national scale and to analyse their relationships with multiple predictor variables.Thus, the model for SOC (0-30 cm depth) prediction accounted for 69% of the variation of this soil property across the whole native forest areas in Argentina with values ranging from 2.18 to 17.35 kg m −2 .The SOC stock was mainly a function of climatic variables related mainly to precipitation (mean annual precipitation and precipitation of the warmest quarter) and temperature (diurnal land surface temperature, seasonality of temperature, maximum temperature of the warmest month, monthly maximum temperature, nocturnal land surface temperature and monthly minimum temperature), reflecting the influence of climate on the ecosystem's forests.AREA = native forest cover (km 2 ), MLL = lower limit (5%) of the mean SOC contents, M = mean SOC contents, MUL = upper limit (5%) of the mean SOC contents, TLL = lower limit (5%) of the total SOC stocks, T = total SOC stocks, and TUL = upper limit (5%) of the total SOC stocks  Temperature influenced SOC by increasing the quantity at lower temperatures (Patagonian forests with 14.45 kg m −2 ) compared to the average SOC stock for the entire country (mean of 6.10 kg m −2 ) and northern forest regions.Some models and observations suggest that forests may behave as a C source in response to increased decomposition of soil organic matter resulting from temperature increases (Davidson et al. 2000).Thus, the temperature sensitivity of decomposing organic matter in soil partly determines how much carbon will be transferred to the atmosphere because of global warming.This is consistent with Li et al. (2020) who reported that soil respiration rates were correlated strongly to air and soil temperatures by evaluating seasonal dynamics in contrasting forests across climate gradients.The negative effect of temperature on SOC observed is supported by studies in other types of ecosystems showing soil organic matter contents decrease with increasing temperature (He et al. 2014) due to increased mineralization rates (Kirschbaum 1995).
Rainfall and other related rainfall variables (seasonality) also influenced SOC by increasing the quantity of precipitation.It has been demonstrated that increased variability in rainfall and soil water content significantly affected SOC in forests (Duarte-Guardia et al. 2019).This relationship could be explained by the increase in primary productivity due to higher soil moisture and higher precipitation values and lower soil temperatures at elevated water supply that enhanced SOC storage (Chapin et al. 2002;Adhikari et al. 2014).However, SOC is probably controlled by the complex interaction of environmental and biotic factors.At the national scale, in this study, patterns of SOC were positively associated with mean annual precipitation and negatively correlated with   mean annual temperature across a diverse range of soils and forest types.This has also been documented in grasslands in North America (Burke et al. 1989) and global dry-lands (Gaitán et al. 2019).

Variable
The surface reflectance of the blue band (0.43-0.45 μm) used to improve sensitivity to chlorophyll represents a proxy of vegetation type through estimations of vegetation canopy moisture and greenness indices (Rahaman et al. 2017).Thus, SOC may vary in different native forest types across the country as a result of the modification of both carbon inputs (e.g.net carbon gain by plants) and losses (e.g.microbial decomposition), as was reported previously for different vegetation types (Post et al. 1982;Yang et al. 2007).Rasse et al. (2005) found that vegetation with higher plant diversity leads to greater ANPP and belowground biomass, which enhances belowground carbon inputs and, thereby, SOC storage.In contrast, Chen et al. (2018) reported that the negative effect of ANPP on SOC in forests and shrublands indicates that the carbon losses associated with microbial respiration may offset the positive effect of belowground biomass on SOC storage.Also, SOC values related to variability in vegetation characteristics may be influenced by topography.For example, in Argentinean Nothofagus forests, evergreen species (e.g. less variability in primary productivity) are more common in lowlands, near lakes, and locations with alluvial soils, while upland forests are generally deciduous (e.g. more variability in primary productivity (Peri et al. 2021a).This is consistent with Martínez Pastur et al. (2022) who reported that SOC values were highest in valleys of the Andes Mountains and southern Tierra del Fuego, ranging from 5.35 to 27.78 kg C m −2 for the whole Patagonia region.
Soil erosion rate influence negatively SOC.Erosion influences the vertical and horizontal distribution patterns of SOC in two ways: redistribution of C within the watershed and loss of C to the atmosphere by dispersing the soil, altering microbiological activity and inducing mineralization (Polyakov and Lal 2004;Fiener et al. 2015).The extent of ecosystems' capacities to reduce erosion depends on many factors related to environmental conditions (amount of precipitation, wind velocity, soil properties, slope, and vegetation characteristics) and pressures (forest management practices and overgrazing).In Patagonia, Peri et al. (2021b) determined a SOC loss by erosion from 85.3 to 250.1 kg C ha −1 year −1 .While topsoil formation takes thousands of years, erosion can disintegrate all the organic matter and nutrients in the topsoil in a few years, dramatically reducing soil productivity.This would also explain the higher SOC recorded in the valley or lowland forests than on the slopes of the Andes Mountains and in southern Tierra del Fuego.Furthermore, forest ecosystems contain a large number of terrestrial carbon stocks globally, with more than half allocated in soils (Liu et al. 2018;Li et al. 2019;Martínez Pastur et al. 2022;Mayer et al. 2020).It was proposed that SOC stocks and biodiversity were closely related (Peri et al. 2019).However, other studies did not find a clear relationship between them (Canedoli et al. 2020).Ecosystem services in natural ecosystems are underpinned by biophysical structures and processes (de Groot et al. 2010), and driven by biodiversity (Balvanera et al. 2006) with positive effects (e.g.regulating and supporting services) (Brockerhoff et al. 2017).Biodiversity is also important because contributes to soil formation, thereby contributing to enhance SOC contents (Canedoli et al. 2020).Many studies showed a positive effect of species and functional richness on SOC stocks (Saha et al. 2009;Dawud et al. 2016;Liu et al. 2018;Li et al. 2019;Magnano et al. 2023), e.g. higher plant diversity promotes higher litter accumulation in natural ecosystems (Steinbeiss et al. 2008;Lange et al. 2015).These plant-soil relationships play an important role in soil C accumulation, where high plant richness results in large litter biomass (Steinbeiss et al. 2008;Leff et al. 2012;Liu et al. 2018;Li et al. 2019).A recent study reported using a prioritization analysis that sites that maximize protection of carbon stocks (living aboveground and belowground biomass, litter, dead wood, and active SOC) had high values for biodiversity conservation (all threatened terrestrial mammals, birds, reptiles, and amphibians, with high threat status and endemism), and potential sources of clean water.By conserving 30% of the territory, it is possible to achieve 80% of the conservation targets for biodiversity and water provisioning, and 50% of the conservation targets for carbon stocks (Frank et al. 2023).
Climate change constitutes one of the greatest threats to forest biodiversity and ecosystem services, where model predictions warn about significant impacts for the next decades (Sala 2000;Seidl et al. 2014;Thom et al. 2019).It was found that climate change impacts will be spatially heterogeneous and nonlinear due to differences in climate, edaphic conditions, and competitive interactions among other factors (Frey et al. 2016;Creutzburg et al. 2017;Thom et al. 2019).Defining this spatial heterogeneity can assist land managers and policymakers to prioritize areas where a reallocation of resources for climate adaptation could be concentrated (Seddon et al. 2016;Thom et al. 2017).Identifying spatial-temporal patterns of land surface temperature and vegetation greenness at the landscape level is important for decision-making in forest management and biodiversity conservation (Silveira et al. 2021), e.g.identification of areas to implement multi-objective adaptive management (e.g.land-sharing strategies) or areas to priories protection efforts (e.g.creation of new natural reserves) (Peri et al. 2022;Rosas et al. 2022Rosas et al. , 2023)).Forests with high seasonality of temperature or inter-annual variability in the phenology may require management actions that improve spatial variability to enhance resilience to biodiversity loss from high interannual variability (Silveira et al. 2021).In contrast, areas where spatial variability is higher, are more adequate for protection because this is where species can persist if inter-annual variability increases due to climate change (Martínez Pastur et al. 2020a;Trew and Maclean 2021;Ye et al. 2023).
In our analyses, we found significant differences among levels of SOC and forest regions in Argentina, both for inter-annual and spatial variability.Similarly, Silveira et al. (2021) found differences for different eco-regions in Argentina and suggest that species living in different areas face dissimilar threats from variability in temperature seasonality and associated changes in the ecosystem phenology.As we found for the different SOC contents and forest regions, they also highlighted the weakly correlation between land surface temperature and phenology.This suggests the importance of mapping both variables when assessing the threats from phenological and seasonal variations (e.g.areas with high inter-annual variability had a strong threat, and low spatial variability meant that resilience is low).If we agree that high SOC content is a good proxy for biodiversity (e.g.Peri et al. 2019), the results are consistent with many other findings for forests around the world (Fischer et al. 2006;Mori et al. 2013;Schmitt et al. 2020;Messier et al. 2022).However, our landscape analyses showed contradictory results.Inter-annual variability in phenology based on land surface temperature followed the described trend (e.g.rainforests presented higher biodiversity conservation concern than dry forests), while the inter-annual variability in phenology based on vegetation greenness showed the contrary.Besides, high resilience was found in the Patagonian forests indicating more stable ecosystem functioning due to climate stability at low temperature regions (Forest et al. 2015;Payette and Frégeau 2019;Barbé et al. 2020).These mismatches in the general trends and the different forests regions can be explained due to the different ecosystem functionality across Argentina (Silveira et al. 2022), forest structure and composition (Silveira et al. 2023), and human-derived drivers of land use change (Martinuzzi et al. 2021), that greatly influenced native forests.
Argentina implemented a specific national legislation to protect native forests (law 26,331/07) to promote sustainable uses and more effective conservation strategies (Martinuzzi et al. 2018;Martínez Pastur et al. 2020b).
Unfortunately, the resulting land use planning incorporates little biodiversity explicitly information (Guida-Johnson and Zuleta 2013;Lorenzo et al. 2018;Silveira et al. 2022).However, these plans are updated every five years, thus, better knowledge about SOC contents (Martínez Pastur et al. 2022), as well as patterns of inter-annual and spatial variability of vegetation greenness and land surface temperature, can assist to identify those areas where many species may have difficulty to persist under more extreme climate events (Silveira et al. 2021).Under the negative effects of climate change on native forests, these tools can be used to identify vulnerable areas where the slow adaptation process implies that adaptive forest management strategies require long lead-in times (Maciver and Wheaton 2005;Thom et al. 2019).We need to promote in situ biodiversity persistence to avoid negative impacts on future provisioning ecosystem services (Bellard et al. 2012;Thom et al. 2017;D'Orangeville et al. 2018).

Conclusions
In the context of climate change mitigation, managing forest ecosystems to conserve existing carbon stocks in soils and to remove carbon from the atmosphere by adding to stocks is an important issue.In the framework of the Kyoto Protocol and REDD+ that aim to conserve ecosystems and their carbon stocks and achieve land degradation neutrality, information derived in the present work from the estimate of SOC in the native forests can be incorporated into the annual National Inventory Report of Argentina.As Argentina strives to find measures to achieve vital climate targets, this work indicates the importance of considering forest soils in improved management practices to increase climate change mitigation under forest management practices.Increasing forest soils' capacity to store carbon and reduce net GHG emissions is crucial for the Argentinean's target to achieve carbon neutrality by 2050.
The inclusion of SOC as biodiversity surrogates allowed to improve conservation planning through positive correlation with remotely sensed data indices that capture the inter-annual and spatial variability of vegetation greenness and land surface temperature.These indices analysed together can be used to identify areas of high biodiversity threat and resilience to phenological and seasonal variability.Thus, results highlight the importance of long-term monitoring to know the processes that determine the magnitude of SOC variation to forecast how it may operate as climate and land use change in the future and to assist forest management proposals to strengthen the resilience of native forests in order to minimize biodiversity loss.

Appendix 2
Predictive environmental covariates used for modelling the SOC stock in forest native of Argentina.Bold shows the 10 covariates selected for the final model after applying the Random-Forest-Recursive Feature Elimination algorithm , including Tierra del Fuego forests (TDF), Andean Patagonian forests (APF), Monte forests (MON), Espinal forests (ESP), Parque Chaqueño forests (PCH), Yungas rainforests (YUN), and Paranaense rainforests (PAR) (see Appendix 1).

Fig. 1
Fig. 1 Soil organic carbon stock (30 cm depth) in native forests of Argentina ) soil carbon contents (SOC) at 30 cm depth and total SOC stocks (Pg C) calculated for the forest regions in Argentina (PF = Patagonian forests including Tierra del Fuego and Andean Patagonian forests, ESP = Espinal forests, MON = Monte forests, PCH = Parque Chaqueño forests, PAR = Paranaense rainforests, YUN = Yungas rainforests)

Fig. 2
Fig. 2 Covariates importance of the top ten variables for soil organic carbon (SOC) stock (0-30 cm depth) prediction based on the Random Forest model.Acronyms are defined in Appendix 2

Table 2
Multiple analyses of variance of biodiversity conservation concern, considering levels of soil organic carbon contents (SOC, L = low, M = medium, MH = medium high, H = high) at the different forest regions (TDF = Tierra del Fuego forests, APF = Andean Patagonian forests, MON = Monte forests, ESP = Espinal forests, PCH = Parque Chaqueño forests, YUN = Yungas rainforests, PAR = Paranaense rainforests) as main factors IAV-LST = inter-annual variability in phenology based on land surface temperature, SV-LST = spatial variability based on land surface temperature, I-LST = integration between IAV-LST and SV-LST, IAV-PHE = inter-annual variability in phenology based on vegetation greenness, SV-PHE = spatial variability based on vegetation greenness, I-PHE = integration between IAV-PHE and SV-PHE, I-GLO = global integration between I-LST and I-PHE F = Fisher test value, p = probability.Different letters showed differences among levels using the Tukey test at p < 0.05

Fig. 3
Fig. 3 (See legend on previous page.) = inter-annual variability in phenology based on land surface temperature, SV-LST = spatial variability based on land surface temperature, I-LST = integration between IAV-LST and SV-LST, IAV-PHE = inter-annual variability in phenology based on vegetation greenness, SV-PHE = spatial variability based on vegetation greenness, I-PHE = integration between IAV-PHE and SV-PHE, I-GLO = global integration between I-LST and I-PHE.F = Fisher test value, p = probability.Different letters showed differences among levels using Tukey test at p < 0.05

Table 1
Mean (kg m −2 Root mean squared error (RMSE) for different numbers of covariates included in the Random Forest (RF) model obtained with Recursive Feature Elimination algorithm.The filled circle refers to the optimal number of covariates select for the final modellingPeri et al.Ecological Processes  (2024)13:1 Simple analyses of variance of biodiversity conservation concern considering different levels of soil organic carbon contents (SOC, L = low, M = medium, MH = medium high, H = high) as main factor for the different forest regions (TDF = Tierra del Fuego forests, APF = Andean Patagonian forests, MON = Monte forests, ESP = Espinal forests, PCH = Parque Chaqueño forests, YUN = Yungas rainforests, PAR = Paranaense rainforests)