Long-Term Carbon Sequestration in Pine Forests under Different Silvicultural and Climatic Regimes in Spain

: Proactive silviculture treatments (e.g., thinning) may increase C sequestration contributing to climate change mitigation, although, there are still questions about this effect in Mediterranean pine forests. The aim of this research was to quantify the storage of biomass and soil organic carbon in Pinus forests along a climatic gradient from North to South of the Iberian Peninsula. Nine experimental Pinus spp trials were selected along a latitudinal gradient from the pre-Pyrenees to southern Spain. At each location, a homogeneous area was used as the operational scale, and three thinning intensity treatments: unthinned or control (C), intermediate thinning (LT, removal of 30–40% of the initial basal area) and heavy thinning (HT, removal of 50–60%) were conducted. Growth per unit area (e.g., expressed as basal area increment-BAI), biomass, and Soil Organic Carbon (SOC) were measured as well as three sets of environmental variables (climate, soil water availability and soil chemical and physical characteristics). One-way ANOVA and Structural Equation Modelling (SEM) were used to study the effect of thinning and environmental variables on C sequestration. Biomass and growth per unit area were higher in the control than in the thinning treatments, although differences were only signiﬁcant for P. halepensis. Radial growth recovered after thinning in all species, but it was faster in the HT treatments. Soil organic carbon (SOC 10 , 0–10 cm depth) was higher in the HT treatments for P. halepensis and P. sylvestris , but not for P. nigra . SEM showed that Pinus stands of the studied species were beneﬁced by HT thinning, recovering their growth quickly. The resulting model explained 72% of the variation in SOC 10 content, and 89% of the variation in silvicultural condition (basal area and density) after thinning. SOC 10 was better related to climate than to silvicultural treatments. On the other hand, soil chemical and physical characteristics did not show signiﬁcant inﬂuence over SOC 10 - Soil water availability was the latent variable with the highest inﬂuence over SOC 10 . This work is a new contribution that shows the need for forest managers to integrate silviculture and C sequestration in Mediterranean pine plantations


