Disentangling the Relationship between Tree Biomass Yield and Tree Diversity in Mediterranean Mixed Forests

: Tree biomass and the diversity relationship in mixed forest have an impact on forest ecosystem services provisions. Tree biomass yield is driven by several aspects such as species identity, site condition, stand density, tree age and tree diversity expressed as species mingling and structural diversity. By comparing diverse degrees of tree mixtures in natural forests, we can gain insight into the ecosystem services provision level and dynamic. Two monitoring sites in the Castilian Northern Plateau (Spain) have been analyzed to disentangle the relationships between biodiversity levels and tree biomass yield. Two permanent one hectare (ha) squared plots were established at Llano de San Marug á n and Valdepoza. In each plot, all individual trees were measured (diameter and height), georeferenced and its species identity deﬁned. Tree species in the two sites were Pinus sylvestris, Pinus nigra, Pinus pinea, Quercus pyrenaica, Quercus ilex, Quercus faginea and Juniperus thurifera . From these datasets, ten diversity indices that fall in three categories (species richness indices, species compositional/mingling indices and vertical structural indices) were used as predictor variables to ﬁt several candidate models. By merging the trees by site (without considering the species identity) selected models include individual tree basal area as an explanatory variable combining by addition or interaction with diversity indices. When species are analyzed independently, structural diversity impacts on biomass yield in combination (additive or multiplicative) with tree size is negative for Pinus nigra and positive for the other species.


Introduction
Mediterranean forests are characterized by a remarkable set of features that make them naturally and aesthetically attractive, but also quite fragile [1]. Mediterranean forestry is focused on a multi-functional approach, providing a wide range of goods and services for society ranging from products with high market value (fuelwood, cork, mushroom, pinecones, etc.) and non-market value ecosystem services (soil and landscape protection, water regulation, biodiversity conservation, carbon dioxide fixation, recreation, aesthetic view, etc.). The latter is more significant than their productive value, especially their significant role for carbon sequestration [2]. One of the notable characteristics of Mediterranean forests is their rich biodiversity, reflected by high genetic variability, exemplified by the large number of tree species in comparison to Nordic forests resulting from the survival of many conifer and broadleaf species during the glacial periods. Long-term exploitation (manipulation) of trees and forestland since ancient times is another feature of Mediterranean forest which has resulted in the dispersion of species such as Pinus pinea, Castanea sativa, species diversity, which is defined by the number of species and abundance of each species living within a certain area [20]. The species coexisting in a certain area are interconnected and dependent on one another for survival. While doing so, they perform important ecosystem functions and offer different ecosystem services for human life and society: provisional services (products obtained from ecosystem: many different type of food, fresh water, etc.); regulating services (the benefits obtained from the regulation of ecosystem processes: air quality and pollination); cultural services (the non-material benefits that people obtain: spiritual enrichment, recreation and aesthetic experiences); supporting service needed to maintain other services (i.e., photosynthesis and nutrient recycling). The provision of ecosystem for such goods and services depends basically on functions performed by living plants [21].
Two main mechanisms explain the reasons that biodiversity influences on productivity: selection effects and complementary effects. Different plant species in a mixture have different physiologies, morphologies and life history traits that might allow them to fully utilize limiting resources in different space and time compared to a monoculture of any species [21]. Complementarity between species traits can lead to over-or under-yielding in mixed forests in relation to monocultures. For instance, tree species that have different root morphologies occupy different soil profiles, which potentially allows them to exploit soil resources from different soil depths. Additionally, differences in shade tolerance or in foliage persistence during the year can lead to complementarity. However, it should be noted that these complementarities occur solely when co-existing species exhibit various forms of niche differentiation that allow them to capture resources in different space or time [22,23]. Another mechanism of productivity that diversity affects is selection effect (sampling effect) which describes species specific-effect on biomass: a greater productivity in more diverse communities is due to the most productive species which become dominant in the community due to competition. The likelihood of becoming a highly productive species increases as diversity increases. Thus, this leads to the increment in the total productivity of the community. Such considerations have led to the general perception of having higher productivity in an area where more plant species co-exist. Therefore, our main objective was to explore the relationship between tree and structural diversity and forest biomass in Mediterranean mixed forests. To achieve that, two intensive sample plots were analyzed in the Spanish northern plateau and all individual trees in each 1-ha plot were recorded for their species, measured and geopositioned. With these data, different biodiversity metrics and forest biomass values were calculated. Several lineal models were fitted, evaluated, and selected to provide an insight into the diversity and biomass relationship.

