domain: sdv sde sdv.bv

text:
Franciska T De Vries
email: franciska.devries@manchester.ac.uk

Rob I Griffiths

Mark Bailey

Hayley Craig

Mariangela Girlanda

Hyun Soon Gweon

Sara Hallin

Aurore Kaisermann

Aidan M Keith

Marina Kretzschmar

Philippe Lemanceau

Erica Lumini

Kelly E Mason

Anna Oliver

Nick Ostle

James I Prosser

Cecile Thion

Bruce Thomson

Richard D Bardgett

Soil bacterial networks are less stable under drought than fungal networks

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

S

oils harbour highly diverse microbial communities that are of crucial importance for soil functioning [START_REF] De Vries | Belowground connections underlying aboveground food production: a framework for optimising ecological connections in the rhizosphere[END_REF][START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF] . A major challenge is to understand how these complex microbial communities respond to and recover from disturbances, such as climate extremes, which are predicted to increase in frequency and intensity with climate change [START_REF]Climate Change 2014: Impacts, Adaptation, and Vulnerability[END_REF] . Many studies have demonstrated that climate extremes, such as drought, can have considerable effects on soil microbial communities, often with consequences for ecosystem processes and plant community dynamics [START_REF] De Vries | Controls on soil microbial community stability under climate change[END_REF][START_REF] Kaisermann | Legacy effects of drought on plant-soil feedbacks and plant-plant interactions[END_REF][START_REF] Meisner | Soil biotic legacy effects of extreme weather events influence plant invasiveness[END_REF] . It has also been shown that different components of the microbial community respond differently to drought, with soil fungi being generally more resistant, but less resilient, than bacteria [START_REF] De Vries | Controls on soil microbial community stability under climate change[END_REF][START_REF] Barnard | Responses of soil bacterial and fungal communities to extreme desiccation and rewetting[END_REF][START_REF] Bapiri | Drying-rewetting cycles affect fungal and bacterial growth differently in an arable soil[END_REF][START_REF] De Vries | Land use alters the resistance and resilience of soil food webs to drought[END_REF] . Moreover, the recovery of fungi and bacteria from drought is differentially governed by plant physiological responses to drought [START_REF] Bardgett | Hierarchical responses of plant-soil interactions to climate change: consequences for the global carbon cycle[END_REF] , for example by reducing the transfer of recently plant-assimilated C to bacteria, but not to fungi [START_REF] Fuchslueger | Experimental drought reduces the transfer of recently fixed plant carbon to soil microbes and alters the bacterial community composition in a mountain meadow[END_REF] . Although these past studies provide important insights into the impacts of climate extremes and other disturbances on soil microbial communities, they have almost exclusively focussed on single properties of soil microbial communities and their functioning, rather than on the multitude of direct and indirect interactions that occur between the networks of microbial taxa that co-exist in soil. This represents an important gap in understanding, given that climate extremes not only have the potential to reorganise networks of interactions between co-existing soil microbial taxa [START_REF] Zhou | Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2[END_REF] , but also the properties of these networks themselves could determine their response to disturbances [START_REF] De Vries | Belowground connections underlying aboveground food production: a framework for optimising ecological connections in the rhizosphere[END_REF][START_REF] De Vries | Land use alters the resistance and resilience of soil food webs to drought[END_REF][START_REF] Thompson | Dispersal governs the reorganization of ecological networks under environmental change[END_REF] .

Evidence is mounting that properties of ecological networks, which might represent interactions between co-existing organisms, can influence the response of communities to environmental change, including climate extremes [START_REF] De Vries | Belowground connections underlying aboveground food production: a framework for optimising ecological connections in the rhizosphere[END_REF][START_REF] De Vries | Land use alters the resistance and resilience of soil food webs to drought[END_REF][START_REF] Kuiper | Food-web stability signals critical transitions in temperate shallow lakes[END_REF] . Theoretical studies, for example, predict that ecological networks that consist of weak interactions are more stable than those with strong interactions [START_REF] Neutel | Stability in real food webs: weak links in long loops[END_REF][START_REF] Coyte | The ecology of the microbiome: networks, competition, and stability[END_REF] , and that compartmentalisation and presence of negative interactions increase the stability of networks under disturbances [START_REF] Coyte | The ecology of the microbiome: networks, competition, and stability[END_REF][START_REF] Rooney | Structural asymmetry and the stability of diverse food webs[END_REF][START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF] . Further, modelling studies show that increasing strength of a few key interactions within a food web can destabilise trophic cascades [START_REF] Kuiper | Food-web stability signals critical transitions in temperate shallow lakes[END_REF] . Thus, communities in which a large proportion of members are connected through positive links are deemed to be unstable; in such communities, members may respond in tandem to environmental fluctuations, resulting in positive feedback and co-oscillation [START_REF] Coyte | The ecology of the microbiome: networks, competition, and stability[END_REF] . Negative links might stabilise co-oscillation in communities and promote stability of networks [START_REF] Coyte | The ecology of the microbiome: networks, competition, and stability[END_REF] . Despite this knowledge, and an increasing use of network analysis in ecology [START_REF] De Vries | Belowground connections underlying aboveground food production: a framework for optimising ecological connections in the rhizosphere[END_REF][START_REF] Morrien | Effects of native and exotic range-expanding plant species on taxonomic and functional composition of nematodes in the soil food web[END_REF][START_REF] Pilosof | The multilayer nature of ecological networks[END_REF][START_REF] Barberan | Using network analysis to explore co-occurrence patterns in soil microbial communities[END_REF] , our understanding of cooccurrences or potential interactions within complex soil microbial communities, and how they respond to and recover from disturbances such as climate extremes, remains scant [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF][START_REF] Zhou | Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2[END_REF] .