Introduction
The burning of fossil fuels is the primary cause of the increasing concentration of CO 2 and other greenhouse gases (GHGs) in the atmosphere [1].Only in 2019, emissions de-rived from fossil-fuel combustion and industrial processes accounted for 52.4 GtCO 2 eq or 57.4 GtCO 2 eq when also including those derived from land-use changes.These figures coincided with the second warmest year in the 140-year record [2].Apart from setting emission reduction targets at the mid-term for country members, international agreements aiming at mitigating the effects of climate change (Kyoto Protocol, Paris Agreements, European Climate Policy) also promote carbon sequestration in the biosphere as a key strategy to be used to offset GHGs emissions.Forest ecosystems have an important role in C cycling at the global and regional levels, contributing to store approximately 80% and 40% of the total aboveground and belowground terrestrial carbon storage, respectively [3,4].Within the world's forests C stock, the soil organic matter constitutes the main pool (383 GtC or 44%, to 1-m depth), followed by the sum of the above and below ground biomass (363 GtC or 42%), deadwood (6 GtC or 8%) and litter (4.3 GtC or 5%) [5].
Forests offer two main options for increasing carbon sequestration, by planting trees in currently non-forested land (afforestation), or by allowing the existing forests to accumulate higher biomass [4,6].In addition, the afforestation of degraded or low productive areas through woody vegetation encroachment has clearly enhanced soil organic carbon [7,8].Afforestation has been suggested not only as an effective way to restore degraded ecosystems, but also as a way to mitigate elevated atmospheric CO 2 , hence contributing towards the reduction of global warming [9].
During the last 50 years, Europe has, on average, multiplied the forest biomass and carbon stock per hectare by 1.75 [10].However, the increases in forested area alone, with 5% for conifer forests and 8% for broadleaved ones since 1950, are not enough to explain the significant increments in forest productivity [10].Recently, the impact of summer drought on dry deposition of tropospheric ozone (O 3 ) has been highlighted as an important factor on forest gas dynamics during severe drought period [11].These scenarios suggest that other factors such as changes in environmental conditions [12,13], and in forest management practices are also greatly affecting the carbon sequestration rates [14].In addition, although the objective of afforestation plans to increase carbon stocks is clear, the contribution of forest management in C sequestration is not always considered given the high sitedependency on local factors [15].For example, it has been established that reducing tree density increases litter decomposition by allowing more radiation to reach the soil surface and higher soil moisture levels due to increased net rainfall [16].However, these findings are not applicable to Mediterranean conditions, and in fact null or opposite effects have been reported [17,18].This situation points out the need to better understand the effects of management practices on semiarid environments such as the Mediterranean, where limited information about ecological effects of management exists as compared to temperate forests [19].
In Spain, pine forests occupy 32.6% of the total forested area, resulting the most extended forest type in the country [20].Given their rapid growth and high adaptability to different types of soil, Pinus species have been profusely used in afforestation, especially in areas with harsh climatic conditions and high erosion rates such as those characterized by arid and semiarid conditions [21].Although thinning interventions are typically applied during rotation of Pinus spp.stands to stimulate the growth of the remaining trees [22][23][24], most of them are currently characterized by high tree density because of low timber profitability.Unmanaged pine forests may hold higher amounts of carbon than managed ones mainly due to the differences related to thinning on living biomass [25,26], increasing normally increase with stand age [27].In contrast, no clear consensus has emerged regarding the influence of forest thinning on soil organic carbon (SOC) on the long-term [28,29].SOC is more resistant to changes promoted by forest management than the living biomass [26,28] with the addition that most studies focused on quantifying its temporal dynamics are carried out under natural or unmanaged conditions [28,30].The balance between soil organic matter input from litterfall and CO 2 release during decomposition determines the amount of SOC in forest soils.The chemical nature of C molecules (labile or stable C), site conditions (climate), and soil parameters (clay content, soil moisture, pH, nutritional status) all influence carbon turnover [16].Forest management has an impact on several of these factors, either directly or indirectly.Detailed reviews on how forest management can influence soil carbon have been published [16,28,[31][32][33][34][35], and modelling estimations on how forest SOC can change over long-time frames in different biomes have been compiled [36].All these reviews have pointed to the complex interactions between ecological and management factors, particularly in regions with multiple growthlimiting factors fluctuating along the year, such as the Mediterranean region.Thus, abiotic environmental variables (including climate and soil) and thinning intensity are the main factors influencing the primary response of managed forests in terms of C sequestration (e.g., tree growth, and SOC) [16,26].This has led to a lack of consensus on the effect of how both abiotic variables and thinning intensity might influence C sequestration, and even how that might differ across species and locations [37][38][39].In fact, some authors have questioned whether C sequestration by plants are widely hold because any gain in carbon storage from faster tree growth will be transitory [40].
Thinning affects SOC by altering the microclimate, water and surface energy balances, litterfall input, fine root density, nutrient budgets, and the composition and activity of soil microbial community [29].Thinning intensity, but also time elapsed from intervention, are key aspects when analyzing carbon sequestration in a same forest type [23,25,26].In this sense, an optimum scheme on to study the effect of silviculture in Mediterranean pine forests would address how carbon sequestration (living biomass and SOC) varies in each of the phases of the management operations along the life span of the stand [41].Since long-term experiments addressing the role of forest management on carbon sequestration in Mediterranean Pinus spp.forests are scarce, complementary approaches combining different chrono-sequences and stands growing under different ecological gradients are needed.Here we describe the approach followed in the framework of the Silvadapt.netresearch network [42].
The emerging paradigm in forestry states that proactive silviculture treatments (e.g., thinning) can be used to face the emerging environmental changes by both increasing C sequestration and contributing to climate change mitigation by reducing the impacts of some stressors (e.g., drought, forest pests or wildfire), while maintaining tree growth rate and improving stand adaptation [23,43].On the other hand, other studies contradict this argument showing that heavier thinning reduces biomass pools and C sequestration rates [37,39].Therefore, there are still open questions about the effects of thinning on C sequestration in Mediterranean pine forests [44].Our initial hypothesis was that despite the currently accepted site-and species-specific responses to thinning, general rules on the effect of thinning on pine plantations in semiarid sites could be developed.Hence, to test this hypothesis, the aim of our research was to quantify the storage of biomass and soil organic carbon in Pinus forests along a climatic gradient from the North to the South of the Iberian Peninsula.With this, we intended to analyze how C reserves vary according to different species and site conditions according to thinning intensities, thus establishing a pattern of C sequestration.Therefore, this type of studies will help to accurately analyze the C sequestration of natural and planted Pinus forests widely used in the Mediterranean Basin and to understand the possible relationships between thinning and C sequestration during forest management.

Study Area
Experimental trials were completed at nine Pinus spp.forests along a latitudinal gradient from the pre-Pyrenees to the southern mountains in Spain (Pyrenees in Navarra; Moncayo in Aragón, Cuenca in Castilla La Mancha, La Hunde in Valencia, Los Cuadros in Murcia, Guadix in Andalucía, and Sª de los Filabres in Andalucía) (Figure 1; Table S1 Supplementary Material).
Forests 2022, 13, x FOR PEER REVIEW 4 of 17 Experimental trials were completed at nine Pinus spp.forests along a latitudinal gradient from the pre-Pyrenees to the southern mountains in Spain (Pyrenees in Navarra; Moncayo in Aragón, Cuenca in Castilla La Mancha, La Hunde in Valencia, Los Cuadros in Murcia, Guadix in Andalucía, and Sª de los Filabres in Andalucía) (Figure 1; Table S1 Supplementary Material).
The locations have elevations ranging from 693 to 1350 m.a.s.l and medium to high slopes.Forest cover at these conifer stands is dominated by either Pinus halepensis Mill., Pinus nigra Arnold., or Pinus sylvestris L. These are planted or regenerated forests, mature in most cases, and located in thermo-to supra-Mediterranean climates.Soils vary from stony to sandy and are characterized by low amount of organic carbon and nutrients (pH 6-8, organic matter content 1-3%, coarse texture and low soil water reserve).Forests have been under silvicultural treatments, and similar thinning schemes have been successfully established in all experimental sites.As a result, homogeneous structures of Mediterranean of Pinus forests have been achieved by well-established standard silvicultural treatments oriented to soil protection and valuable timber production.

Experimental Design and Field Data
Experimental sites were selected based on the presence of similar thinning schemes.In this regard, thinning intensity was the silvicultural treatment evaluated and the common factor across sites (Table 1 and Table S1 Supplementary Material).Besides site-specific attributes, thinning intensity at each site was based on the requirements of the existing Pinus forest and used for the primary aim of eliminating overtopped, small-sized, dead, or supressed trees to stimulate future growth, with uniform spacing being a secondary priority.At each selected location, a homogeneous area was used as the operational scale, and thinning treatments were randomly assigned and carried out in plots of 20 × 30 The locations have elevations ranging from 693 to 1350 m.a.s.l and medium to high slopes.Forest cover at these conifer stands is dominated by either Pinus halepensis Mill., Pinus nigra Arnold., or Pinus sylvestris L. These are planted or regenerated forests, mature in most cases, and located in thermo-to supra-Mediterranean climates.Soils vary from stony to sandy and are characterized by low amount of organic carbon and nutrients (pH 6-8, organic matter content 1-3%, coarse texture and low soil water reserve).Forests have been under silvicultural treatments, and similar thinning schemes have been successfully established in all experimental sites.As a result, homogeneous structures of Mediterranean of Pinus forests have been achieved by well-established standard silvicultural treatments oriented to soil protection and valuable timber production.

Experimental Design and Field Data
Experimental sites were selected based on the presence of similar thinning schemes.In this regard, thinning intensity was the silvicultural treatment evaluated and the common factor across sites (Table 1 and Table S1 Supplementary Material).Besides site-specific attributes, thinning intensity at each site was based on the requirements of the existing Pinus forest and used for the primary aim of eliminating overtopped, small-sized, dead, or supressed trees to stimulate future growth, with uniform spacing being a secondary priority.At each selected location, a homogeneous area was used as the operational scale, and thinning treatments were randomly assigned and carried out in plots of 20 × 30 m with a 15-m-wide buffer strip around each plot, considering the similarity of the canopy structural parameters, and with analogous microclimatic and edaphic characters at each plot.This experimental design allowed a comparative analysis of three thinning intensity treatments: unthinned or control (C), light thinning (LT) with removal of 30-40% of the initial basal area, and heavy thinning (HT) with removal of 50-60% of the initial basal area in replicate blocks (Table 1).Although the timeframe of thinning differs among locations (Table S1 Supplementary Material), all measurements used in this study were always obtained after 10-12 years the thinning treatment was carried out to make the data comparable.All the biomass with commercial value was extracted (e.g., stems and large branches), but small residues (small branches and part of the crown) were left in situ, and the remaining trees were tagged.Before the thinning, diameter at breast height (dbh, cm, measured at 1.3 m above ground level, caliper Haglöf, Sweden, Table 1) and tree height (m, Vertex III hypsometer, Haglöf Sweden) were measured.

Growth and Biomass Estimation
Dendrochronological samples were also obtained for the studied sites.Tree-ring widths were calculated and mean basal area increment (BAI) series for each age class were obtained by averaging the annual values of the measured series using the equation [1]: where R is the radius of the tree and t is the year of tree-ring formation.We focused on BAI trends for the 20-year (BAI 20 , cm 2 ) and 10-year (BAI, cm 2 ) periods after thinning in each plot (Table 1).
Forest aboveground biomass was assessed by using allometric equations for Pinus sp. in Spain [45].We considered that the percentage of carbon in all aboveground dry biomass for Pinus is 50.0%.These equations estimate biomass for several fractions of the tree: stem, thick branches (dbh > 7 cm), medium sized branches (2 cm < dbh < 7 cm), twigs and leaves, and roots (Table S2 Supplementary Material).

Soils Samples and SOC Estimation
For each thinning treatment tested, at least nine soil samples were obtained at each site and treatment, and soil extractions were obtained at fixed-equal intervals (0-10 and 10-20 cm depth) below the forest floor using a steel corer (8 cm diameter and 60 cm length) (Table 1) and transported undisturbed to the laboratory.In the laboratory, soil samples were carefully placed to be air-dried, and coarse (≥2 mm) and fine (<2 mm) particles were separated with a sieve.The resultant fine-earth fraction was ground to pass through a 0.5 mm mesh and stored for SOC determination.In the sieving process, all particulate organic matter (rootlets, leaves, seeds, and other plant material) was manually extracted.The gravel (>2 mm) was weighed and stored separately.All the soil samples were processed separately.Organic carbon content of the fine fraction was calculated using the Walkey-Black method via wet oxidation and bulk density (BD) layered in the soil was calculated following Equations ( 2)-( 4) [46] (2)-( 4) in g cm −3 : BD = 100/((OM)/(0.244)+ (100 Soil organic matter percentage (OM, %) was obtained using 1.64 constant for mineral bulk density [47] and soil organic carbon stock (SOC-S) was consequently obtained in Mg ha −1 following [48]: where D = thickness of the analyzed layer (cm) and SOC40-S = sum of the SOC in the first 40 cm from the soil surface obtained for each horizon.

Environmental Variables
We used three sets of environmental variables for building response models: annual climatic variables (temperature, precipitation, and hydric deficit), topographic variables (slope, aspect), and edaphic conditions (soil chemical, physical and water availability characteristics) (Table 2).Climatic data were downloaded from the climatic data sets developed by [49] based on Spanish Meteorological Agency (AEMET) stations.To harmonize these data sets, we interpolated them to the same 1-km resolution (30 arc-seconds, ≈1 km 2 cell size) using the R package meteoland [50].From the final dataset, we calculated annual and seasonal time series of climatic variables (Table S2 Supplementary Material).Soil data were obtained from each plot based on soil geodatabases developed within different research projects (Table S1, Supplementary Material).The soil geodatabases contain information of chemical (% organic carbon, pH and % carbonates) and physical measurements (% sand, % silt and % clay, soil depth, bulk density and % gravel volume) recorded for each soil profile.Digital Elevation Model (DEM) with a spatial resolution of 5 m was downloaded from the National Geographic Institute (IGN) website, and 'Surface Analysis' function (ARC GIS version 10.1, [51]) was used to generate slope and aspect values.Prior to the analysis, we used the Pearson correlation coefficient to assess for any collinearity issues among the explanatory variables [52].Variables with a pair-wise correlation of less than 0.6 were chosen (Figure S1, Supplementary Material).We chose the variables with the most widespread use in the literature and the clearest biological relevance in relation to silviculture response from the sets of strongly correlated variables.

Statistical Analysis
The variables were examined for normality using the Shapiro-Wilk normality test and the homogeneity of the variances using the Levene test.The non-normal distributed variables were normalized through different types of transformation (logarithmic and exponential).The non-parametric Kruskal-Wallis's test was applied to those variables that could not be normalized.All data are represented as means ± standard error (SE).Once the basic requirements have been met, the effects of thinning with respect to the silvicultural, SOC and growth variables were tested using 1-way ANOVA for each of the three levels of thinning intensity.A significant value of F in the ANOVA was followed by mean values comparisons using Tukey Post-hoc at p < 0.05 [53].This procedure was also applied for the comparison between SOC with the main dasometric variables of the studied forests.A two-way ANOVA was carried out to test possible interactions between thinning intensities or treatments and species in the SOC variable.
To model multivariate relationships between response variables (e.g., biomass, SOC and growth) and latent constructs (e.g., silvicultural and environmental variables), we used Structural Equation Modeling (SEM) [54].The full set of explanatory variables were used in previous Pearson's correlation analyses and the selected subset of those variables without collinearity problems was used in SEM (Figure S1 Supplementary Material).Due to the high number of interpolated environmental variables used, the dataset was reduced to five different latent variables, each one describing a group of environmental characteristics with common origin, nature or influence (Soil water availability-SWAV; climatic conditions-CLIM; soil chemical and physical characteristics-SCHR; forest stand structure-silvicultural characteristics, SILV; treatment characteristics-TREAT).The working hypothesis was based on considering the climatic and treatment variables as fixed, and then evaluating the direct and indirect interactions of these two variables on the rest, identifying SOC as the endogenous (dependent) variable (Figure 2).For each grouping, the most representative variables were selected using Principal Component Analysis (PCA), except for treatment where all the possible interactions, including a linear combination of the two variables (Treatment classification-C, LT and HT-the exact thinning intensity) was tested to represent treatments (Table 3).Only variables with communality over 0.6 for the retained Principal Components (eigenvalues over 1), and without collinearity problems were considered.To ensure representativity of latent variables they were only considered for the SEM analysis when the first component resulting of the linear combination of all the chosen variables explained over the 70% of For each grouping, the most representative variables were selected using Principal Component Analysis (PCA), except for treatment where all the possible interactions, including a linear combination of the two variables (Treatment classification-C, LT and HT-the exact thinning intensity) was tested to represent treatments (Table 3).Only variables with Forests 2022, 13, 450 9 of 17 communality over 0.6 for the retained Principal Components (eigenvalues over 1), and without collinearity problems were considered.To ensure representativity of latent variables they were only considered for the SEM analysis when the first component resulting of the linear combination of all the chosen variables explained over the 70% of the total variance of the subset.All correspondence analysis were used as initial hypotheses for the SEM approach.To evaluate the variability explained by the model, R 2 values were calculated for the silvicultural variables and SOC.To determine the best model, different combinations of interactions between the four latent variables and the dependent variable were tested, including the saturated model, and the best one was chosen testing differences in the Aikake Information Criteria trough ANOVA analysis [55].The statistical analysis was carried out using the lavaan package [56] and were represented using semPlot [57] of R software version 3.4.1.[58].In addition, the relationships between SOC and the main influential variables identified in SEM were analysed through generalized linear models to see their direct influence considering each species separately.
Growth trends of the three pine species were similar across all thinning intensities, with a notable increase in growth after thinning (Table 1).Except for P. sylvestris, the control treatments had the lowest annual growth rates (BAI 20 ).BAI recovered after thinning at both thinning intensities, but BAI increased faster in the heavy thinning treatments.For P. halepensis, P. nigra, and P. sylvestris, BAI after heavy thinning was 4.20, 5.60, and 7.29 cm 2 year −1 , respectively; and higher than light thinning (3.73, 4.95, and 7.07 cm 2 year −1 , respectively); however, differences were not significant in any instance.

Structural Equations Model
The saturated model did not fit well data structure (χ 2 User = 180.5;df = 68; p < 0.001), but after eliminating interaction between soil water availability (SWAV) and treatment (TREAT), the results improved significantly (χ 2 User = 41.2;df = 70; p = 0.99) (Figure 2).The model showed significant convergence of predictions (χ 2 Model = 2803; df = 91; p < 0.01) with high values of CFI (1) and TLI (1.01), and good values of RMSEA (≈ 0) and SRMR (<0.06).Moreover, model user test showed significant fit of the user's model with data structure (χ 2 User = 41.2;df = 70; p = 0.99).The resulting model explained 72% of the variation in SOC10 content, and 89% of the variation in silvicultural condition after thinning (Figure 3).In the resulting function model, it is possible to notice that SILV variables were more sensitive to climate and soil water content than CARB.Among the fixed characteristics, climate (CLIM) presented high negative influence in stand silvicultural characteristics (SILV), and negative significant indirect effect over SOC 10 (z = 2.27, p < 0.05) but without significant direct effect.Higher values of evapotranspiration and lower values of temperature and hydric deficit were related with low values of stand characteristics after thinning (lower Gpos and Npos) (Figure S2, Supplementary Material).Treatments (TREAT) influenced positively in SOC 10 , meaning that lower intensities of thinning were related with high values of SOC 10 (Figure S2, Supplementary Material).
Regarding site-dependent characteristics, the latent variable soil characteristics (SCHR) did not show significant influence over SOC 10 or silvicultural response of the stand in the best model.Water availability in the soil (SWAV) was the latent variable with the highest positive influence over the dependent variable SOC 10 .Among the parameters selected for the construction of the latent variable, all of them except the total number of days under soil moisture stress, presented similar positive influence.On the other hand, climate and soil water availability did not present direct relationship but strongly covariate in the same sense for all the studied sites.
The species was significantly determinant in the results of SOC (F = 23.92,p < 0.001) but without interaction with the treatment (Two-way ANOVA Treatment:Species interaction; p > 0.05).When SOC 10 was analysed for each species regarding the most influential characteristics identified in the SEM, there were significant differences in the trend between species (Figure 4).A significant increasing logarithmic relationship was found for Pinus sylvestris between SOC 10 and Npos (F = 56.4;r 2 = 0.45; p < 0.001) and SOC 10 and Gpos (F = 48.38;r 2 = 0.65; p < 0.001).Pinus halepensis showed a significant logarithmic adjustment with Gpos (F = 7.83; r 2 = 0.23, p < 0.01) and marginally significant with Npos (F = 3.83; r 2 = 0.13; p < 0.1), decreasing from lower values of x-axis to higher ones, a significant positive relationship between thinning intensity and SOC 10 (F = 8.74; r 2 = 0.25; p < 0.01), and a significant negative differences regarding treatment (F = 4.24; p < 0.05), with lower values of SOC 10 in heavy thinning respect the control.Pinus nigra did not show significant adjustment with any variable, for logarithmic or linear regressions.

Discussion
The study of thinning treatments across a range of pine species and site conditions aids in the implementation of forest management for pine forest adaptation to climate change.Previous research in Pinus forests in Spain has combined dendrochronology, stable isotopes and growth productivity to quantify the impacts of thinning [23,28,[59][60][61][62][63][64][65].In accordance with those studies, silvicultural treatments provide evidence that intense thinning has the effect of reducing forest biomass but increasing tree growth per unit area (e.g., expressed as BAI).However, when also considering SOC, thinning only showed a significant positive effect on C sequestration for P. halepensis, but negative for P. sylvestris and P. nigra [26,28,66].Our data suggest that the regulation of tree density improves tree growth performance in P. halepensis, P. nigra and P. sylvestris stands, although climate variables showed to have a stronger relationship with SOC10 than silvicultural treatments.

Effects of Thinning on Radial Growth
Growth and biomass variation related to thinning treatments have been one of the silvicultural principles.Thinning has shown a positive influence on tree growth per unit area in a variety of Pinus species across Spain (e.g., [25,59,64]).Heavy thinning (HT) resulted in the highest growth rates for all species, albeit the change in post-thinning BAI was greater in P. nigra and P. sylvestris than in P. halepensis, which only showed a modest

Discussion
The study of thinning treatments across a range of pine species and site conditions aids in the implementation of forest management for pine forest adaptation to climate change.Previous research in Pinus forests in Spain has combined dendrochronology, stable isotopes and growth productivity to quantify the impacts of thinning [23,28,[59][60][61][62][63][64][65].In accordance with those studies, silvicultural treatments provide evidence that intense thinning has the effect of reducing forest biomass but increasing tree growth per unit area (e.g., expressed as BAI).However, when also considering SOC, thinning only showed a significant positive effect on C sequestration for P. halepensis, but negative for P. sylvestris and P. nigra [26,28,66].Our data suggest that the regulation of tree density improves tree growth performance in P. halepensis, P. nigra and P. sylvestris stands, although climate variables showed to have a stronger relationship with SOC 10 than silvicultural treatments.

Effects of Thinning on Radial Growth
Growth and biomass variation related to thinning treatments have been one of the silvicultural principles.Thinning has shown a positive influence on tree growth per unit area in a variety of Pinus species across Spain (e.g., [25,59,64]).Heavy thinning (HT) resulted in the highest growth rates for all species, albeit the change in post-thinning BAI was greater in P. nigra and P. sylvestris than in P. halepensis, which only showed a modest rise in BAI when compared to BAI 20 .Previous research focused on the thinning effects in different pine species (e.g., [59,67,68]) found growth per unit area increment in the three Pinus species.P. nigra and P. sylvestris plots were located in more mesic sites, thus a possible higher tree water availability led by thinning was translated into a greater overall assimilation capacity, agreeing with the result of the SEM analysis, which showed high influence of soil water availability after thinning.Hence the post-thinning remaining trees at these mesic sites were more able to take advantage of the reduced competition pressure for resources than those trees in more arid sites [69].
In relative values, the better response of P. halepensis and P. nigra to heavy thinning respect to P. sylvestris might be due to higher drought tolerance capacity in highdensity stands [61,63].P. sylvestris, a more drought-sensitive species, showed the highest growth increment when using light thinning (LT) which is in concordance with previous studies [70].Low growth rates in dense and unthinned plots, on the other hand, indicate less drought resilience in high density stands, making these forests more susceptible to drought-induced dieback and mortality [65,71].

Thinning as a Tool to Optimize C Sequestration
Previous studies have shown contradictory effects of thinning on the increase in biomass and SOC in thinned plots, and thus, this silvicultural practice does not appear to have a clear effect on C sequestration [72].Our results show that thinning reduced aboveground biomass, but this was not translated into a significant decrease in SOC.In fact, for P. halepensis the highest value of SOC was observed in the heavy thinning for both soil depths of 10 cm and 20 cm.Previous research has found that numerous pine species respond similarly to thinning treatments in terms of C sequestration related to a more conservative growth strategy for these species at high densities [73].
SEM showed that thinning influenced positively SOC10, which was also positively affected by silvicultural condition after thinning (density and basal area).Those results are consistent with previous studies showing the positive effects of the remaining trees in thinned stands on C sequestration.It is possible that this response was linked to the presence of logging residues with increased nutrient availability for growth and leaf formation [74][75][76].We observed that in the thinned stands, growth per unit area (e.g., BAI) was higher, implying that thinning increases resource availability [77].This effect was stronger in P. halepensis, the most drought-tolerant species, possibly because of a hydraulic adaptation that resulted in higher total water usage efficiency than in P. nigra and P. sylvestris [78][79][80].
There was also a negative correlation between climate variables (CLIM) and postthinning silvicultural condition (SILV), and indirectly over SOC.This relationship was, on the other hand, positive with soil water availability (SWAV).The negative relationship between SOC and climate variables through post-thinning condition suggests that increasing stress over time (higher hydric deficit and annual mean temperature) may limit growth.This could be linked to tree-level reactions to reduced competition, particularly enhanced water availability, which is the most significant constraint to forest productivity in Mediterranean ecosystems.As a result, changes in water availability induced by reduced tree competition in dense stands may have a significant influence in these managed stands' ability to store carbon.Due to more favourable gas-exchange and photosynthetic rates, trees with limited competition are more sensitive to water availability.Reduced competition boosts soil water content, soil-to-canopy hydraulic conductance, stomatal conductance, and photosynthetic rates in the near term, although it is unclear whether these effects persist over time [63].Thinning may also raise water demand and increase transpiration of the remaining trees due to their larger leaf area, and stimulate understory vegetation development, also limiting soil water availability due to competition [69].These long-term implications of thinning were not completely covered in this research, although the different response between the thinned and control stands appears to extend over time, corroborating the notion that the remaining trees have more growth resources (see [74] for a review).
This work brings together the results from different types of Mediterranean pine sites and species.Thinning has shown a species-related effect on total carbon sequestration (total biomass and soil organic carbon) but thinning improves the state of the remaining trees.The results are more promising for P. halepensis plantations due to a significant increase in SOC with thinning.Although results are not conclusive, this study shows that forest management can affect the SOC, even if it is subordinate to the soil water availability.In Mediterranean sites, it is common to observe contradictory results, since there are many limitations that interact at the same time (aridity, but also high and low temperatures, poor soil due to ancient practices, etc.).

Conclusions
One of the most serious issues confronting Pinus forests in southern Europe and other parts of the Mediterranean Basin is the lack of silvicultural treatments that could reduce climate-related mortality processes.Thinning could assist to minimize the effects of climate change on Mediterranean drought-sensitive Pinus species by boosting the availability of nutrients and water to the remaining trees.As a result, stand structure management may be able to modify C sequestration capacity, and various authors have stressed the importance of thinning on the C cycle.We add to this subject by providing a new example on the effects of silviculture on C sequestration.Thinning had a positive effect on growth per unit area in three ecologically different pine species across a climate gradient in the Iberian Peninsula.Our findings demonstrate that high-intensity thinning helped Pinus stands to increase their growth quickly, although climate variables showed to have a stronger relationship with SOC 10 than silvicultural treatments.Soil parameters, on the other hand, had no significant influence on SOC 10 , with the latent variable soil water availability having the greatest influence on this variable.This study is a new addition that highlight the need for forest managers to integrate silviculture with carbon sequestration practices in Mediterranean pine forests.Despite its limited effect, thinning is the most important silvicultural technique for promoting steady adaptation of forest structure to climate change, especially creating low or very low dense forests, while decreasing tree mortality processes that cause longterm damage and broad social alarm.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10.3390/f13030450/s1, Figure S1.Correlation matrix computing the correlation coefficients between environmental and silvicultural variables with response indicators of interest.Features with high potential are identified by color legend.See variables description in Table 2. Figure S2.Structural Equation Model path with the standardized estimators of each variable (boxes) in the nodes (circles).Arrow labels show the standardized effects.Table S1.Site descriptions and characteristics of selected stands.Values are means ± SD.Table S2.Selected biomass models from SUR fitting statistics for softwood species [45].Table S3.Environmental data used to predict the effect of thinning on soil C sequestration (SOC10).

Forests 2022 , 17 Figure 2 .
Figure 2. Conceptual model for relationships between environmental variables, site variables and carbon content.Regular arrows show direct effects, and dashed line indicates indirect effects.

Figure 2 .
Figure 2. Conceptual model for relationships between environmental variables, site variables and carbon content.Regular arrows show direct effects, and dashed line indicates indirect effects.

Figure 3 .
Figure 3. Structural Equation Model path describing the effects of multiple predictors in carbon stock and dasometric variables after treatments.Only significant effects or covariances were showed.Nodes represent different categories of predictors and endogenous variables: CLIM = Climatic variables; SWAV = Soil water availability; TREAT = Treatment variables; SILV = Silvicultural status of forest after treatment; CARB = SOC10.The ellipsoid of CARB does not represent a latent

Figure 3 .
Figure 3. Structural Equation Model path describing the effects of multiple predictors in carbon stock and dasometric variables after treatments.Only significant effects or covariances were showed.represent different categories of predictors and endogenous variables: CLIM = Climatic variables; SWAV = Soil water availability; TREAT = Treatment variables; SILV = Silvicultural status of forest after treatment; CARB = SOC 10 .The ellipsoid of CARB does not represent a latent variable.Boxes represent the R 2 values for the explained variance of dependent variables.Red arrows mean negative effect.Green arrows mean positive effect.Double-headed pivot arrows mean covariance between groups of variables.Arrow labels show standardized effects magnitude.χ 2 Model represents the parameters for the model fit, and χ 2 User represents the parameters for the comparison between variances of the model and the dataset.

Table 2 .
Environmental data used to predict the effect of thinning on soil C sequestration (SOC 10 ) after the PCA selection for the construction of latent variables.

Table 3 .
Selection of exogenous variables for the SEM analysis, carried out through Principal Component Analysis (PCA).#comp.: Number of components retained for the PCA optimum solution without rotation (eigenvalue above 1).Cum.Var.: Percentage of the variance of the dataset explained by the retained components (eigenvalue above 1).SWAV = Soil water availability; CLIM = Climatic characteristics; SCHR = Soil characteristics, SILV = Silvicultural characteristics after treatment.See Table2for variables description.