Tree—Open Grassland Structure and Composition Drive Greenhouse Gas Exchange in Holm Oak Meadows of the Iberian Peninsula

: Iberian holm oak meadows are savannah-like ecosystems that result from traditional silvo-pastoral practices. However, such traditional uses are declining, driving changes in the typical tree—open grassland structure of these systems. Yet, there are no studies integrating the whole ecosystem—including the arboreal and the herbaceous layer—as drivers of greenhouse gas (GHG: CO 2 , CH 4 and N 2 O) dynamics. Here, we aimed at integrating the inﬂuence of tree canopies and interactions among plant functional types (PFT: grasses, forbs, and legumes) of the herbaceous layer as GHG exchange drivers. For that purpose, we performed chamber-based GHG surveys in plots dominated by representative canopy types of Iberian holm oak meadows, including Quercus species and Pinus pinea stands, the last a common tree plantation replacing traditional stands, and unraveled GHG drivers through a diversity-interaction model approach. Our results show the tree–open grassland structure, especially drove CO 2 and N 2 O ﬂuxes, with higher emissions under the canopy than in the open grassland. Emissions under P. pinea canopies are higher than those under Quercus species. In addition, the inclusion of diversity and compositional terms of the herbaceous layer improve the explained variability, with legumes enhancing CO 2 uptake and N 2 O emissions. Changes in the tree cover and tree species composition, in combination with changes in the structure and composition of the herbaceous layer, will imply deep changes in the GHG exchange of Iberian holm oak meadows. These results may provide some guidelines to perform better management strategies of this vast but vulnerable ecosystem.


