Can trees buffer the impact of climate change on pasture production and digestibility of Mediterranean dehesas?

Can trees buffer the impact of climate change on pasture production and digestibility of Mediterranean dehesas? Maria Dolores Hidalgo-Galvez a,b, ⁎, Karim Barkaoui c,d , Florence Volaire e , Luis Matías f , Jesús Cambrollé f , Pilar Fernández-Rebollo g , Maria Dolores Carbonero h , Ignacio Manuel Pérez-Ramos a a Institute of Natural Resources and Agrobiology of Sevilla (IRNAS-CSIC), 10 Reina Mercedes Avenue, 41012 Seville, Spain b Integrated Biology Doctoral Program, University of Seville, 6 Reina Mercedes Avenue, 41012 Seville, Spain c CIRAD, UMR ABSys, F-34398 Montpellier, France d ABSys, University of Montpellier, CIHEAM-IAMM, CIRAD, INRAE, Institut Agro, Montpellier, France e Centre d'Écologie Fontionnelle et Évolutive de Montpellier (CEFE-CNRS), 1919 Route de Mende, 34293 Montpellier cedex 5, France f Department of Plant Biology and Ecology, Faculty of Biology, University of Seville, 6 Reina Mercedes Avenue, 41012 Seville, Spain g Department of Forestry Engineering ETSIAM, School of Agricultural and Forestry Engineering ETSIAM, University of Córdoba, 14071 Córdoba, Spain h Department of Agricultural Production, Institute of Agricultural and Fishing Research and Education (IFAPA), km. 15, El Viso Road, 14270 Hinojosa del Duque, Córdoba, Spain HIGHLIGHTS •Habitat type had a significant influence on biomass productivity and digestibility. •Experimental warming increased by 33% biomass productivity, but reduced its quality. •Reduced rainfall decreased by 13% biomass productivity. •Changes in plant functional composition affected productivity and digestibility. •Tree canopy played a buffering role in mitigating the impact of climate change. GRAPHICAL ABSTRACT ABSTRACTARTICLE INFO Editor: Elena Paoletti Sustainability and functioning of silvopastoral ecosystems are being threatenedby the forecasted warmer and drierenvironments in the Mediterranean region. Scattered trees of these ecosystems could potentially mitigate the impact of climate change on herbaceous plant community but this issue has not yet tested experimentally. We carried out a field manipulative experiment of increased temperature (+2–3 °C) using Open Top Chambers and rainfall reduction (30%) through rain-exclusion shelters to evaluate how net primary productivity and digestibility respond to climate change over three consecutive years, and to test whether scattered trees could buffer the effects of higher aridity in Mediterranean dehesas. First, we observed that herbaceous communities located beneath tree canopy were less productive (351 g/m 2 ) than in open grassland (493 g/m 2 ) but had a higher digestibility (44% and 41%, respectively), likely promoted by tree shade and the higher soil fertility of this habitat. Second, both habitats responded similarly to climate change in terms of net primary productivity, with a 33% increase under warming and a 13% decrease under reduced rainfall. In contrast, biomass digestibility decreased under increased temperatures (−7.5%), since warming enhanced the fiber andlignin content and decreased thecrude protein content of aerial biomass. This warming-induced effect on biomass digestibility only occurred in open grasslands, suggesting a buffering role of trees in mitigating the impact of climate change. Third, warming did not only affect these ecosystem processes in a direct way but also indirectly via Keywords: Aridity Pasture quality Net primary productivity Silvopastoral ecosystem Warming Water stress Science of the Total Environment 835 (2022) 155535 ⁎Corresponding author at: Institute of Natural Resources and Agrobiology of Sevilla (IRNAS-CSIC), 10 Reina Mercedes Avenue, 41012 Seville, Spain. E-mail addresses: [email protected] (M.D. Hidalgo-Galvez), [email protected] (K. Barkaoui), fl[email protected] (F. Volaire), [email protected] (L. Matías), [email protected] (J. Cambrollé), [email protected] (P. Fernández-Rebollo), [email protected] (M.D. Carbonero), [email protected] (I.M. Pérez-Ramos). http://dx.doi.org/10.1016/j.scitotenv.2022.155535 Received 9 February 2022; Received in revised form 11 April 2022; Accepted 22 April 2022 Available online 27 April 2022 0048-9697/© 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv changes in plant functional composition. Our findings suggest that climate change will alter both the quantity and quality ofpasture production, with expectedwarmerconditions increasing net primary productivity but at the expense of reducing digestibility. This negative effect of warming on digestibility might be mitigated by scattered trees, highlighting the importance of implementing strategies and suitable management to control tree density in these ecosystems. 1. Introduction Sustainability and functioning of silvopastoral ecosystems are being threatened by climate change (Campos et al., 2013). Thus, the warmer and drier conditions forecasted for the Mediterranean region might impact net primary productivity and biomass quality, which represent major services in this type of ecosystems (Fay, 2009;Martin et al., 2014), potentially influencing other ecosystem services such as cattle production (Dumont et al., 2015;Lee et al., 2017;Rojas-Downing et al., 2017). Previous studies have detected a differential impact of climate change depending on the identity and intensity of the climatic stressor (Ma et al., 2017;Wilcox et al., 2017). On the one hand, moderate increases in temperature tend to advance plant maturity and increase photosynthetic performance, resulting in improved nutrient use efficiency and greater plant biomass (Martinez et al., 2014;Song et al., 2019;Viciedo et al., 2019). Nevertheless, warming reduces the period when plant biomass shows optimal nutritional quality, increasing the fiber and lignin contents and reducing its digestibility (Dumont et al., 2015;Lee et al., 2017;Habermann et al., 2019b). In contrast, when temperature overpass an optimal threshold, excessive warming may negatively affect seed germination, photosynthesis and respiration, thus decreasing net primary productivity and quality (Fahad et al., 2017; D'Orangeville et al., 2018). On the other hand, a moderate water deficit slows plant maturation, although if it does not cause severe leaf loss, biomass quality can be maintained or even improved because it delays stem development, reduces the lignification and increases the digestibility (Buxton, 1996;Reddy et al., 2003). Nevertheless, severe droughts usually reduce net primary productivity and quality because water deficit constrains the growth and development of new organs, hastens leaf senescence, alters fiber and sugar content, and relocates different nutrients and carbohydrates from leaves to roots, decreasing thus biomass digestibility (Nandintsetseg and Shinoda, 2013;Dumont et al., 2015;Li et al., 2018).Despite the unquestionable value of these studies, understanding the isolated role of these two abiotic stressors only represents a part of the overall picture since the effect of one of them may be exacerbated or mitigated by the other. For instance, warming can trigger higher soil mineralization and higher nitrogen absorption by plants, while drought can reduce plant nutrient uptake (Matías et al., 2011;Song et al., 2019). Although the two abiotic stressors associated with climate change may impact net primary productivity and biomass quality differently, the combined effect of both factors remains poorly known, particularly in Mediterranean areas (Volaire et al., 2014). Therefore, an accurate prediction of climate change projections on ecosystem functioning requires the evaluation of the joint and interactive effects of temperature and precipitation. The impact of climate change on these ecosystem properties (pasture productivity and quality) is expected to be spatially heterogeneous due to the mosaic of abiotic and biotic conditions caused by the typical scattered trees that characterize many silvopastoral ecosystems. Thus, trees reduce abiotic stress by attenuating air and ground temperature, decreasing plant evapotranspiration and reducing water stress for herbaceous plants growing under their canopies (Moreno, 2008;Abraham et al., 2014; Gargaglione et al., 2014;Gomes da Silva et al., 2021). In addition, trees favor hydraulic lift from the wettest soil layers (Ludwig et al., 2003, 2004a, 2004b;Yu and D'Odorico, 2015) and improve water infiltration (Joffre and Rambal, 1988;Ellison et al., 2017).Besides, treelitter decomposition provides nutrients and organic matter, which improves biomass quality under tree canopy (Ludwig et al., 2008;Cubera et al., 2009;Barneset al., 2011). All these benefits of scattered trees on net primary productivity and biomass quality indicate that they could mitigate the impact of climate change on ecosystem functioning (Fay et al., 2011;Tramblay et al., 2020; IPCC, 2021). However, this issue has barely been experimentally tested, particularly in the Mediterranean region. Trees might also attenuate the impact of climate change on ecosystem functioning indirectly via changes in the plant functional composition and diversity of pastures. The herbaceous species established under tree canopies usually develop morphological and physiological traits that allow them to toleratetheexcessshadeandthestrongcompetition caused by trees, for example by having high values of height, leaf area and SLA that allow them to increase light uptake in this habitat where light is more limited (Scholes and Archer, 1997;Ludwig et al., 2004b;Dalke et al., 2018). This increases the heterogeneity of plant communities at the landscape level and contributes to spatial changes in community functional composition, diversity of pastures and, ultimately, in ecosystem functioning (Bruno et al., 2003;Hisano et al., 2018). However, the influence of tree canopy on different ecosystem properties seems to be highly context-dependent, with variable responses as a function of the characteristics of the dominant herbaceous species and the local climate (Maestre et al., 2009;Rivest et al., 2013). Therefore, it is essential to analyze the indirect effects of the forecasted rise in temperature and water stress on ecosystem functioning based on the interactions established between trees and the co-existing herbaceous species. In Southwestern Spain, silvopastoral ecosystems termed ‘dehesas’cover around 3.1 million ha (Moreno and Pulido, 2009). This ecosystem with scattered trees –mainly Quercus sp. –provide pasture for extensive cattle breeding (Moreno, 2008;Santamaría et al., 2014;López-Díaz et al., 2015; Iglesias et al., 2016). Due to their essential economic relevance and their high susceptibility to changes in climate, management and habitat structure (Moreno and Pulido, 2009), it is timely to investigate the potential role of trees to mitigate the impact of climatechangeonecosystemfunctioning.To addressthisissue,wecarriedoutafield experiment simulating rainfall reduction and increased temperature as predicted by climate change models during three consecutive years. The main objective was to evaluate how these silvopastoral ecosystems will respond to future climatic scenarios, and how scattered trees could modify net primary productivity and biomass digestibility and modulate these potential climate-induced changes. Specifically, we addressed the following questions: (i) Does climate change have an impact on net primary productivity and biomass digestibility? (ii) Could we differentiate the direct from the indirect effects of climate change on these two ecosystem properties?; and (iii) Would trees be able to buffer the impact that climate change has on net primary productivity and biomass digestibility? We hypothesize that net primary productivity and biomass digestibility will be altered under future scenarios of climate change. We expect to find stronger effects of warming than rainfall reductiononthevariablesanalyzeddue to the recognized tolerance of Mediterranean plant species to survive long periods of drought. We hypothesize that warmer conditions will cause an increase in net primary productivity at the expense of its digestibility. Moreover, we hypothesize that scattered trees will buffer the impact of abiotic stress (especially thermal stress) on net primary productivity and biomass digestibility. To our knowledge, this is the first study that provides an overall picture about the impact of climate change (including warming and drought) on pasture functioning in silvopastoral ecosystems, and evaluates the potential buffering role of scattered trees under experimental field conditions (but see Moreno, 2008, focused on watering effects). 2. Material and methods 2.1. Experimental design The study was carried out at the Pedroches Valley (Córdoba, Southwestern Spain; 38° 22′50.64″N, 4° 45′27.69″W). The Mediterranean M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 2 climate of this area is characterized by cold, wet winters and hot, dry summers, showing a large intraand interannual variability. The annual precipitation is 419 mm/yr and the mean annual temperature is 15.4 °C, being January the coldest month with an average temperature of 5.9 °C, and July the warmest month with an average temperature of 26.5 °C (IFAPA meteorological station, Hinojosa del Duque; records over the period 2009–2019; web reference at References section). The florais characterized by a dense herbaceous layer (≥80% of cover) dominated by a diverse community of annual native species such as Geranium dissectum L., Hordeum murinum L. subsp. leporinum (Link) Arcangeli, Calendula arvensis L. subsp. arvensis and Echium plantagineum L. The mean species richness is 9.6 ± 0.3 species/m 2 , and they coexist with scattered oak trees [Quercus ilex L. subsp. ballota (Desf.) Samp., with around 20% of canopy cover]. Tree density is 14.5 ± 1.3 trees/ha. At this location, 36 experimental plots of 4 × 6 m were randomly established, placing half of them under the tree canopy and the other half in open grasslands (more than 10 m apart from the edgeof tree canopies). The minimum distance between trees was 20 m. In each plot, four climatic treatments were set up to simulate future climate forecasts for the XXI century (2040–2070) in the Mediterranean area (SRES A-2 model by the IPCC, 2021). We delimited an area of 2.5 × 2.5 m to evaluate the impact of current climate conditions (‘control’, C, hereafter) on ecosystem functioning. We installed Open Top Chambers (OTC, Marion et al., 1997) to simulate the predicted annual temperature increase of 2–3°C(‘warming’treatment, W, hereafter). They were built of methacrylate sheets with a hexagonal design and sloping sides of 40 × 50 × 32 cm (projected area in the soil of 0.65 m 2 ) without UV-Filter to avoid modifying the light spectrum and to allow wavelength transmission between 280 and 750 nm (Faberplast, Madrid). A ‘drought’treatment (D, hereafter) was implemented through the use of rain-exclusion shelters (2.5 × 2.5 × 1.5 m height) following the design of Matías et al. (2012). These shelters intercept 30% of the precipitation thanks to six methacrylate gutters of 0.14 m wide, and inclined at 20° (Yahdjian and Sala, 2002). Finally, to assess the impact of the combination of both sources of abiotic stress (‘warming + drought’treatment, W + D, hereafter), OTCs were installed under the rain-exclusion shelters. This experimental design resulted in 144 experimental units (18 plots × 2 habitat types × 4 climatic treatments). To avoid interference with large herbivores, experimental plots were enclosed with cattle fences. The experimental design outline and main objectives of this paper are detailed in Fig. 1. 2.2. Environmental characterization of the experimental plots 2.2.1. Soil conditions Soil texture was analyzed once in all 36 plots at the beginning of the experiment (2017) using samples from the first 10 cm soil layer (topsoil). These samples were collected by using a 3 cm diameter auger. Holes made with the auger were not filled afterwards to avoid using soil from an area situated outside our experimental units, which could alter the microbial community of the sampling points. Soil fertility was assessed each Fig. 1. Experimental design and the main variablesanalyzed in thestudy. We monitored 36 plots equally divided into two habitat types(‘under tree’and ‘open grassland’)and four climatic treatments(control, warming, drought and warming + drought) installed ineach one. We analyzed the direct, indirect and buffering effects of (i) abiotic factors such astree shade (estimated bytree leaf areaindex), soiltexture and fertility (ii)thermal and waterstress indices and (iii)community functionalcomposition(e.g.proportion in grasses, forbs and legumes) and structure (e.g. functional traits such as plant height, SLA or LNC) on net primary productivity and biomass digestibility. M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 3 year of study (2017–2019) using a composite top-soil sample collected at four different locations in each of the 144 experimental units.In the laboratory, samples were air-dried, crushed, and sieved, and the <2mmfraction was analyzed for standard chemical properties (Sparks, 1996). Eight chemical propertieswere determined: acidity (measured in a 1:2.5 soil:water suspension), total organic carbon (through the oxidation of organic matter with H 2 SO 4 and K 2 Cr 2 O 7 ), available P (following Watanabe and Olsen, 1965), NH 4 + and NO 3 − content (extracted with KCl 2 N and determined by distillation-titration), and available Ca + ,Mg + and K + (extracted with NH 4 CH 3 CO 2 1 N and determined by atomic absorption spectroscopy). 2.2.2. Microclimate 2.2.2.1. Tree cover. Tree cover was estimated in each of the 144 experimental units by measuring leaf area index [(LAI = leaf area (m 2 ) / ground area (m 2 )] with hemisphericalcanopyphotographs. They were taken in the middle of each of the 144 experimental units in spring 2018, before sunrise or after sunset. We used a horizontally-leveled digital camera (Coolpix 4500; Nikon, Tokyo, Japan) with a fish-eye lens of 180° field of view (FCE8; Nikon) situated at 0.5 m above ground level. Images were analyzed using Hemiview Canopy Analysis software version 2.1 (1999; Delta-TDevices Ltd., Cambridge, UK). 2.2.2.2. Thermal stress index. Air temperature was registered at ground level hourly with a resolution of 0.1 °C using HOBOMX2201 data loggers (Onset Computer, Bourne, MA, USA) in one sixth of the total experimental units (n = 24). A sensor was installed in each of the 24 randomly-chosen experimental units. According to some authors, plants generally suffer stress above 30 °C. From this temperature, photosynthesis experiences a rapid decline (Saini and Aspinall, 1982), which causes a reduction in biomass production (Benavides et al., 2009;Fahad et al., 2017). In addition, seed sterility appears (Saini and Aspinall, 1982) and root development is prevented (Ferris et al., 1998;Huang et al., 2012) when temperatures overpass 30 °C. Taking this information into account, cumulated thermal stress was calculated for mean temperature during January–July (covering the growth period of herbaceous plants in the study area) as in Eq. (1): Thermal stress ¼Accumulated temperature Tmean 30ðÞ Monitoring days January JulyðÞ (1) where T mean is the mean temperature, and monitoring days represent the number of days of the period sampled. 2.2.2.3. Water stress index. In half of the experimental units (n = 72), volumetric water content (% VWC) in the soil was periodically (i.e. weekly in spring and monthly in the rest of the year) assessed using a PR2 humidity probe (Delta-T Devices Ltd., Cambridge, UK) using permanent access tubes of 40 cm depth. It records moisture values at intervals of 10 cm. A permanent tube was installed in each of the 72 randomly-chosen experimental units. In order to determine the amount of water available for plants, the total transpirable soil water (% TTSW) was calculated for the average of the four depths as in Eq. (2): TTSW ¼VWCmax VWCmin (2) where VWC max and VWC min represent the maximum and minimum volumetric water content in the soil (%) for each experimental unit during January–July period, respectively. TTSW shows the ability of plant communities to uptake water under particular soil conditions because it is based on the observed SWC dynamics, assuming that the minimum SWC in summer reflects the limit of plant water uptake. Assuming that 30% of the TTSW is the threshold under which plants suffer from water stress, i.e. when their accessible soil water reserve is almost empty (e.g. Barkaoui et al., 2017), the water stress index was calculated during January–July using the following formula (Eq. (3)): Water stress index ¼VWCstress threshold VWCi TTSW x 0,3(3) where VWC i is the volumetric water content in each monitoring data, and VWC stress threshold shows the result of the sum of 30% TTSW and VWC min . 2.3. Plant measurements 2.3.1. Species composition Plant species abundance and composition of herbaceous species were determined yearly from 2017 to 2019 during the peak of vegetative growth. Two annual censuses were carried out between March and May to record the maximum number of species, including the entire phenological spectrum. Four 21 × 21 cm quadrats (divided into 9 squares of 7 × 7 cm) were randomly placed in each experimental unit. Species frequencies were calculated in each of the 144 experimental units from the number of squares where each plant species was present.Species of each experimental unit were classified into three groups: ‘grasses’,‘leguminous’and ‘forbs’. 2.3.2. Trait measurements The most dominant species (i.e. those making up at least the 70% of the total frequency of the community) were selected in each of the 144 experimental units. Five key functional traits were measured during the peak vegetation growth in spring of 2017 and 2018: plant height (cm), leaf area [LA (cm 2 )], leaf dry-matter content [LDMC = leaf dry weight (mg)/ hydrated leaf weight (g)], specific leaf area [SLA = leaf surface (cm 2 )/ leaf dry weight (g)], leaf nitrogen content [LNC (mg/g)] and leaf carbon content [LCC (mg/g)]. LNC and LCC were analyzed by the Analysis Service of the Institute of Natural Resources and Agrobiology of Seville. Plant height was measured in five individuals per species and experimental unit, whereas the five leaf traits were quantified in 10 individuals per species and experimental treatment combination (habitat type and climatic treatment) following standardized protocols (Cornelissen et al., 2003). Leaf size was quantified using an image analysis software (Image Pro-plus 4.5; Media Cybernetic Inc., Rockville, MD, USA). Traits were measured at the species level during two successive years, with contrasting climate conditions. Thus, 2017 was a typically dry year (275.4 mm), while 2018 was particularly wet (533.6 mm). During 2019, which was also dry (228.5 mm), no functional trait measurements were taken, but they were estimated from the average traits quantified in the two previous years. Previously, we evaluated the interannual variability of species trait values between 2017 and 2018, with the aim of identifying which traits were the most variable and possibly less accurate to estimate for 2019 (Fig. A.1). Although the analyzed traits varied from one year to another (especially LA, LNC and LCC),no clear bias was identified since points were distributed around the 1:1 line (Fig. A.1). 2.3.3. Functional composition of plant communities We employed the formula provided by Garnier et al. (2004) to calculate community weighted means (CWMs, hereafter) for each of the traits mentioned above (Eq. 4): CWM ¼∑n i¼1pitraiti(4) where p i is the relative contribution of the species i to the community, n is the number of most abundant species, and trait i is the trait value of the species i. 2.4. Net primary productivity and biomass quality 2.4.1. Net primary productivity At the end of the vegetative cycle (late June) of 2017, 2018 and 2019, aerial biomass produced in each of the 144 experimental units was M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 4 collected using 50 × 50 cm adjacent quadrats. Vegetation was cut at ground level, trying to collect the maximum possible amount of aerial biomass. Collected samples were cleaned of non-herbaceous material, and dry biomass was determined using a precision balance after drying it in a forced-air oven at 60 °C for 48 h. 2.4.2. Biomass digestibility Samples of dry biomass were finely ground (<2 mm sieve) in an IKA mill. Enzymatic digestibility of organic matter (hereafter, digestibility) was estimated using the near-infrared reflectance spectroscopy technique (Vis-NIRS). Samples were scanned with the LabSpec 5000 spectrophotometer (350–2500 nm; ASD Inc., Boulder, Colorado, USA) using IndicoPro6.0 software (ASD Inc., Boulder, CO, USA). Four replicates of each biomass sample were scanned (each being an average of 50 internal scans). The final sample spectrumwas obtainedby averaging the four replicates. Statistic of Vis-NIRS equation used to predict digestibility is presented as an appendix (Table A.1). 2.5. Statistical analyses We conducted a multiple factor analysis (MFA) to analyze the multivariate correlations among the abiotic and biotic factors quantified in this study and toidentify the major axes of variation among the 144 experimental units. First, we analyzed differences between habitat types and their interaction with climatic treatments on thermal and water stress indices and CWMs using ANOVAs, followed by Tukey post-hoc tests. Second, we tested the individual and combined effects on net primary productivity and biomass digestibility of the following fixed factors: habitat type, climatic treatment, year, proportion of grasses, CWMs of plant height, LDMC, SLA, and LNC. For this purpose, we used Linear Mixed-effect Models (LMM), including the plot as a random factor, which was coded to affect only intercepts (not slopes). We compared a ‘full model’based on the combination of all factors to a series of ‘reduced models’in which factors were deleted one by one according to their significance (elimination of the least significant) until reaching the ‘null model’. To accomplish this, we used the AICc (Akaike Information Criterion) with a greater penalty for an extra parameter, as recommended by Burnham and Anderson (2002). Models were considered different when ΔAICc <2. The best model withthe lowest AICc was evaluated using maximum likelihood estimation. The proportion of variance explained by the fixed factors alone and by both the fixed and random factors was estimated by R 2 m(‘marginal R 2 ’) and by R 2 c(‘conditional R 2 ’), respectively (Nakagawa and Schielzeth, 2013). We derived the effect size of each factor using the estimates of the best model. Previously, normality of these variables was tested using Kolmogorov-Smirnov test, and variables were transformed in the cases in which it was necessary. Third, we then evaluated the relationships between stress indices and net primary productivity and biomass digestibility using bivariate linear regressions.Statistical analyseswere performed with STATISTICA® Software and within the R environment v.4.0.0 using tidyverse (for data manipulation and visualization), cowplot (for multiple plots compilation), Hmisc (for graphical representation), ggpubr (for color and shapes), FactoMineR (for MFA computation), lme4 and lmerTest (for running linear mixed models), AICcmodavg (for computing AICc and comparing models), car and pbkrtest (for testing the estimated parameters, p-value)andMuMIn (for determining R 2 mandR 2 c) packages (R Core Group, 2020). 3. Results 3.1. Descriptive analysis of the abiotic environment and plant community functional structure 3.1.1. Differences between habitat types The 144 experimental units captured a large heterogeneity in both abiotic and biotic conditions, primarily promoted by differences associated to the two habitat types considered in our study, as summarized by the MFA (Fig. 2). The MFA explained 54.80% of the total variability, with 33.40% accounted by dimension 1 (Dim-1, hereafter) and 21.40% by dimension 2 Fig. 2. Plot of the first and second axes of the multiple factor analysis (MFA) relating plant community functional structure, soil characteristics [soil fertility (organic matter, OM;C;N;P;NO 3 − ;NH 4 + ;Ca + ;Mg + ;K + ), and soil texture (coarse-sand; fine-sand; clay; silt)] and tree shade (represented as leaf area index, LAI). Plant community functional structure was characterized in each of the 144 experimental units according to the following information: the proportion of different plant functional types (grasses, forbs and legumes), and the possession ofsome key functional traits [plant height (Height); leaf dry-matter content (LDMC); leaf area (LA); specific leafarea (SLA); leaf nitrogencontent (LNC); leaf carbon content (LCC)] (panel a). Projection of the 36 plots on the plane defined by the two main MFA dimensions, with emphasis on differences due to habitat types (open grassland versus under tree). Three years of sampling (2017, 2018 and 2019) were represented (panel b). M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 5 (Dim-2, hereafter). Dim-1 was determined mainly by the level of tree shade (LAI) and, to a lesser extent by soil organic matter (OM), the proportion of grasses and some leaf morphological traits such as LA and SLA. Dim-2 was mainly determined by the two leaf chemical traits (LNC and LCC) and soil texture (Fig. 2a). Factor loadings are detailed in Fig. A.2. Habitat types differed strongly along Dim-1, with plots located under trees presenting higher OM and K + content inthe soil, and their plant communities having a higher proportion of grasses, higher SLA and LNC than those located in open grasslands (Fig. 2b). Microclimatic conditions also differed in the two habitat types. On the one hand, the annual mean air temperature was slightly lower under trees than in the open grasslands (15.5 ± 0.1 vs. 16.2 ± 0.1 °C). Thus, habitat located in open grasslands presented higher thermal stress than those situated under tree canopies, as evidenced by the negative relationship between LAI and thermal stress (Fig. A.3a; p<0.001, R 2 = 0.66). No significant differences in water stress were detectedwhen comparing both habitat types (Fig. A.3b; p= 0.54, R 2 =−0.003). On the other hand, soil moisture was slightly lower under trees (p= 0.08) when using an average of the 4 seasons of theyear. These differences between habitat typesin soil moisture were more pronounced in winter (p<0.001) and spring (p<0.001) than in the other two seasons (p=0.84insummerandp= 0.78 in autumn) (Table A.2). 3.1.2. Differences among climatic treatments As expected, the climatic treatments exerted a significant influence on microclimatic conditions. First, the two treatments using Open Top Chambers (W and W + D)significantly increasedair temperature for C or D treatments by about 2 °C (Fig. A.4). Consequently, these plots subjected to temperature increase exhibited the highest thermal stress values, but differences were only significant in open grasslands (Table A.2). Second, spring moisture in the top-soil was significantly higher in the C and W treatments than in the D and W + D treatments, but these differences were only significant in open grasslands (Table A.2 and Fig. A5). Water stress, however, did not differ significantly among climatic treatments (Table A.2). Our experimental climatic treatments also induced significant changes in species composition and plant community functional structure. First, we detected a higher abundance of grasses in those units subjected to increased temperature. Second, plants subjected to warmer conditions were taller compared to those growing under C and D conditions, although these differences were only significant in open grasslands (Table A.3). 3.1.3. Differences among years The three sampling years (2017, 2018 and 2019) did not differ significantly in thermal stress but did in water stress accumulated over the growing season (Table A.4). 2017 and 2019 were particularly dry (especially in spring), whereas 2018 was extremely rainy compared with the yearly mean precipitation registered for the last 20 years at the study area. Theseamongyear differences in precipitation resulted in significant differences in SVWC and water stress. Specifically, 2019 presented higher values of accumulated water stress than 2018 and 2017 (Table A.4). 3.2. Factors affecting net primary productivity and biomass digestibility Net primary productivity and biomass digestibility were affected by most of the fixed factors considered in this study, as evidenced by the inclusion of habitat type, climatic treatment, year and some variables related with plant community functional structure in the best linear mixed models (i.e. those with the lowest AICcvalues). The models presented a high goodness of fit in the two response-variables, with R 2 c values ranging from 0.33 to 0.47. In general, the accuracy of the models only improved slightly when including the random effects of the experimental plots, as denoted by the small differences between R 2 candR 2 m(Table 1). 3.2.1. Differences between habitat types Significant differences in net primary productivity and biomass digestibility were found between habitat types (Tables 2 and 3), with lower values of net primary productivity and higher values of biomass digestibility under tree canopies than in open grasslands (p<0.001; 351 g/m 2 vs. 493 g/m 2 for net primary productivityand 44% vs.41% for biomass digestibility, respectively; Fig. 3). 3.2.2. Differences among climatic treatments Climatic treatments influenced both net primary productivity and biomass digestibility, but with relevant differences depending on the response variable and the habitat type (Tables 2 and 3). In general, plant communities experimentally subjected to warmer conditions had higher net primary productivity (an increase of 33%; p<0.001) (Fig. 3a) but lower biomass digestibility (a decrease of 7.5%; p= 0.01) (Fig. 3b). Both habitat types exhibited a similar trend, with enhanced net primary productivity in the W treatment (Fig. 3a). However, reduced biomass digestibility under warming conditions was significant in plots located in open grasslands but not in those under tree canopy (Fig.3b). Conversely, reducedrainfall caused a significant decrease of net primary productivity by 13% (p<0.001) with no significant impact on biomass digestibility (p= 0.45). This above-reported variation in net primary productivity and biomass digestibility can be related to differences between habitat types and climatic treatments in abiotic stress. On the one hand, higher thermal stress Table 1 Results from the linear mixed-effects models used to analyze the effects of climate change on the response variables considered in this study (net primary productivity and biomass digestibility). Habitat type (under tree and open grassland) and sampling year (2017, 2018 and 2019) were also included as fixed factors into the models to control the effects mediated by (habitat-associated) microclimatic and (inter-year) climatic differences, respectively. The best models were selected using the corrected Akaike Information Criterion (AICc). Differences between the best and the null model (ΔAICc) are indicated. Variability explained by fixed factors (R 2 marginal, R 2 m), by both fixed and random factors together (R 2 conditional, R 2 c). Response variable Best model ΔAICc R 2 mR 2 c Productivity Habitat type × Climatic treatment + Plant height + LNC + Grasses + Year + (1|Plot) 101.42 0.28 0.47 Digestibility Habitat type × Climatic treatment + Grasses + Year + (1|Plot) 86.43 0.27 0.33 Table 2 Results from the analysis of variance evaluating the influence of fixed factors (habitat type, climatic treatment, plant height, LNC, grasses, year and the interaction habitat type × climatic treatment) on net primary productivity. Factor SS d.f. F p Habitat type 615049 1 18.55 <0.001 Climatic treatment 2249541 3 22.62 <0.001 Plant height 157965 1 4.76 0.03 LNC 139633 1 4.21 0.04 Grasses 229839 1 6.93 0.01 Year 1315215 2 19.84 <0.001 Habitat type × Climatic treatment 94572 3 0.95 0.42 Table 3 Results from the analysis of variance evaluating the influence of fixed factors (habitat type, climatic treatment, grasses, year and the interaction habitat type × climatic treatment) on biomass digestibility. Factor SS d.f. F p Habitat type 1075.73 1 42.40 <0.001 Climatic treatment 294.49 3 3.87 0.01 Grasses 826.50 1 32.58 <0.001 Year 375.48 2 7.40 <0.001 Habitat type × Climatic treatment 85.13 3 1.12 0.34 M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 6 caused an increase innet primary productivity (p<0.001, R 2 =0.17)anda decrease in biomass digestibility (p<0.001, R 2 =0.19)(Fig. 4aandb,respectively). On the other hand, higher water stress caused a decrease in both net primary productivity (p=0.08,R 2 = 0.03) and biomass digestibility (p= 0.002, R 2 = 0.11) (Fig. 4c and d, respectively). 3.2.3. Differences among years Net primary productivity and biomass digestibility varied among the three sampling years, with 2019 presenting a significant reduction in these variables compared to 2017 and 2018 (p<0.001). Specifically, net primary productivity decreased from 459.5 g/m 2 in the most productive year (2017) to 356.5 g/m 2 in 2019 (22.4% reduction). Biomass digestibility decreased by 9% from 2018 (year with the highest values) to 2019. 3.2.4. Influence of plant community functional structure Higher plant height and higher LNC were associated with greater net primary productivity (Fig. 5a; p=0.04,R 2 =0.04;and5b;p=0.71, R 2 =−0.01; respectively). A higher proportion of grasses was slightly and positively related with net primary productivity (Fig. 5c; p= 0.29, R 2 = 0.001). By contrast, biomass digestibility decreased strongly as the proportion Fig. 3. Effects of climatic treatments (control, warming, drought and warming + drought) on net primary productivity (g/m 2 ) (panel a) and biomass digestibility (%) (panel b) for the two habitat types considered in this study (under tree and open grassland). The mean of the three sampling years of (2017, 2018 and 2019) is represented. Letters (lowercase for habitat located under tree and uppercase for open grassland) are the result of factorial ANOVA, and subsequent post-hoc Tukey analyseswereperformedto compare the four climatic treatments within each habitat. Fig. 4. Response of net primary productivity (g/m 2 ) (panels a and c) and biomass digestibility (%) (panels b and d) to thermal stress (°C) and accumulated water stress (%), respectively. Open symbolsrepresentthe climatic treatments inopen grassland habitat,whereas closed symbolsrepresent the climatic treatments under tree habitat. The line is the result of the linear regression for the whole data set, whose equation is also indicated. Three years of sampling (2017, 2018 and 2019) are represented. M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 7 of grasses increased in the communities (Fig. 5d; p<0.001, R 2 =0.14). Detailed significance values are shown in Tables 2 and 3. 4. Discussion We present a pioneering work on the buffering effect of trees on climate change and their impact on net primary productivity and biomass digestibility in semi-arid silvopastoral ecosystems of Southern Spain. We show how increased temperatures and reduced rainfall affect these ecosystem properties through direct (impact of abiotic stress) and indirect (changes in plant functional structure) effects. We discuss the role of trees in modifying aridity impacts on ecosystem properties of dehesa ecosystems. 4.1. Climate change had an impact on net primary productivity and biomass digestibility Our experimental study of climate change simulation indicated that future climatic conditions could alter net primary productivity (both in quantity and quality), with contrasting effects of the two studied abiotic stressors. 4.1.1. Direct effects of climate change (abiotic factors) On the one hand, increased temperaturefavored net primary productivity in the studied dehesas. These findings are in agreement with the results found by previous studies (Wu et al., 2011;de Assis et al., 2016;Yu et al., 2019), who showed that higher temperatures increased ecosystem respiration, causing an increase in plant productivity, as a consequence of an enhanced metabolism and growth (Lawlor, 1993). On the other hand, reduced rainfall caused a significant decrease in net primary productivity, probably because a lower soil water availability impeded seedling establishment and reduced plant growth (Espigares and Peco, 1993; Dimitrakopoulos and Bemmerzouk, 2003). Moreover, a decrease in soil water content tends to reduce nutrient availability to plants because it limits soil microbial processes and ecosystem CO 2 fluxes, which translates into a decrease in plant productivity (De Dato et al., 2008;Sardans et al., 2008;Wu et al., 2011). This result is in agreement with previous works that reported a negative impact of water stress on vegetation diversity and net primary productivity (Fu et al., 2013, 2018;Fu and Shen, 2016; Wang et al., 2013;Sternberg et al., 2015). As a consequence of these contrasting effectsof both climatic stressors, the amount of net primary productivity reached intermediate values in the combined W + D treatment,likely as a result of offsetting the positive effect of warming with the negative effect of drought. This lack of combined effects of the two climatic stressors contrasts with the results of other recent studies in regions with low rainfall (e.g. Suzuki et al., 2014;Zandalinas et al., 2018), where a moderate warming exacerbated the negative impacts of drought on plant growth. Conversely to productivity, biomass digestibility was significantly lower under experimental warming, revealing that an increase in temperature negatively affects biomass quality (Čop et al., 2009;Bloor et al., 2010; Habermann et al., 2021). According to Wilson (1982), temperature is the most important environmental factor influencing biomass quality, reporting that for every 8 °C that temperature increases, digestibility of grasses decreases ~0.5%. Plant species, especially annuals, tend to advance their reproductive period in response to rising temperature, reaching earlier plant maturity and senescence (Whittington et al., 2015;Valencia et al., 2016;Moore and Lauenroth, 2017). In general, plant maturation increases fiber and lignin content and reduces crude protein content and digestibility, decreasing thus the nutritional value of plant biomass (Neel et al., 2008;Nordheim-Viken et al., 2009;Habermann et al., 2019b). In contrast, reduced rainfall did not significantly affect biomass digestibility. This lack of drought effects on biomass digestibility could be a Fig. 5. Influence of plant functional traits [expressed in terms of CWMs of plant height (cm) and LNC (mg/g) (panels a and b, respectively)] and plant species composition [expressed in terms of percentageof grasses (panels c and d)] on net primary productivity (g/m 2 ) and biomass digestibility (%).Open symbols represent those locatedin open grassland, whereas closed symbols represent the experimental units located under tree habitat. The line is the result of the linear regression for the whole data set, whose equation is also indicated. Three years of sampling (2017, 2018 and 2019) are represented. M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 8 consequence of the lower productivity derived from reduced rainfall, which implies a reduction of the total transpiration surface and, therefore, water availability per plant biomass may be similar independently from water availability. Our results suggest that the predicted temperature increase will decrease biomass digestibilityin thedehesas of Southern Spain, withpotential negative consequences for the cattle that feed on (lower nutritional intake) and for the farmers who obtaineconomic benefits from them (lower quality of milk, meat, etc.) (Pearce et al., 2010;Alothman et al., 2019). Additionally to the significant influence of climatic treatments on net primary productivity and biomass digestibility, we also detected a large variation between years, likely promoted by the inter-annual differences in water stress. Specifically, the lowest values of net primary productivity appeared in the driest year (2019), whereas the highest values ofpasture digestibility were collected in the rainiest year (2018). This result is in agreement with those found in previous studies documenting greater inter-annual variation of plant production than that promoted by experimental climatic treatments (Yu et al., 2019). 4.1.2. Indirect effects of climate change Climate change might also indirectly affect the ecosystem processes quantified in our study via significant changes in plant species composition and functional trait distribution within plant communities (Cowles et al., 2016; Napier et al., 2016;Yang et al., 2017;Li et al., 2018). On the one hand, our results revealed that grasses have a higher abundance under an experimental increased temperature (Shi et al., 2016). This enrichment in grasses under warmer conditions could partially explain the higher productivity and the lower biomass digestibility detected in the warming treatment. Similarly, Vázquez de Aldana et al. (2008) showed that an increase in precipitation over spring favored the presence of legumes, which provides a pasture with higher protein content and digestibility, while in drier springs the proportion of grasses increased and therefore the lignin content was higher. On the other hand, plant communities exhibited higher plant heights in the warming treatment and the lowest ones under to reduced rainfall. This could explain the enhanced productivity with warming detected in our study given the high correlation detected between this response variable and plant height at the community level. Previous studies have also found that plant height and productivity increase with higher temperature (Debouk et al., 2015;Zhang et al., 2015;Habermann et al., 2019a, 2019c) as a result of an acceleration of biochemical reactions (Copeland, 2000). 4.2. Tree canopy had an influence on net primary productivity and biomass digestibility We observed different trends of net primary productivity and biomass digestibility according to the habitat type. Thus, plant communities located beneathtreecanopywerelessproductivethanthoselocatedinopengrasslands but exhibited biomass of higher quality. The reduced productivity under tree canopy could be explained by the microclimatic alterations induced by tree shade. As stated by previous studies (Benavides et al., 2009;Hussain et al., 2009;Pang et al., 2019),thelevelofshadeanditsdurationaffectthephysiological processes of plants, decreasing plant carbohydrate fabrication and net dry matter production. In contrast, other studies showed a positive effect of tree canopy on plant growth attributed to the reduced incoming solar irradiation and the temperature attenuation provided by the tree canopy (Frost and McDougald, 1989;Belsky, 1994;Moreno, 2008). Some studies revealed that tree enhance water availability by hydraulic lift from the wettest soil layers and water infiltration (Joffre and Rambal, 1988;Ludwig et al., 2003, 2004a, 2004b), whereas other studies found the opposite pattern due to the large competition for resources between trees and herbaceous plants (Fay et al., 2003;Köchy et al., 2008). All thesefindings suggest thatthe impact of tree canopy on herbaceous productivity is highly context-dependent, being strongly influenced by the environmental characteristics of the ecosystem. Our results indicated that between-habitat type differences in productivity were due to an ameliorating effect of the tree on thermal stress since similar values of water stress were detected in both habitat types. Therefore, these results suggest that trees exercise a buffer effect on temperature but not on available water, potentially reducing the impact of climate change on the ecosystem processes quantified in this study (as explained in detail below). Similarly, Mordelet and Menaut (1995) comparatively analyzed ten studies where the enhancement or suppression of biomass production by tree canopies was unrelated to annual rainfall. Interestingly, we found the opposite trend for pasture digestibility, with higher values under tree habitat compared with those located in open grasslands. This enhanced pasture quality in more shaded microsites could be explained by the higher levels of soil fertility found under tree canopies, which provide soils large amounts of organic matter through leaf fall (Treydte et al., 2007, 2008;Ludwig et al., 2008;Cubera et al., 2009;Barnes et al., 2011). In addition, these microsites are usually enriched in nutrients due to the foraging patterns of large herbivores, which spend more time in this habitat to take shelter from the sun and consume acorns (Belsky et al., 1989;Belsky, 1994;Bardgett et al., 1998;Dijkstra et al., 2006). 4.3. Trees buffer the impact of climate change on biomass digestibility but not on net primary productivity Our findings revealed that tree canopy could buffer the impact of climate change regarding only biomass digestibility. Both habitat types exhibited similar patterns in response to climate manipulation, with enhanced productivity under warming and decreased productivity under reduced rainfall. However, the decrease of biomass digestibility under warming conditions was only significant in open grasslands but not under tree canopy. This buffering effect of the tree canopy was likely promoted, at least partially, by an ameliorating effect of tree shade on thermal stress, which has been identified in our study as the main abiotic stressor that reduced biomass digestibility. Under these more favourable conditions, a slight increase of temperature was not likely enough to cause significant changes in the phenological cycle of the plant communities inhabiting this habitat, thus reducing the impact of warming on biomass digestibility. It has been proven that tree shade has a buffering effect on thermal stress and excess irradiance, with potential benefits on overall crop productivity and quality (Lin, 2007; Lott et al., 2009;Bayala et al., 2014;Sida et al., 2018). Tree buffering effects on biomass digestibility could be also a consequence of the alterations of plant species composition and functional structure induced by tree shade. We detected a marked functional trait selection according to the microclimatic conditions generated under the tree canopies. On the one hand, we observed that plant communities located under trees developed higher SLA values, likely as an effective strategy to increase light uptake in this habitat where light is limited (Steinger et al., 2003;Curt et al., 2005). On the other hand, these plant communities growing beneath tree canopy exhibited a higher content in nitrogen, likely as a result of the mineralization of excess organic matter and nutrients through leaf fall (Gallardo, 2003;Cubera et al., 2009;Barnes et al., 2011) and the higher accumulation of herbivore droppings (Belsky et al., 1989;Kendall et al., 2006;Tucker et al., 2008). These functional attributes induced by tree canopy, which often confer higher nutritional value on plant biomass, probably attenuated the impact of simulated climate change on biomass digestibility. Further studies including a large array of functional traits are necessary to assess whether plant communities growing under tree canopies contain plant species inherently more adapted to climate change than those dominant in open grasslands. 5. Conclusions This study provides an analysis of the major potentially influencing factors that can affect net primary productivity and biomass digestibility under different climatic scenarios. Although we detected an expected negativeeffect of reduced rainfall on net primaryproductivity, the strongest influence was exerted by warming. Higher temperatures caused an increase in net primary productivity and a decrease in biomass digestibility, M.D. Hidalgo-Galvez et al. Science of the Total Environment 835 (2022) 155535 9