Materials and Methods
Two one ha squared plots (100 × 100 m) were established in 2015. Each plot is part of the regional marteloscopes network and are located in Valladolid and Palencia provinces (northern plateau in Castile and Leon region) with the latitude and longitude coordinates 41 • Figure 1). Both marteloscopes are in public forests managed by the regional forest service.
Marteloscopes are silvicultural training sites that, in our case (1 ha squared plots), are divided into sixteen permanent subplots, and hereafter referred to as quadrats, of 25 × 25 m length. Wooden stakes were placed with GPS subcentimetric to easily access the plot and its quadrants. Tree positions are determined with headings (in degrees) and distances, and they are identified with a number and measured (diameter and height). Trees at the boundary of the experimental plot are included in the analysis when the stem and at least 50% of the crown are within the plot (Figure 1). Different management scenarios and their consequences can be discussed in marteloscopes, thus supporting and improving the decision-making capacities in a variety of target groups including forest practitioners, forest owners, decision makers and scientists and students from different sectors. Within each quadrat, the locations of all trees were recorded in geographic coordinates and their species were identified and corresponding diameter at breast height (dbh, in cm) and, total height (in m). Both marteloscopes are dominated by pines in terms of basal area and biomass ( Figure 2) with oaks and juniper as associated species but the mixture degree in terms of trees per ha is more evenly distributed. The diameter at breast height (in cm) and total height (in m) for tree with dbh greater than 5 cm. To avoid interferences from past management practices, the sites were established in stands where no silvicultural treatments had been applied during the previous ten years.
50% of the crown are within the plot (Figure 1). Different management scenarios and their consequences can be discussed in marteloscopes, thus supporting and improving the decision-making capacities in a variety of target groups including forest practitioners, forest owners, decision makers and scientists and students from different sectors. Within each quadrat, the locations of all trees were recorded in geographic coordinates and their species were identified and corresponding diameter at breast height (dbh, in cm) and, total height (in m). Both marteloscopes are dominated by pines in terms of basal area and biomass ( Figure 2) with oaks and juniper as associated species but the mixture degree in terms of trees per ha is more evenly distributed. The diameter at breast height (in cm) and total height (in m) for tree with dbh greater than 5 cm. To avoid interferences from past management practices, the sites were established in stands where no silvicultural treatments had been applied during the previous ten years.  The Llano de San Marugán marteloscope was installed in the municipality of Portillo (Valladolid). The mixed forest is characterized by a plantation of Pinus pinea with resprouting of Quercus (Q. faginea and Q. ilex) and natural regeneration of Juniperus thurifera ( Figure 2). The mean annual temperature is 12 • C and the mean annual precipitation is 500 mm, falling irregularly throughout the year, with a minimum in the summer with a mean precipitation of 64 mm ( Figure 1). Fog is frequent in the long, cold winter, while summers are dry and hot with average temperatures around 30 • C. The altitude of the experimental site is 853 m. Köppen classification is defined as warm-summer Mediterranean climate (Csb). Following the World Reference Base for Soil Resources classification (WRB), the soil belongs to the category of leptosols, which are dominant and associated with regosols [25]. Leptosols comprise very thin soils over continuous rock and soils that are extremely rich in coarse fragments, while regosols are weakly developed mineral soils in unconsolidated materials that do not have a mollic or umbric horizon [26].
Valdepoza marteloscope installed in the municipality of Saldaña (Palencia), presents the species Pinus sylvestris and Pinus nigra from reforestation during the 1950s and Quercus pyrenaica as coppice ( Figure 2). The mean annual temperature is 10 • C and the mean annual precipitation is 640 mm, with a mean precipitation of 92 mm in summer ( Figure 1). The altitude of the experimental site is 1027 m. Köppen classification is defined as temperate oceanic climate (Cfb). Following WRB, the soil belongs to category of cambisols [25], which are soils with at least an incipient subsurface soil formation [26]. The Llano de San Marugán marteloscope was installed in the municipality of Portillo (Valladolid). The mixed forest is characterized by a plantation of Pinus pinea with resprouting of Quercus (Q. faginea and Q. ilex) and natural regeneration of Juniperus thurifera ( Figure 2). The mean annual temperature is 12 °C and the mean annual precipitation is 500 mm, falling irregularly throughout the year, with a minimum in the summer with a mean precipitation of 64 mm ( Figure 1). Fog is frequent in the long, cold winter, while summers are dry and hot with average temperatures around 30 °C. The altitude of the experimental site is 853 m. Köppen classification is defined as warmsummer Mediterranean climate (Csb). Following the World Reference Base for Soil Resources classification (WRB), the soil belongs to the category of leptosols, which are dominant and associated with regosols [25]. Leptosols comprise very thin soils over continuous rock and soils that are extremely rich in coarse fragments, while regosols are weakly developed mineral soils in unconsolidated materials that do not have a mollic or umbric horizon [26].  [24]) representing the percentage of species mixture in terms of biomass, basal area, and number of trees per hectare. The representation in the axis is the identity of the species within the marteloscope: Pinus pinea, Juniperus thurifera, Quercus (mixture of Quercus ilex and Quercus faginea), Pinus sylvestris, Quercus pyrenaica and Pinus nigra. Color of each of the points (quadrats) corresponds to biomass (Mg/ha), basal area (m2/ha) and number of trees per hectare, according to the legend.