Introduction
Holm oak meadows, also called dehesas in Spain and montados in Portugal, are seminatural savannah-like agroecosystems that result from the thinning of the Mediterranean forest, in which an herbaceous and an arboreal (mostly Quercus species) layer coexist.They are one of the largest agroforestry systems in Europe [1], covering 3.5-4 million ha, mostly along the South-West of the Iberian Peninsula [2], and are also present in other world regions with Mediterranean climate, mainly in California [3][4][5].In addition, Mediterranean savannahs, with trees from other taxonomic groups, can be found in South Africa, south-western Australia, and central Chile [4]; and savannah-like ecosystems, in a broad sense, can be found worldwide, including mesquite grasslands in the United States, savanna ecosystems in Africa, and tropical savannas in Brazil (called cerrado).All of them, being ecosystems of high cultural, economic, and ecological value, are able to provide higher ecosystem services than open grasslands or field crops growing under similar conditions [6].
Holm oak meadows, in particular, have traditionally provided a wide variety of goods and services, including pasture for livestock, acorns, timber, and cork.However, such traditional uses that have shaped holm oak meadows into a matrix of trees and open grassland are changing, with the consequent implications that this may have on ecosystem functioning.Extensive grazing is declining towards intensive farming; plantations of fast-growing trees, mostly Eucalyptus and Pinus species; shrub encroachment due to land abandonment; and there is a worrying lack of tree regeneration [7,8].
Hence, although the canopy influence has been described to some extent on soil [9][10][11][12], water fluxes [13], evapotranspiration [14], vegetation structure [15][16][17], and composition [4,18,19]; many fewer studies have assessed the influence of the tree-open grassland structure-on greenhouse gas (GHG: CO 2 , CH 4 , and N 2 O) exchange.In this sense, studies conducted in Iberian holm oak meadows [17,[20][21][22][23][24], as well as on other savannah-like ecosystems [25][26][27][28], have highlighted the relevance of the ecosystem structure as drivers of CO 2 fluxes.However, some authors have described an enhanced soil respiration under the canopy compared to the open grassland (i.e., Tang and Baldocchi, 2005 [24]; Uribe et al., 2015 [20]), related to a higher soil C and N content, despite lower soil temperature, while others reported increased CO 2 exchange rates in the open grassland, due to higher herbaceous biomass and light availability, as main drivers of CO 2 uptake, and due to higher soil temperature, as main driver of CO 2 release [17].
On the other hand, the only previous studies addressing CH 4 and N 2 O exchange in Iberian holm oak meadows were conducted by Shvaleva et al., 2014 [29] and 2015 [30], which related CH 4 and N 2 O emissions to soil water content, but the effect of soil water content (positive or negative) was dependent on the canopy [29].In addition, CH 4 and N 2 O assessments from other savannah-like ecosystems have identified the environment under the canopy as a possible source of N 2 O emissions [27,31].Yet, CH 4 and N 2 O assessments from savannah-like ecosystems, especially in Mediterranean environments, are still very scarce [32], and there is a lack of studies integrating the whole ecosystem structure and composition, combining the arboreal and the herbaceous layer (both profitable resources), and assessing how both interact to drive GHG exchange.Accordingly, it must be considered that the vegetation of the herbaceous layer of holm oak meadows is highly diverse [33], and sorting the wide variety of species of the herbaceous layer into plant functional types (PFT) may provide a mechanistic link between diversity and ecosystem functioning [34].Indeed, most common (PFT) of the herbaceous layer, including grasses, non-legume forbs (hereafter "forbs"), and legume forbs (hereafter "legumes"), have different nitrogen (N) and light (and, therefore, CO 2 ) use and acquisition strategies, which may result in different responses to tree canopies in their GHG exchange.In addition, most diversity-GHG studies have focused on the effect of specific or functional diversity on GHG dynamics [35,36], but much fewer have disentangled identity and interaction effects among PFT [37].
In the present study, we aimed at integrating the influence of trees and the herbaceous layer structure and composition as GHG exchange drivers.In particular, we aimed at (1) assessing the canopy effect under representative canopy types of Iberian holm oak meadows, including traditional Quercus species and Pinus pinea L. stands, the last a common tree plantation replacing traditional canopies; and (2) unraveling the influence of the main PFT of the herbaceous layer, using a diversity-interaction modeling approach.

Study Sites and Sampling Design
Field work was carried out in spring (05/04/2016-10/04/2016) and autumn (13/12/-17/12/2016), coinciding with the most productive moments of the system, to capture seasonal variability of the studied variables and effects that may be season dependent.Study plots were distributed in two locations in the South-West of the Iberian Peninsula: Doñana Natural Park (DN, 37 • 15 34" N, 6 • 19 55" W, 30 m a.s.l.) and Sierra Morena mountains (SM, 37 • 39 50" N, 5 • 56 20" W, 296 m a.s.l.).Both locations have Mediterranean climate regime with warm, dry summers, and mild winters [38].Mean annual temperature in DN is 18.1 • C and in SM is 16.8 • C, and mean annual precipitation in DN is 543 mm and in SM is 648 mm.Grassland in both locations is dominated by herbaceous annual species, including grasses, non-legume forbs (hereafter "forbs"), and legume forbs (hereafter "legumes").Both locations are extensively grazed at similar stocking rates: DN grazed by goat and cattle (0.40 livestock units (LSU) ha −1 ), and SM by cattle and Iberian pigs (0.36 LSU ha −1 ).
To characterize soil properties of the study plots, soil samples were extracted and analyzed according to standard methods (Table 1).Texture in the SM-ilex plot varied from sandy clay loam (0-40 cm depth) to clay (40-80 cm depth).All DN soils had a sandy clay loam texture, except the deeper layer of the DN-suber plot, which was sandy loam.The SM-ilex plot had a slightly acid pH and DN plots had a neutral-basic pH.Organic C was very low in all the plots (Table 1), although the value in the first 30 cm of the DN-pinea plot was markedly above the average.Total N was in general also quite low (Table 1).

Plot
Tree individuals of the UC treatment of each species were selected with a similar diameter at breast height (DBH, Q. ilex 0.43 ± 0.03, Q. suber 0.63 ± 0.03 and P. pinea 0.57 ± 0.06 m).In addition, sampling points of the UC treatment were always placed at 1 m distance from the selected tree trunk, and sampling points of the OG treatment were placed at a minimal distance of 3 m from the selected tree, clearly outside the canopy, to minimize the influence of stem water flow.Sampling points were systematically placed following the north orientation with respect to the tree trunks.
For each treatment level, we selected 3-4 replicates, totaling 73 sampling points.In the DN-mixed plot, we discriminated between both Quercus species (Q.suber and Q. ilex) to establish sampling points.However, preliminary comparative analysis in the DN-mixed plot on environmental and vegetation characteristics under the canopy of both Quercus species indicated no relevant differences.DN-mixed plot results are then always presented, combining both tree species.
At each sampling point, we hammered a metal collar, which was necessary for measuring GHG exchange (Section 2.2), and that defined the area in which vegetation (Section 2.3) and soil (Section 2.4) were sampled.

Greenhouse Gas Exchange (GHG) Measurements
To measure GHG exchange, including CO 2 , CH 4 , and N 2 O, we used a portable gasexchange system, consisting of a cylindrical chamber (volume = 0.019 m 3 ), connected to a photoacoustic spectroscopy gas analyzer (PAS, INNOVA 1412, LumaSense Technologies, Denmark; see further system set-up details in Debouk et al., 2018).PAS nominal detection limit is 5, 0.24, and 0.03 ppm for CO 2 , CH 4 , and N 2 O, respectively [42].PAS was calibrated prior to field campaigns by the vendor in the customary way [43].
The cylindrical chamber was secured to the ground by fixing it to metal collars (h = 8 cm, diameter = 25 cm), previously hammered into the ground (3 cm deep) at least 24 h before each measurement, to limit disturbances to the soil-vegetation system.In addition, we inspected any possible disturbance on the soil that could have caused an increase on GHG emissions, but no disturbance was detected, neither during the current survey nor in previous PAS surveys [44].
Flux measurements were done placing the chamber consecutively over the metal collars for five minutes.Measurements were done over intact vegetation, at light and dark (by covering the chamber) conditions.Afterwards, the vegetation inside the metal collar was harvested, and bare soil measurements at dark conditions were also performed, at least two hours after harvesting.As a result, 219 flux measurements (considering all treatment levels and replicates-73 sampling points-×3 different conditions) were recorded.
In the case of CO 2 , resulting fluxes measured over vegetation at light conditions can be approximated as the understory net ecosystem CO 2 exchange (NEE); over vegetation at dark conditions can be approximated as ecosystem respiration (R eco ); and, on bare soil at dark conditions, can be approximated as soil respiration (R soil ).In the case of CH 4 and N 2 O, we did not detect significant differences among measuring conditions.Therefore, we present the average value considering all three measuring conditions [44].
Flux measurements, flux calculation, and data quality checks were done according to Debouk et al. (2018).This included a fitting goodness assessment based on the R Adj 2 value (fluxes below a R Adj 2 of 0.8 for CO 2 , and below 0.2 for CH 4 and N 2 O, were excluded), and filtering fluxes below flux detection limit, calculated as the standard deviation of the ambient concentration over the measuring time.Negative values refer to the flux from the atmosphere to the biosphere and positive values correspond to the flux from the biosphere to the atmosphere, according to the micrometeorological sign convention [45].

Vegetation Sampling
After GHG measurements were made, we harvested the vegetation at ground level at each sampling point.Thereafter, in the laboratory, we separated the vegetation into aboveground biomass (AGB) and litter (dead plant material detached from the herbaceous vegetation and tree leaves on soil surface).All fractions are summarized in Figure S2.In addition, we separated the AGB into PFT (forbs, grasses, and legumes, summarized in Figure S3).Vegetation was oven dried at 60 • C until constant weight.

Belowground Biomass Sampling and Soil Water Content Determination
Two soil cores of 9 cm 2 surface and 0-10 cm depth were extracted at each sampling point.In the laboratory, one of the cores was washed and filtered with a 0.2 mm pore size strainer to obtain belowground biomass (BGB, summarized in Figure S2).The second core was used for SWC determination by gravimetric method, as the difference between fresh and dry soil weight.Both, BGB and soil samples were oven dried at 60 • C until constant weight.

Data Analysis: Greenhouse Gas Exchange Modeling
To assess main GHG (CO 2 , CH 4 and N 2 O) drivers, especially focusing on the influence of trees and the herbaceous layer structure and composition, we run a diversity-interaction modeling [46,47].The modeling compares a null model, in which the response variable is not affected by plant diversity and/or composition, to models that address diversity and composition at different levels.In our study, we compared the null model (Equation ( 1)), in which the corresponding GHG depended only on treatment variables, including plot, season, and canopy; environmental variables, including PAR (µmol photons m −2 s −1 ), temperature (T, • C), and SWC (fraction); and structural components of the herbaceous layer, including AGB, litter, and BGB (g DW m −2 ): Null model to models that included PFT composition and diversity of the herbaceous layer in different ways: (a) Identity model, which includes PFT identity effects, meaning the biomass proportion of each PFT (Equation ( 2)), where P indicates the proportion of the given PFT and the sub-index F indicates forbs, G grasses, and L legumes, respectively: Identity model (b) Average interaction model, which includes PFT identity effects, plus evenness [46] as an average interaction term (Equation ( 3)).Evenness (E) calculated as E = 2PFT (PFT−1) PFT ∑ i<j P i P j , where PFT is the number of PFT present in the community, and P i the relative abundance of the PFTs.Evenness lies between 0 for mono-PFT plots, and 1 when all PFT are equally represented.Average interaction model (c) Specific interaction model, which includes specific interactions between PFT in addition to the identity effects (Equation ( 4)): Specific interaction model The models were run without intercept to test the effect of all three PFT at the same time [46].In addition, we modeled all microclimatic sampling conditions (PAR, T s , and SWC, summarized in Table S1) and vegetation fractions (AGB, litter, BGB, and PFT, summarized in Table S2), as function of plot, season, and canopy.To compare and determine the final model in each case, we used several methods, including Akaike information criteria (AIC), the adjusted determination coefficient R 2 (R 2 adj ), and model comparison by analysis of variance (ANOVA) to test significant differences between models, using the "anova" function in R. The ANOVA comparison between models tested whether the reduction in the residual sum of squares was statistically significant [48].The most explicative and parsimonious model of each GHG is shown and discussed.The modeling of microclimatic sampling conditions and vegetation fractions was used to interpret and discuss our results.

Results
Tree canopies drove microclimatic conditions (Figure S1 and Table S1), as well as the structure (Figure S2) and composition (Figure S3 and Table S2) of the herbaceous layer.Factors, in combination with season, interacted among them to drive GHG fluxes.PAR and T s decreased under the canopy compared to the open grassland, being this difference less marked in autumn than in spring (Figure S1 and Table S1).The structure and composition of the herbaceous layer was also dependent on the presence of tree canopies.The green fraction (AGB) slightly decreased under the canopy compared to the open grassland in all DN plots (Figure S2 and Table S2).Reduction that was mostly due to a change in the PFT composition, with the biomass of forbs and legumes decreasing under the canopy compared to the open grassland (Figure S3 and Table S2).On the contrary, such change in the AGB between the under the canopy and the open grassland was not noticeable in the SM-ilex plot (Figure S2 and Table S2).The litter fraction was markedly higher under the canopy than in the open grassland, although this difference almost disappeared in autumn (Figure S2 and Table S2).

CO 2 Exchange
CO 2 net uptake was similar in spring and autumn (NEE, Figure 1A), while CO 2 emissions (R eco and R soil , Figure 1A) were lower in autumn than in spring (season effect, Table 2).The SM-ilex plot showed the highest net CO 2 uptake rates (NEE, Figure 1A) and the highest R eco rates in the open grassland (Figure 1A), while DN-pinea was the plot with the highest CO 2 emissions under the canopy (R eco and R soil , Figure 1A).Thus, the canopy itself had a strong influence over CO2 fluxes.Generally, NEE under the canopy was dominated by CO2 emissions (canopy effect, t = 5.58, p < 0.001, Table 2),  Thus, the canopy itself had a strong influence over CO 2 fluxes.Generally, NEE under the canopy was dominated by CO 2 emissions (canopy effect, t = 5.58, p < 0.001, Table 2), NEE being strongly driven by PAR (PAR effect, t = −4.58,p < 0.001, Table 2).The exception was the SM-ilex plot, where there was net CO 2 uptake under the canopy, instead of emissions (NEE, Figure 1A).
On the other hand, R eco decreased under the canopy compared to the open grassland in spring but increased under the canopy in autumn (season x canopy effect, t = 4.73, p < 0.001, Table 2).R soil increased under the canopy, especially in autumn (season x canopy effect, t = 2.40, p = 0.02, Table 2), R soil being enhanced by T a (T a effect, t = 3.06, p = 0.003, Table 2) and BGB (BGB effect, t = 2.21, p = 0.03, Table 2).
Moreover, PFT composition drove significantly NEE.The most parsimonious and explanatory model was a specific interaction model (Equation ( 4)).Legumes enhanced net CO2 uptake (NEE, legumes effect, t = −3.97,p < 0.001, Table 2 and Figure 2), which slightly decreased with the addition of forbs (forbs × legumes effect, t = 3.84, p < 0.001, Table 2) and grasses (grasses × legumes effect, t = 4.07, p < 0.001, Table 2) in the mixture.2).F indicates forbs, G grasses, and L legumes.Predicted NEE modeled with a maximum proportion of legumes of 50%, according to the legumes proportion observed in the field (Figure S3).Negative NEE values indicate net CO2 uptake, and positive values indicate CO2 emissions.Environmental conditions set as the mean of the given treatment level.
N2O fluxes were uptake dominated in spring, except in the DN-pinea plot, where there were emissions both under the canopy and in the open grassland (Figure 3).Under the canopy N2O fluxes tend to increase (canopy effect t = 2.68, p = 0.01, Table 3), especially in autumn (Figure 3).N2O emissions also increased with litter (litter effect, t = 2.84, p = 0.006, Table 3).Moreover, PFT composition of the herbaceous layer drove significantly  2).F indicates forbs, G grasses, and L legumes.Predicted NEE modeled with a maximum proportion of legumes of 50%, according to the legumes proportion observed in the field (Figure S3).Negative NEE values indicate net CO 2 uptake, and positive values indicate CO 2 emissions.Environmental conditions set as the mean of the given treatment level.

