Tree-size heterogeneity modulates the forest age-dependent carbon density in biomass and top soil stocks on Mediterranean woodlands

Carbon accumulation in terrestrial ecosystems is inherent to the vegetation development and ageing process. Primary productivity synthetize biomass which is constantly incorporated to soil. Vegetation community composition, and other ecological drivers, are known to mediate biomass production. However, links between forest developmental stage and ecological drivers of carbon stocks are unexplored. We address this topic under the prediction that species-rich and uneven-sized forests can improve carbon storage potential in biomass and topsoil fraction across its development. The study was carried in forest stands growing under Mediterranean conditions in Central Spain. Carbon content in both above- and below-ground tree biomass and in topsoil organic matter (0–40 cm) was measured in 30 sampling plots of variable size (900–3000 m2). We also assessed Shannon species diversity index, Gini tree-size inequality index and forest developmental stage using dendrochronological procedures to derive the mean age of the oldest trees. First-order interaction terms between diversity factors and forest age were regressed against carbon density in compartment-independent regressions. Forest-age and tree-size heterogeneity coupling was the main factor driving carbon density of both compartments. The interaction showed that woodlands maximize density in aged forests composed by uneven-sized trees. Models gave not support to consider species diversity as a mediator of carbon stocks. Our results shed light on how tree-size heterogeneity can regulate the temporal dimension of forest ageing to rise the carbon storage potential. Mature forests in semi-arid environments cannot store carbon due to their intrinsic ontogeny, they need to grow structurally diverse.