Here, we experimentally investigated how prolonged summer drought impacts soil fungal and bacterial networks, and whether this response is consistent with the properties of these networks. We also aimed to test the role of plant communities in the response of bacterial and fungal networks and communities to drought, and what the implications are for soil functioning. Because soil bacterial communities are less resistant, but more resilient (i.e., show stronger fluctuations over time), to drought than fungal communities, we expected bacterial networks to display more destabilising properties. First, we expected correlations between bacterial OTUs to be stronger overall than those between fungal OTUs, and that a larger proportion of these correlations are negative in fungal networks. We tested these assumptions using all possible correlations in bacterial and fungal networks separately, thus including weak and non-significant interactions. Second, we expected a larger proportion of bacterial than fungal OTUs to be significantly correlated in co-occurrence networks, and these networks to be more connected for bacteria than for fungi. Third, we expected networks of positive correlations to have lower modularity for bacteria than for fungi. We tested these two latter assumptions using co-occurrence networks that only included significant, positive correlations. Positive correlations within fungal and bacterial communities can represent a range of interactions or simply indicate that they are responding in the same way to a change in environmental conditions [START_REF] Barberan | Using network analysis to explore co-occurrence patterns in soil microbial communities[END_REF][START_REF] Freilich | Species co-occurrence networks: can they reveal trophic and non-trophic interactions in ecological communities?[END_REF] . Thus, although caution is needed in interpreting them [START_REF] Freilich | Species co-occurrence networks: can they reveal trophic and non-trophic interactions in ecological communities?[END_REF][START_REF] Berry | Deciphering microbial interactions and detecting keystone species with co-occurrence networks[END_REF] , co-occurrence networks can reveal information on cooscillation of microbial taxa [START_REF] Barberan | Using network analysis to explore co-occurrence patterns in soil microbial communities[END_REF] and the stability of communities. Fourth, we expected drought-induced changes in plant growth and community composition to be associated with the trajectory of recovery of microbial networks and communities, and this link to be stronger for fungal communities, with these being the first users of belowground plant C inputs [START_REF] Ballhausen | The sapro-rhizosphere: carbon flow from saprotrophic fungi into fungus-feeding bacteria[END_REF] . Finally, given evidence for links between soil microbial communities and biogeochemical cycles [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF] , we expected drought-induced changes in both fungal and bacterial communities and networks to influence soil functioning.

The above expectations were tested using a field-based mesocosm experiment consisting of plant communities varying in relative abundance of four common grassland species with contrasting life history strategies: a fast-growing, resource exploitative grass (Dactylis glomerata) and herb (Rumex acetosa), and a slowgrowing, resource conservative grass (Anthoxanthum odoratum) and herb (Leontodon hispidus). Plant communities were dominated by one of these species and varied in evenness, such that each plant species dominated a low and medium evenness community in which the three other species had equal abundances. In addition, a high evenness treatment was included in which all species had equal abundances (see Methods and Supplementary Table 1). These treatments were intended to create a gradient of drought-induced shifts in plant communities in response to a summer drought that was imposed in the second growing season. We then quantified responses of plant communities, soil bacterial and fungal communities and their overall and co-occurrence networks, and fluxes of greenhouse gases (CO 2 and N 2 O; products of plant and microbial respiration and the microbial processes nitrification and denitrification) over time in response to drought.

Our results show that bacterial and fungal networks differ significantly in key properties that might inform on their stability under disturbance: bacterial networks show more destabilising properties than fungal networks, and these properties are stimulated by drought. We also show that a drought-induced shift to dominance of the fast-growing grass has long-lasting direct and indirect legacy effects on bacterial networks and communities, with the potential to reinforce changes in plant community composition and affect aboveground and belowground responses to future disturbances.

Results

Bacterial and fungal communities. Both fungal and bacterial communities were strongly, but contrastingly, affected by drought (Fig. 1). Fungal richness and evenness increased during the drought period, but rapidly recovered to control levels 1 week after rewetting. In contrast, drought decreased bacterial richness and evenness; this effect was strongest 1 week after rewetting and persisted for 2 months. We used Bray-Curtis similarity as a measure of resilience of bacterial and fungal communities. When similarities between drought and control communities were significantly lower than those before drought, communities were assumed to be affected by drought and had not recovered. When similarities did not differ significantly, the communities were assumed to be unaffected or had recovered [START_REF] Shade | Resistance, resilience and recovery: aquatic bacterial dynamics after water column disturbance[END_REF] . We found that similarities between bacterial drought and control communities were significantly lower than they were before drought at all postdrought samplings, indicating an ongoing effect of drought on those communities. In contrast, similarities between fungal drought and control communities were only lower than before drought at the end of the drought, and weakly so (Fig. 1). In addition, we also found that a larger proportion of fungal OTUs increased in relative abundance in response to drought than did bacterial OTUs (Fig. 2a and Supplementary Fig. 1).

We identified individual bacterial and fungal indicator OTUs that were both high in abundance and responded strongly in their abundance to drought (Methods, Supplementary Fig. 1). We also distinguished between drought-sensitive and drought-tolerant indicators, which respectively decreased and increased in response to drought. Both drought-tolerant and droughtsensitive bacterial indicators were mainly found within the phylum Verrucomicrobia and the class Alphaproteobacteria, and additional bacterial drought-sensitive indicators belonged to the phylum Actinobacteria (Supplementary Data 1). Fungal droughttolerant indicators belonged to the phyla Ascomycota and Glomeromycota (Supplementary Data 2), whereas all fungal drought-sensitive indicators belonged to the phylum Zygomycota, specifically of the family Mortierellaceae. 1 Bacterial and fungal community response to drought. The duration of the drought is indicated as a shaded area in b-d and f-h. NMDS of bacterial (a) and fungal (e) community composition shows that drought significantly affected bacterial and fungal communities at all sampling dates following drought (ADONIS). Bacterial (b, c) and fungal (f, g) richness and evenness were strongly affected by drought, but this response differed over time (Repeated measures ANOVA sampling × Drought interaction F 3,153 = 10.4, P < 0.001 and F 3,195 = 14.8, P < 0.001 for richness of bacteria and fungi, respectively; Sampling × Drought interaction F 3,153 = 23.1, P < 0.001 and F 3,195 = 23.1, P < 0.001 for evenness of bacteria and fungi, respectively). The similarity between drought and control communities changed over time for both bacteria and fungi, but bacterial drought and control communities were still less similar than before drought at late recovery (ANOVA sampling effect F 3,80 = 40.9, P < 0.001 and F 3,103 = 2.8, P = 0.04 for bacteria and fungi, respectively). In b-d and f-h, dots represent individual observations and lines indicate means ± 1 SE. In d, h, symbols with the same letter are not statistically different Bacterial and fungal networks. When considering all correlations, the links between bacterial OTUs were consistently stronger than those between fungal OTUs (Fig. 3, ANOVA; F 1,18,778,692 = 5577, P < 0.0001), supporting our expectation that potential interactions in bacterial networks are stronger than those of fungal networks. We also found that fungal networks consistently had fewer negative correlations than bacterial networks (Fig. 3, two-sided χ 2 -test of proportions; χ 2 > 500 and P < 0.0001 for all fungal-bacterial control network pairs). When considering only significant correlations (ρ > 0.6, P < 0.01), again we found that fungal networks contained fewer negative correlations than bacterial networks (two-sided χ 2 -test of proportions; χ 2 > 10 and P < 0.001 for all fungal-bacterial control network pairs). Drought reduced the proportion of negative correlations in bacterial networks at the end of the drought, and in fungal networks after 2 months of recovery (two-sided χ 2 -test of proportions; χ 2 = 9.6, P = 0.002 and χ 2 = 9.5, P = 0.002, respectively).

