1 Running Head 1 Plant and microbes affect multifunctionality 2 3 Cascading effects from plants to soil microorganisms explain how plant species 4 richness and simulated climate change affect soil multifunctionality 5 6 7 Enrique Valencia1,2*, Nicolas Gross1,3,4, José L. Quero5, Carlos P. Carmona6, Victoria 8 Ochoa1, Beatriz Gozalo1, Manuel Delgado-Baquerizo1,7, Kenneth Dumack8,9, Kelly 9 Hamonts7, Brajesh K. Singh7,10, Michael Bonkowski8,9 & Fernando T. Maestre1 10 11 1Departamento de Biología y Geología, Física y Química Inorgánica, Escuela Superior de 12 Ciencias Experimentales y Tecnología, Universidad Rey Juan Carlos, C/ Tulipán s/n, 13 28933 Móstoles, Spain. 14 2Department of Botany, University of South Bohemia, Na Zlate stoce 1, 370 05 Ceske 15 Budejovice, Czech Republic. 16 3INRA, USC1339 Chizé (CEBC), F-79360, Villiers en Bois, France. 17 4Centre d’étude biologique de Chizé, CNRS - Université La Rochelle (UMR 7372), F-18 79360, Villiers en Bois, France 19 5Departamento de Ingeniería Forestal. Escuela Técnica Superior de Ingeniería Agronómica 20 y de Montes Universidad de Córdoba Campus de Rabanales Crta. N-IV km. 396 C.P. 14071 21 Córdoba, Spain. 22 6Institute of Ecology and Earth Sciences, University of Tartu, Lai 40, 51005, Tartu, Estonia 23 7Hawkesbury Institute for the Environment, University of Western Sydney, Penrith 2751, 24 New South Wales, Australia. 25 8Universität zu Köln, Zoologisches Institut, Terrestrische Ökologie. 26 2 9Cluster of Excellence on Plant Sciences (CEPLAS), University of Cologne. 27 10Global Centre for Land-Based Innovation, Western Sydney University, Penrith, NSW 28 2751, Australia 29 30 *Corresponding author: 31 Enrique Valencia Gómez: 32 Tel: +346650133 33 Email: valencia.gomez.e@gmail.com 34 mailto:valencia.gomez.e@gmail.com 3 Abstract 35 Despite their importance, how plant communities and soil microorganisms interact to 36 determine the capacity of ecosystems to provide multiple functions simultaneously 37 (multifunctionality) under climate change is poorly known. We conducted a common 38 garden experiment using grassland species to evaluate how plant functional structure and 39 soil microbial (bacteria and protists) diversity and abundance regulate soil 40 multifunctionality responses to joint changes in plant species richness (1, 3 and 6 species) 41 and simulated climate change (3°C warming and 35% rainfall reduction). The effects of 42 species richness and climate on soil multifunctionality were indirectly driven via changes 43 in plant functional structure and their relationships with the abundance and diversity of soil 44 bacteria and protists. More specifically, warming selected for the larger and most 45 productive plant species, increasing the average size within communities and leading to 46 reductions in functional plant diversity. These changes increased the total abundance of 47 bacteria that, in turn, increased that of protists, ultimately promoting soil multifunctionality. 48 Our work suggests that cascading effects between plant functional traits and the abundance 49 of multitrophic soil organisms largely regulate the response of soil multifunctionality to 50 simulated climate change, and ultimately provides novel experimental insights into the 51 mechanisms underlying the effects of biodiversity and climate change on ecosystem 52 functioning. 53 54 Keywords: bacteria, biodiversity, climate change, ecosystem functioning, environmental 55 filtering, nutrient cycles, protist, species richness. 56 4 Introduction 57 A large body of literature suggests that alterations in climate, such as warming and shifts 58 in rainfall regimes, largely impact plant communities, and the multiple ecosystem functions 59 that depend on them (Peñuelas et al., 2013). Similarly, most recent studies suggest that the 60 diversity and abundance of soil microbes are also highly vulnerable to climate change 61 (Castro, Classen, Austin, Norby, & Schadt, 2010; Maestre et al., 2015). Given the strong 62 links between plant and microbial diversity and the provision of multiple ecosystem 63 functions simultaneously (multifunctionality; Delgado-Baquerizo, Maestre, Reich, Jeffries, 64 et al., 2016; Maestre et al., 2012), alterations in the diversity of multiple trophic levels 65 (microbial and plant communities) could largely regulate the cascading effects of climate 66 change on ecosystem functioning (Jing et al., 2015). However, the indirect role of 67 multitrophic diversity in driving multifunctionality responses to climate change has not 68 been deciphered yet. Previous studies have rarely considered the indirect effect of climate 69 change on multifunctionality driven by changes in soil microbes and plant traits (Jing et al., 70 2015), nor the joint effects of simulated climate change and plant diversity on ecosystem 71 functioning (but see Kardol, Cregger, Campany, & Classen, 2010). Furthermore, most 72 studies exploring the impacts of global change drivers on ecosystem functioning have 73 focused on a single trophic level (e.g. plants or particular soil microorganisms), ignoring 74 the trophic interactions across multiple aboveground and belowground organisms taking 75 place in natural ecosystems. These gaps in our knowledge limit our understanding of the 76 direct and indirect impacts of changes in climate and biodiversity on terrestrial ecosystems, 77 and thus our ability to accurately predict the ecological consequences of global change. 78 We know that the effects of plant diversity on ecosystem functioning are largely 79 driven by variations in plant species composition and functional traits, such as height and 80 specific leaf area (Díaz et al., 2016). Moreover, the functional structure of plant 81 5 communities can, either directly (e.g., via litter quality and inputs) or indirectly (e.g., via 82 changes in abiotic factors), explain variation in soil microorganisms (de Vries et al., 2012). 83 For instance, Maestre et al. (2009) showed how the increase in height associated with shrub 84 encroachment in grasslands, and the modifications in environmental conditions associated 85 to it (e.g. reductions in soil temperature under plant canopies), affected the biomass of 86 fungi, actinomycetes and other bacteria. Plant functional traits are highly sensitive to 87 climate, with important consequences for multiple ecosystem properties (Butler et al., 88 2017; Valencia et al., 2015). The traits of the dominant species and the functional diversity 89 of plant communities may determine how plant species respond to climate change, possibly 90 via intraspecific shifts in trait values, changes in species relative abundances, and/or species 91 turnover. For example, climate change could promote species with traits that confer a better 92 adaptation to warmer conditions, such as thicker leaves (more stress-tolerant plant 93 strategies) and smaller plant size (Westoby, Falster, Moles, Vesk, & Wright, 2002), and 94 these changes would cascade to affect soil communities and ecosystem functioning (Díaz 95 et al., 2007; Grime, 1998). 96 Within the soil food web, protists occupy a key position, as they feed on bacteria 97 and yeasts, linking the flow of energy and the cycling of nutrients to higher trophic levels 98 (Bonkowski, 2004). These organisms also promote defence mechanisms and/or suppress 99 some pathogens in bacterial communities (Jousset, 2012), and their activity could explain 100 an important part of nitrogen mineralization (Trap, Bonkowski, Plassard, Villenave, & 101 Blanchart, 2016), among other processes. Despite their importance for multiple ecosystem 102 process and the control they exert on bacterial populations, protists are not routinely 103 considered in either plant-soil or global change biology studies (Geisen et al., 2017). 104 Therefore, to gain a more mechanistic understand of the ecological impacts of climate 105 6 change we must consider the ‘indirect effects’ of climate change on ecosystem functioning 106 via changes in plant functional traits and soil microorganisms, including protists. 107 To the best of our knowledge, no previous study has explicitly evaluated, using an 108 experimental approach, whether the simultaneous impacts of biodiversity and climate 109 change on soil multifunctionality are mediated by trophic level interactions across plants 110 and soil microorganisms. We aimed to do so with a full factorial mesocosm experiment 111 involving the use of open top chambers (OTCs) to increase temperature, rainfall shelters to 112 alter precipitation, and manipulated plant species richness to evaluate its direct and indirect 113 (via changes in plant functional structure and soil microbes) effects on soil 114 multifunctionality (Figure S1). First, we assessed how climate change (warming and 115 rainfall reduction) and initial plant species richness affect plant functional structure and soil 116 microorganisms. Second, we used a variance partitioning analysis to quantify how plants, 117 bacteria and bacterivorous protists alter the impacts of simulated climate change and initial 118 species richness on soil multifunctionality. Finally, we used a multi-scale path analysis (d-119 sep, Shipley, 2013) to quantify the direct and indirect effects of these different trophic levels 120 on soil multifunctionality. We expect that: (1) warming and rainfall reduction will increase 121 plant mortality, reducing its species richness (Klein, Harte, & Zhao, 2004) and promoting 122 the dominance of plants with more stress-tolerant strategies, ultimately altering the 123 diversity and abundance of soil microbes (e.g., via rizosphere interactions); (2) altered 124 microbial abundance and/or diversity linked to climate change and shifted plant attributes 125 will regulate soil multifunctionality; (3) mesocosms with the highest diversity will be 126 positively correlated with the highest soil multifunctionality (Maestre et al., 2012); (4) the 127 proportion of soil multifunctionality variance explained will substantially increase when 128 trophic level complexity is considered (Jing et al., 2015); and (5) protists –rarely studied in 129 this context- are expected to be a key intermediary of the impacts of global change on soil 130 7 multifunctionality via their role in regulating bacterial populations (Trap et al., 2016; 131 Weidner, Latz, Agaras, Valverde, & Jousset, 2017). 132 133 Materials and Methods 134 Experimental design 135 We conducted a mesocosm experiment in the Climate Change Outdoor Laboratory 136 (CCOL), established at the facilities of Rey Juan Carlos University (URJC, Móstoles, 137 Spain: 40°20’37’’N, 3°52’00’’W, 650 m a.s.l.; Figure S2), between April 2011 and 138 September 2013. Mean annual temperature and precipitation for this site are 16.6°C and 139 362 mm. 140 The experiment was designed as a fully factorial design, with three treatments: 141 initial plant richness (one, three and six species), warming (ambient and 3ºC increase) and 142 rainfall reduction (control and ~35% reduction in total annual rainfall). We selected 27 143 herbaceous species typical of grasslands from semi-arid areas in central Spain (García-144 Palacios et al., 2010) for our experiment. In each mesocosm, the species were randomly 145 selected (García-Palacios, Maestre, & Gallardo, 2011). The selected species (Table S1) 146 belong to three main functional groups (grasses, nitrogen-fixing legumes and forbs) 147 differing in traits that, such as specific leaf area and maximum plant height (Díaz et al., 148 2016), are potentially relevant to ecosystem functioning, such as biomass production, 149 resource use and N-fixation ability (Table S1; McLaren & Turkington, 2010). Because of 150 these differences, the selected species are expected to show contrasted strategies regarding 151 their responses to climate change and their impacts on ecosystem functioning (Suding et 152 al., 2008). 153 Monocultures and mixtures were established in plastic pots (depth 38 cm, internal 154 diameter 28 cm, volume 0.023 m3). The base of the pots was filled with 3 cm of expanded 155 8 clay for drainage, and then we added 32 cm of natural soil (sand content: 73.5 %, silt 156 content: 18.5 %, clay content: 8.0 %). In April 2011, we randomly sowed seeds of each 157 species within each mesocosm obtained from a commercial supplier (Intersemillas Ltd, 158 Valencia, Spain), reaching a final density of 194 individuals m-2. After seeding, all the pots 159 were buried and levelled to the ground (Figure S2b). After this burial, we added 500 mL of 160 a soil microbial inoculum to all mesocosms to recreate realistic soil microbial communities, 161 as described in Maestre et al. (2005). We irrigated all the mesocosms with 1000 mL of 162 water three times per week during the first 6 weeks of the experiment to improve seed 163 establishment. We also watered once a week in July and August 2011 to ensure the survival 164 of the plants during summer before the installation of the warming and rainfall exclusion 165 treatments, which took place in December 2011 (once all the mesocosms had an established 166 population). Individuals that failed to survive in our mesocosms once these treatments were 167 installed were not replaced. We regularly removed weeds throughout the experiment. 168 The warming treatment aimed to simulate climatic predictions for central Spain by 169 the end of the 21st century, i.e., an increase in annual temperature ranging between 2.6 - 170 4.8°C (de Castro, Martín-Vide, & Alonso, 2005; RCP8.5 scenario, IPPC, 2013). To achieve 171 this temperature increase, we used open top chambers (OTCs), described in detail in 172 Appendix 1. The rainfall exclusion treatment consisted on passive rainfall shelters, which 173 did not modify the frequency of rainfall events but reduced the total amount of rainfall 174 reaching the soil surface (35.4% ± 2.2 on average; means ± SE; n = 25 rain events). 175 Additional details about these shelters can be found in Appendix 1. Rainfall reduction 176 values obtained with them are consistent with predictions from climatic models in central 177 Spain, which forecast reductions between 10% and 50% in the total amount of rainfall 178 received during spring and fall (de Castro, Martín-Vide, & Alonso, 2005). 179 9 The effects of the OTCs and rainfall shelters on air temperature and humidity were 180 monitored using automated sensors (HOBO U23 Pro v.2 Temp/RH, Onset Corporation, 181 Bourne, MA, USA). Ten replicates per richness level (one, three or six species) were 182 maintained in all possible warming × rainfall exclusion combination (plant assemblages), 183 resulting in a total of 120 mesocosms. We eliminated a monoculture of Lotus corniculatus 184 from all climate change drivers (four mesocosms), because that species experienced a 185 severe mortality just after the end of irrigation. 186 187 Experimental measurements and harvest 188 The experiment was harvested in September 2013. We measured plant functional traits, 189 following standard protocols (Cornelissen et al., 2003) at the end of the second growing 190 season (June 2013), just after the peak of vegetation growth and before summer drought to 191 sample mature leaves and to avoid damage of leaf materials. The following traits were 192 measured on one individual per species and mesocosm: vegetative height (VH, cm), lateral 193 spread (LS, cm²), leaf area (LA, cm²), leaf length (LL, cm), leaf width (LW, cm) and leaf 194 thickness (LT, mm), phenology, measured using a phenology index (1 = no reproductive 195 stem; 2 = reproductive stem starting to grow; 3 = flowering; 4 = flower fading; 5 = fruit 196 present; and 6 = fruit absent and senescence of the reproductive stem); and specific leaf 197 area (SLA, cm² g-1). These traits were selected because they are related to water use 198 efficiency, competitive ability, light interception, water stress tolerance, leaf-level carbon 199 gain strategies, plant relative growth rate and resource capture and utilization (Westoby et 200 al., 2002; Wright et al., 2004) and are important traits in determining plant community 201 responses to climate in drylands (Gross et al., 2013). Additionally, we visually estimated 202 the total cover (cm²) per species in each mesocosm. For species with compound leaves, 203 foliar traits were measured in one leaflet per individual. 204 10 In September 2013, we collected three soil cores (0–7.5 cm depth) from each 205 mesocosm, which were bulked and homogenized to obtain a composite sample per 206 mesocosm. Soil samples were sieved (2 mm mesh) in the laboratory and separated into 207 three fractions. The first fraction was air-dried for one month for measurement of soil 208 variables related to nutrient cycles, the second fraction was stored at 4°C for estimations of 209 protist diversity and total abundance, and the last fraction was frozen at -20º C for 210 quantification of the abundance and diversity of bacterial communities. 211 212 Molecular analyses 213 Soil DNA was extracted from defrosted soil using the Powersoil DNA Isolation Kit (Mo 214 Bio Laboratories, Carlsbad, CA, USA), according to the manufacturer’s instructions. The 215 abundance of bacteria was estimated using quantitative PCRs on an ABI 7300 Real-Time 216 PCR (Applied Biosystems, Foster City, CA, USA), and the bacterial diversity (Shannon) 217 was obtained from 16s amplicon sequencing using the Illumina MiSeq platform, as 218 described in detail in Appendix 1. The abundance (total number of protists·g-1 of soil) and 219 diversity (Shannon´s index) of bacterivorous protists was determined using a modified 220 version of the liquid aliquot method, as described in detail in Appendix 1. Both protist 221 variables represent the feeding pressure on bacteria, since the method estimates the viable 222 encysted protists that accumulated in soil over time due to high bacterial production rates. 223 224 Assessing soil multifunctionality 225 We measured seven variables related to carbon (organic C, β-glucosidase activity), nitrogen 226 (total N, ammonium, nitrate) and phosphorus (available inorganic P and phosphatase 227 activity) storage and cycling. These variables constitute a good proxy for nutrient cycling, 228 biological productivity, and buildup of nutrient pools, and are important determinants of 229 11 ecosystem functioning. They also have been extensively used as proxies of ecosystem 230 functioning in many different ecosystem types (Jax, 2010). These variables were measured 231 in the laboratory as described in detail in Appendix 1. We estimated soil multifunctionality 232 from all soil variables measured using the multifunctionality index of Maestre et al. (2012). 233 Raw data for each soil variable were first normalized using a log10-transformation. 234 Following this, the Z scores of the seven soil variables were averaged to obtain soil 235 multifunctionality. This index is statistically robust (Maestre et al., 2012) and provides an 236 intuitive way and easily interpretable measure to assess changes in multifunctionality, as 237 the higher the values for the different ecosystem functions we measured, the higher the 238 multifunctionality (Byrnes et al., 2014). We acknowledge that using an a priori 239 standardized average may not allow to discriminate when all functions are performing at 240 similar levels from situations when one function could be strongly outperforming the 241 others. However, all correlations between the soil variables measured were positive (Table 242 S2), suggesting that there are not potential trade-offs between the surrogates of ecosystem 243 functioning evaluated. We also acknowledge that if two variables are highly correlated, the 244 inclusion of both provides redundancy, albeit the interpretation of soil multifunctionality 245 values is much more simple. However, among our soil variables only one out of the 25 246 correlations had had a r value higher than 0.7, suggesting that our dataset did not contain 247 high redundancy among the variables used to estimate soil multifunctionality (Table S2). 248 Finally, the comparison between our estimates of soil multifunctionality and other 249 alternative approaches revealed that all indices were strongly related, even when a multiple-250 threshold soil multifunctionality approach (Byrnes et al., 2014) was used (Figure S3). 251 Indeed the relationships between soil multifunctionality and their predictors were similar 252 when averaging soil multifunctionality (Figure S4) or multiple-threshold soil 253 multifunctionality (Figure S5) were used, hence only the former are presented in the main 254 12 text. For instance, CW-VH, CW-SLA, bacterial abundance and protist abundance all had 255 positive relationships with averaging soil multifunctionality, but the opposite was found for 256 bacterial diversity (Figure S4). The same relationships were found across all thresholds 257 when using the multiple-threshold approach (Figure S5). Additionally, we calculated 258 functions related to carbon, nitrogen and phosphorus cycles by averaging the Z scores of 259 the variables involved in each of these cycles. 260 261 Evaluating functional traits variation among plant species 262 We considered the species pool sampled in the experiment to obtain an ordination of 263 functional traits. The trait variables obtained per species and mesocosm were normalized 264 using a log transformation to reduce the impact of species with extreme trait values. We 265 then performed a principal component analysis (PCA) with Varimax rotation on trait data 266 to identify the main axes of specialization and covariation among traits. We used a PCA 267 for these analyses because different traits are associated with different niche axes 268 (Butterfield & Suding, 2013), and to avoid selecting highly correlated traits for further 269 analyses (Table S3). We identified two independent PCA components, which together 270 explained around 59% of the total variance in the data (Table S4). To represent each axis 271 of variation, we selected one trait for each PCA component as a functional marker of 272 specialization (Gross et al., 2013). The selected traits (VH and SLA) are related with two 273 main plant ecological strategies (Díaz et al., 2016). VH reflects water use efficiency and/or 274 competitive ability, whereas SLA is related with light interception and water stress 275 tolerance (Westoby et al., 2002). 276 To assess the effect of simulated climate change drivers and initial species richness 277 on plant community structure we quantified community-weighed functional traits (CWT 278 hereafter, Violle et al., 2007) and functional diversity (FD, called functional dispersion in 279 13 Laliberté & Legendre, 2010) using the trait values measured at the end of the experiment. 280 CWT indicates the “average trait value” for the selected traits, and was calculated using the 281 traits measured in each mesocosm weighted by the relative abundance (total cover) of each 282 species within the mesocosm. FD quantifies the variance of the community trait distribution 283 weighted by the relative abundance of each species within the mesocosm, and reflects the 284 functional differences between the species present in a given community. We focused on 285 VH and SLA for CWT and in all of the traits measured for FD. At the community level, 286 changes in CWT reflect a change in the trait values of the dominant species. Higher FD 287 values suggest higher niche complementarity among the species forming the community 288 (Maire et al., 2012), enhancing the use of resources, and in turn ecosystem functioning 289 (Petchey & Gaston, 2006). This index is one of the few indices that can be used when there 290 are monocultures, and it is widely used (e.g., Castro-Díez, Pauchard, Traveset, & Vilà, 291 2016; Valencia et al., 2015). 292 293 Statistical analyses 294 We first evaluated how the functional structure of plant assemblages (CWT of selected 295 traits and FD) and both the abundance and diversity of soil microorganisms (bacteria and 296 protists) responded to climate change drivers (as a combination of warming and rainfall 297 reduction treatments, resulting in 4 levels: control, warming, rainfall reduction, and 298 warming + rainfall reduction) and initial plant species richness using Linear Mixed Models 299 (LMM). In all models, we used climate change drivers, plant species richness and their 300 interaction as fixed factors, while plant assemblages (30, ten per species richness level) 301 were considered as a random factor (Kuebbing, Maynard, & Bradford, 2018). We log-302 transformed bacterial and protist abundance, and Box–Cox transformed community-303 weighed vegetative height (CW-VH), FD, bacterial diversity and protist diversity data prior 304 14 to LMM (and subsequent) analyses to improve their normality. A Tukey post-hoc analysis 305 was used to evaluate the differences within levels of the factors evaluated (climate change 306 drivers and plant species richness). 307 We tested the direct effects of climate change drivers and initial plant species 308 richness on soil multifunctionality and C, P and N nutrient cycles, using a LMM that 309 considered plant assemblages as a random factor. Note that we removed the interaction 310 between climate change drivers and initial plant species richness as a predictor in our 311 models because it did not explain additional variation. We conducted a Tukey post-hoc 312 analysis to evaluate the differences among levels of climate change drivers or plant species 313 richness, when significant responses were found. We then performed a variance 314 partitioning analysis with all predictors (full model, without model simplification), climate 315 change drivers, initial plant species richness, the functional structure of plant assemblages 316 (CWT of selected traits and FD) and the attributes of microbial communities (abundance 317 and diversity of bacteria and protists and interactions between the abundance and diversity 318 of both bacteria and protists). We conducted this analysis to quantify the relative 319 importance (unique portion of the variance) of each factor for soil multifunctionality and 320 C, P and N nutrient cycles. Finally, we assessed different hypotheses to evaluate the 321 responses of soil multifunctionality (and nutrient cycles) to different trophic levels. These 322 included: model 1) climate change drivers and initial richness; model 2) previous variables 323 and plant functional structure; model 3) previous variables and abundance and diversity of 324 bacteria; model 4) previous variables and abundance and diversity of protists; and model 325 5) previous variables and interactions between the abundance and diversity of bacteria and 326 protists. Then, we fitted all five models, and dropped terms that did not improve model fit, 327 based on the Akaike information criterion (AIC) (Akaike, 1973), until the best model 328 remained. 329 15 In addition to LMM analyses, we conducted a confirmatory path analysis using a d-330 sep approach (Shipley, 2013) to test direct and indirect causal relationships between climate 331 change drivers, plant species richness, the functional structure of plant assemblages, the 332 attributes of microbial communities and soil multifunctionality (Figure S1, Appendix 1, 333 Table S6). This approach offers a flexible way to test causal relationships between variables 334 in path analyses by relaxing some important limitations of standard structural equation 335 models, including non-normal data distribution, non-linear relationships between variables 336 and small sample sizes (Grace, 2006). This method is based on an acyclic graph that 337 summarizes the hypothetical relationships between variables to be tested using the C 338 statistic (Shipley, 2013). We selected the most appropriate predictors for soil 339 multifunctionality using the stepAIC procedure described above. Plant assemblages were 340 considered as random factors in all the models, except on those including soil 341 multifunctionality as response variable. Finally, we used standardised path coefficients to 342 measure the direct and indirect effects of predictors (Grace & Bollen, 2005). We also 343 evaluated the effect of all these predictors on C, P and N nutrient cycles, instead of soil 344 multifunctionality. All statistical analyses were conducted with R (R Core Team, 2016), 345 using different R packages (see Appendix 1). Data are available on Figshare (Valencia et 346 al., 2018). 347 348 Results 349 Our warming treatment increased air temperatures and reduced air relative humidity by 2.9 350 °C and 0.9% on average, respectively (Figure S6). Rainfall shelters decreased the amount 351 of rainfall reaching the soil by an average 35%, varying between a minimum of 9% and 352 maximum of 52% depending on the event (Figure S6). This treatment promoted a very 353 slight increase in air temperatures of 0.1 ºC compared to the control, and a reduction in the 354 16 air relative humidity of 1.2% (Figure S6). Finally, the combination of rainfall shelters and 355 warming treatment promoted an average increase of 3.2 ºC in air temperature and an 356 average reduction of 2.2% in air relative humidity (Figure S6). 357 358 Effects of climate change and plant species richness on ecosystem structure and functioning 359 Climate change drivers and plant species richness largely affected the functional structure 360 of plant assemblages (Figure 1, Table S7). We found a significant interaction between plant 361 species richness and climate change drivers when evaluating CW-VH data (Table S7). We 362 observed that plant species richness was positively correlated with high CW-VH and high 363 FD in each of the climate change drivers evaluated (Figure 1). The analyses of the effects 364 of climate change drivers at each richness level showed that warming promoted a 365 significant increase of CW-VH, but only in mesocosms with high richness (χ² = 14.85; P = 366 0.002). However, warming tended to decrease of CW-VH in mesocosms with low and 367 medium species richness (not significant relationship). Warming decreased FD regardless 368 of rainfall reduction, but this response was only significant in mesocosms with intermediate 369 richness (χ² = 11.10; P = 0.011). CW-SLA did not respond to the experimental treatments 370 evaluated (Figure 1). 371 Rainfall reduction decreased the abundance of protists (Table S7). A significant 372 plant species richness × climate change drivers interaction was found when analysing 373 protist diversity data (Table S7). Thus, the combination of warming and rainfall reduction 374 increased protist diversity in mesocosms with low richness (χ² = 11.17; P = 0.011), but 375 decreased protist diversity in mesocosms with high richness (χ² = 10.06; P = 0.018). There 376 were no significant differences between bacterial abundance and diversity in any of the 377 treatments evaluated (Figure 1, Table S7). 378 17 Mesocoms with more plant species had increased soil multifunctionality and C 379 cycling values (χ² = 7.43; P = 0.024 and χ² = 12.08; P = 0.002, respectively, see model 1 in 380 Tables 1 and S8, Figure 2). Climate change drivers significantly affected P cycling (χ² = 381 7.83; P = 0.050, Figure 2); a positive direct effect of warming plus rainfall reduction 382 compared to control was observed when analysing this variable (Figure 2, Tukey post-hoc, 383 P = 0.057). 384 385 Direct and indirect effects of climate change drivers, plants and soil microorganisms on 386 soil multifunctionality 387 The quantification of the unique portion of the variance, using all predictors, showed that 388 CW-VH, the abundance and diversity of bacteria and plant species richness were the 389 variables with higher effects on soil multifunctionality (Figure 3). The sum of their unique 390 effects explained the 37% of total variation in soil multifunctionality. Results obtained for 391 variables from the C, N and P nutrient cycles were similar, but a higher contribution of 392 CW-SLA and climate change drivers on the P cycle and protist abundance on the C cycle 393 was observed (Figure 3). 394 We evaluated the responses of soil multifunctionality to the addition of different 395 trophic levels. Remarkably, the best model included the combined effects of all trophic 396 levels studied but without the interaction between bacterial and protist communities (model 397 4; Table 1). Additionally, for the nutrient cycles the best models included the combined 398 effects of all trophic levels studied (model 4; Table S8a and S8c), except for the N cycle 399 where the best model did not include protists (model 3; Table S8b). In the case of C and N 400 cycling, models with the interaction between bacterial and protist abundances and diversity 401 had almost the same AIC and higher R2 than models without interaction (Table S8a and 402 S8b). 403 18 The model accepted in the d-sep regarding direct and indirect effects on soil 404 multifunctionality explained 68% of the total variation (Figure 4 and Table 1). Climate 405 change drivers did not have a direct effect on soil multifunctionality (Figure 4 and Table 406 1). Our path analysis revealed that initial plant species richness had a positive direct effect 407 on soil multifunctionality (Figure 4 and Table 1). However, these initial treatments had 408 multiple indirect effects on soil multifunctionality via altering plant functional structure 409 and soil microorganism community (Figure 4). For instance, climate change drivers and 410 plant species richness affect CW-VH and FD, but higher CW-VH values also had a negative 411 effect on FD. High CW-VH and low FD values were associated with higher soil 412 multifunctionality. By contrast, both climate change drivers and plant species richness did 413 not affect CW-SLA values, but increasing CW-SLA values increased soil 414 multifunctionality (Figure 4). The effect of plant species richness and climate change 415 drivers on bacterial communities was mediated by plant functional structure (Figure 4). 416 Increases in bacterial abundance were linked to high CW-VH and low FD (plant functional 417 structure), which in turn increased soil multifunctionality. Bacterial diversity had an 418 opposite response to CW-VH and FD, and was associated with a negative effect on soil 419 multifunctionality. Increases in protist abundance and diversity were mainly favoured by 420 increases in bacterial abundance and diversity, respectively, and had a positive impact on 421 soil multifunctionality. Results from C, N and P nutrient cycles were similar, but it is 422 interesting highlight the significant positive effect of protist abundance on variables from 423 the P cycle (Figure S7). 424 425 Discussion 426 Our study provides experimental evidence that the effects of diversity and simulated 427 climate change on soil multifunctionality are largely mediated by changes in the attributes 428 19 of plant, bacteria and protist communities. We found that warming impacted soil 429 multifunctionility through a complex multitrophic cascading effect. For instance, in plant 430 species-rich communities, warming selected for larger, and more productive, plants that 431 increased the abundance of soil microorganisms, ultimately enhancing soil 432 multifunctionality. The joint consideration of a wide variety of factors beyond the climate 433 and biodiversity treatments employed, including functional structure of plant assemblages 434 and the abundance and diversity of soil bacteria and protists, improved significantly the 435 proportion of explained variance on soil multifunctionality (Figure 4), and provided novel 436 insights into the cascading mechanisms underlying the effects of key global change drivers 437 on this variable. 438 439 Climate change and plant species richness alter the functional structure of plant 440 assemblages 441 Higher initial plant species richness resulted in communities with higher community-442 weighed vegetative height (CW-VH) and FD (Figures 1 and 2). This response could be 443 related to initial “sampling effect” (Tilman, Lehman, & Thomson, 1997), i.e. treatments 444 with higher species richness had higher probabilities of having keystone species. 445 Alternatively, this response may have occurs during the experiment and may reflect an 446 effect of biotic interactions (Gross et al., 2013). Since tall species are better competitors 447 than short species for space, the latter can be outcompeted. Both mechanisms are related 448 with CWT and are not exclusive. The relationship between CW-VH and richness at the 449 beginning of the experiment (equal abundances) was not significant (F = 1.57; P = 0.225), 450 suggesting that initial sampling effect is not the mechanism underlying the results observed. 451 The decrease in FD observed with warming (Figures 1 and 4) agrees with our first 452 hypothesis: climate change reduces plant diversity (Klein et al., 2004). Such result is not 453 20 surprising, but provides empirical evidence that some species could be wiped out under 454 climate change scenarios (see also Reich, 2009; Isbell et al., 2015). Also, warming had 455 contrasting effects depending on the species richness level considered. Warming promoted 456 a decrease, although not significant, of CW-VH in mesocosms with one and three species. 457 This result agrees with our expectations, since we hypothesized that temperature increase 458 and rainfall reduction would promote conservative plant strategies and smaller plant size 459 because plants with these traits are better adapted to arid or semiarid environments 460 (Westoby et al., 2002), as we had in the experimental area. However, the positive 461 significant relationship between CW-VH and warming in high plant species richness was 462 contrary to our expectations. In our study, warming favored the large and most productive 463 species, and also reduced FD, driven by the disappearance of the smaller species from the 464 mesocosms. Overall, warming resulted in higher average height and lower variability in 465 trait values at the community level. This finding indicates that warming may first shift plant 466 communities by selecting larger species, that are able to better exploit resources when 467 available and outcompete more tolerant species in mixtures. In natural communities, this 468 process may leave available niches that could be occupied by other species, but our weeding 469 treatment avoided this possibility. 470 These results highlight the importance of environmental filtering (Maire et al., 471 2012), as warming select species with a similar set of traits, suggesting that this mechanism 472 was the main driver of plant functional structure. Our findings advance our undestanding 473 on how grassland species might respond to ongoing warming. Future studies should also 474 consider the role of warming in regulating these large species in earlier developmental 475 stages (i.e. germination and seedling stage), as well as responses of other functional groups 476 that, such as non-herbaceous perennial species or annual species, may differentially 477 respond to climate change. 478 21 479 Cascading effects of plant species richness and climate change on soil microorganisms via 480 functional structure of plant assemblages 481 Warming and rainfall reduction are known to affect the abundance of soil microorganisms 482 through changes in soil temperature and moisture conditions (Castro, Classen, Austin, 483 Norby, & Schadt, 2010; Maestre et al., 2015). However, these treatments did not alter 484 bacterial abundance or diversity, and only slightly affected (i.e. explained a low amount of 485 variance) that of protists (Figures 1 and 2). Similarly, the abundance and diversity of soil 486 bacteria and protists did not respond to initial plant species richness (Figure 1). However, 487 we observed that plant species richness and climate change drivers affected bacterial and 488 protist communities via changes in the functional structure of plant communities (Figure 489 4). Mesocosms with higher plant community size (CW-VH) had higher bacterial 490 abundance, but lower bacterial diversity. Increases in CW-VH are correlated with increases 491 in aboveground plant biomass (Valencia, Méndez, Saavedra, & Maestre, 2016), which in 492 turn increases root biomass and root exudation (Eisenhauer et al., 2017; Mellado-Vázquez 493 et al., 2016). The bacterial communities on plant roots are mainly composed of copiotrophic 494 taxa with reduced diversity compared to bulk soil (Bulgarelli, Schlaeppi, Spaepen, van 495 Themaat, & Schulze-Lefert, 2013). Also, CW-VH had a positive significant relationship 496 with organic carbon (χ² = 32.61; P = <0.001), possibly because of increased amounts of 497 litter material. High ammounts of organic matter are related with high total bacterial 498 abundance and diversity (Delgado-Baquerizo, Maestre, Reich, Trivedi, et al., 2016; 499 Maestre et al., 2015). However, other studies showed a positive relationship with total 500 abundance of bacterial communities, but negative with diversity of soil bacteria, at least in 501 soils with enough carbon content (Delgado-Baquerizo et al., 2017). This agrees with our 502 results and, together with the negative relationship between bacterial abundance and 503 22 bacterial diversity observed, suggests an increase in the dominance of some bacterial 504 species/groups when CW-VH is higher, which overcompensates the decrease in the 505 abundance of other species. We propose that this process could be related to the competitive 506 exclusion principle, i.e. increases in the abundance of some species results in decreases or 507 even extinction of others (Eldridge et al., 2017). Additionally, increases in FD augmented 508 bacterial diversity, as found in other studies (Carney & Matson, 2005; Hendriks et al., 509 2013). Higher plant FD may promote new niche spaces, increasing microclimate and 510 habitat variability, allowing the survival of more types of species. Finally, all these 511 cascading effects of plant community on bacterial communities also affected protists, as 512 increases in bacterial abundance and diversity were related to increases in protist abundance 513 and diversity, respectively (Figure 3). This result supports the tight coupling of protists to 514 their bacterial food source (Bonkowski, 2004). 515 516 Direct and indirect effects of plant species richness and climate change on soil 517 multifunctionality 518 We found a direct effect of initial plant species richness, but not of climate change drivers, 519 on soil multifunctionality (Table 1, Figure 3). Our results agree with studies finding 520 enhanced soil functioning with increasing plant species richness in drylands (Maestre et al., 521 2012), and provide insights on the mechanisms underlying the interactive effects of 522 multiple trophic levels on soil multifunctionality. We found that the impact of plant species 523 richness and climate change drivers on soil multifunctionality were mostly indirect and 524 mediated by their effects on plant functional traits and soil microorganisms (Figure 4). As 525 found previously by field studies (Valencia et al., 2015), soil multifunctionality was driven 526 by changes in functional traits in our experiment. More specifically, plant communities 527 with larger plants, higher SLA, early flowering phenology and short leaves enhanced soil 528 23 multifunctionality (Figure 4, Table S3 and S4). These traits are associated to a fast-growing 529 strategy (Reich, 2014, and also CSR scheme, Grime, 1977), characterized by rapid 530 acquisition of water and other resources, as well as by the production of a dense litter layer 531 (Quero, Gómez-Aparicio, Zamora, & Maestre, 2008), which in turn reduces losses in soil 532 moisture due to evaporation (Holmgren, Gómez-Aparicio, Quero, & Valladares, 2012). 533 These large species have higher water demands, but also obtain more resources than their 534 neighbours, as they have extensive root systems, and under warming conditions they may 535 outcompete smaller plant species. This process enhances the biomass and productivity of 536 the mesocosms, but as we previously commented, decreases FD, especially under warming 537 conditions. This result is contrary to our expectations because plant FD should enhance 538 resource distribution between species and efficiency of resource use. Our results highlight 539 the importance of the functional identity of dominant species over niche complementarity 540 (Mokany, Ash, & Roxburgh, 2008), and how they affect the whole community among the 541 different trophic levels. 542 Together, our results suggest that communities dominated by fast-growing plants, 543 characterized by high SLA, and indirect cascading effects on soil microbial trophic 544 networks might result in a higher soil multifunctionality. For instance, fast-growing plants 545 promote bacterial production, via root exudation. This increases the release of N of bacteria 546 consumed by protist grazers (microbial loop; Clarholm, 1985) and potentially the release 547 of nutrients through priming of soil organic matter (Bonkowski & Clarholm, 2012), 548 improving uptake of N by mycorrhiza (Koller, Rodriguez, Robin, Scheu, & Bonkowski, 549 2013) and further increasing plant growth. Also, soil drying could periodically kill 550 significant amounts of microbial biomass, leading to larger nutrient pulses after rewetting. 551 In this context, large and fast-growing plants could take faster advantage of these nutrients 552 (e.g. dead microbial biomass). 553 24 Increases in plant FD and bacterial diversity were negatively related to soil 554 multifunctionality. Also, a higher abundance of soil bacteria and protists lead to higher soil 555 multifunctionality, regardless of their diversity. Then, the cascading effect between plant 556 and soil microorganisms and their feedbacks may promote higher sizes (biomass or 557 abundance) for the entire plant and soil microorganism community, ultimately promoting 558 soil multifunctionality. Finally, interactions between bacterial and protist communities did 559 not affect soil multifunctionality, but their inclusion increased the explained variance of N 560 cycle variables (even when the best selected model did not include them, Table S8b). 561 Despite these effects, we also detected a positive direct impact of warming and rainfall 562 reductions on P cycling variables (Figures 2 and S7). This positive effect can be related to 563 an increase in the activity of decomposers and microbes (Wolters et al., 2000). 564 Our results provide a strong evidence that climate change impacts on soil 565 multifunctionality are indirectly driven via changes in biotic attributes such as plant growth 566 rates and microbial abundance. Changes in plant functional structure promoted by warming 567 affected the microbial community and finally soil multifunctionality. Hence, grasslands 568 dominated by fast-growing plants (tall species with high SLA) favor a few groups of 569 bacteria, but these explode in numbers and have high turnover due to exudation and 570 constant protist predation, and finally could show higher levels of soil multifunctionality. 571 Higher bacterial and protist abundance positively impacted soil multifunctionality, 572 accelerating the decomposition of litter materials and increasing nutrient contents. 573 However, increases in plant and bacterial diversity promoted decreases in soil 574 multifunctionality. Overall, our results suggest that the trait values of the dominant plant 575 species, through their influence on other trophic groups, were the main driver of soil 576 multifunctionality in our experiment, and highlight the interest of evaluating multiple 577 25 trophic levels to increase our ability to explain complex ecosystem responses to climate 578 change. 579 580 Acknowledgements 581 This research was funded by the European Research Council under the European 582 Community's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 583 242658 (BIOCOM) and by the Spanish Ministry of Economy and Competitiveness 584 (BIOMOD project, Ref CGL2013-44661-R). EV is supported by the project (GACR 16-585 15012S) funded by the Czech Science Foundation and by the 2017 program for attracting 586 and retaining talent of Comunidad de Madrid (n° 2017-T2/AMB-5406). FTM 587 acknowledges support from the European Research Council (BIODESERT project, ERC 588 Grant agreement n° 647038). NG was supported by the AgreenSkills+ fellowship 589 programme which has received funding from the EU’s Seventh Framework Programme 590 under grant agreement N° FP7-609398 (AgreenSkills+ contract). CPC is supported by the 591 Estonian Research Council (project MOBJD13), the Estonian Ministry of Education and 592 Research (IUT20-29), and the European Union through the European Regional 593 Development Fund (Centre of Excellence EcolChange). M.D-B. acknowledge support 594 from the Marie Sklodowska-Curie Actions of the Horizon 2020 Framework Programme 595 H2020-MSCA-IF-2016 under REA grant agreement n° 702057. B.K.S. and M.D-B. work 596 on this topic is supported by Australian Research Council (DP 170104634). 597 26 References 598 Akaike, H. (1973). Information theory and an extension of the maximum likelihood 599 principle. In Second international symposium on information theory, ed. B. N. 600 Petrov and F. Csaki. Minnesota Studies in the Philosophy of Science, Budapest: 601 Akademinai Kiado. 602 Bonkowski, M. (2004). Protozoa and plant growth: The microbial loop in soil revisited. 603 New Phytologist, 162(3), 617–631. https://doi.org/10.1111/j.1469-604 8137.2004.01066.x 605 Bonkowski, M., & Clarholm, M. (2012). Stimulation of plant growth through interactions 606 of bacteria and protozoa: testing the auxiliary microbial loop hypothesis. Acta 607 Protozoologica, 51(3), 237–247. https://doi.org/10.4467/16890027AP.12.019.0765 608 Bulgarelli, D., Schlaeppi, K., Spaepen, S., van Themaat, E. V. L., & Schulze-Lefert, P. 609 (2013). Structure and functions of the bacterial microbiota of plants. Annual Review 610 of Plant Biology, 64(1), 807–838. https://doi.org/10.1146/annurev-arplant-050312-611 120106 612 Butler, E. E., Datta, A., Flores-Moreno, H., Chen, M., Wythers, K. R., Fazayeli, F., … 613 Reich, P. B. (2017). Mapping local and global variability in plant trait distributions. 614 Proceedings of the National Academy of Sciences of the United States of America, 615 114(51), E10937–E10946. https://doi.org/10.1073/pnas.1708984114 616 Butterfield, B. J., & Suding, K. N. (2013). Single-trait functional indices outperform 617 multi-trait indices in linking environmental gradients and ecosystem services in a 618 complex landscape. Journal of Ecology, 101(1), 9–17. https://doi.org/10.1111/1365-619 2745.12013 620 Byrnes, J. E. K., Gamfeldt, L., Isbell, F., Lefcheck, J. S., Griffin, J. N., Hector, A., … 621 27 Emmett Duffy, J. (2014). Investigating the relationship between biodiversity and 622 ecosystem multifunctionality: Challenges and solutions. Methods in Ecology and 623 Evolution, 5(2), 111–124. https://doi.org/10.1111/2041-210X.12143 624 Carney, K. M., & Matson, P. A. (2005). Plant communities, soil microorganisms, and soil 625 carbon cycling: does altering the world belowground matter to ecosystem 626 functioning? Ecosystems, 8(8), 928–940. https://doi.org/10.1007/s10021-005-0047-0 627 Castro-Díez, P., Pauchard, A., Traveset, A., & Vilà, M. (2016). Linking the impacts of 628 plant invasion on community functional structure and ecosystem properties. Journal 629 of Vegetation Science, 27(6), 1233–1242. https://doi.org/10.1111/jvs.12429 630 Castro, H. F., Classen, A. T., Austin, E. E., Norby, R. J., & Schadt, C. W. (2010). Soil 631 microbial community responses to multiple experimental climate change drivers. 632 Applied and Environmental Microbiology, 76(4), 999–1007. 633 https://doi.org/10.1128/AEM.02874-09 634 Clarholm, M. (1985). Interactions of bacteria, protozoa and plants leading to 635 mineralization of soil nitrogen. Soil Biology and Biochemistry, 17(2), 181–187. 636 https://doi.org/10.1016/0038-0717(85)90113-0 637 Cornelissen, J. H. C., Lavorel, S., Garnier, E., Diaz, S., Buchmann, N., Gurvich, D. E., … 638 Pooter, H. (2003). A handbook of protocols for standardised and easy measurement 639 of plant functional traits worldwide. Australian Journal of Botany, 51(4), 335–380. 640 https://doi.org/10.1071/BT02124 641 de Castro, M., Martín-Vide, J., & Alonso, S. (2005). El clima de España: pasado, 642 presente y escenarios de clima para el siglo XXI. Impactos del cambio climático en 643 España. 3n Moreno J.M. (Ed), Ministerio Medio Ambiente, Madrid, Spain. 644 de Vries, F. T., Manning, P., Tallowin, J. R. B., Mortimer, S. R., Pilgrim, E. S., Harrison, 645 28 K. A., … Bardgett, R. D. (2012). Abiotic drivers and plant traits explain landscape-646 scale patterns in soil microbial communities. Ecology Letters, 15(11), 1230–1239. 647 https://doi.org/10.1111/j.1461-0248.2012.01844.x 648 Delgado-Baquerizo, M., Maestre, F. T., Reich, P. B., Jeffries, T. C., Gaitan, J. J., Encinar, 649 D., … Singh, B. K. (2016). Microbial diversity drives multifunctionality in terrestrial 650 ecosystems. Nature Communications, 7, 1–8. https://doi.org/10.1038/ncomms10541 651 Delgado-Baquerizo, M., Maestre, F. T., Reich, P. B., Trivedi, P., Osanai, Y., Liu, Y.-R., 652 … Singh, B. K. (2016). Carbon content and climate variability drive global soil 653 bacterial diversity patterns. Ecological Monographs, 86(3), 373–390. 654 https://doi.org/10.1002/ecm.1216 655 Delgado-Baquerizo, M., Reich, P. B., Khachane, A. N., Campbell, C. D., Thomas, N., 656 Freitag, T. E., … Singh, B. K. (2017). It is elemental: soil nutrient stoichiometry 657 drives bacterial diversity. Environmental Microbiology, 19(3), 1176–1188. 658 https://doi.org/10.1111/1462-2920.13642 659 Díaz, S., Kattge, J., Cornelissen, J. H. C., Wright, I. J., Lavorel, S., Dray, S., … Gorné, L. 660 D. (2016). The global spectrum of plant form and function. Nature, 529(7585), 167–661 171. https://doi.org/10.1038/nature16489 662 Díaz, S., Lavorel, S., de Bello, F., Quétier, F., Grigulis, K., & Robson, T. M. (2007). 663 Incorporating plant functional diversity effects in ecosystem service assessments. 664 Proceedings of the National Academy of Sciences of the United States of America, 665 104(52), 20684–20689. https://doi.org/10.1073/pnas.0704716104 666 Eisenhauer, N., Lanoue, A., Strecker, T., Scheu, S., Steinauer, K., Thakur, M. P., & 667 Mommer, L. (2017). Root biomass and exudates link plant diversity with soil 668 bacterial and fungal biomass. Scientific Reports, 7, 44641. 669 29 https://doi.org/10.1038/srep44641 670 Eldridge, D. J., Delgado-Baquerizo, M., Travers, S. K., Val, J., Oliver, I., Hamonts, K., & 671 Singh, B. K. (2017). Competition drives the response of soil microbial diversity to 672 increased grazing by vertebrate herbivores. Ecology, 98(7), 1922–1931. 673 https://doi.org/10.1002/ecy.1879 674 García-Palacios, P., Maestre, F. T., & Gallardo, A. (2011). Soil nutrient heterogeneity 675 modulates ecosystem responses to changes in the identity and richness of plant 676 functional groups. Journal of Ecology, 99(2), 551-562. 677 https://doi.org/10.1111/j.1365-2745.2010.01765.xGarcía-Palacios, P., Soliveres, S., 678 Maestre, F. T., Escudero, A., Castillo-Monroy, A. P., & Valladares, F. (2010). 679 Dominant plant species modulate responses to hydroseeding, irrigation and 680 fertilization during the restoration of semiarid motorway slopes. Ecological 681 Engineering, 36(10), 1290–1298. https://doi.org/10.1016/J.ECOLENG.2010.06.005 682 Geisen, S., Mitchell, E. A. D., Wilkinson, D. M., Adl, S., Bonkowski, M., Brown, M. W., 683 … Lara, E. (2017). Soil protistology rebooted: 30 fundamental questions to start 684 with. Soil Biology and Biochemistry, 111, 94–103. 685 https://doi.org/10.1016/j.soilbio.2017.04.001 686 Grace, J. B. (2006). Structural equation modeling and natural systems. Cambridge 687 University Press, Cambridge, UK. 688 Grace, J. B., & Bollen, K. A. (2005). Interpreting the Results from Multiple Regression 689 and Structural Equation Models. Bulletin of the Ecological Society of America, 690 86(4), 283–295. https://doi.org/10.2307/bullecosociamer.86.4.283 691 Grime, J. P. (1977). Evidence for the existence of three primary strategies in plants and its 692 relevance to ecological and evolutionary theory. The American Naturalist, 111(982), 693 30 1169–1194. https://doi.org/10.1086/283244 694 Grime, J. P. (1998). Benefits of plant diversity to ecosystems: Immediate, filter and 695 founder effects. Journal of Ecology, 86(6), 902–910. https://doi.org/10.1046/j.1365-696 2745.1998.00306.x 697 Gross, N., Börger, L., Soriano-Morales, S. I., Le Bagousse-Pinguet, Y., Quero, J. L., 698 García-Gómez, M., … Maestre, F. T. (2013). Uncovering multiscale effects of 699 aridity and biotic interactions on the functional structure of Mediterranean 700 shrublands. Journal of Ecology, 101(3), 637–649. https://doi.org/10.1111/1365-701 2745.12063 702 Hendriks, M., Mommer, L., de Caluwe, H., Smit-Tiekstra, A. E., van der Putten, W. H., 703 & de Kroon, H. (2013). Independent variations of plant and soil mixtures reveal soil 704 feedback effects on plant community overyielding. Journal of Ecology, 101(2), 287–705 297. https://doi.org/10.1111/1365-2745.12032 706 Holmgren, M., Gómez-Aparicio, L., Quero, J. L., & Valladares, F. (2012). Non-linear 707 effects of drought under shade: reconciling physiological and ecological models in 708 plant communities. Oecologia, 169(2), 293–305. https://doi.org/10.1007/s00442-709 011-2196-5 710 Intergovernmental Panel on Climate Change (IPPC). (2013). Fifth Assessment Report: 711 Climate Change 2013: The Physical Science Basis. Cambridge University Press, 712 Cambridge, UK. 713 Isbell, F., Craven, D., Connolly, J., Loreau, M., Schmid, B., Beierkuhnlein, C., … 714 Eisenhauer, N. (2015). Biodiversity increases the resistance of ecosystem 715 productivity to climate extremes. Nature, 526(7574), 574–577. 716 https://doi.org/10.1038/nature15374 717 31 Jax, K. (2010). Ecosystem functioning. Cambridge University Press, Cambridge, UK. 718 Jing, X., Sanders, N. J., Shi, Y., Chu, H., Classen, A. T., Zhao, K., … He, J.-S. (2015). 719 The links between ecosystem multifunctionality and above- and belowground 720 biodiversity are mediated by climate. Nature Communications, 6(1), 8159. 721 https://doi.org/10.1038/ncomms9159 722 Jousset, A. (2012). Ecological and evolutive implications of bacterial defences against 723 predators. Environmental Microbiology, 14(8), 1830–1843. 724 https://doi.org/10.1111/j.1462-2920.2011.02627.x 725 Kardol, P., Cregger, M. A., Campany, C. E., & Classen, A. T. (2010). Soil ecosystem 726 functioning under climate change: plant species and community effects. Ecology, 727 91(3), 767–781. https://doi.org/10.1890/09-0135.1 728 Klein, J. A., Harte, J., & Zhao, X.-Q. (2004). Experimental warming causes large and 729 rapid species loss, dampened by simulated grazing, on the Tibetan Plateau. Ecology 730 Letters, 7(12), 1170–1179. https://doi.org/10.1111/j.1461-0248.2004.00677.x 731 Koller, R., Rodriguez, A., Robin, C., Scheu, S., & Bonkowski, M. (2013). Protozoa 732 enhance foraging efficiency of arbuscular mycorrhizal fungi for mineral nitrogen 733 from organic matter in soil to the benefit of host plants. New Phytologist, 199(1), 734 203–211. https://doi.org/10.1111/nph.12249 735 Kuebbing, S. E., Maynard, D. S., & Bradford, M. A. (2018). Linking functional diversity 736 and ecosystem processes: A framework for using functional diversity metrics to 737 predict the ecosystem impact of functionally unique species. Journal of Ecology, 738 106(2), 687–698. https://doi.org/10.1111/1365-2745.12835 739 Laliberté, E., & Legendre, P. (2010). A distance-based framework for measuring 740 functional diversity from multiple traits. Ecology, 91(1), 299–305. 741 32 https://doi.org/10.1890/08-2244.1 742 Maestre, F. T., Bowker, M. A., Puche, M. D., Belén Hinojosa, M., Martínez, I., García-743 Palacios, P., … Escudero, A. (2009). Shrub encroachment can reverse desertification 744 in semi-arid Mediterranean grasslands. Ecology Letters, 12(9), 930–941. 745 https://doi.org/10.1111/j.1461-0248.2009.01352.x 746 Maestre, F. T., Bradford, M. A., & Reynolds, J. F. (2005). Soil nutrient heterogeneity 747 interacts with elevated CO2 and nutrient availability to determine species and 748 assemblage responses in a model grassland community. New Phytologist, 168(3), 749 637–650. https://doi.org/10.1111/j.1469-8137.2005.01547.x 750 Maestre, F. T., Delgado-Baquerizo, M., Jeffries, T. C., Eldridge, D. J., Ochoa, V., Gozalo, 751 B., … Singh, B. K. (2015). Increasing aridity reduces soil microbial diversity and 752 abundance in global drylands. Proceedings of the National Academy of Sciences, 753 112(51), 201516684. https://doi.org/10.1073/pnas.1516684112 754 Maestre, F. T., Quero, J. L., Gotelli, N. J., Escudero, A., Ochoa, V., Delgado-Baquerizo, 755 M., … Zaady, E. (2012). Plant species richness and ecosystem multifunctionality in 756 global drylands. Science, 335(6065), 214–218. 757 https://doi.org/10.1126/science.1215442 758 Maire, V., Gross, N., Börger, L., Proulx, R., Wirth, C., Pontes, L. da S., … Louault, F. 759 (2012). Habitat filtering and niche differentiation jointly explain species relative 760 abundance within grassland communities along fertility and disturbance gradients. 761 New Phytologist, 196(2), 497–509. https://doi.org/10.1111/j.1469-762 8137.2012.04287.x 763 McLaren, J. R., & Turkington, R. (2010). Ecosystem properties determined by plant 764 functional group identity. Journal of Ecology, 98(2), 459–469. 765 33 https://doi.org/10.1111/j.1365-2745.2009.01630.x 766 Mellado-Vázquez, P. G., Lange, M., Bachmann, D., Gockele, A., Karlowsky, S., Milcu, 767 A., … Gleixner, G. (2016). Plant diversity generates enhanced soil microbial access 768 to recently photosynthesized carbon in the rhizosphere. Soil Biology and 769 Biochemistry, 94, 122–132. https://doi.org/10.1016/J.SOILBIO.2015.11.012 770 Mokany, K., Ash, J., & Roxburgh, S. (2008). Functional identity is more important than 771 diversity in influencing ecosystem processes in a temperate native grassland. Journal 772 of Ecology, 96(5), 884–893. https://doi.org/10.1111/j.1365-2745.2008.01395.x 773 Peñuelas, J., Sardans, J., Estiarte, M., Ogaya, R., Carnicer, J., Coll, M., … Jump, A. S. 774 (2013). Evidence of current impact of climate change on life: A walk from genes to 775 the biosphere. Global Change Biology, 19(8), 2303–2338. 776 https://doi.org/10.1111/gcb.12143 777 Petchey, O. L., & Gaston, K. J. (2006). Functional diversity: back to basics and looking 778 forward. Ecology Letters, 9(6), 741–758. https://doi.org/10.1111/j.1461-779 0248.2006.00924.x 780 Quero, J. L., Gómez-Aparicio, L., Zamora, R., & Maestre, F. T. (2008). Shifts in the 781 regeneration niche of an endangered tree (Acer opalus ssp. granatense) during 782 ontogeny: Using an ecological concept for application. Basic and Applied Ecology, 783 9(6), 635–644. https://doi.org/10.1016/J.BAAE.2007.06.012 784 R Core Team. (2016). R: A language and environment for statistical computing. R 785 Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. 786 Retrieved from https://www.r-project.org/ 787 Reich, P. B. (2009). Elevated CO2 reduces losses of plant diversity caused by nitrogen 788 deposition. Science (New York, N.Y.), 326(5958), 1399–402. 789 34 https://doi.org/10.1126/science.1178820 790 Reich, P. B. (2014). The world-wide “fast-slow” plant economics spectrum: A traits 791 manifesto. Journal of Ecology, 102(2), 275–301. https://doi.org/10.1111/1365-792 2745.12211 793 Shipley, B. (2013). The AIC model selection method applied to path analytic models 794 compared using a d-separation test. Ecology, 94(3), 560–564. 795 https://doi.org/10.1890/12-0976.1 796 Suding, K. N., Lavorel, S., Chapin, F. S., Cornelissen, J. H. C., Díaz, S., Garnier, E., … 797 Navas, M. L. (2008). Scaling environmental change through the community-level: A 798 trait-based response-and-effect framework for plants. Global Change Biology, 14(5), 799 1125–1140. https://doi.org/10.1111/j.1365-2486.2008.01557.x 800 Tilman, D., Lehman, C. L., & Thomson, K. T. (1997). Plant diversity and ecosystem 801 productivity: theoretical considerations. Proceedings of the National Academy of 802 Sciences of the United States of America, 94(5), 1857–61. 803 https://doi.org/10.1073/PNAS.94.5.1857 804 Trap, J., Bonkowski, M., Plassard, C., Villenave, C., & Blanchart, E. (2016). Ecological 805 importance of soil bacterivores for ecosystem functions. Plant and Soil, 398(1–2), 1–806 24. https://doi.org/10.1007/s11104-015-2671-6 807 Valencia, E., Gross, N., Quero, J. L., Carmona, C. P., Ochoa, V., Gozalo, B., … Maestre, 808 F. T. (2018). Data from “Cascading effects from plants to soil microorganisms 809 explain how plant species richness and simulated climate change affect soil 810 multifunctionality. Figshare. https://doi.org/10.6084/m9.figshare.6231395 811 Valencia, E., Maestre, F. T., Le Bagousse-Pinguet, Y., Quero, J. L., Tamme, R., Börger, 812 L., … Gross, N. (2015). Functional diversity enhances the resistance of ecosystem 813 35 multifunctionality to aridity in Mediterranean drylands. New Phytologist, 206(2), 814 660–671. https://doi.org/10.1111/nph.13268 815 Valencia, E., Méndez, M., Saavedra, N., & Maestre, F. T. (2016). Plant size and leaf area 816 influence phenological and reproductive responses to warming in semiarid 817 Mediterranean species. Perspectives in Plant Ecology, Evolution and Systematics, 818 21, 31–40. https://doi.org/10.1016/j.ppees.2016.05.003 819 Violle, C., Navas, M. L., Vile, D., Kazakou, E., Fortunel, C., Hummel, I., & Garnier, E. 820 (2007). Let the concept of trait be functional! Oikos, 116(5), 882–892. 821 https://doi.org/10.1111/j.2007.0030-1299.15559.x 822 Weidner, S., Latz, E., Agaras, B., Valverde, C., & Jousset, A. (2017). Protozoa stimulate 823 the plant beneficial activity of rhizospheric pseudomonads. Plant and Soil, 410(1–2), 824 509–515. https://doi.org/10.1007/s11104-016-3094-8 825 Westoby, M., Falster, D. S., Moles, A. T., Vesk, P. A., & Wright, I. J. (2002). Plant 826 ecological strategies: some leading dimensions of variation between species. Annual 827 Review of Ecology and Systematics, 33(1), 125–159. 828 https://doi.org/10.1146/annurev.ecolsys.33.010802.150452 829 Wolters, V., Silver, W. L., Bignell, D. E., Coleman, D. C., Lavelle, P., Van Der Putten, 830 W. H., … Van Veen, J. A. (2000). Effects of global changes on above- and 831 belowground biodiversity in terrestrial ecosystems: implications for ecosystem 832 functioning. BioScience, 50(12), 1089–1098. https://doi.org/10.1641/0006-833 3568(2000)050[1089:eogcoa]2.0.co;2 834 Wright, I. J., Reich, P. B., Westoby, M., Ackerly, D. D., Baruch, Z., Bongers, F., … 835 Villar, R. (2004). The worldwide leaf economics spectrum. Nature, 428(6985), 821–836 827. https://doi.org/10.1038/nature02403 837 36 838 839 37 Table 1. Results of the stepwise procedure to evaluate the effects of our experimental treatments, the functional structure of plant assemblages and 840 microbial abundance and diversity on soil multifunctionality. We tested the responses of soil multifunctionality to the following models: 1) climate 841 change drivers (warming and rainfall reduction, CCD) and initial richness, 2) plant functional structure, 3) bacterial abundance and diversity, 4) 842 protist diversity and abundance, and 5) interactions between the abundance and diversity of bacteria and protists. Grey-filled cells indicate the 843 variables that were not included in the models. CW-VH: Community Weighted mean of vegetative size, CW-SLA: Community Weighted mean 844 of specific leaf area, Est: direction of relationship; DF: degrees of freedom. 845 model 1 model 2 model 3 model 4 model 5 AIC -10.11 -75.37 -108.05 -112.12 -112.12 Adjusted R2 0.14 0.528 0.657 0.676 0.676 DF Est Fvalue Pvalue Est Fvalue Pvalue Est Fvalue Pvalue Est Fvalue Pvalue Est Fvalue Pvalue CCD 3 - - - - - - - - - - Richness 2 7.9 0.001 2.7 0.072 5.5 0.005 5.8 0.004 5.8 0.004 CW-VH 1 0.6 59.7 <0.001 0.3 10 0.002 0.3 11.5 0.001 0.3 11.5 0.001 CW-SLA 1 0.1 4.4 0.039 0.2 8.1 0.005 0.2 13.8 <0.001 0.2 13.8 <0.001 Functional diversity 1 -0.3 8.7 0.004 -0.2 5.6 0.02 -0.2 7.2 0.008 -0.2 7.2 0.008 Bacterial abundance (Ba) 1 0.3 15.7 <0.001 0.3 12.9 0.001 0.3 12.9 0.001 Bacterial diversity (Bd) 1 -0.3 16.4 <0.001 -0.3 16.8 <0.001 -0.3 16.8 <0.001 Protist abundance (Pa) 1 0.1 3 0.085 0.1 3 0.085 Protist diversity (Pd) 1 0.1 2.1 0.146 0.1 2.1 0.146 Ba:Pa 1 - - - Bd:Pd 1 - - - 846 38 List of Figures: 847 FIGURE 1. Box plot showing the effects of climate change drivers (warming and rainfall 848 reduction) and plant species richness (one, three, and six plant species) on plant functional 849 structure (a, b, c), soil bacterial abundance (d) and diversity (e), and soil protist abundance 850 (f) and diversity (g). Significant differences among each treatment combination are 851 indicated by lowercase letters (Tukey post-hoc, P < 0.05). n = 10. CW-VH: Community 852 Weighted mean of vegetative size, and CW-SLA: Community Weighted mean of specific 853 leaf area. 854 855 FIGURE 2. Box plot showing the effects of climate change drivers (warming and rainfall 856 reduction) and plant species richness (one, three, and six plant species) on soil 857 multifunctionality (a) and similar indices calculated with nitrogen (b), phosphorus (c), and 858 carbon (d) cycling variables (n = 10). 859 860 FIGURE 3. Variance components showing the unique portion of variation (percentage of 861 total R²) explained by each predictor on soil multifunctionality and variables from carbon, 862 nitrogen and phosphorus cycling. The predictors used are CCD (climate change drivers: 863 warming and rainfall reduction), plant species richness, community weighted mean of 864 vegetative size (CW-VH) and specific leaf area (CW-SLA), plant functional diversity, 865 abundance and diversity of bacterial and protist communities and the interactions between 866 them. 867 868 FIGURE 4. Relationships between climate change drivers, plant species richness, the 869 functional structure of plants, abundance and diversity of soil microorganisms and soil 870 multifunctionality. The width of each arrow is proportional to the standardized path 871 39 coefficients. Solid and dashed lines represent positive and negative effects, respectively. 872 CW-VH: Community Weighted Mean of vegetative size, CW-SLA: Community Weighted 873 Mean of specific leaf area, and FD: Functional diversity. Significance levels are as follows: 874 * p < 0.05 and ** p < 0.01. 875