Introduction
Forests store vast amounts of carbon above-and belowground in compartments subject to ecological drivers that are presumed stationary. Forest age and diversity (functional dispersion, species diversity and structural heterogeneity included) are among the factors most frequently studied in driving forest productivity (Chun et al. 2020;Dănescu et al. 2016;Silva Pedro et al. 2017). However, little is known about whether these variables control simultaneously the carbon densities of two dissimilar but connected compartments in semi-arid environments (Godlee et al. 2021): i.e., tree and root biomass carbon (AGBC + BGBC) and topsoil organic carbon (SOC). Moreover, evidence of a possible coupled effect of diversity variables and forest age on carbon stocks of forest biomass and soils is scarce (Zeller and Pretzsch 2019). Identifying any possible interconnections could lead to new insight on how to approach sustainable management practices focused on forest carbon sequestration across the forest developmental stages.
Plants, particularly those producing woody tissues, have a prominent role in the carbon sequestration, where trees excel at this function owing to their capacity to store large volumes of biomass during centuries. As an initial step in understanding the carbon fixation, forest productivity is one of the most commonly studied ecosystem processes either per se or by using proxies like cumulated above-ground biomass or tree radial growth. In this sense, factors such as climate (Bastos et al. 2016;Zscheischler et al. 2014) and biotic factors like tree size (Stephenson et al. 2014), tree age (Pregitzer and Euskirchen 2004) and/or stand dynamics (Molina-Valero et al. 2021) are well known to impact vegetation and forest productivity. Additionally, tree species diversity has also been proposed to foster tree growth (Forrester and Bauhus 2016;Liang et al. 2016;del Río and Sterba 2009;Ruiz-Benito et al. 2014). Mechanisms to explain the causal relationship between diversity and productivity are mainly grouped under two main ideas mutually compatible (Carrol and Nisbet 2015): the complementarity effect and the selection/sampling effect. While the former postulates that a more efficient use of resources can be achieved through niche partitioning, the latter supports the idea that more species increase the probability of finding certain species more efficient at resource exploitation and usage. Nevertheless, in forest communities the diversity-productivity coupling can also be interpreted as a masked influence of structure heterogeneity, since tree size is often dependant on tree species identity (Forrester 2019). Consequently, several authors have identified the positive effects of tree size heterogeneity on productivity (Dănescu et al. 2016;Madrigal-González et al. 2016;Ruiz-Benito et al. 2014;Silva Pedro et al. 2017), as well as negative effects in some instances (Bourdier et al. 2016;Yáñez et al. 2017). Biomass in forests is mainly accumulated in woody tissues as the result of tree productivity which in turn is influenced by multiple drivers described above. Nevertheless, it still remains unclear whether the forest age-carbon biomass relationship is mediated by different diversity metrics (Aponte et al. 2020;Zeller and Pretzsch 2019).
According to the main ideas introduced above, SOC is partially related to the same factors affecting biomass production because a large part of the SOC derives from the plant residues (Liang et al. 2019;Mayer et al. 2020). In forests, vegetal inputs are heavily represented by tree aerial parts (leaves, small and big branches, flowers and fruits, bark litter, dead stems…) but also by underground processes like the continuous detachment of fine roots that rise carbon density according forest age (Pregitzer and Euskirchen 2004). This phenomenon produces a vertical gradient of biomass in which SOC concentration decreases according to soil depth (Jobbágy and Jackson 2000), which suggests that root vertical distribution may play an important role in the distribution of carbon in soils. Due to this complex configuration, it is not easy to clarify the extent at which drivers of soil carbon stocks operates. For instance, in grassland communities, the positive effect of plant diversity on carbon density driven by root inputs was found to be contingent on soil microbial activity (Lange et al. 2015). This mechanism may explain as well how species diversity can also promote directly both below-ground biomass and carbon stocks (Chen et al. 20,218). Other evidence suggests that the absence of direct diversity effects that may be true only if this influence is mediated by environmental conditions (Felipe-Lucia et al. 2018). The lack of consensus is more evident if considering a possible joint influence of tree diversity and temporal dimensions like forest age over SOC density (Conti and Díaz 2013;Jonard et al. 2017;Vesterdal et al. 2013). In secondary forests, tree species diversity and structure heterogeneity are usually dependent on human activities that also have an impact on carbon stocks. Besides, the effect of forest management is difficult to quantify and compare (Wiesmeier et al. 2019). A rule may be that soils of unmanaged forests hold more carbon than the managed counterparts (Schulp et al. 2008). To turn this situation, the study of species diversity and forest structure factors may help to clarify a debate that could help to carry out a successful land management focused on carbon sink.
The aim of this study was to examine quantitatively whether the influence of diversity drivers of above-and below-ground carbon density is contingent on forest ageing process. Specifically, we applied linear models to evaluate differences on how species diversity-forest age and tree size heterogeneity-forest age couplings drive stocks of AGBC + BGBC and topsoil SOC. We hypothesize that (i) carbon stored in compartments is dependent on forest age, tree species diversity and tree size heterogeneity. We assume to find differences in the size-effect of diversity to control carbon stocks depending on the compartment considered (woody biomass or top soil). The temporal dimension represented by forest age should appear equally relevant for both compartments bearing in mind that older stands may hold larger carbon stocks than younger forests. Also, we propose that increased size heterogeneity and species diversity could help to retain more carbon in the standing biomass. Moreover, such influence should also be evident on SOC, owing to a more efficient use of space and a more intense filling belowground due to niche segregation in root systems (Brassard et al. 2011;Aponte et al. 2020). In addition, the relationship of carbon density with forest age may appear modulated by diversity. In this regard we expect that (ii) carbon stored will be comparatively higher in diverse and uneven-sized forest plots. This control will be represented by the firstorder interaction between forest age and tree species diversity and tree size heterogeneity and between species diversity and tree size heterogeneity.