CH 4 and N 2 O Exchange
CH 4 emissions were lower in autumn than in spring (season effect, t = −9.68,p < 0.001, Table 3) and were enhanced by T s (T s effect, t = 2.72, p = 0.008, Table 3) and SWC (marginally significant SWC effect, t = 1.83, p = 0.07, Table 3).N 2 O fluxes were uptake dominated in spring, except in the DN-pinea plot, where there were emissions both under the canopy and in the open grassland (Figure 3).Under the canopy N 2 O fluxes tend to increase (canopy effect t = 2.68, p = 0.01, Table 3), especially in autumn (Figure 3).N 2 O emissions also increased with litter (litter effect, t = 2.84, p = 0.006, Table 3).Moreover, PFT composition of the herbaceous layer drove significantly N 2 O fluxes, being the most parsimonious and explanatory model an identity model (Equation ( 2)), with the presence of legumes significantly enhancing N 2 O emissions (legumes effect, t = 2.49, p = 0.02, Table 3 and Figure 3).2)), with the presence of legumes significantly enhancing N2O emissions (legumes effect, t = 2.49, p = 0.02, Table 3 and Figure 3).

CO 2 Exchange Drivers
Our results highlight the relevance of tree canopies as drivers of CO 2 exchange, with NEE being uptake dominated in the open grassland in all the study plots, while, under the canopy, there were emissions in all DN plots (Figure 1A).These results support the well-known effect of light as driver of CO 2 uptake, as well as agree with previous CO 2 assessments in Mediterranean holm oak meadows [17,26,49] and other savannahlike ecosystems [25,28,50,51] that identified an enhancement of CO 2 emissions under the canopy compared to the open grassland.However, note that, in the SM-ilex, there was net CO 2 uptake under the canopy instead of emissions (Figure 1A), which interestingly suggests a differential canopy effect between locations (DN vs. SM).Thus, although the light reduction under the canopy in the SM-ilex plot was in the same range of magnitude than in all DN plots (Figure S1), the vegetation in the SM-ilex plot took up CO 2 at similar rates, both under the canopy and in the open grassland (Figure 1A).This might be related to the aboveground biomass (AGB) in the SM-ilex, which was similar under the canopy and in the open grassland, in comparison to a certain AGB reduction observed under the canopy in all DN plots (Figure S2 and Table S2).In line with this differential canopy effect between locations (SM vs. DN), Ibañez (2019) [52] reported that forbs from the same study plots presented more 13 C-depleted tissues under the canopy than in the open grassland in all DN plots, while this depletion was not noticeable in the SM-ilex plot.A fact that indicate that the vegetation was photosynthesizing at similar rates under the canopy and in the open grassland in the SM-ilex plot, while the photosynthetic rate differed in the DN plots.This suggests that, in the SM-ilex plot, the environment created under the canopy does not differ so much to the open grassland, in opposition to the strong differences found in DN, and emphasizes the relevance of the canopy effect in more environmentally constrained holm oak meadows, as is the case of DN compared to SM.
The inclusion of PFT composition and diversity improved NEE's explained variability and highlights the relevance of specific PFT identity and interaction effects (Table 2 and Figure 2), which were even more explicative than the AGB.Legumes were key NEE drivers, enhancing net CO 2 uptake (Table 2 and Figure 2), in agreement with previous studies that reported higher photosynthetic capacity of legumes compared to grasses and forbs [53][54][55][56], resulting in higher net CO 2 uptake [57].In addition, legumes have the ability to transfer symbiotic N to other species, increasing photosynthetic rates of the overall community [58,59], and, at the same time, the acquisition of symbiotic N by legumes can be stimulated by the presence of grasses [60].Accordingly, N fixed by legumes could be an important source of soil fertility [61] and provide an advantage for photosynthesis, especially in N limited systems, as might be the case of our plots, with very low soil N content (Table 1).
The tree-open grassland structure also drove ecosystem respiration dynamics (R eco , Figure 1A), but this effect was dependent on season and the magnitude of R eco components (R eco = R autotrophic + R heterotrophic ).In spring, the lower R eco rates under the canopy compared to the open grassland were the result of a lower temperature (T s effect, Table 2), which is known to be an important R eco driver [17], especially enhancing the heterotrophic component of R eco [62].In addition, enhanced R eco rates in spring and in the open grassland could also be indirectly related to increased gross CO 2 uptake rates, increasing markedly R autotrophic rates [63,64].In autumn conversely, the increased respiration rates (R eco and R soil ) recorded under the canopy compared to the open grassland (canopy × season, Table 2) were probably related to litter decomposition processes, increasing the heterotrophic component of R eco .The big amount of litter that was present under the canopy in spring was no longer present in autumn (Figure S2), suggesting that litter was already incorporated into the soil, and that it was probably a mineralization process.Respiratory fluxes, R eco and R soil , were also both influenced by the structure of the herbaceous layer, enhanced by the belowground biomass (BGB effect, Table 2), which has been directly linked to auto and heterotrophic respiration [65][66][67].
Moreover, our data also interestingly suggested that respiration rates (R eco and R soil ) under P. pinea canopies were higher than under Quercus species (Figure 1A, spring).This might be because the canopy of P. pinea is more open than that of Quercus species, with lower leaf area index [68], resulting in the lowest buffering of PAR (Figure S1A), and the lowest buffering of T s among all DN plots (Figure S1B) under the canopy, conditions that could be enhancing respiration rates.In addition, litter of P. pinea has been reported to be richer in carbon I content than that of Quercus species [69,70], which could be increasing soil organic C at the plot level (Table 1), and the carbon, nitrogen ratio in the soil under the canopy, as reported by Ibañez (2019) [52] from the same study plots.This, leading to enhanced R heterotrophic rates, since soil organic C has been positively related to CO 2 emissions [71]; in addition, microorganisms respire more C when decomposing C-rich and/or N-poor substrates [72].