We then used these significant correlations to construct cooccurrence networks, consisting only of positive correlations. Bacterial co-occurrence networks were larger, more connected and less modular than fungal networks. A larger proportion of bacterial OTUs was included in co-occurrence networks than of fungal OTUs (Fig. 2b). The larger bacterial co-occurrence networks also had a higher node-normalised degree (the number of connections a node has standardised by the total number of connections in the network [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] ) and betweenness (the number of paths through a node [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] ) than fungal networks (Figs. 4, 5), while a lower clustering coefficient [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] indicated a marginally lower modularity in bacterial networks than in fungal networks (ANOVA; F 1,8 = 3.9, P = 0.08, Supplementary Table 2).

Membership of bacterial networks was more constant over time than that of fungal networks, indicating consistent correlations between taxa (Supplementary Fig. 2). In addition, drought increased the connectedness and centrality of nodes in bacterial networks, while it decreased these properties in fungal networks (Fig. 5). When we analysed the combined co-occurrence networks, including both fungi and bacteria, networks were dominated by bacterial nodes and showed similar dynamics to bacteria-only networks (Supplementary Fig. 3 and Supplementary Table 2).

We found that indicator OTUs, regardless of whether they were drought-sensitive or drought-tolerant, were more connected than non-indicator OTUs in bacterial networks (Fig. 5), while in fungal networks drought-sensitive indicators were more central and connected than non-indicators. In particular, highly central in bacterial drought networks were drought-tolerant taxa of the genera DA101 and Candidatus xiphinematobacter (Verrucomicrobia) and Rhodoplanes (Alphaproteobacteria) (Supplementary Data 1). We also found that the centrality (betweenness) and connectedness (normalised degree) of bacterial network OTUs was positively related to their relative abundance, but only in drought networks (regression; R 2 = 0.016, P = 0.0002, df = 884 for betweenness in drought late recovery network; R 2 = 0.037, P = 0.0002, df = 367; R 2 = 0.010, P = 0.002, df = 884, for normalised degree in drought early and late recovery networks). In contrast, these properties were strongly positively related to OTU relative abundance in fungal control networks. Before the simulated drought, connectedness was predicted by OTU abundance in both fungal control and drought networks (regression; R 2 = 0.37 and P < 0.0001, df = 86 and R 2 = 0.22 and P < 0.0001, df = 111, respectively), and both centrality and connectedness were predicted by relative abundance in the late recovery control network (regression; R 2 = 0.49 and P < 0.0001, and R 2 = 0.24 and P < 0.0001, df = 101 for betweenness and normalised degree, respectively). For most bacterial networks, there was a negative relationship between OTU relative abundance and the significance (P value) of their response to drought. This indicates that the most abundant bacterial OTUs showed the strongest response to drought (regression; R 2 > 0.01 and P < 0.001 in end of drought and late recovery control and drought networks; R 2 = 0.024 and P = 0.003 in early recovery drought network); these relationships were absent in fungal networks.

Links between plant and microbial community properties.

Drought caused a strong shift in plant community composition (Fig. 6a-d). Biomass of the fast-growing grass D. glomerata increased significantly relative to other plant species, and this increase was responsible for a divergence between control and drought plant communities (Fig. 6e,f). We assessed the extent to which this increased D. glomerata biomass was associated with the relative abundances of bacterial and fungal network OTUs during recovery from drought. We found more significant correlations of D. glomerata biomass with OTUs in bacterial networks, and these correlations were stronger than those with fungal network OTUs (Fig. 7). Moreover, in bacterial drought networks, the strength of the correlations with D. glomerata increased when nodes were more central and connected and might thus be potential keystone taxa [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] (Fig. 7). We also found that the resilience of the bacterial community to drought (calculated as pairwise Bray-Curtis similarities between control and drought communities) was linked to the resilience of the plant community, but this relationship was not significant for fungal communities (Fig. 6g,h).

We further inferred the potential mechanisms through which plant community change affected fungal and bacterial community composition as well as functional guilds involved in denitrification and N 2 O reduction by constructing a structural equation model (SEM) for the final, late recovery sampling (Supplementary Fig. 4 and Supplementary Note 1). We found that the increase in D. glomerata biomass was directly associated with bacterial community composition. The increase in D. glomerata was also indirectly associated with fungal and bacterial community composition through reducing soil moisture content, which was caused by an increase in aboveground biomass and thus evapotranspiration (Fig. 8). The abundances of nir and nosZ genes, used as proxies for denitrifiers and nitrous oxide reducers, respectively, were not predicted by changes in plant community composition. Surprisingly, soil inorganic and organic dissolvable nitrogen, and dissolved organic carbon were not associated with bacterial or fungal community composition, despite our observation that the D. glomerata-associated increase in aboveground biomass reduced soil inorganic N availability (regression; R 2 = 0.09, P = 0.010).

Links with ecosystem functioning. We constructed a SEM to assess the direct and indirect effects of drought on soil functioning (Supplementary Fig. 5 and Supplementary Note 2). We used a multigroup modelling approach to assess how our hypothesised relationships between plant community composition, microbial communities, and CO 2 and N 2 O fluxes changed during drought and recovery [START_REF] De Vries | Land use alters the resistance and resilience of soil food webs to drought[END_REF] . Drought strongly and directly affected both fungal and bacterial community composition at the end of the drought through changing soil moisture, but plant community composition only affected bacterial community composition directly (Fig. 9). Neither fungal nor bacterial community composition were linked to ecosystem respiration, but bacterial community composition was linked to the abundance of nir genes. However, 1 week after rewetting, bacterial but not fungal community composition was directly associated with both ecosystem respiration and N 2 O production, and indirectly with N 2 O production through its link with the abundance of nosZ genes. As expected, the abundance of nir genes was linked to increased N 2 O production, while nosZ gene abundance linked to reduced N 2 O production at early recovery. At the final sampling, 2 months after ending the drought, this association of nir and nosZ gene abundances ratio with N 2 O fluxes disappeared (Fig. 9), and all rates of ecosystem processes had returned to control levels (Supplementary Fig. 6) irrespective of drought, nir and nosZ gene abundances increased during the duration of the experiment (Supplementary Fig 7). The N 2 O reducers were dominated by nosZ clade I indicating that complete denitrifiers rather than nondenitrifying N 2 O reducers were promoted in the experiment [START_REF] Graf | Intergenomic comparisons highlight modularity of the denitrification pathway and underpin the importance of community structure for N2O emissions[END_REF] .

