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

agronomy Article Tree—Open Grassland Structure and Composition Drive Greenhouse Gas Exchange in Holm Oak Meadows of the Iberian Peninsula Mercedes Ibañez 1,2,* , María JoséLeiva 3, Cristina Chocarro 4, Salvador Aljazairi 2,Àngela Ribas 5,6 and Maria-Teresa Sebastià1,2   Citation: Ibañez, M.; Leiva, M.J.; Chocarro, C.; Aljazairi, S.; Ribas, À.; Sebastià, M.-T. Tree—Open Grassland Structure and Composition Drive Greenhouse Gas Exchange in Holm Oak Meadows of the Iberian Peninsula. Agronomy 2021,11, 50. https://doi.org/10.3390/ agronomy11010050 Received: 4 November 2020 Accepted: 24 December 2020 Published: 28 December 2020 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2020 by the authors. LicenseeMDPI, Basel, Switzerland. This articleisanopenaccessarticledistributed under the terms and conditions of the CreativeCommonsAttribution(CCBY) license(https://creativecommons.org/ licenses/by/4.0/). 1Laboratory of Functional Ecology and Global Change, Forest Sciences Centre of Catalonia (CTFC), 25280 Solsona, Spain; [email protected] 2GAMES Group & Department HBJ, School of Agrifood and Forestry Science and Engineering (ETSEA), University of Lleida (UdL), 25198 Lleida, Spain; [email protected] 3Plant Biology and Ecology Department, University of Sevilla (US), 41080 Sevilla, Spain; [email protected] 4 Department of Crop and Forest Science and Agrotecnio-Center, School of Agrifood and Forestry Science and Engineering (ETSEA), University of Lleida (UdL), 25280 Lleida, Spain; [email protected] 5Centre for Ecological Research and Forestry Applications (CREAF), 08193 Bellaterra, Spain; [email protected] 6Departamento de Biología Animal, de Biología Vegetal y de Ecología (BABVE), Edifici C, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain *Correspondence: mer[email protected] Abstract: 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 influence 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 fluxes, 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. Keywords: canopy; CH 4 ; CO 2 ; dehesas; diversity-interaction model; N 2 O; plant functional types (PFT) 1. 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 – 5 ]. In addition, Mediterranean savannahs, with trees from other taxonomic groups, can be found in South Agronomy 2021,11, 50. https://doi.org/10.3390/agronomy11010050 https://www.mdpi.com/journal/agronomy Agronomy 2021,11, 50 2 of 16 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 – 12 ], water fluxes [ 13 ], evapotranspiration [ 14 ], vegetation structure [ 15 – 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 – 24 ], as well as on other savannah-like ecosystems [ 25 – 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 CO2release [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. Agronomy 2021,11, 50 3 of 16 2. Materials and Methods 2.1. 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 0 34” N, 6 ◦ 19 0 55” W, 30 m a.s.l.) and Sierra Morena mountains (SM, 37 ◦ 39 0 50” N, 5 ◦ 56 0 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). Table 1. Soil characteristics per plot and depth. Soil analysis performed according to: pH [ 39 ], organic C [ 40 ], total N (Elemental analyzer—measuring the total amounts of nitrogen, carbon and Sulfur (CNS methodology)—Trumac, LECO Corporation, St. Joseph, MI, USA) and texture [41]. Plot SM-ilex DN-mixed DN-suber DN-pinea Depth (cm) 0–40 40–80 0–30 30–60 0–30 30–60 0–30 30–60 pH 6.7 7.2 7.9 7.8 7.5 7.8 7.0 7.4 Organic C (%) 0.80 0.30 0.32 0.34 0.60 0.02 1.52 0.51 Total N (%) 0.85 0.60 0.15 0.04 0.06 0.05 0.20 0.11 Clay (%) 18 30 10 11 13 4 16 16 Silt (%) 29 28 25 15 18 9 22 21 Sand (%) 54 42 65 74 69 87 62 63 Study plots were selected according to their tree composition, with representative canopy types of Iberian holm oak meadows. One pure Q. ilex stand, in the SM location (SM- ilex), and one pure Q. suber stand in the DN location (DN-suber), both the most abundant stands; one Q. ilex and Q. suber mixed stand (DN-mixed), the next most abundant; and a pure Pinus pinea L. stand (DN-pinea), a common tree plantation replacing traditional canopies [ 7 ]. Plot tree densities (trees ha −1 ) were 34 ± 1 in SM-ilex, 26 ± 1 in DN-mixed, 26 ± 4 in DN-suer, and 48 ± 6 DN-pinea. Study treatments were, therefore, established according to plot (SM-ilex, DN-mixed, DN-suber, and DN-pinea), season (spring and autumn), and canopy (open grassland, OG, and under the canopy, UC). 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 Agronomy 2021,11, 50 4 of 16 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. 2.2. 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 Adj2 value (fluxes below a R Adj2 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]. Ancillary meteorological variables, including photosynthetically active radiation (PAR) outside the chamber (AccuPAR model LP-80 PAR/LAI ceptometer, Decagon Devices, Inc., Washington, D.C., USA); air temperature (T a ) inside and outside the chamber (multi-logger thermometer, TMD-56, Amprove, Washington, D.C., USA); soil temperature (T s ) 1–10 cm; and soil water content (SWC, gravimetric method, Section 2.4) were recorded to calculate and model GHG fluxes, as well as to characterize microclimatic sampling conditions (summarized in Figure S1). 2.3. 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. Agronomy 2021,11, 50 5 of 16 2.4. 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. 2.5. 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 GHG =βPlotPlot +βSeasonSeason ×βCanopyCanopy +βTT+βSWCSWC +βPARPAR +βAGBAGB +βLitterLitter +βBGBBGB +ε(1) 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 Pindicates the proportion of the given PFT and the sub-index Findicates forbs, Ggrasses, and Llegumes, respectively: Identity model GHG =Null model +βFPF+βGPG+βLPL+ε(2) (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 PiPj , where PFT is the number of PFT present in the community, and Pi 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 GHG =Null model +βFPF+βGPG+βLPL+evenness +ε(3) (c) Specific interaction model, which includes specific interactions between PFT in addition to the identity effects (Equation (4)): Specific interaction model GHG =Null model +βFPF+βGPG+βLPL+βFGPFPG+βFLPFPL+βGLPGPL+βFGLPFPGPL+ε(4) 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 2adj ), 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. Agronomy 2021,11, 50 6 of 16 3. 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). 3.1. CO2Exchange 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 CO2emissions under the canopy (Reco and Rsoil, Figure 1A). Agronomy2021,11,507of17    Figure1.Greenhousegasexchangeperplot,season,andcanopy:opengrassland(OG)andunderthecanopy(UC).(A) NetecosystemCO2exchange(NEE,green),ecosystemrespiration(Reco,garnet),andsoilrespiration(Rsoil,brown);(B)CH4 (purple);and(C)N2O(yellow)exchange.Mean±standarderror. 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), 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CO2uptakeunderthecanopy,insteadofemissions(NEE,Figure1A). Ontheotherhand,Reco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).Rsoilincreasedunderthecanopy,especiallyinautumn(seasonxcanopy effect,t=2.40,p=0.02,Table2),RsoilbeingenhancedbyTa(Taeffect,t=3.06,p=0.003, Table2)andBGB(BGBeffect,t=2.21,p=0.03,Table2).  Figure 1. Greenhouse gas exchange per plot, season, and canopy: open grassland (OG) and under the canopy (UC). ( A ) Net ecosystem CO 2 exchange (NEE, green), ecosystem respiration (R eco , garnet), and soil respiration (R soil , brown); ( B ) CH 4 (purple); and (C) N2O (yellow) exchange. Mean ±standard error. Agronomy 2021,11, 50 7 of 16 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. CO2Flux (µmol m−2s−1) NEE Reco Rsoil Specific Interaction Model Null Model Null Model Par. SE t p-Value Par. SE t p-Value Par. SE t p-Value Plot SM-ilex −6 3 − 1.81 0.07 8.9 0.6 14.60 <0.001 0 1 0.35 0.7 Plot DN-mixed −2 3 − 0.54 0.6 9.0 0.5 18.82 <0.001 1 1 0.55 0.6 Plot DN-suber −2 3 − 0.57 0.6 8.6 0.6 13.74 <0.001 2 1 1.33 0.2 Plot DN-pinea −1 3 − 0.38 0.7 9.2 0.6 14.69 <0.001 1 1 0.65 0.5 Season −6.2 0.6 − 9.88 <0.001 −1.3 0.5 − 2.56 0.01 Canopy 4 0.7 5.58 <0.001 −1.8 0.6 − 2.98 0.004 1.1 0.5 2.10 0.04 Season x canopy 4.3 0.9 4.73 <0.001 1.8 0.7 2.40 0.02 PAR (µmol photons m−2s−1) − 0.003 0.001 − 4.58 <0.001 Ta(◦C) 0.15 0.05 3.06 0.003 BGB (g DW m−2) 0.0015 0.0007 2.21 0.03 Forbs (fraction) 0 3 0.15 0.9 Grasses (fraction) −1 3 − 0.22 <0.001 Legumes (fraction) −67 17 − 3.97 <0.001 Forbs × legumes (fraction) 97 25 3.84 <0.001 Grasses × legumes (fraction) 96 24 4.07 <0.001 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 CO 2 uptake (NEE, legumes effect, t = − 3.97, p< 0.001, Table 2and 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. Agronomy 2021,11, 50 8 of 16 Agronomy2021,11,509of17   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.  Figure2.Predictednetecosystemexchange(NEE,μmolCO2m−2s−1)accordingtothespecificinteractionmodel(Table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. 3.2.CH4andN2OExchange CH4emissionswerelowerinautumnthaninspring(seasoneffect,t=−9.68,p<0.001, Table3)andwereenhancedbyTs(Tseffect,t=2.72,p=0.008,Table3)andSWC(marginallysignificantSWCeffect,t=1.83,p=0.07,Table3). 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 Figure 2. Predicted net ecosystem exchange (NEE, µ mol CO 2 m −2 s −1 ) according to the specific interaction model (Table 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. 3.2. CH4and N2O 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). Agronomy 2021,11, 50 9 of 16 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. CH4(nmol m−2s−1) N2O (nmol m−2s−1) Null Model Identity Model Par. SE t p-Value Par. SE t p-Value Plot SM-ilex −24 81 −0.30 0.8 Plot DN-mixed −83 81 −1.03 0.3 Plot DN-suber −22 85 −0.26 0.8 Plot DN-pinea −28 85 −0.33 0.7 Season spring −13 5 −2.73 0.008 Season autumn −150 15 −9.68 <0.001 −11 5 −2.30 0.03 Canopy 2.1 0.8 2.68 0.010 Ts(◦C) 13 5 2.72 0.008 0.4 0.2 2.44 0.02 SWC (fraction) 1.4 0.8 1.83 0.07 −0.05 0.03 −2.00 0.05 Litter (g DW m−2)−0.1 0.1 −1.89 0.06 0.007 0.003 2.84 0.006 BGB (g DW m−2)0.002 0.001 1.75 0.09 Forbs (fraction) 3 3 1.02 0.3 Grasses (fraction) 2 3 0.47 0.6 Legumes (fraction) 10 4 2.49 0.02 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 3and Figure 3). Agronomy2021,11,5010of17   N2Ofluxes,beingthemostparsimoniousandexplanatorymodelanidentitymodel(Equation(2)),withthepresenceoflegumessignificantlyenhancingN2Oemissions(legumes effect,t=2.49,p=0.02,Table3andFigure3).  Figure3.PredictedN2Oexchange(nmolm−2s−1)accordingtotheidentitymodel(Table3).Findicatesforbs,Ggrasses,andLlegumes.PredictedN2Omodeledwithamaximumproportionof legumesof50%,accordingtothelegumesproportionobservedinthefield(FigureS3).Negative N2Ovaluesindicatenetuptake,andpositivevaluesindicateN2Oemissions.Environmentalconditionssetasthemeanofthegiventreatmentlevel. Table3.CH4andN2Oexchangediversity‐interactionmodelresults.CH4andN2Oexchangeas functionofplot,season,canopy,soiltemperature(Ts),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.  CH4(nmolm−2s−1)N2O(nmolm−2s−1) NullModelIdentityModel Par.SEtp‐ValuePar.SEtp‐Value PlotSM‐ilex−2481−0.300.8    PlotDN‐mixed−8381−1.030.3    PlotDN‐suber−2285−0.260.8    PlotDN‐pinea−2885−0.330.7    Seasonspring   −135−2.730.008 Seasonautumn−15015−9.68<0.001−115−2.300.03 Canopy2.10.82.680.010 Ts(°C)1352.720.0080.40.22.440.02 SWC(fraction)1.40.81.830.07−0.050.03−2.000.05 Litter(gDWm−2)−0.10.1−1.890.060.0070.0032.840.006 BGB(gDWm−2)   0.0020.0011.750.09 Forbs(fraction)    331.020.3 Grasses(fraction)    230.470.6 Legumes(fraction)    1042.490.02 4.Discussion 4.1.CO2ExchangeDrivers OurresultshighlighttherelevanceoftreecanopiesasdriversofCO2exchange,with NEEbeinguptakedominatedintheopengrasslandinallthestudyplots,while,under Figure 3. Predicted N 2 O exchange (nmol m −2 s −1 ) according to the identity model (Table 3). F indicates forbs, G grasses, and L legumes. Predicted N 2 O modeled with a maximum proportion of legumes of 50%, according to the legumes proportion observed in the field (Figure S3). Negative N 2 O values indicate net uptake, and positive values indicate N 2 O emissions. Environmental conditions set as the mean of the given treatment level. Agronomy 2021,11, 50 16 of 16 74. Barnard, R.; Leadley, P.W.; Hungate, B.A. Global change, nitrification, and denitrification: A review. Global Biogeochem. Cycles 2005,19, 1–13. [CrossRef] 75. Baggs, E.M.; Rees, R.M.; Smith, K.A.; Vinten, A.J.A. Nitrous oxide emission from soils after incorporating crop crop residues. Soil Use Manag. 2000,16, 82–87. [CrossRef] 76. Lin, S.; Iqbal, J.; Hu, R.; Cai, J.; Shaaban, M.; Chen, X. Nitrous Oxide Emissions from Yellow Brown Soil as Affected by Incorporation of Crop Residues With Different Carbon-to-Nitrogen Ratios: A Case Study in Central China. Arch. Environ. Contam. Toxicol. 2013,65, 183–192. [CrossRef] 77. Shaaban, M.; Peng, Q.; Hu, R.; Lin, S.; Zhao, J. Soil Nitrous oxide and Carbon dioxide emissions following incorporation of above- and below-ground biomass of green bean. Int. J. Environ. Sci. Technol. 2016,13, 179–186. [CrossRef]