CH 4 and N 2 O Exchange Drivers
Although there are some CH 4 and N 2 O assessments from savannah-like ecosystems (i.e., Herman et al., 2003 [31]; McLain et al., 2008 [51]; McLain and Martens, 2006 [27]), this type of surveys from Mediterranean holm oak meadows are very scarce [32], and our study provides one of the first GHG datasets of the Iberian Peninsula, unraveling the mechanisms behind these fluxes.CH 4 emissions were driven by seasonality rather than by the holm oak meadow structure and composition.CH 4 shifted from high emissions in spring, to low emissions or uptake in autumn (Figure 1B), probably as result of several factors combined.First, CH 4 emissions in spring (Figure 1B) were enhanced by temperature (T s effect, Table 3), which could be favoring methanogenic activity and methane diffusion [32].And second, SWC enhanced CH 4 emissions (although marginally significant, Table 2), which, in combination with drying-rewetting cycles, typical in spring, could also be enhancing CH 4 soil emissions [73].
On the other hand, we reported N 2 O emissions under the canopy in autumn, which could be the result of relatively high soil temperatures and soil moisture (T s and SWC effect, Table 3, and see the range of these variables in Figure S1), in combination with organic matter mineralization (see litter effect, Table 3, and the decrease in the litter fraction in autumn compared to spring, Figure S2).Those factors are known to enhance nitrification and denitrification processes, resulting in N 2 O emissions [74].In agreement, other studies have also reported higher N 2 O emissions under the canopy in combination with high temperatures and soil water moisture [27,51], as well as higher N 2 O emissions linked to a higher N availability and higher N-cycling rates under the canopy than in the open grassland [31].Interestingly, the DN-pinea plot performed again higher GHG fluxes, with enhanced N 2 O emissions in spring (Figure 1C) compared to plots dominated by Quercus species.The lower buffering effect of P. pinea canopy on PAR (Figure S1A) and T s (Figure S1B), and the higher soil organic C of the plot (Table 1), combined with higher soil N content under the canopy [52], could be enhancing N 2 O emissions.In agreement, Liu et al. (2014) [71] reported a positive relationship between soil organic C, soil N, and N 2 O emissions.However, the role of different tree species modifying soil properties, and this, in turn driving N 2 O dynamics, needs further investigation.
Moreover, N 2 O exchange was influenced by the structure and composition of the herbaceous layer, increasing N 2 O emissions with litter and the presence of legumes in the community (Table 3).The influence of litter on N 2 O exchange at field conditions is not well understood, but there are some experiments assessing the effect of cover crops on N 2 O exchange that may provide some understanding [75][76][77].For instance, an experimental study by Shaaban et al. (2016) [77] reported increasing N 2 O emissions with the addition of litter onto soil surface.The authors related the input of organic matter to stimulated microbial activity and denitrification, resulting in N 2 O production [77].In addition, N 2 O emissions have been negatively correlated to the C/N ratio of the substrate [77].Findings that agree with our results of legumes enhancing N 2 O emissions (Table 3 and Figure 3).

