Thinning Effect of C Sequestration along an Elevation Gradient of Mediterranean Pinus spp. Plantations

Forests are key elements in mitigating the effects of climate change due to the fact of their carbon sequestration capacity. Forest management can be oriented to optimise the carbon sequestration capacity of forest stands, in line with other productive objectives and the generation of ecosystem services. This research aimed to determine whether thinning treatments have a positive influence on the growth patterns of some of the main Mediterranean pine species and, therefore, on their Carbon (C) fixation capacity, both in terms of living biomass and soil organic carbon. The results obtained show that C sequestration capacity (biomass and SOC) increased at higher thinning intensities due to the induced alterations in tree growth patterns. We observed almost a 1.5-fold increase in P. nigra and P. sylvestris, respectively, and over a two-fold increase in P. pinaster under heavy thinning treatments; SOC stocks were affected by the intensity of the thinning treatments. These results can contribute to improving silvicultural practices aimed at C sequestration in forest plantations located in dry areas of the Mediterranean.


Introduction
Forests ecosystems are key elements in the mitigation of global warming through their carbon-fixing capacity [1]. Over the last decades, it has been witnessed how politicians are keen on controlling and reducing excessive anthropogenic activities such as the conversion of forests into croplands and urban areas. The Intergovernmental Panel on Climate Change's report [2] addressed the response of countries to reduce the constant increase in atmospheric CO 2 and, consequently, global warming alterations. Terrestrial carbon sinks are considered by the Kyoto Protocol (KP) to offset greenhouse gas emissions. By doing this, an opportunity is provided to count net changes and develop evidence that soil carbon stocks could be affected by management practices [3]. Recent concerning underestimations regarding global CO 2 emission sources have been brought to attention by policymakers. These measurements are critically pointing at global C sequestration to assess the impact of forest management on climate change mitigation (EASAC). A sharp increase in forest harvesting has been found in 26 European countries with potential climate change impacts from forest carbon sequestration and other ecosystem services containing C impacts associated with increased forest harvesting [4].
Silvicultural treatments, normally used for ecological or economic purposes, aim to maximise forest resources, although biodiversity and biological interactions might be indirectly improved [5]. Notwithstanding, timber production, or the process of managing stands seeking woody outputs, is the main goal when seeking an economic asset. A solution combining both purposes would provide the possibility of improving the carbon storage of forests through new forest management practices [6]. Forest management, especially harvesting, could shift soil C balance [7]. Thinning and harvesting techniques are the most influential measures to affect the amount of carbon sequestered [8]. Thinning generally changes the microclimate in pine forests, producing an impact on the amount of litter in the soil, which is temporarily lowered in heavily thinned stands. Thus, thinning contributes to boosting the decomposition rate of forest floor C, since evapotranspiration and soil C pools decrease, affecting the understory herbs as well as soil productivity in nutrient availability [3,[9][10][11]. Numerous studies have described the benefits of thinning on increasing tree growth rate and forest productivity due to the reduction in competition for resources [12][13][14].
Soil organic carbon (SOC) present in the forest, known as a natural form of energy storage, is a major contributor to overall soil health and climate change [15], and it could be a determinant barrier to controlling excess anthropogenic CO 2 emissions. It is also known that from the total carbon accumulated in boreal forests, approximately 30% is stored in soil [16]. SOC concentrations are low in arid regions [1], and due to the impact of climate change, the rise in temperatures will generally lead to a decrease in SOC of 28% in humid, 20% in sub-humid, and approximately 15% in arid zones, which are the objects of this study [17,18]. It is widely recognised that forest soil carbon is a key factor in the global carbon cycle [6]. The understanding of all links in stocks and fluxes in the global C balance is extremely important for categorising the sources and sinks of C and to developing mitigation strategies for climate change. These mitigation strategies require increasing the C sequestration capacity of forest ecosystems and designing adaptive forest management strategies [19,20]. Mediterranean arid zones are considered vulnerable in the context of climate change [21]. Pine plantations established at the rear edge of species distribution are under a notorious growing threat, becoming drier and vulnerable to drought and wildfires as well as limiting its growing resources over time, producing mortality that affects forest ecosystem dynamics such as the carbon balance. Particularly, in dry areas in the Mediterranean Basin, an increase in arid conditions were observed to be linked to a more frequent diminution of precipitation triggering forest dieback [22]. In this sense, adult trees are seeing an increased mortality rate because of summer drought and limited silvicultural practices [23]. The continuous debilitation of forest vigour promotes an increase in carbon emissions [24]. In this paper, we studied the effect of different thinning intensities on C pools (biomass and soil) and sequestration for three Mediterranean pine species (i.e., Pinus pinaster L., P. nigra Arnold. and P. sylvestris L.). These practices can certainly lead to long-term C sequestration in Mediterranean forest soils. With this in mind, we aimed to obtain promising answers for the following hypotheses: (a) Do thinning treatments positively alter tree growth patterns in Mediterranean pine species and, therefore, the biomass C dynamic? (b) Are biomass and soil organic carbon storage that are modulated by thinning treatments beneficial to climate sensitivity? (c) Are we then managing our forests efficiently? The results of this work can contribute to improving the silvicultural practices oriented to the sequestration of C in forest plantations located in dry Mediterranean areas.

Study Area
The study area was located at Sª de Los Filabres, southeastern mountains of the Andalusia region, Almeria, Spain ( Figure 1): 37 • 22 N, 2 • 50 W. The terrain, extending across 45,000 ha, possess a notable altitudinal distribution that varies between 750 and 2168 m.a.s.l.; thus, a wide range of bioclimes are present [25]. The mean annual precipitation is near 320 mm, and the mean annual temperature oscillates between 7 and 16 • C, being appertained within the Mediterranean semi-arid climate region ( Figure S1). Geologically, the area has a pH between 6.5 and 7 and is characterised with xerorthents regosols and limestone substrate [26], mainly developed on schists and quartzites enriched with abundant slates and shales with steep slopes greater than 35% as the predominant topography. Most of the forest stands in the area come from reforestation undertaken toward the third quarter of the last century, with the most frequently used species being logically, the area has a pH between 6.5 and 7 and is characterised with xerorthents regosols and limestone substrate [26], mainly developed on schists and quartzites enriched with abundant slates and shales with steep slopes greater than 35% as the predominant topography. Most of the forest stands in the area come from reforestation undertaken toward the third quarter of the last century, with the most frequently used species being Pinus pinaster, P. nigra, P. sylvestris and P. halepensis. Some natural areas of native pine stands can also be found over the mentioned area.

Experimental Design and Field Data
In June 2008, a thinning experiment was established. The experimental design was based on a randomised complete block setup with three blocks and three treatments for each of the studied species [27]. The chosen layout of the plots and the blocks' location were placed next to each other with a 15 m wide buffer strip around the 1600 m 2 area related to each plot to minimise microclimatic, site and edaphic variations across the experimental design. Three intensities of thinning were applied:heavy thinning (T60) with the removal of 60% of the initial basal area (BA), moderate thinning (T30) removing 30% of the initial BA, and control (C) or no thinning area (Table 1). Thinning treatments were carried out with the main objective of removing overgrown, undersized, dying or suppressed trees to promote future development under natural conditions and uniform densities and spacing. Thinning debris, such as slash and stumps, were removed from the experimental plots. All trees with a diameter at breast height (1.3 m above ground level) larger than 10 cm were measured, labelled and had their dasometric variables measured, that is: diameter at breast height (dbh, cm) using a calliper (Haglöf Mantax, Långsele, Sweden), total height (H, m) using a Vertex III hypsometer (Haglöf, Sweden), stand density (N, stems/ha −1 ) and basal area (G, m 2 ha −1 ) ( Table 1).

Experimental Design and Field Data
In June 2008, a thinning experiment was established. The experimental design was based on a randomised complete block setup with three blocks and three treatments for each of the studied species [27]. The chosen layout of the plots and the blocks' location were placed next to each other with a 15 m wide buffer strip around the 1600 m 2 area related to each plot to minimise microclimatic, site and edaphic variations across the experimental design. Three intensities of thinning were applied:heavy thinning (T60) with the removal of 60% of the initial basal area (BA), moderate thinning (T30) removing 30% of the initial BA, and control (C) or no thinning area (Table 1). Thinning treatments were carried out with the main objective of removing overgrown, undersized, dying or suppressed trees to promote future development under natural conditions and uniform densities and spacing. Thinning debris, such as slash and stumps, were removed from the experimental plots. All trees with a diameter at breast height (1.3 m above ground level) larger than 10 cm were measured, labelled and had their dasometric variables measured, that is: diameter at breast height (dbh, cm) using a calliper (Haglöf Mantax, Långsele, Sweden), total height (H, m) using a Vertex III hypsometer (Haglöf, Sweden), stand density (N, stems/ha −1 ) and basal area (G, m 2 ha −1 ) ( Table 1). Table 1. Characterisation of the Pinus spp. of dasometric variables per treatment plot at Sª de Los Filabres (Almería, Spain). Treat: thinning intensities of stem density removed-T0 (0%), T30 (30%) and T60 (60%); Age: years of the plantation; N: stand density (stems ha −1 ); HO: dominant height (m); DBHm: mean diameter at breast height (cm); G: basal area (m 2 ha −1 ). Values are the means ± standard error (in brackets). Above and belowground biomass (Mg ha −1 ) were obtained using the allometric models developed for softwood species [28], as the compute of the compartment was calculated using biomass Equations (1)-(4) ( Table S1). At this point, a 0.5 standard coefficient was applied to assess the biomass C content [29].

Soil Sampling and Analysis
In June 2016, an extensive field survey took place to select the soil sampling locations in which both orographic and physiographic conditions were considered. Soil sampling was performed by establishing a 12-point soil sampling grid inside the 40 × 40 m 2 experimental plot for each of the proposed treatments. Points were taken with a 5 m margin to the border and a 15 m distance between them and the central samples. Hereafter, systematic distribution of transects lines according to the maximum slope was represented. Soil extractions were conducted and sampled together by layers at four fixed-equal intervals (0-10, 10-20, 20-30 and 30-40 cm) below the forest floor. These were taken using a steel corer (8 cm diameter and 60 cm length) reaching up to a 40 cm depth. Following the chosen extraction procedure, 48 soil soundings were priory projected in each plot. Nevertheless, in some cases, the point location slightly shifted due to the rocky formation of the soil, where it was not feasible to extract the entire sample length. There was a litter layer on these sites; however, this was not sampled since the C present on it produces an overestimation in the final output.
Upon receiving the soil extractions in the laboratory, they were carefully placed to be air dried. Once finished, coarse (≥2 mm) and fine (<2 mm) particles were separated with a sieve. Coarse material was then weighted and stored, and the resultant fine particles fraction was sieved through a 0.5 mm mesh and kept for SOC determination, drawing out all particulate organic matter. At this point, three replicates of each measurement were processed for each of the soil samples. The Walkey-Black method via wet oxidation [30] was applied over the fine fraction to obtain the organic carbon content, and the bulk density (BD) layered in the soil was calculated following [31] Soil organic matter percentage (SOM) was obtained using a 1.64 constant for mineral bulk density [32], and soil organic carbon stock (SOC-S) was consequently obtained in Mg ha −1 following [33]: SOC-S = SOC * BD * D Forests 2021, 12, 1583 5 of 12 where 1.724 is the Van Bemmelen constant factor, D is the thickness of the analysed layer (cm) and SOC40-S is the sum of the SOC in the first 40 cm from the soil's surface obtained for each horizon. The plots' establishment, soil samples' position, tree layer's structure as well as tree positions were obtained using Field-Map technology (IFER-MMS, Field-Map Technology 6 July 2017, and https://www.field-map.com). This technology permits gauging at the level of a single tree measurement through to a level of research or inventory plot for both mapping and dendrometric measurements.

Dendrochronological Sampling and Analysis
As a complementary tool to assess the impact of thinning on tree growth, in June 2016, 15 dominant trees in each plot with a dbh larger than 15 cm were cored with a Pressler increment borer at 1.3 m above ground level. Samples were extracted in pairs from each tree following the direction perpendicular to the maximum slope [34], separated by at least 90º. Upon arrival of the individually labelled samples at the laboratory, cores were air dried at room temperature before preparing them with sandpaper using a progressively finer grain to allow for clear tree-ring exposure for cross-dating analysis [35]. Annual increments were scanned at a 1600 dpi resolution and measured to the nearest 0.01 mm accuracy with a LINTAB device (Rinntech, Heidelberg, Germany), where tree-ring width series were obtained and subsequently evaluated using the computer program COFECHA [36], which adds a high degree of confidence that tree-ring series have been, first, measured accurately and, secondly, cross-dated correctly. Hereafter, with the objective of determining the best growth response for each thinning treatment, tree-ring width data were transformed into the basal area increment (BAI, cm 2 ) using the following formula: where R corresponds to the tree radius, and t is the year of tree-ring formation.
Residual chronologies of the ring-width index (RWI) were complementarily obtained once the influence of long-term biological trends on radial growth was discarded, due to the increasing tree size and age, by applying a double-detrending of the tree ring series upon fitting spline and negative exponential curves to each series, provided by the dplR package in R software.

Statistical Analysis
The impact of forest thinning treatments at the plot level on C sequestration was statistically analysed. The normality of the variables (Shapiro-Wilk test, n ≤ 50, p < 0.05) and homoscedasticity (Levene's test, p < 0.05) were studied. Variables that did not fit the normal distribution were transformed by the square root function, since the exponential function did not fit the criteria. The comparison of growth between thinning treatments and tree species was performed using one-way analysis of variance (ANOVA, p < 0.05 for variables with normal distribution and homoscedasticity, penalizing by Bonferroni's procedure (p < 0.047). Tamhane tests (p < 0.05) were applied when the variables fit a normal distribution but were heteroscedastic. SOC was analysed according to soil layers (0-10, 10-20, 20-30 and 30-40 cm) as well as for the whole soil depth (0-40 cm). To quantify the trends in radial growth, dendrochronological results were also statistically assessed for the last 20 years' growth frame as well as pre-and post-thinning year treatments upon applying the same analysis test. An additional array of dendrochronological parameters was statistically derived to provide greater consistency to our analysis, considering the entire growth period and the growth achieved pre-and post-thinning year for the studied species across treatments. Mean sensitivity (MS), a measure of the relative difference of indexed tree-ring width from one year to the next, mean among trees correlation (Rbar) and derived chronology reliability through expressed population signal (EPS) were calculated. All statistical analyses were performed using "R" version 3.3.1 [37].

Results
Our results showed no significant differences in the total carbon stock accumulated in the above and belowground biomass. In P. pinaster, C content increased from 33.2 Mg of C ha −1 in the control plots to 49.6 Mg of C ha −1 in moderate thinning and 73.1 Mg of C ha −1 in heavy thinning 8 years after thinning (F = 1.40, p = 0.316). The same trend was observed in P. nigra (29.0, 33.0 and 36.4 Mg of C ha −1 , respectively, F = 0.07, p = 0.925) and P. sylvestris (38.5, 43.0 and 51.5 Mg of C ha −1 , respectively, F = 0.27, p = 0.765) ( Table 2). Table 2. Combination of the structural characteristics and biomass C stock of Pinus spp. according to thinning intensities of basal area percentage removed (control or unthinned plots 0%; moderate thinning 30%; heavy thinning 60%) and quantification of soil organic carbon by depth horizon per soil sampling plot at Sª de Los Filabres (Almeria, Spain). Values are the means ± standard error (in brackets), and superscripts (a, b) indicate when pairwise comparisons were significantly different.

Species
Control Moderate Thinning Heavy Thinning

Biomass Fraction C in Biomass (Mg C ha −1 )
Pinus pinaster Total SOC in the top 40 cm of mineral soil showed different patterns across species and treatments, ranging from 18.5 Mg C ha −1 in the control up to 28.1 Mg C ha −1 in moderate thinning but reducing to 15.5 Mg C ha −1 in heavy thinning for P. pinaster (F = 6.39, p = 0.700). On the contrary, SOC decreased in the control treatment (27.6 Mg C ha −1 ) to moderate (23.4 Mg C ha −1 ) and heavy (22.6 Mg C ha −1 ) thinning treatments (F = 1.89, mboxemphp = 0.158). With considerably higher values, SOC also decreased in the control (37.0 Mg C ha −1 ) to moderate (31.6 Mg C ha −1 ) and slightly higher in heavy (34.0 Mg C ha −1 ) thinning treatments (F = 1.23, p = 0.297) ( Figure S2, Tables S2 and S3). BAI values, considering the last 20 years of growth, were significantly affected by thinning intensity. For P. pinaster, BAI20 increased from 4.27 cm 2 in the control to 5.33 cm 2 in moderate thinning and 5.23 cm 2 in heavy thinning (F = 1.431, p > 0.243). Similar trends were observed for P. sylvestris (3.57, 2.93 and 5.50 cm 2 , respectively, F = 4.691, p < 0.010), and for P. nigra (2.49, 2.79 and 2.99 cm 2 , respectively), although in this case differences were not significant (F = 0.510, p > 0.602). Growth rates before and after thinning year (2008) showed no significant differences in P. pinaster (4.10 vs. 5 Table 3 and Table S4). Table 3. Dendrochronological statistics of sampled Pinus pinaster, Pinus nigra and Pinus sylvestris for the three treatments (C, control or unthinned plots; MT, moderate thinning; HT, heavy thinning). TRW: mean tree-ring width; BAI20: mean basal area increment over the last 20 years; BAIpreT: mean basal area increment before thinning year; BAIpostT: mean basal area increment after thinning year (all in mm); MS: mean sensitivity; Rbar, mean among trees correlation; EPS, expressed population signal. Values are the means ± standard errors (in brackets), and superscripts (a, b) indicate when pairwise comparisons were significantly different. The hyphen symbol indicates no thinning treatment in the control plots.

Discussion
Through this work deployment, we have brought to the fore the importance of modulating forest structure after analysing experimental data to extricate the fluctuations in forest growth dynamics, biomass generation and carbon storage affections, which were generated upon varying forest structures via thinning treatment intensities for three of the most important pine species across the Mediterranean Basin. These results further strengthen our confidence, demonstrating how methodical forest management practices can certainly modulate carbon sequestration and forest health over time.
Differences found in the amount of C in the live tree biomass components (above and below ground) were previously demonstrated, where biomass C stock was significantly higher in unthinned plots [8,[38][39][40][41] in P. halepensis, P. pinaster, P. nigra and P. sylvestris populations at a non-approximate location to one another. In contrast, biomass content was assessed considering proximity among species at its harshest edafo-climatic conditions, with a positive increase in C stocks across thinning treatments. Even though the basal area diminished with biomass loss due to the thinning operations, the trees chosen to remain in the thinned plots were those dominants in the stand; thus, the mean diameter at breast height considerably increased due to the reduced competition for resources. Consequently, an exponential trend in C stocks was observed where the thinning was more intense, resulting in an almost 1.5-fold increase in P. nigra and P. sylvestris, respectively, and over a two-fold increase in P. pinaster under heavy thinning treatments (Table 2). Hence, in line with our first hypothesis, thinning treatments altered tree growth patterns in response to mass adaptation when competition for essential resources is diminished, allowing those thinned stands to produce greater trees upon internally fixing C via live biomass generation.
Forest management treatments, such as thinning and harvesting, modify soil C dynamics, resulting in the variation of soil conditions. However, these variations through thinning treatments across soil horizons fluctuated depending on the species, location, parent material and whether the harvesting biomass removed was left on site after thinning operations in addition to its historical management. In the Mediterranean area, soil C stocks are usually lower in content than those located in less arid conditions [42]. Surveyed stands were situated in shallow soils with a low or non-existing humus surface horizon directly overlying the weathering rock in a landscape with active erosion processes, associated with soil degradation and the reduction of valuable land to desertification due to the strong interconnection between biophysical and anthropogenic components caused by agricultural abandonment over the past decades. Hence, looking at our soil data, the tested thinning intensities showed how the soil C stock trend suffered a reduction from the top horizon layer up to 40 cm depth, which is in accordance with [43], except in the unthinned plots where P. nigra was present. P. pinaster significantly affected the C stock in moderate thinning plots, up to the first 30 cm. In addition, the experimental thinning treatments did not greatly increment across intensities and species as reported in previous findings [44,45] due to the harvest management residue impact, except for P. pinaster, where a significant difference was notable in the moderate treatment. This carbon recovery 8 years after harvesting was the principal source of fluctuation [7].
Forest rotation lengths were designed to maintain timber production regularly as well as to sequester carbon through above-and below-ground biomass generation. However, the ability to prove the efficiency of such management practices relies on how we manage our forest plantations.
Looking at the historical climatic characterisation, including mean annual temperature, spring precipitation and annual SPEI ( Figure S1), and supporting this with tree growth in terms of basal area increment (Figure 2), and ring width index ( Figure S3) in the years characterised by severe droughts before (1995 and 2005) and after thinning treatments (2012 and 2015), a decline in radial growth was induced [22] and, therefore, in C sequestration also, indicating how biomass and soil organic carbon storage are modulation benefits climate sensitivity, which is in concordance with our second hypothesis.
Such differences found in terms of C sequestration in unmanaged to managed forests outline the relevant importance of combating forest decline against mitigation effects. As previously reported [40,46], special care should be taken in efficiently managing our forests by applying silvicultural treatments that ensure sustainable ecosystem regulation.

Conclusions
This research work highlighted the importance of forest management in modulating the capacity of forest stands as carbon sinks. By managing the intensities of the thinning treatments within the designed silvicultural itineraries, the carbon sequestration capacity of P. pinaster, P. nigra and P. sylvestris forest stands in Mediterranean environments can be increased, and it can also contribute to the improvement of forest health over time. Thus, carbon sequestration capacity over time tends to increase markedly at higher thinning intensities due to the induced alterations in tree growth patterns. Remnant trees showed a positive response when competition for resources decreased, allowing for a higher generation of live biomass and increasing SOC stocks at stand level and, thus, a higher carbon sequestration capacity of the forest stand.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12111583/s1, Figure S1: (Above) Historic climate characterisation of the study considering mean annual temperature (T); (Below) spring precipitation (P) and the annual standardised precipitation evapotranspiration index (SPEI) for the 1950-2016 period (r 2 and P correspond to annual temperature and the SPEI statistics); Figure S2: Depth distribution of SOC (Mg C ha −1 ) by species and grade of thinning. Control or unthinned plots (0% of BA removed), moderate plots (30% of BA removed) and heavy plots (60% of BA removed) at Sª de Los Filabres (Almeria, Spain); Figure S3: Chronologies of growth index (RWI) for Pinus pinaster, Pinus nigra and Pinus sylvestris. The vertical, dashed lines correspond to the 2008 thinning year. Treatments: C-control (0% BA removed); MT-moderate thinning (30% BA removed), HT-heavy thinning (60% BA removed); Table S1: Selected biomass models for Pinus spp. Ws: biomass weight of the stem fraction; Wb7: biomass weight of the thick branch fraction with a diameter greater than 7 cm; Wb2-7: biomass weight of the medium branch fraction with a diameter between 2 and 7 cm; Wb2+n: biomass weight of the thin branch fraction with a diameter smaller than 2 cm, with needles; Wr: biomass weight of the below-ground fraction (all in kg); d = dbh (cm); h = tree height (m). Equations were indexed by compartments for stem biomass (1), thick branches (2), medium branches (3), thin branches + needles (4), and roots (5); Table S2: ANOVA comparison of species with SOC distribution across horizons and total depth at Sª de Los Filabres (Almeria, Spain). Significant differences indicated by p < 0.05 and no significant differences by p > 0.05; Table S3: ANOVA multiple comparisons of total SOC for the studied species among treatments. Control (0%), moderate thinning (30%) and heavy thinning (60%) at Sª de Los Filabres (Almeria, Spain). Significant differences indicated by p < 0.05 and no significant differences by p > 0.05; Table S4: ANOVA comparison and Post Hoc analysis of standard BAIn, BAI pre & post-thinning year against treatments (C, control; MT, moderate thiunnig; HT, heavy thinning at Sª de Los Filabres (Almeria, Spain). Significant differences (p < 0.05) and no significant differences (p > 0.05) (p < 0.05, * p < 0.01). Funding: This research was collaboratively funded by the LIFE FOREST CO 2 "Assessment of Forest-Carbon Sinks and Promotion of Compensation Systems as Tools for Climate Change Mitigation" (European Commission LIFE14 CCM/ES/001271), ESPECTRAMED (CGL2017-86161-R) and ISO-Pine (UCO-1265298) projects. The information included reflects only the opinion of the authors, and the European Commission/Agency is not responsible for any use that may be made of the information contained in this scientific paper.

Data Availability Statement:
The data that support the findings of this study are available on request from the last author, R.M.N.-C. The data are not publicly available due to the fact of project restrictions.