Biomass Estimation
Biomass in the different components of the tree, e.g., stem and branches (thick, medium and thin with needles or leaves) and roots were calculated from diameter at breast height (dbh) and total height (h) using existing allometric equations for softwood species [27] and for the hardwood species [28]. From stem diameter at dbh and h, the total biomass value is estimated by the sum of the biomass of the stem, thick branch (diameter larger than 7 cm), medium branch (diameter between 2 and 7 cm) and thin branch (diameter smaller than 2 cm) with needles or leaves. Biomass data were transformed from kilograms to tons. Tree component biomass values were computed for individual trees within each quadrat and summed up to derive a summary of tree biomass for each species and marteloscope (Tables 1 and 2).
Average variable values are presented for stand level: B-total above-ground biomass (kg/ha); d-mean breast height diameter (m); h-mean total height (m), V-mean tree volume (m 3 /ha), G-total basal area (m 2 /ha), and mean tree value of: Sm-Simpson index; Sn-Shannon index; E-evenness index; D-Berger-Parker index, Mi-mingling index; MS-spatial diversity status; W-uniform angle index; S-segregation index; A-vertical profile index; TH-height differentiation index.

Tree Species Diversity Estimation
A forest is a 3-dimensional system whose structure is a key element in ecosystem function and biological diversity, by regulating resources related to forest function (light, water, and soil nutrients supply, capture, and use), intra-and inter-specific interactions [2,29], regeneration patterns, consequent self-thinning and past and present disturbance events [30,31]. Stand structural diversity leads to increases in species richness and contributes to forest stability and integrity [32]. Stand structural diversity combines the concepts of species richness (diversity), species composition (mixture), and spatial diversity (tree positioning) and size differentiation [33]. Accordingly, three distinct types of stand structural indices and methods have frequently been purposed in the preceding literature for explaining the influence of stand structural diversity on the productivity and function of forest stand: (I) species richness-Simpson index [34], Shannon index [35], Berger-Parker index [36] and evenness index [37]; (II) species composition indices-mingling index [38], spatial diversity status [39] and segregation index [40]; (III) tree distributional indices including horizontal and vertical patterns and size differentiation, uniform angle index [41], vertical species profile [42], and height differentiation index [43]. Since forest structure is determined in three dimensions, it is appropriate to analyze the effect of tree diversity on biomass by the metrics that can fully address the 3-dimensionality of mixed forest structures.
The diversity indices used in this study were classified into three categories: (I) species richness (Simpson index, Sm; Shannon index, Sn; evenness index, E; Berger-Parker index, D), (II) species composition (mingling, Mi; spatial diversity status, MS; segregation index, S) and (III) vertical (vertical profile index, A; height differentiation index, TH) and horizontal distribution pattern (uniform angle index, W).The basic idea behind the diversity indices is to obtain a quantitative estimate of biological variability that can be used to compare biological entities, composed of discrete components, in space or in time [44]. In this study, we extend this definition to also include structural variability in terms of relative size and spatial organization of trees. Average variable values are presented for tree species level: N-number of trees, B-total above-ground biomass (kg/ha); d-mean tree diameter at breast height (cm); h-mean total height (m); v-mean volume (m3/per tree); g-mean basal area (m2/per tree) and tree mean  The Simpson diversity index (Sm) [34] accounts for the number of species and the abundance of each species. The index represents the probability that two individuals randomly selected from a sample belong to different species; therefore, this index increases with species diversity. The Shannon index (Sn) [35] considers the abundance and evenness of the species. When all species are equally common, the value increases until the logarithm of the number of species in the sample, ln (S). On the contrary, if all abundance is concentrated in one species, and the other species are very rare, its value reduces to 0. Species evenness (E) [37] expresses how equally the individuals in the community are distributed over the different species. Low values indicate that one or a few species dominate, and high values indicate that relatively equal numbers of individuals belong to each species. The Berger-Parker index (D) [36] is a measure of the numerical importance of the most abundant species in the population. The reciprocal of the index is frequently used, so that an increase in number explains a higher diversity and thus a reduced dominance.
Species spatial mingling index (Mi) describes the proportion of neighbor trees which belong to different species identity as the reference tree [38]. A high degree of mingling means that trees are surrounded by different species. However, the number of different species in the structure unit was not considered, and this was a shortcoming of this index. This limitation has fulfilled the spatial diversity status index (MS) [2], which considers not only the spatial mingling, but also the number of tree species [39]. The reference tree of a frequent species is more likely to have the neighbors of the same species, reflecting a low index value. On the contrary, rare species have less probability to have the same neighbor species, resulting in high value of this index. Thus, spatial diversity status index is considered as a sensitive index for rare species. Segregation index (S) developed by Pielou [40] describes the degree of intermingling of two species groups [45] based on the nearest-neighbor method. Values greater than 0 and below 0 indicate a trend towards segregation and association, respectively, while an independent distribution of species will reach values close to 0 [2].
The vertical profile index (A) [42] calculation is based on the diversity index of Shannon [35]. 'A' considers both proportion of the species within a stand and the presence of each species in different height zones [46] and rises as the vertical profile increases in heterogeneity. The height differentiation index (TH) [43] measures the variability in the vertical size between the reference tree and neighboring trees, reflecting high differentiation when reaching values close to 1. Increasing values in the uniform angle index (W) indicate a transition from regular to random and then to clumped spatial pattern and define the degree of spatial dispersion of nearest neighbors around the reference tree based on angles [41].