We did not find any relationship between bacterial and archaeal ammonia oxidiser gene abundances and N 2 O emissions (see Supplementary Fig. 8 for drought effects on amoA gene abundances).

Discussion

Our findings provide insight into the immediate and delayed response of belowground networks to drought and highlight the role of plant community composition in governing the dynamics of this response. We show that soil bacterial and fungal networks have different properties and respond differently to drought, with drought having a much stronger impact on bacterial than on fungal networks. We also found that drought-induced changes in plant communities had long-lasting associations with bacterial networks and communities and strongly governed their recovery, but much less so for fungal networks and communities. Bacterial co-occurrence networks were characterised by properties that indicate low stability under disturbance, such as high connectivity and centrality, and low modularity, while fungal networks had properties that suggest higher stability. In bacterial networks, abundant, drought responsive indicator OTUs were highly central and connected, regardless of whether they were drought-sensitive or drought-tolerant (Fig. 5), suggesting that they might drive the observed drought-induced changes in bacterial networks [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF][START_REF] Shi | The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages[END_REF][START_REF] Steele | Marine bacterial, archaeal and protistan association networks reveal ecological linkages[END_REF] . The central indicators of drought in the bacterial networks, Rhodoplanes and DA101, are highly abundant in soils in general, although little is known about their ecology [START_REF] Brewer | Genome reduction in an abundant and ubiquitous soil bacterium 'Candidatus Udaeobacter copiosus[END_REF] . However, their drought tolerance might be a result of their ability to maintain dominance under a wide range of environmental conditions. The most dominant bacterial taxa were the strongest responders to drought (negative correlation between abundance and P value of drought response), which is in contrast to previous studies [START_REF] Zhao | Altered precipitation seasonality impacts the dominant fungal but rare bacterial taxa in subtropical forest soils[END_REF][START_REF] Sayer | Links between soil microbial communities and plant traits in a species-rich grassland under long-term climate change[END_REF] , and we found that these dominant taxa drive network structure (positive correlations between connectedness and betweenness, and abundance). Thus, the most abundant OTUs were driving bacterial, but not fungal, network reorganisation in response to drought [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF][START_REF] Shi | The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages[END_REF][START_REF] Steele | Marine bacterial, archaeal and protistan association networks reveal ecological linkages[END_REF] . While the relationships between bacterial phylogeny and function are complex [START_REF] Philippot | The ecological coherence of high bacterial taxonomic ranks[END_REF][START_REF] Martiny | Phylogenetic conservatism of functional traits in microorganisms[END_REF] , shifts in the abundance of indicator taxa might inform on the stability of these networks, and consequently on the response of soil bacterial communities to drought.