Management Implications
Holm oak meadows are productive ecosystems, and any management strategy must take that into account.In this regard, our results suggest that an increase in the tree cover may reduce the net CO 2 sink capacity of the understory, while also reduce forage production and quality.Under the canopy, there is a reduction on net CO 2 uptake rates compared to the open grassland, via (a) a direct canopy effect (especially in DN, Table 2), with the corresponding implications that this may have on forage production; and (b) an indirect effect through changes in the herbaceous layer composition.The latter related to a decrease in the presence of legumes under the canopy (Figure S3), which were enhancing net CO 2 uptake rates (NEE, Table 2 and Figure 2).With the implications that this may have on forage quality.Therefore, it is highly advisable to preserve open grassland spaces to maximize net CO 2 uptake and preserve forage provision.In addition, a change in the tree species composition, shifting from Quercus species stands to Pinus species plantations may increase CO 2 (Figure 1A) and N 2 O (Figure 1B) emissions, being highly advisable to preserve traditional Quercus, stands to minimize CO 2 and N 2 O emissions.

Conclusions
Our study provides insight into Iberian holm oak meadows functioning, integrating the arboreal and the herbaceous layer structure and composition as GHG drivers.Tree canopies, especially, drove CO 2 and N 2 O fluxes, increasing emissions under P. pinea canopy compared to that of Quercus species.The inclusion of the herbaceous layer composition and diversity terms improved explained variability, legumes enhancing net CO 2 uptake and N 2 O emissions.Thus, PFT identity and interaction effects were even more explanatory than some structural components (i.e., aboveground biomass).Changes in the tree cover and species will imply significant changes in the GHG exchange of Iberian holm oak meadows mediated by changes in the structure and composition of the herbaceous layer.This may provide some keys to improve ecosystem services provision and guarantee the preservation of this vast but vulnerable ecosystem.S1.Microclimatic conditions linear modeling results, as function of plot, season, and canopy.Photosynthetically active radiation (PAR), soil temperature (T s ) and soil water content (SWC).Plot with SM-ilex as reference level, season with spring as reference level, and canopy with open grassland (OG) as reference level.Parameter estimates (Par.), standard error (SE), t and p-value.Table S2.Herbaceous layer structure and composition linear modeling results, as function of plot, season, and canopy.Aboveground biomass (AGB), litter and plant functional types (PFT) biomass: forbs, grasses, and legumes.Season with spring as reference level, plot with SM-ilex as reference level, and canopy with open grassland (OG) as reference level.Parameter estimates (Par.), standard error (SE), t and p-value.BGB was quite variable among treatments (Figure S2), and neither season, plot or canopy explained its variability.Linear modeling on BGB is not shown.