Statistical Analysis
Different models representing different plausible hypothesis were fitted to explain how diversity indices impact on individual tree biomass. The three general model structures tested were as follows in Equations (1)- (3): where x i and y i represent different explanatory variables and β i represents the estimated parameters. Parameters were estimated by ordinary least squares (OLS) multiple linear regression. When needed, the models were transformed by logarithms. Basic assumptions of multiple linear regression (linearity, normality of the distributions of the residuals, homoscedasticity and independence of variables) were tested. The models were defined for explanatory use and not for prediction. Eleven predictor variables were tested: four species richness indices (Sm, Sn, D, E), three species composition indices (Mi, MS, S), four spatial distribution indices (A, TH, W) and the individual basal area per tree (Gi, m 2 ). Different variables combinations, allowing interactions between variables, were tested (Table 3) leading to 623 alternative models for each site by species and for the whole tree datasets (without considering species identity) for each marteloscope. Table 3. Models tested including its general structure and the predictor variables.

# Fitted Models
Predictor Variable Multivariate models with 3 predictors (188 alternative models) 10. To select the best model for dominant tree species and for each marteloscope, five criteria were employed: (I) significance of variables (in the ANOVA analysis, an effect is considered to be significant when its coefficients have a probability less than or equal to the significant probability (p < 0.001), (II) biological meaning of parameters signs and values, (III) graphical distribution of the residuals (normal Q-Q plot and predicted versus adjusted plot), (IV) the adjusted R squared and (V) Akaike information criteria (AIC). In addition, we checked collinearities between the variables (biodiversity metrics and basal area) [47] where higher collinearities were found in both sites for the Simpson index and the Shannon index, Mi and MS and G and B (see Supplementary Materials, Figure S1). All statistical modeling was carried out with R software version 3.6.2 and tidyverse package [48].

Results
Relationships among biodiversity metrics and forest biomass were broadly examined at stand and species level across the two mixed forest sites (the best three significant models per each site are presented at Table 4). In the selected models, the response variable is the natural logarithm of the biomass in megagrams per tree (Ln B) while species spatial mingling index (Mi), segregation index (S), uniform angle index (W) and height differentiation index (TH), modified in some cases by Gi (individual tree basal area), are the explanatory variables. These models selected according to the five criteria, are shown in Table 4 using both the Akaike information criteria (AIC) and the adjusted R 2 to determine which model was the most parsimonious while simultaneously explaining most of the information in the data. The significance level is p < 0.001. Selected models for the whole plots (not species identity considered) are codified by site and the number (i.e, for Llano de San Marugán are M1, M2 and M3 while for the Valdepoza site are V1, V2 and V3). All the selected models were checked against deviations that prevent the use of multiple linear regressions (see Figure S2 for details). Selected models for Llano de San Marugán (Table 4), three linear models have been highlighted (structure as shown in Equation (3)): M1 including as explanatory variables Gi and W, M2 including as explanatory variables Gi and Mi and M3 including as explanatory variables Gi and TH. Between these three models, M1 was found to be the best model for explaining tree biomass in this site, showing the positive influence of the interaction of individual tree basal area (Gi) and the spatial distribution pattern with the uniform angle index (W) combined (transformed by logarithms) in one variable (parameter equal to 0.936). For the Valdepoza site, the three selected models are V1 (structure as shown in eq 1) including as explanatory variables G and Mi, V2 (structure as shown in eq 1) including as explanatory variables Gi and S and V3 (structure as shown in eq 3) including as explanatory variables Gi combined with Mi. In this case, V2 was found to be the best model for understanding the variables' relationships in Valdepoza, indicating a positive relationship with tree size expressed by Gi (parameter equal to 1.057) and a negative relationship with the degree of intermingling with other species expressed by the segregation index (S) (parameter equal to −2.398).
The species analysis (Table 5 for Llano de San Marugán and Table 6 for Valdepoza sites) indicates that biomass (transformed by natural logarithms) is closely related with individual tree size (expressed by Gi) for all the species but modified by diversity indices (mainly vertical and horizontal structure and compositional but also diversity expressed by Simpson index). When the relationship between variables was sufficiently strong, three alternative models were selected. For Pinus pinea, the relation was so weak that no model was selected. For the other species in the Llano de San Marugán site, the selected models (Table 5) included, beside the tree size (Gi) individually or modifying the other variables, the height differentiation index, TH, (Quercus faginea), species segregation, S, (Quercus ilex) and uniform angle index, W, expressing the horizontal distribution pattern (Juniperus thurifera). We can observe that structural diversity indexes (species composition, i.e S index and distribution pattern, i.e., W and TH) together with tree basal area (Gi) were the strongest independent predictors of forest biomass in Llano de San Marugán and for the species Quercus faginea (parameter equal to 0.360), Quercus ilex (parameter equal to 1.570) and Juniperus thurifera (parameter equal to 0.582), whereas relationships for Pinus pinea were comparatively weak in terms of fulfilling the criteria of the distribution of the residuals (predicted versus adjusted plot). Distributions of the residuals of the selected models are exposed in Figure S3.    The more humid site of Valdepoza (Table 6) shows the main results for Pinus sylvestris, Pinus nigra and Quercus pyrenaica. In this plot, for the pines, four models were selected as candidates because all of them were very close in terms of normality assessment but the evaluation criteria allowed us to select the best one. The selected models for the Valdepoza site always included tree size (Gi) in combination with structural diversity indices (height structure, TH; for Quercus pyrenaica (parameter equal to 0.433) and segregation index, S, for Pinus sylvestris (parameters equal to 1.145 for Gi and to −0.225 for S; this trend is the same as that shown by the model selected by the Valdepoza site) or the Simpson index for Pinus nigra (parameter equal to 6.830 for Gi and to −5.190 for Sm). All the selected models were checked against deviations that prevent the use of multiple linear regressions (see Figures S4-S6 for details). Only the models for Quercus pyrenaica show a slight heteroskedasticity ( Figure S4); therefore, those models must be analyzed with caution.

Discussion
Our main findings supported the importance of structural and spatial diversity (see Table 4 for models not considering specie identity), in combination with tree identity and size (see selected models in Tables 5 and 6 when species identity is considered) to assess biomass yield in Mediterranean mixed forests. Individual basal area, as a tree size proxy, is included in all the models both by species and without considering species identity. Our results agree with the findings from Bohn and Huth [30], who examined the influence of species diversity and forest structure on aboveground biomass over a broad range of forest stands and found out a positive relation between forest structural diversity and forest productivity. Focusing on absorptive root biomass the contrary has been reported in a wide geographical transect covering boreal, hemiboreal, mountainous and thermophilous dediduous forests in Europe [49]. Aponte et al. [50] found that forest structure diversity is the key variable to explain carbon storage in Australian forests while species-and functional-diversity indices hold a weak impact. Interestingly, Hime and Puettman [51] found no influence of tree mixture on biomass while the understory diversity was positively associated with biomass. Previously, our team analyzing pine mixed forests in northern Spain [4] also found the positive impact of species combination on aboveground forest yield. Size heterogeneity enables bigger trees to obtain a greater amount of a certain resource and use them more efficiently than small trees [29], while species complementarity due to differential traits leads to a more efficient resource use [52,53]. These facts are reflected in the models for our drier site (Llano de San Marugán) where the mixture of Pinus pinea, Quercus faginea, Quercus ilex and Juniperus thurifera create a multilayered open canopy strata where the top layer is occupied by Pinus pinea, enabling the light-demanding species (pine) to capture more light and grow in a superior way than the other cohabitant species and become dominant in the stand. As competition (negative interaction) and facilitation (positive interaction) occurs simultaneously, through direct field observation we are only able to detect the net effect of interactions. That means that cohabitant species can simultaneously facilitate and compete, but we only can observe the joint (positive or negative) outcome. This limitation of direct observation results must be considered when interspecific and intraspecific tree interaction is analyzed. Additionally, facilitation can represent a null or negative effect on the benefactor tree [54] contributing to the negative outcome that can be observed in some observational studies.
The positive relation between productivity (biomass accumulation) and tree diversity largely depends on presenting highly productive species in multispecific communities [23]. In our case, this role of highly productive species is played by pine (Pinus pinea in Llano de San Marugán and Pinus sylvestris and Pinus nigra in Valdepoza). It is important to recall the moderate productivity of the analyzed pine species here; the concept of highly productive species must be considered by comparison with the oak species. Thus, the outcome obtained is lower than when more contrasting species are present in the mixture. In both sites, pine species also facilitate the development of a low layer of oaks. In areas where pine density is lower, oaks can dominate at the microsite level because pines may facilitate oak development by increasing seed protections and enhancing the habitat to facilitate seedling establishment. Although there is an opposition between suitable conditions for recruitment and suitable conditions for further oak growth and development after seedlings and saplings stages [55], the open conditions in the mixed pinewoods are adequate sites for the development of oak patches that enhance the mixture by including species with differential traits. As root development has not been considered, it is possible that oaks trees' belowground biomass accumulation overtakes the pine trees' values. Limited research on mixed forest root systems has been conducted (with exceptions such as [49]) so further studies are needed to gain insight into the effect of mixing on belowground tree biomass.
Moreover, the biomass of individual trees in a given stand is not only a reflection of diversity (species richness, composition and structure) but also various internal and external factors such as age, stand density, site productivity, competition at the tree level, climate, soil (texture, moisture content), geographical location, and length of growing season [56,57] or tree genetics which might not be reflected by the metrics considered in this analysis.

Conclusions
After thoughtful examination of 623 models (developed to represent alternative hypotheses) with 10 predictor tree diversity indices (grouped as proxies of species richness, species composition and size distribution pattern) plus a tree size proxy (individual basal area), we found out that there is a relation between tree biomass and diversity, although this relation varies among the analyzed species. As expected, all the selected models include individual tree basal area as explanatory variables combining by addition or interaction with diversity indices. When the data of each site is analyzed without considering the species identity, the spatial dispersion results in a positive outcome in terms of biomass in our drier site (Llano de San Marugán) exposing the importance of the effect site selection in the dominance of the two similar species (Quercus faginea and Quercus ilex) and thus more productive. On the mesic site (Valdepoza) the spatial segregation shows a negative effect on yield that can be explained by the cohabitation of two pines with similar traits (Pinus sylvestris and Pinus nigra). Structural diversity shows a positive impact on biomass yield in combination (additive or multiplicative) with tree size. Only in Pinus nigra is structural diversity (segregation index) negatively related with biomass. Operational forestry by selection thinning can modify structural and specific diversity impacting on biomass yield and carbon sequestration. Local studies, where a careful analysis of the species composition, structural and spatial arrangement jointly with site conditions is performed, are needed to orientate operational forestry to obtain diverse ecosystem services under multifunctional forest management.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12070848/s1. Figure S1. Correlation matrix using corrplot package in R [47]. The correlation number is presented at significance level 0.05 for both experimental sites. Figure S2. Residual graphs are referred in the first column to Llano de San Marugán (models M1, M2 and M3), and in the second column to Valdepoza (models V1, V2 and V3). On the left we examine the predicted versus residual plot with the predicted values on the x-axis and the residuals on the y-axis. On the right, quantile-quantile plot (q-q plot). Notice the x-axis for the theoretical quantiles and the y-axis for sample quantiles. Figure S3. Residual graphs are referred in the first column to Quercus faginea (Qf1, Qf 2. Qf 3), the second column to Quercus ilex (Qi1, Qi 2. Qi 3), and the third to Juniperus thurifera (Jt1, Jt 2. Jt 3). On the left we examine the predicted versus residual plot with the predicted values on the x-axis and the residuals on the y-axis. On the right, quantile-quantile plot (q-q plot). Notice the x-axis for the theoretical quantiles and the y-axis for sample quantiles. Figure S4. Residual graphs of Quercus pyrenaica models (Qp1, Qp2. Qp3). On the left, we examine the predicted versus residual plot with the predicted values on the x-axis and the residuals on the y-axis. On the right, quantile-quantile plot (q-q plot). Notice the x-axis for the theoretical quantiles and the y-axis for sample quantiles. Figure S5. Residual graphs of Pinus sylvestris models (Ps1, Ps2. Ps3, Ps4). On the left we examine the predicted versus residual plot with the predicted values on the x-axis and the residuals on the y-axis. On the right, quantile-quantile plot (q-q plot). Notice the x-axis for the theoretical quantiles and the y-axis for sample quantiles. Figure S6. Residual graphs of Pinus nigra models (Pn1, Pn2. Pn3, Pn4). On the left we examine the predicted versus residual plot with the predicted values on the x-axis and the residuals on the y-axis. On the right, quantile-quantile plot (q-q plot). Notice the x-axis for the theoretical quantiles and the y-axis for sample quantiles. Funding: This paper has been funded partly by the projects PCIN-2017-027 (Mixed species forest management. Lowering risk, increasing resilience, REFORM: REsilience of FORest Mixtures funded by the Spanish Ministry of Science under the EU ERA-NET Sumforest: Sustainable forest management, multifunctional forestry, European Forest Policy) and COMFOR-SUDOE: Integrated and intelligent management of complex forests and mixed-species plantations in Southwest Europe (Project Code: SOE4/PA/E1012) through the European Regional Development Fund (ERDF) Narangarav Dugarsuren was funded by the EJMD MEDFOR (Erasmus +) The APC was funded by Universidad de Valladolid trough COMFOR-SUDOE funds.

Conflicts of Interest:
The authors declare no conflict of interest.