Drought caused a long-lasting shift in plant community composition, and this shift was strongly associated with bacterial Fig. 5 Node connectedness and centrality of bacterial and fungal networks. In bacterial networks, drought increased connectedness (normalised degree) and centrality (betweenness) at post-drought samplings (a, b Sampling × Drought interaction F 3,5973 = 67.3, P < 0.001 and F 3,5973 = 107.5, P < 0.001, respectively), while in fungal networks, drought did either not affect, or decrease these properties (e, f Sampling × Drought interaction F 3,833 = 21.6, P < 0.001 and F 3,833 = 1.4, P = 0.234 for normalised degree and betweenness, respectively). In bacterial networks, the normalised degree of droughtsensitive and drought-tolerant indicator OTUs was higher than those of non-indicator OTUs (c F 2,5977 = 12.5, P < 0.001), but betweenness was not affected (d). In fungal networks, drought-sensitive OTUs had a higher normalised degree and betweenness than non-indicator OTUs (g, h F 2,834 = 10. communities and networks. Fungal communities have been found to be more responsive than bacterial communities to vegetation change [START_REF] Dassen | Differential responses of soil bacteria, fungi, archaea and protists to plant species richness and plant functional group identity[END_REF] , and there is accumulating evidence that fungi are the first consumers of belowground inputs of plantderived C [START_REF] Ballhausen | The sapro-rhizosphere: carbon flow from saprotrophic fungi into fungus-feeding bacteria[END_REF][START_REF] De Deyn | Rapid transfer of photosynthetic carbon through the plant-soil system in differently managed species-rich grasslands[END_REF] . However, our findings point to a novel mechanism by which vegetation change can affect soil microbial community reorganisation via changes in soil moisture content. Drought increased the abundance of the fast-growing grass D. glomerata [START_REF] De Vries | Grassland species root response to drought: consequences for soil carbon and nitrogen availability[END_REF] , resulting in increased total aboveground biomass, which in turn resulted in a prolonged reduction in soil moisture due to higher evapotranspiration. This increase in D. glomerata is likely explained by the tolerance of this species to drought and its ability to capitalise on the flush of N that occurred on rewetting, thereby enforcing its dominance relative to other species [START_REF] Hoover | Resistance and resilience of a grassland ecosystem to climate extremes[END_REF] . Our SEM (Fig. 8) showed that while plant community induced changes in soil moisture affected both bacterial and fungal community composition 2 months after the drought had ended, only bacterial community composition was directly affected by the biomass of D. glomerata. This corresponds with our observation that D. glomerata biomass was linked to central and connected OTUs in bacterial networks, but not in fungal networks. In addition, both the immediate and delayed reductions in moisture affected bacterial communities more than they did fungal communities. Thus, drought can have long-lasting legacy effects on soil microbial communities by promoting the dominance of a fastgrowing grass. Such responses could potentially be common given that increases in abundance of fast-growing grasses are often reported following drought, contributing to the maintenance of grassland productivity [START_REF] Hoover | Resistance and resilience of a grassland ecosystem to climate extremes[END_REF][START_REF] Dreesen | Summer heat and drought extremes trigger unexpected changes in productivity of a temperate annual/biannual plant community[END_REF][START_REF] Hofer | Yield of temperate forage grassland species is either largely resistant or resilient to experimental summer drought[END_REF] . While these vegetation-mediated legacy effects on soil microbial communities did not result in long-term-altered soil functioning, they were shown in a related study to feedback positively on growth of D. glomerata 5 , thereby potentially reinforcing its dominance. As such, changes in soil microbial communities resulting from drought have the potential to have legacy effects on plant community composition, in this case triggering a positive feedback on D. glomerata abundance [START_REF] Kaisermann | Legacy effects of drought on plant-soil feedbacks and plant-plant interactions[END_REF] , potentially contributing to long-term-altered plant community composition.

Not only were bacterial communities affected more by drought, both directly and indirectly through vegetation change, but bacterial communities also showed stronger links to CO 2 and N 2 O fluxes than did fungal communities during drought recovery. As expected, bacterial communities predicted nosZ gene abundances [START_REF] Graf | Intergenomic comparisons highlight modularity of the denitrification pathway and underpin the importance of community structure for N2O emissions[END_REF] , and both N 2 O emissions and abundances of genes indicative of denitrifier genes (nir and nosZ clade I) were sharply increased 1 week after rewetting. These dynamics, together with the low variance of N 2 O explained in our SEM, suggest that the flush in inorganic N availability after rewetting (which is not included in the SEM) caused growth of denitrifiers [START_REF] Leitner | Linking NO and N2O emission pulses with the mobilization of mineral and organic N upon rewetting dry soils[END_REF] that did not use the final step in the pathway under N replete conditions. Thereby, N 2 O emission rates increased, but once available N content decreased, the denitrifiers likely used the entire pathway and terminated with nitrogen gas. This would explain why gene abundances remained high at the final sampling, while N 2 O production returned to control levels. Our observations are in contrast with previous studies that show a reduction or no response of denitrifier genes in response to drought [START_REF] Gschwendtner | Climate change induces shifts in abundance and activity pattern of bacteria and archaea catalyzing major transformation steps in nitrogen turnover in a soil from a mid-European Beech Forest[END_REF][START_REF] Hartmann | Effects of drought and N-fertilization on N cycling in two grassland soils[END_REF] . Moreover, we show that the indirect effects of drought can increase both denitrifier and ammonia-oxidiser gene abundances, potentially 'priming' N cycling communities to respond to increased N availability after a second drought.

Our findings have important implications for understanding how complex soil microbial communities respond to climate extremes and other disturbances. While the exact nature of inferred bacterial and fungal interactions remains unknown, our data show that drought promoted destabilising properties of bacterial co-occurrence networks. We also reveal an important association between drought-induced shifts in plant community composition and changes in bacterial communities, bacterial networks and the relative abundance of denitrification and nitrous oxide reduction genes long after the drought has ended. These findings suggest that bacterial communities might not be as resilient as previously thought, and that changes in vegetation resulting from drought could have long-lasting effects on soil bacterial communities, potentially influencing plant community composition and the ability of aboveground and belowground communities to withstand future disturbances [START_REF] Kaisermann | Legacy effects of drought on plant-soil feedbacks and plant-plant interactions[END_REF][START_REF] Meisner | Soil biotic legacy effects of extreme weather events influence plant invasiveness[END_REF][START_REF] De Vries | Land use alters the resistance and resilience of soil food webs to drought[END_REF][START_REF] Kroel-Dulay | Increased sensitivity to climate change in disturbed ecosystems[END_REF] .

Methods

Experimental set-up and soil and greenhouse gas analyses. Topsoil (silt loam of the Brickfield 2 association (Avis & Harrop, 1983), pH 6.2, C and N content 3.13 and 0.25 g kg -1 , respectively) was collected from mesotrophic grassland at Hazelrigg Field Station, Lancaster University, UK (54°1′N, 2°46′W, 94 m a.s.l), where the experiment was conducted. Average yearly rainfall at the site is 1092 mm, average yearly temperature is 9.3 °C, based on daily site measurements since 1966. After homogenisation and removal of stones and roots, 72 pots of 42L (38 × 38 cm, 40 cm depth) were filled with soil in May 2012 and planted with nine different plant communities consisting of 36 individuals, arranged in a randomised block design. Each community consisted of four plant species (Anthoxanthum odoratum, Dactylis glomerata, Leontodon hispidus and Rumex acetosa, germinated from seed in glasshouse conditions) and plant communities differed in evenness (low, medium and high) and dominant plant species (each species alternating in the low and medium evenness treatment; the high evenness treatment had equal abundances of all species; see Supplementary Table 2). Plant communities were left for two growing seasons, and during the second growing season, an extreme drought designed to simulate a 100-year drought event was imposed in a splitblock design by placing transparent rain covers over the pots from May to July 2013, following Bloor and Bardgett [START_REF] Bloor | Stability of above-ground and below-ground processes to extreme drought in model grassland ecosystems: Interactions with plant species diversity and soil nitrogen availability[END_REF] . All plant community treatment and drought combinations were replicated four times. Local weather data (1967-2008) were used to fit a Gumbel I distribution to the annual extremes of drought duration for the local growing period. During the experimental drought, 75.6 mm of rainfall was excluded in June, 84.9 in July and 131.1 mm in August (for daily rainfall and temperature data during the experimental period, see Supplementary Fig. 12). All pots, and control treatments during the drought, were supplemented with 700 ml of additional water (equivalent to 0.5 mm of rainfall) every 2 days in dry weeks. Aboveground biomass was cut to 4 cm two times between plant community establishment and start of the drought to simulate grazing (end of growing season 2012 and 1 week before starting the drought in 2013), and the biomass of the four individual species was dried (70 °C) and weighed for each pot. Samples were taken to quantify microbial properties and gaseous fluxes immediately before starting the drought, at the end of the drought, and 1 week and 2 months after ending the drought. At each sampling date, five randomly distributed soil samples per pot were taken (1 cm diameter, 0-10 cm depth), pooled, sieved (2 mm) and kept at -80 °C until processing. CO 2 fluxes (ecosystem respiration and net ecosystem exchange) were measured using a portable Infra Red Gas Analyser (EGM4, PP Systems) using a blacked-out and transparent chamber (to prevent and allow photosynthesis), respectively. Photosynthesis rates were calculated as net ecosystem exchange minus respiration rates, with negative values representing an uptake of CO 2 . N 2 O fluxes were measured using a static chamber approach as in Ward et al. [START_REF] Ward | Warming effects on greenhouse gas fluxes in peatlands are modulated by vegetation composition[END_REF] . Briefly, dark, airtight chambers fitted with a septum for gas collection were placed over each mesocosm, and headspace was sampled manually every 10 min for 30 min (four time points) using a syringe (10 ml) fitted with a needle. Samples were injected into 3 ml vacuum exetainers (Labco, Lampeter, UK) and analysed by gas chromatography using Autosystem XL GCs (Perkin Elmer, Waltham, MA, USA). Fluxes were adjusted for field temperature when sampling, headspace volume and chamber area, and calculated by linear regression using all time points sampled. At the final sampling (September 2013), a full sampling of plant and soil properties was done. Aboveground vegetation was cut to 4 cm, split per species, dried (70 °C) and weighed. Three soil cores (3 cm diam, 10 cm depth) were pooled per pot, homogenised, sieved (4 mm) and used for analysis of water soluble soil carbon (dissolved organic C-DOC) and nitrogen (NO 3 -, NH 4 + , dissolved organic N-DON) pools. For dissolved organic C and N, 5 g of soil was extracted using 35 ml of Milli-Q water, extracts were filtered (0.45 µm), and analysed using a 5000A TOC analyser (Shimadzu, Japan) (DOC) and a AA3 HR Auto Analyzer (Seal Analytical, UK) (DON). For inorganic N, 5 g of soil was extracted using 25 ml of 1 M KCl, filtered through Whatman no 1, and analysed on a AA3 HR Auto Analyzer (Seal Analytical, UK). One core (3 cm diam, full depth of pot) was used for determination of total root biomass and root trait analysis. Structural root traits were analysed using WinRhizo® root analysis software (Regent Instruments Inc., Canada) and an Epson flatbed scanner. After analysis, roots were blotted dry, weighed, dried at 70 °C and re-weighed to calculate root dry matter content (RDMC). Specific root length (SRL) was calculated by dividing the dry biomass by the total root length (cm g -1 ); root tissue density (RTD) was calculated by dividing the weight of the dry biomass by the root volume (g cm -3 ).

Quantitative PCR of nir and nosZ genes. DNA was extracted from 0.3 g of soil using the MoBIO PowerSoil-htp 96-Well DNA Isolation kit (Carlsbad, CA) according to the manufacturer's protocols. The DNA quality was checked by agarose gel electrophoresis and prior to performing quantitative PCR, the quantity was measured using a Qubit fluorimeter (Invitrogen, Carlsbad, CA). To determine the genetic potential for denitrification and N 2 O reduction, the genes nirS and nirK, coding for the cytochrome-like nitrite reductase and the copper-dependent nitrite reductase in denitrifiers, and nosZI and nosZII, coding for the nitrous oxide reductase from clade I and II, were quantified by quantitative real-time PCR based on SYBR green detection using a Biorad CFX Connect Real-Time System (Biorad) and according to the amplification protocols in Supplementary Table 3. Each reaction had a volume of 15 µl using iQ™ SYBR® Green Supermix (Bio-Rad Laboratories, Inc.), 0.05% bovine serum albumin, 0.5 µM of each primer except 0.25 µM for the nirK primers, and 5 ng DNA. Standard curves for each gene were obtained by serial dilutions of linearised plasmids containing fragments of the respective genes, which were obtained from pure cultures. Standard curves were linear (R 2 = 0.990) in the range used. The amplifications were verified by melting curve analyses and agarose gel electrophoreses, and non-template controls resulted in negligible values. Prior to quantification, potential inhibition of the PCR reactions was checked for each sample by amplifying a known amount of the pGEM-T plasmid (Promega) with the plasmid-specific T7 and SP6 primers when added to the DNA extracts or non-template controls. No inhibition of the amplification reactions was detected with the amount of DNA used.

Quantitative PCR of archaeal and bacterial amoA genes. Archaeal and bacterial amoA gene abundances were determined by real-time PCR in a RealPlex2 Mastercycler (Eppendorf, Stevenage, UK) using QuantiFast® SYBR® Green PCR Master Mix (Qiagen, Crawley, UK), as described by Thion and Prosser [START_REF] Thion | Differential response of nonadapted ammoniaoxidising archaea and bacteria to drying-rewetting stress[END_REF] . Briefly, 20 µl final volume reaction comprised of 0.2 mg ml -1 bovine serum albumin (BSA), 1.5 µM of both primers, 10 µl of QuantiFastTM qPCR master mix (Qiagen, Crawley, UK) and 2 µl DNA template, diluted to reach 5-10 ng µl -1 . Crena-moA23f/CrenamoA616r 48 and amoA-1F/amoA-2R 49 primers were used to amplify AOA and AOB amoA genes, respectively (Supplementary Table 3). Two complete standard dilution series, prepared as described by Thion and Prosser 47 , from 10 1 to 10 7 AOA or AOB amoA gene copies, were run for every real-time PCR reaction. Amplification efficiency ranged from 83% to 94% and from 92% to 101% for AOB and AOA, respectively, and only reactions with r 2 values ≥0.98 were considered.

16S and ITS amplicon sequencing. Amplicon libraries were constructed according to the dual indexing strategy of Kozich et al. [START_REF] Kozich | Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform[END_REF] , with each primer consisting of the appropriate Illumina adaptor, 8-nt index sequence, a 10-nt pad sequence, a 2-nt linker and the gene-specific primer. For bacteria, the V3-V4 hypervariable regions of the 16S rRNA gene were targeted using primers 341F 51 and 806R [START_REF] Caporaso | Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample[END_REF] ; for fungi, the ITS2 region was amplified using fITS7f and ITS4r primer sequences described by Ihrmark et al. [START_REF] Ihrmark | New primers to amplify the fungal ITS2 region-evaluation by 454-sequencing of artificial and natural communities[END_REF] . Amplicons were generated using a high-fidelity DNA polymerase (Q5 Taq, New England Biolabs) and pooled. PCR was conducted on 20 ng of template DNA employing an initial denaturation of 30 s at 95 °C, followed by (25 for 16S and 30 cycles for ITS) 30 s at 95 °C, 30 s at 52 °C and 2 min at 72 °C. A final extension of 10 min at 72 °C was also included to complete the reaction. Amplicon sizes were determined using an Agilent 2200 TapeStation system (~550 bp: 16S; ~350-425: ITS) and libraries normalised using SequalPrep Normalisation Plate Kit (Thermo Fisher Scientific). Library concentration was calculated using a SYBR green quantitative PCR (qPCR) assay with primers specific to the Illumina adaptors (Kappa, Anachem). Libraries were sequenced at a concentration of 5.4 pM with a 0.6 pM addition of an Illumina-generated PhiX control library. Sequencing runs, generating 2 × 300 bp reads, were performed on an Illumina MiSeq using V3 chemistry. The read 1 (R1), read 2 (R2) and index sequencing primers used were also 16S specific: R1 = sequence of the combined pad, linker and 341F; R2 = sequence of the combined pad, linker and 806R; I = reverse compliment of the R2 primer.

Bioinformatic and statistical analyses of sequencing data. Sequenced pairedend reads were joined using PEAR [START_REF] Zhang | PEAR: a fast and accurate Illumina Paired-End reAd mergeR[END_REF] , quality filtered using FASTX tools (hannonlab.cshl.edu), length filtered with the minimum length of 300 bpsm, presence of PhiX and adaptors were checked and removed with BBTools (jgi.doe.gov/data-andtools/bbtools/), and chimeras were identified and removed with VSEARCH_-UCHIME_REF [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF] using Greengenes Release 13_5 (at 97%) [START_REF] Desantis | a chimera-checked 16S rRNA gene database and workbench compatible with ARB[END_REF] . Singletons were removed and the resulting sequences were clustered into operational taxonomic units (OTUs) with VSEARCH_CLUSTER 55 at 97% sequence identity [START_REF] Tindall | Notes on the characterization of prokaryote strains for taxonomic purposes[END_REF] . Representative sequences for each OTU were taxonomically assigned by RDP Classifier with the bootstrap threshold of 0.8 or greater [START_REF] Wang | Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy[END_REF] using the Greengenes Release 13_5 (full) [START_REF] Desantis | a chimera-checked 16S rRNA gene database and workbench compatible with ARB[END_REF] as the reference. Unless stated otherwise, default parameters were used for the steps listed. The fungal ITS sequences were analysed using PIPITS [START_REF] Gweon | PIPITS: an automated pipeline for analyses of fungal internal transcribed spacer sequences from the Illumina sequencing platform[END_REF] with default parameters as outlined in the citation. Briefly, this involved quality filtering and 97% clustering of the ITS2 region as indicated above for the 16S processing, using the UNITE database for chimera removal and taxonomic identification of representative OTUs. Subsequent analyses of taxon abundances were conducted in R using principally the vegan library (https://github.com/vegandevs/vegan) for rarefaction and ordination. Both bacterial and fungal OTU abundance tables were resampled to a minimum of 4000 reads per sample, and samples with a Shannon diversity higher than 6 and 3.2 for 16S and ITS, respectively, were removed prior to further analyses.

We generated non-metric multidimensional scaling (NMDS) plots to visualise and determine the effect of drought on fungal and bacterial community composition using the functions decostand and adonis in the R library vegan. We also performed principal components analysis (PCOA) based on Bray-Curtis dissimilarities on both communities using the function cmdscale; PCOA scores were used as proxies for community composition in subsequent structural equation modelling. In addition, we used Bray-Curtis similarities between control and drought fungal and bacterial communities as a measure of community resistance and resilience to drought. To identify OTUs significantly associated with the drought treatments at each time point, we used indicator species analysis as implemented within the R library labdsv (http://ecology.msu.montana.edu/labdsv/R). The indval score for each gene is the product of the relative frequency and relative average abundance within each treatment, and significance was calculated through random reassignment of groups (1000 permutations). In the circle plots and subsequent network analyses we focus on only the abundant indicator OTUs which are significant (P < 0.05) and present at >1% abundance. Full tables of indicator scores are provided in supplementary material. For calculating indicator taxa and generating circle plots all taxa with fewer than ten reads across all samples were also excluded to ease presentation.

We analysed fungal and bacterial networks for each sampling and drought and control treatment separately. Thus, each network was based on 36 communities, but only OTUs that occurred in at least 8 communities were included in the analysis, as in Shi et al. [START_REF] Shi | The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages[END_REF] . All network analyses were done in R, using the package igraph, as in Williams et al. [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] ; we adapted the code available at https://github.com/ ryanjw/co-occurrence. Interactions consisted of Spearman's rank correlations and co-occurrence networks were constructed using only significant correlations (P < 0.01 as in Barberan et al. [START_REF] Barberan | Using network analysis to explore co-occurrence patterns in soil microbial communities[END_REF] ) of ρ > 0.6; this cutoff was chosen to include a range of interactions strengths (not only strong interactions). Random networks were constructed as in Williams et al. [START_REF] Williams | Demonstrating microbial cooccurrence pattern analyses within and between ecosystems[END_REF] . Subsequently, we detected network modules using the edge.betweenness.community function and calculated clustering coefficients using the transitivity function, analysed the centrality of network nodes using the betweenness function, and analysed the connectedness of network nodes using the degree function. All bacterial and fungal networks were significantly more clustered than random networks. Networks were visualised in the R library igraph.

General statistical analyses. All data were checked for normality and logtransformed if necessary. Gas fluxes, N cycling gene data and measures derived from sequencing data were analysed using repeated measures ANOVA taking into account the split-plot design (error term with block, drought and mesocosm), soil nutrient data were analysed using ANOVA (error term with block and drought). All analyses were performed in R version 3.3.2 [START_REF]R: A language and environment for statistical computing[END_REF] .

Structural equation modelling. We constructed our a priori models based on current knowledge on plant-microbe-functioning interactions and tested whether the data fit these models using the standard and the multigroup modelling approach in the R library lavaan. We used model modification indices and stepwise removal of non-significant relationships as in De Vries and Bardgett [START_REF] De Vries | Plant community controls on short-term ecosystem nitrogen retention[END_REF] . We used a minimum set of parameters to assess model fit, including root mean square error of approximation (RMSEA), and comparative fit index (CFI). We used PCOA axis 1 scores as a proxy for fungal and bacterial community composition, PCA axis 1 scores as a proxy for plant community composition, and the ratio between the sum of nirS and nirK and the sum of nosZI and nosZII gene abundances to indicate relative changes in genetic potential for the soil functions denitrification and N 2 O reduction.

  Fig.2Relative abundance and proportion of total bacterial and fungal OTUs included in co-occurrence networks. A significantly larger proportion of fungal OTUs than bacterial OTUs increased under drought at the end of the drought and early recovery (a χ 2 = 122, df = 1, P < 0.001 and χ 2 = 23.2, df = 1, P < 0.001, respectively). A significantly larger proportion of bacterial OTUs was included in bacterial co-occurrence networks than of fungal OTUs in fungal co-occurrence networks (b χ 2 -tests for all drought (D) and control (C) treatments over time were significant). Asterisks indicate the significance of the difference between proportions in bacterial and fungal communities (*P < 0.05, **P < 0.01, ***P < 0.001)

  Fig.1Bacterial and fungal community response to drought. The duration of the drought is indicated as a shaded area in b-d and f-h. NMDS of bacterial (a) and fungal (e) community composition shows that drought significantly affected bacterial and fungal communities at all sampling dates following drought (ADONIS). Bacterial (b, c) and fungal (f, g) richness and evenness were strongly affected by drought, but this response differed over time (Repeated measures ANOVA sampling × Drought interaction F 3,153 = 10.4, P < 0.001 and F 3,195 = 14.8, P < 0.001 for richness of bacteria and fungi, respectively; Sampling × Drought interaction F 3,153 = 23.1, P < 0.001 and F 3,195 = 23.1, P < 0.001 for evenness of bacteria and fungi, respectively). The similarity between drought and control communities changed over time for both bacteria and fungi, but bacterial drought and control communities were still less similar than before drought at late recovery (ANOVA sampling effect F 3,80 = 40.9, P < 0.001 and F 3,103 = 2.8, P = 0.04 for bacteria and fungi, respectively). In b-d and f-h, dots represent individual observations and lines indicate means ± 1 SE. In d, h, symbols with the same letter are not statistically different

Fig. 3 7 4

 37 Fig. 3 Frequency distributions of correlations in fungal and bacterial networks. Frequency distributions of the sign and strength of all correlations in fungal and bacterial networks (a), and of the absolute strength of all correlations in fungal and bacterial networks (b), in drought and control treatments over time (top to bottom). Correlations in bacterial networks are light grey, correlations in fungal networks are dark grey

Fig. 6 Fig. 7 7 8

 677 Fig.6Drought effects on plant community composition and the relationship between changes in plant community composition and microbial community composition. Biomass of the four species in our established plant communities in the pre-drought season (a), just before drought (b) and 2 months after ending the drought (c, d), split for control (c) and drought communities (d). Plant community treatments are on the x-axis, with three different evenness levels and each species dominating within each evenness level. Total community biomass increased marginally in response to drought (ANOVA F 1,48 = 273, P = 0.08); D. glomerata (Dg) biomass increased strongly under drought, except in Anthoxanthum odoratum (Ao)-dominated communities (ANOVA dominant species × Drought interaction F 14,48 = 4.52, P = 0.004). Ao Anthoxanthum odoratum, Dg Dactylis glomerata, Lh Leontodon hispidus, Ra Rumex acetosa. PCA-biplot of plant community composition (e) shows that PC-axes 1 and 2 scores were significantly affected by drought at the late recovery sampling (ANOVA F 1,48 = 273, P < 0.001 and F 1,48 = 37.2, P = 0.008, respectively). Resilience of plant community composition (similarity between drought and control, measured as Bray-Curtis distances) explained by the relative change in D. glomerata biomass in response to drought (f). With a larger drought-induced increase in D. glomerata biomass, droughted plant communities were less similar to control communities (P < 0.001). The resilience of bacterial community composition (g) was positively explained by the resilience of plant community composition (R 2 = 0.26, P = 0.008), but this relationship was not significant for fungal community resilience (h R 2 = 0.06, P = 0.182). In a-d, bars represent means ± 1 SE (n = 4, SE for total biomass only); in f-h, dots represent single observations, with shaded areas indicating 95% confidence intervals

Fig. 8

 8 Fig. 8 Structural equation model of relationships between D. glomerata biomass, plant community attributes and microbial properties, at the final, late recovery sampling. For each significant relationship, standardised coefficient and P values are given alongside arrow, and arrow weights are proportional to standardised coefficients. The model fit the data well: P = 0.345, CFI = 0.995, P RMSEA < = 0.05 = 0.439

Fig. 9

 9 Fig. 9 Structural equation model of relationships between plant community attributes, microbial properties and ecosystem functioning for each sampling. For each significant relationship, standardised coefficient and P = values are given alongside arrow, and arrow weights are proportional to standardised coefficients. Our multigroup model fit the data well: P = 0.309, CFI = 0.994, P RMSEA < = 0.05 = 0.574

  ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-05516-7 10 NATURE COMMUNICATIONS | (2018) 9:3033 | DOI: 10.1038/s41467-018-05516-7 | www.nature.com/naturecommunications

2

  Relative abundance and proportion of total bacterial and fungal OTUs included in co-occurrence networks. A significantly larger proportion of fungal OTUs than bacterial OTUs increased under drought at the end of the drought and early recovery (a χ 2 = 122, df = 1, P < 0.001 and χ 2 = 23.2, df = 1,

	a		b
	Before drought	End of drought
	Early recovery	Drought R 2 = 0.09, P < 0.001 Late recovery	Bacterial community richness			Bacterial community evenness				Bray-Curtis similarity between
	Drought R 2 = 0.12, P < 0.001 Drought R 2 = 0.09, P < 0.001
	Control	Drought	Before drought	End of drought	Early recovery	Late recovery	Before drought	End of drought	Early recovery	Late recovery	Before drought	End of drought	Early recovery	Late recovery
	e		f
	Before drought	End of drought
	Early recovery	Drought R 2 = 0.20, P < 0.001 Dominant R 2 = 0.06, P = 0.046 Late recovery	Fungal community richness			Fungal community evenness				Bray-Curtis similarity between	drought and control communities
	Drought R 2 = 0.06, P < 0.001	Drought R 2 = 0.05, P < 0.001
	Dominant R 2 = 0.07, P = 0.011	Dominant R 2 = 0.08, P = 0.002
			Before	End of	Early	Late	Before	End of	Early	Late	Before	End of	Early	Late
			drought	drought	recovery	recovery	drought	drought	recovery	recovery	drought	drought	recovery	recovery
	P < 0.001, respectively). A significantly larger proportion of bacterial OTUs was included in bacterial co-occurrence networks than of fungal OTUs in fungal
	co-occurrence networks (b χ 2 -tests for all drought (D) and control (C) treatments over time were significant). Asterisks indicate the significance of the
	difference between proportions in bacterial and fungal communities (*P < 0.05, **P < 0.01, ***P < 0.001)

Acknowledgements

This study was funded as part of the European project EcoFINDERS (FP7-264465). F.T. d.V. is supported by a BBSRC David Phillips Fellowship (BB/L02456X/1). We thank Caley Brown, Melanie Hartley, Neil Mullinger, Helen Quirk, Deborah Ashworth and Victor van Velzen for help in the lab and field.

Code availability. Code for network analysis was adapted from, and is available at, https://github.com/ryanjw/co-occurrence. Code for indicator analysis is available on request.

Data availability. OTU tables and metadata are available from figshare (https:// doi.org/10.6084/m9.figshare.6548999.v1). Raw sequence data have been deposited at the EBI European Nucleotide Archive under study accessions ERP109485 (fungal amplicons) and ERP109472 (16S amplicons).

Author contributions

Additional information

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-018-05516-7.

Competing interests:

The authors declare no competing interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.