Figure 2 .
Figure 2. Predicted net ecosystem exchange (NEE, μmol CO2 m −2 s −1 ) according to the specific interaction model (Table2).F indicates forbs, G grasses, and L legumes.Predicted NEE modeled with a maximum proportion of legumes of 50%, according to the legumes proportion observed in the field (FigureS3).Negative NEE values indicate net CO2 uptake, and positive values indicate CO2 emissions.Environmental conditions set as the mean of the given treatment level.

Figure 2 .
Figure 2. Predicted net ecosystem exchange (NEE, µmol CO 2 m −2 s −1 ) according to the specific interaction model (Table2).F indicates forbs, G grasses, and L legumes.Predicted NEE modeled with a maximum proportion of legumes of 50%, according to the legumes proportion observed in the field (FigureS3).Negative NEE values indicate net CO 2 uptake, and positive values indicate CO 2 emissions.Environmental conditions set as the mean of the given treatment level.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2073-4395/11/1/50/s1, Figure S1.Microclimatic sampling conditions per plot, season, and canopy: open grassland (OG) and under the canopy (UC).(A) Photosynthetically active radiation (PAR); (B) soil temperature (T s ) and; (C) soil water content (SWC).Boxplot's midline indicates the median; upper and lower limits of the box indicate the third and first quartile; whiskers extend up to 1.5 times the interquartile range from the top/bottom of the respective box, and points represent data beyond the whiskers.Figure S2.Mean ± SE of aboveground biomass (AGB), litter and belowground biomass (BGB) per plot, season, and canopy: open grassland (OG) and under the canopy (UC).
Figure S3.Mean ± SE of plant functional types (PFT: forbs, grasses, and legumes) per plot, season, and canopy: open grassland (OG) and under the canopy (UC).Table

Table 2 .
CO 2 exchange diversity-interaction model results.Net ecosystem exchange (NEE), ecosystem respiration (R eco ), and soil respiration (R soil ) as function of plot, season, canopy, photosynthetically active radiation (PAR), air temperature (T a ), aboveground biomass (AGB), belowground biomass (BGB), and plant functional types (forbs, grasses, and legumes) proportions.Season with spring as reference level, and canopy with open grassland (OG) as reference level.Parameter estimates (Par.), standard error (SE), t, and p-value.

Table 3 .
CH 4 and N 2 O exchange diversity-interaction model results.CH 4 and N 2 O exchange as function of plot, season, canopy, soil temperature (T s ), soil water content (SWC), litter, belowground biomass (BGB), and plant functional types (forbs, grasses, and legumes) proportions.Season with spring as reference level, and canopy with open grassland (OG) as reference level.Parameter estimates (Par.), standard error (SE), t, and p-value.