Study area
We established the study area in the northern slopes of the west subdivision of the Central System Mountains, in central Spain. The geographic area selected comprised two mountain ranges: the Sierra de Francia-Quilamas and the Sierra de Gredos. Mean annual precipitation and temperature values for these regions are 1200 mm and 12.5 °C in the Francia-Quilamas range and 500 mm and 11 °C in the Sierra de Gredos range (Spanish State Meteorological Agency, AEMET). Both areas present the characteristic Mediterranean summer drought, and humid and cold weather conditions in winter. The regional geologic and lithologic configuration is mainly dominated by siliceous substrates. Granites are found in the hillsides of both Gredos and Francia-Quilamas, but metamorphic rocks can also be found in the lowlands of Francia-Quilamas (Martínez-Graña et al. 2017, 2019, i.e., slates. The most common soil types in which trees establish are leptosols and regosols in abrupt relief areas and cambisols in the more gradual slopes.

Sampling design
We selected 30 forest stands (14 in the Francia-Quilamas subsector and 16 in the Gredos subsector, Fig. 1) where silvicultural practices were not detected when sampling plots were defined. One plot was established in each stand for the following soil and dendrochronological data collection. Plot size was variable in the areas, ranging between 900 and 3000 m 2 , because of differences in tree density between stands. All field work, sampling collection and soil analysis were conducted between 2012 and 2019. Plot elevation varied from 625 to 1355 m a.s.l. in the Francia-Quilamas subsector and between 1040 and 1440 m a.s.l. in the Gredos subsector (Appendix S1: Table S1). Due to the distinguishable pattern formed in the geographic location of forest stands ( Fig. 1), we grouped plots into two differentiated clusters. More information about clustering can be found on the "Statistical analysis" section.

Dendrochronological and tree data collection
We followed standard dendrochronological procedures to estimate plant age. Within each sampling plot, all adult woody plants were identified at the species level and tagged. If alive, they were sampled by extracting a core at 0.3-0.4 m above ground using Pressler increment borers (Häglof, Sweden) with the intention of registering as many rings as possible (Veblen 1992). In addition, Diameter at Breast Height (DBH, measured at 1.3 m above ground) was recorded for all living or dead adult woody plants. We categorized adult as those individuals presenting a DBH ≥ 5.0 cm and height ≥ 2.00 m. Cores were sanded until the tree rings were clearly visible and then measured using scanned images together with CooRecorder 7.6 software (Cybis, Sweden). Then, we proceeded with visual and statistical cross-dating using Cofecha (Holmes 1986). If a core did not reach the pith, the number of missing rings was estimated using the geometric procedure described by Duncan (1989).

Forest characteristic data
As a measure of forest structure heterogeneity and species diversity, we selected popular indices used in the fields of ecology and/or forest management Fig. 1 Map showing the location of the sampled plots. Cluster encompasses the plots of Sierra de Francia-Quilamas range (cluster 1) and Sierra de Gredos range (cluster 2). For more information about clustering, please see "Statistical analy-sis" section. Data regarding plots such as altitude, geographic location and forest characteristics is available in Appendix S1; Table S1 (Appendix S1: Table S1). For this purpose, DBH tree individual measurements were converted to Basal Area (BA), assuming circular stems, to proceed with.
The Shannon diversity index (Shannon and Weaver 1949) combines two components, the number of species and the relative abundance of each species. The Shannon index was computed using the proportion of BA relative to a species regarding the total BA of a sampling plot. On the other hand, the Gini coefficient is a measure of inequality among values of a distribution (Gini 1912) and have been proposed as a measure of forest structure diversity (Lexerød and Eid 2006). Theoretically, the Gini coefficient ranges between 0 and 1. Thus, 0 depicts a perfect equality, i.e., all trees within a plot present the same BA, while 1 postulates the highest possible inequality (one tree represent all the BA of the plot). Gini coefficient, and hence forest structure heterogeneity, was calculated as follows: Where BA is the basal area of the tree of rank j, and j is the rank of a tree in ascending order from 1 (the smaller tree) to the total number of trees n.
Tree age was derived as the number of rings. Specifically, forest age was considered as the average number of rings in the three oldest trees living in each plot (Di Filippo et al. 2017) (Appendix S1: Table S1).

Soil data collection
Soil data were recorded at two different depths, from the ground surface (0-20 cm) to the deeper upper soil layer (20-40 cm). This sampling procedure was adopted to record the possible variability of SOC caused by the vertical gradient (Jobbágy and Jackson 2000). All subsamples were obtained at random points within a sampling plot.
On the one hand, to assess organic matter content in each plot, we took four subsamples at 0-20 cm and two subsamples at 20-40 cm. Subsamples were pooled in a single sample per depth and plot for chemical soil analysis. Afterwards, soil samples were processed by grinding them and then passed through a 2-mm sieve to remove the larger particles. Organic matter content was dated by wet oxidation using a dichromate-sulphuric acid mixture (Walkley 1947). Residual dichromate was back titrated using (1) Forest structure heterogeneity Gini coefficient ferrous sulphate. The difference in the amount of FeSO 4 added compared with a blank titration determined the amount of easy oxidizable organic carbon (Walkley 1947).
On the other hand, aiming to quantify soil bulk density, we extracted core subsamples with cylindrical cores (279.9 cm 3 volume). 2-3 subsamples were extracted at 0-20 cm and 1-2 subsamples at 20-40 cm. Subsamples were later oven dried for 48 h and then weighed to the nearest 0.01 g. Bulk density values were averaged by depth and plot.

Carbon density estimations
Regarding woody vegetation, biomass was computed for all standing living and dead individual trees by applying allometric equations derived from experimental designs on DBH measures. A different allometric equation was applied for each species and tree component (aerial or root) (see Table 1 for equation's references). In cases where this was not possible, equations belonging to equivalent species and/or regions, inspired by Gil et al. (2011), were used. Hereafter, tree biomass was considered as the sum of aerial and root components. Biomass was transformed to carbon using conversion factors (0.500 for Pinus spp., 0.475 for Quercus pyrenaica and 0.4735 for the rest of species) according to Gil et al. (2011). The summation of all biomass carbon within a plot established the AGBC + BGBC compartment.
We quantified top soil SOC, assuming that 50 per cent of soil organic matter is carbon (Pribyl 2010), by relating the organic matter content and bulk densities obtained in the previous soil analyses. This calculation was carried out independently for each sampling plot and soil depth. After, top soil SOC was considered as the sum of the SOC contents at two depths (0-20 and 20-40 cm) within each plot.
Finally, carbon stocks were expressed as carbon densities as megagrams per hectare (Mg ha − 1 ).

Statistical analysis
To assess the influence of forest characteristics on carbon density we fitted linear models with the lm function in the R environment (R Core Team 2020). We analysed the effects of forest age, forest structure heterogeneity and species diversity on carbon densities of two compartments independently: AGBC + BGBC and top soil SOC. We also considered an additional categorical explanatory variable to take care of possible statistical independence violations derived from the spatial proximity between plots (Duncan et al. 1961) (clusters 1 and 2, see Fig. 1 and Appendix S1: Table S1).
Specifically, we built models to analyse carbon density as a linear function of plot clustering, forest age,  (Zuur et al. 2010).

Carbon stocks of tree biomass and top soils
Stands presented a total mean carbon value of 214.22 ± 16.86 Mg ha − 1 (mean ± SE) computed as the sum of the two SOC and AGBC + BGBC carbon compartments (Fig. 2). The AGBC + BGBC compartment showed a value of 122.55 ± 12.20 Mg ha − 1 while the top soil SOC content was slightly lower with 91.67 ± 6.76 Mg ha − 1 . The carbon fraction in dead trees was generally low or simply absent in many of the sampling plots (Fig. 2).

Effects of forest characteristics in carbon stocks
The models fitted showed the interaction between forest age and forest structure heterogeneity had a significant effect on the carbon stocks in both compartments, but this relationship was stronger for AGBC + BGBC (Table 2; AGBC + BGBC: t = 4.96, p < 0.001; SOC: t = 2.37, p = 0.027). In the case of AGBC + BGBC, the combined effect of forest age and structure heterogeneity resulted in the estimated coefficient with the strongest influence on carbon density, while in top soil SOC this interaction had a lesser effect ( Table 2). The sign of this interaction was primarily positive in both cases in the sense that higher carbon densities coincided with the highest values of both forest age and structure heterogeneity (Table 2; Fig. 3). Specifically, forest age had a null effect over the carbon fraction when structure heterogeneity was below-than-average as seen in both compartments (Fig. 3); however, the effect of age turned positive as the Gini coefficient increased. An absence of significant differences between the two spatial clusters indicated an absence of spatial autocorrelation. Furthermore, VIF values showed a lack of multicollinearity among predictors since all explanatory variables resulted in VIF < 3.00 units (Appendix S2). In addition, to test for a possible masked relationship between forest heterogeneity structure and age predictor variables, we fitted a linear model with these particular variables which resulted in a non-significant relationship (p = 0.956, Appendix S2).
All terms including species diversity either as a main effect or through the interaction effect with forest age or with structure heterogeneity were non-significant (Table 2), what suggests the Shannon index is not an effective driver of carbon stocks.   . 3 Graphical representation of the predicted carbon density of forest biomass (AGBC + BGBC) and topsoil (SOC) computed based on the interaction effect between forest age (average age of the three oldest trees) and forest structure heterogeneity (Gini coefficient of BA) Table 2). To illustrate this effect the Gini index have been converted to a categorical variable of three levels (-1, 0 and 1 z-scores) to simplify the representation of the interaction term. Negative Gini z-scores describe less-than-average heterogeneous forests regarding tree sizes and positive values the opposite. Lines represent predicted carbon densities and points the observed values

Discussion
Our results evidence how the temporal effect of forest ageing process on the two analysed carbon stocks operate contingent on tree size heterogeneity on semi-arid Mediterranean woodlands (Table 2). At a first glance, this effect is mainly positive, indicating that the largest carbon pools take place in mature uneven-sized forest stands. Also, our results indicate tree species diversity does not influence carbon density, and those factors driving forest productivity may not always mediate carbon storage. Interestingly, our results highlight that mature forests cannot maximize their carbon stocks in a structurally homogeneous environment.
Focusing on AGBC + BGBC, our results suggest species diversity does not influence the benefits of forest age as an effective driver of carbon density. Previous studies have pointed to greater biomass productivity in diverse forests in terms of composition and/or structure (Dănescu et al. 2016;del Río and Sterba 2009;Ma et al. 2020;Ruiz-Peinado et al. 2017). However, this relationship was found to be contingent on the developmental stage of the forest (Bongers et al. 2021;Lexerød and Eid 2006;Zeller and Pretzsch 2019). Accordingly, other research has reported that stand structure is one of the most relevant factors driving biomass carbon density (Ruiz-Benito et al. 2014;Vayreda et al. 2012;Aponte et al. 2020). Complementing our results, Felipe-Lucia et al. (2018) found that stands with high vertical stand heterogeneity harbor larger standing tree biomass pools. By contrast, the effect of tree species diversity was not observed in our work despite the fact that it could be expected that forests rich in species would lead to forests rich in structure (Hakkenberg et al. 2016;Mensah et al. 2020;Godlee et al. 2021). Nonetheless, due to insufficient taxonomic diversity among the sampled stands, or due to diversity-productivity relationships not always been successfully established (Felipe-Lucia et al. 2018;Zeller and Pretzsch 2019), our results help to solve this issue.
The mechanisms governing the storage of carbon with forest ageing, although obvious, have been addressed with subtly different results. Theoretically, stand carbon stocks increase conform forests experiencing ageing processes (Lee et al. 2016). This implies that tree biomass productivity decreases when forests communities reach late-successional stages where ageing processes gain more relative importance over the other (Collati et al. 2019;Curtis and Gough 2018;Pregitzer and Euskirchen 2004). Therefore, carbon density should not increase considerably after some stage-developmental related tipping point. Yet, other authors suggested that the carbon accumulation rate at the tree level increases exponentially and continuously with tree size (Stephenson et al. 2014). Assuming a finite threshold of carbon density in forests, these two a priori divergent affirmations may evidence that in late-successional stages some individual trees can store carbon persistently in detriment of the death or suppression of another smaller trees, what should foster forest structure heterogeneity. Our fitted models indicate that, in our study region, AGBC + BGBC density is still experiencing a linear pattern even though seven plots are, on average, older than 100 years old (Appendix 1: Table S1). Whether the effect of age is linear or not, the role of forest ageing is not sufficient to explain per se carbon densities of AGBC + BGBC since the structure heterogeneity had a conditional effect over forest age. This could be interpreted as the conditioning of structural diversity on the carbon stocks, which may help these ecosystems to maximize carbon density in earlier successional stages (Fig. 3a). Thus, our results suggest a prevalence of the complementarity effect probably driven by uneven-sized trees rather than tree species that becomes more important according to ageing processes. Our results partially disagree with the recent work of Molina-Valero et al. (2021) who found that old forests in the Iberian Peninsula can gather biomass despite their degree of naturalness. In this work, such indicators of naturalness were assessed through chronological and retrospective tree-ring measurements (Di Filippo et al. 2017) that in some cases such indicators were lower in uneven-aged forests. Aside from the methodological and bioregional differences, an explanation for such differences may rely on the lack of interaction considered between naturalness and age measurements. Thus, the need to carry out further research including such interactions is still required due to the relatively younger ages recorded in our forest stands.
Uneven-sized forests may be configured by more than one canopy layer, and as such, the largest trees could develop larger root systems to access water located in the deeper soil layers (Forrester 2019). Within this conjecture, structurally diverse stands could present higher stocks of litterfall and root litter thanks to an optimization of the space occupied by leaves and root systems above and below ground. This would in turn promote higher carbon densities in the forest floor and upper soil layers (Crow et al. 2009;Vesterdal et al. 2013). Changes in forest structure, whether caused by natural dynamics and/or management activities, should not only impact carbon storage but should also modify processes driving carbon sink such as climate, tree species, floor climatic conditions and biotic activity and possible interactions between them (Lal 2005).
In consonance with the previous results regarding AGBC + BGBC, the forest age and structure heterogeneity coupling ( Table 2) was supported as a mediator of top soil SOC density variability, although this influence was slightly inferior in soils than in aboveand below-ground biomass compartment (Table 2; Fig. 3). These similarities between the fitted models with tree biomass and top soil likely rely on the idea that tree biomass is the main driver of top soil SOC in forests Martens et al. 2004;Mayer et al. 2020; but see Lal 2005). In this sense, the role of forest age has been widely addressed (Dybala et al. 2019;Jonard et al. 2017;Ma et al. 2020;Pregitzer and Euskirchen 2004;Nitsch et al. 2018, Smith 2004); however, the role of diversity has been documented in different forms. For instance, Conti and Díaz (2013) report that the diversity of plant functional traits has a positive influence over SOC. Other research has pointed to tree species identity, rather than species diversity, as being an effective driver of SOC density (Andivia et al. 2016;Dawud et al. 2016;Nitsch et al. 2018;Osei et al. 2021). These references both agree on how conifer forests tend to accumulate more carbon in soils due to the existence of thicker organic matter layers in the forest top floor formed by needle litter depositions. In addition, the decomposition rate of pine needles is slower than in broadleaf species, as the former present higher proportions of recalcitrant compounds (Fernández-Alonso et al. 2018). Other processes emerge as being significant when deeper soil horizons are considered, since carbon allocation is driven by different factors such as root activity and SOC transportation (Nitsch et al. 2018;Rumpel and Kögel-Knabner 2011). For instance, root litter decomposition is slower in the deepest layers as a result of less microbial activity, which leads to root litter carbon in soils having longer residence times (Pries et al. 2018). Moreover, fine root biomass tends to be high in mature forests that are diverse in species or can increase according to the proportional presence of conifers (Finér et al. 2017;Ma and Chen 2017). These findings may equally pertain to structurally diverse forest where larger trees are forced to develop profound roots as a consequence of sharing the soil with other trees (Aponte et al. 2020).
Short-term harvesting, clearcutting, dead wood raking, litter removal, species identity and tree size heterogeneity are all factors that have been identified as being relevant for controlling the size of soil and biomass carbon stocks (Aponte et al. 2020;Gimmi et al. 2013;Powers et al. 2011;Segura et al. 2017) and can be employed in forest silviculture prescriptions focused on carbon sequestration (Jandl et al. 2007). In this regard, sustainable practices such as the lengthening of rotation periods or the application of thinning to keep uneven-sized forests would be highly desired in plantations and managed forests to maximize carbon density. Thus, we underscore the need to promote tree size heterogeneity to enhance carbon stocks as mature forests are unable to efficiently store and densify carbon solely through ageing processes.

Conclusion
Given the primary role played by forests as reservoirs of terrestrial carbon, in the face of global change, we set out to quantify the effect of forest diversity drivers on the forest-age and carbon density relationship on two forest compartments. Our results provide correlational evidence that supports the idea that uneven-sized forests can significantly exceed the carbon stored on even-aged forests throughout the canopy ageing process. A more efficient use of the physical environment by trees both above and below ground will help forests to maximize carbon density during early successional stages. On the other hand, tree ontogeny is not sufficient to enhance the ecosystem carbon density, as mature forests must also develop into structurally diverse ecosystems. These finding supplement with previous studies that have shown how biomass productivity of forests is stimulated in structurally diverse stands. Bearing in mind the opportunities and benefits of sustainable forest management for the mitigation of rising atmospheric CO 2 concentration, it is worth noting the advantages of turning plantations to uneven-sized forests that can store carbon more efficiently.