Long-Term Drought and Warming Alter Soil Bacterial and Fungal Communities in an Upland Heathland

The response of soil microbial communities to a changing climate will impact global biogeochemical cycles, potentially leading to positive and negative feedbacks. However, our understanding of how soil microbial communities respond to climate change and the implications of these changes for future soil function is limited. Here, we assess the response of soil bacterial and fungal communities to long-term experimental climate change in a heathland organo-mineral soil. We analysed microbial communities using Illumina sequencing of the 16S rRNA gene and ITS2 region at two depths, from plots undergoing 4 and 18 years of in situ summer drought or warming. We also assessed the colonisation of Calluna vulgaris roots by ericoid and dark septate endophytic (DSE) fungi using microscopy after 16 years of climate treatment. We found significant changes in both the bacterial and fungal communities in response to drought and warming, likely mediated by changes in soil pH and electrical conductivity. Changes in the microbial communities were more pronounced after a longer period of climate manipulation. Additionally, the subsoil communities of the long-term warmed plots became similar to the topsoil. Ericoid mycorrhizal colonisation decreased with depth while DSEs increased; however, these trends with depth were removed by warming. We largely ascribe the observed changes in microbial communities to shifts in plant cover and subsequent feedback on soil physicochemical properties, especially pH. Our results demonstrate the importance of considering changes in soil microbial responses to climate change across different soil depths and after extended periods of time.


ABSTRACT
The response of soil microbial communities to a changing climate will impact global biogeochemical cycles, potentially leading to positive and negative feedbacks. However, our understanding of how soil microbial communities respond to climate change and the implications of these changes for future soil function is limited. Here, we assess the response of soil bacterial and fungal communities to long-term experimental climate change in a heathland organo-mineral soil. We analysed microbial communities using Illumina sequencing of the 16S rRNA gene and ITS2 region at two depths, from plots undergoing 4 and 18 years of in situ summer drought or warming. We also assessed the coloni-sation of Calluna vulgaris roots by ericoid and dark septate endophytic (DSE) fungi using microscopy after 16 years of climate treatment. We found significant changes in both the bacterial and fungal communities in response to drought and warming, likely mediated by changes in soil pH and electrical conductivity. Changes in the microbial communities were more pronounced after a longer period of climate manipulation. Additionally, the subsoil communities of the long-term warmed plots became similar to the topsoil. Ericoid mycorrhizal colonisation decreased with depth while DSEs increased; however, these trends with depth were removed by warming. We largely ascribe the observed changes in microbial communities to shifts in plant cover and subsequent feedback on soil physicochemical properties, especially pH. Our results demonstrate the importance of considering changes in soil microbial responses to climate change across different soil depths and after extended periods of time.

INTRODUCTION
Climate change has and will continue to fundamentally alter global ecosystem functioning. Understanding how ecosystems will respond to future change is essential for societal adaptation as well as mitigation. Soils are a source of major uncertainty in future earth predictions, as we still do not know in sufficient detail how soil biogeochemical cycling, hydrology or biology will change in response to climate change and how these changes will alter climatic feedbacks (Bradford and others 2016). The potential for soils to accentuate or mitigate future climate change is in large part due to the vast quantities of carbon that is currently stored in the soil system ( 2400 Pg, three times as much as is in the atmosphere) (Batjes 1996; Stocker and others 2013; Bossio and others 2020). In order to predict how soil systems will change in the future, we have to take into account how climate alters both the soil physicochemical environment and changes in the soil biota (for example, Robinson and others 2019).
Soil microbial communities have responded to climate change manipulations across various experimental systems (Cavicchioli and others 2019;Jansson and Hofmockel 2020). The response of microbial communities to the different climate change-related stressors is dependent on the type of climatic stress, the ecosystem type and the identity of the microbial communities. There have been recent suggestions that microbial responses to drought are phylogenetically conserved (Amend and others 2016), though this has not been evidenced by another analysis of multiple independent global studies (Oliverio and others 2016). In general, warming has a stronger impact on fungi than bacteria (García-Palacios and others 2015), and stronger impacts on microbial abundances in colder regions (Chen and others 2015). Metaanalyses of the impact of altering precipitation on microbial abundance and composition found that the impact on biomass depended on the climate, with high precipitation areas being more respon-sive to drought and low precipitation areas being more responsive to increased precipitation others 2017, 2018;Zhou and others 2018). Drought has been found to result in increased fungal dominance in a temperate heath (Haugwitz and others 2014), and in grassland ecosystems fungi showed higher resilience to drought than bacteria (de Vries and others 2018).
The response of soil microbial communities to climate change has implications for various essential ecosystem functions, including the provision of nutrients to plants through mycorrhizal associations. The response of mycorrhizal associations to experimental drought or warming varies according to the specific type of mycorrhiza (Olsrud and others 2009;Binet and others 2017). There is some evidence that changes in mycorrhizal associations in response to climate change may be more related to changes in plant composition than to changes in mycorrhizal interactions (Rudgers and others 2014). Studies across a variety of climates and ecosystem types have found that altering precipitation can impact the extracellular enzyme activities within soils (Ren and others 2017), resulting in impacts on soil and root respiration and associated soil carbon loss (Crowther and others 2016;Ren and others 2018). The impact of drought upon soil respiration has been found to be dependent on the local climate, with high precipitation areas being more responsive (Reinsch and others 2017). Drought has been found to affect the microbial community impact on litter decomposition across various studies (Allison and others 2013; Martiny and others 2017; Santonja and others 2017; Tó th and others 2017). The legacy of global change persists within the microbial community for several years and impacts their ability to carry out key functions in response to new or altered environments (Martiny and others 2017).
The increasing awareness of the impact of legacy effects upon the ability of an ecosystem to respond to future change makes the use of long-term ecological experiments increasingly important. Longterm climate change has been found to impact plant communities (Fridley and others 2011; Andresen and others 2016), soil respiration (Crowther and others 2016; Domínguez and others 2017b), hydrological behaviour (Robinson and others 2016), soil mesofauna (Petersen 2011;Holmstrup and others 2013), soil microbial communities (Rousk and others 2013;Sayer and others 2017) and whole ecosystem feedbacks (Abbasi and others 2020). Importantly, many of these impacts emerge only after years of experimental treatment (for example, Andresen and others 2016), indicating how essential long-term experiments are for evaluating future climate change. Within this study, we will look at the climate manipulation experiment in the Clocaenog Forest, NE Wales, UK, which has imposed summer drought and warming treatments over an organo-mineral heathland soil since 1999. There has been a rise in soil respiration in the treatments compared to the controls (Reinsch and others 2017), which has not been matched by changes in the plant community which has shown only a slight shift to greater moss cover in the warming plots (Krö el-Dulay and others 2015). Previous work on this site has found that drought could be impacting the summer fungal community only (Toberman and others 2008), and other studies show limited impact of the warming or drought treatments upon microbial biomass, extracellular enzyme activity, microbial community as measured by phospholipid fatty acid analysis and microbial growth rates (Rousk and others 2013;Domínguez and others 2017a). Within this study, our aims were: to characterise the bacterial and fungal communities after four versus eighteen years of persistent summer drought and warming; to see if the microbial communities are altered by the legacy of repeated drought and warming; to assess whether the microbial response can be clearly demarcated into a bacterial versus a fungal response through identifying co-occurrence patterns; to compare changes in the bulk community composition to changes in mycorrhizal associations; and to see if these changes were associated with changes in the soil physicochemical environment. Our hypotheses were: 1. The bacterial and fungal communities both show significantly different composition within the plots exposed to long-term summer drought and warming compared to control. 2. The impact of climate treatments on microbial community would be greater after eighteen versus four years of treatment, and would also be greater in the topsoil compared to the subsoil. 3. Fungal communities would show stronger responses to the treatments than bacterial communities, with limited interactions across the Kingdoms identified through network analysis. 4. The root fungal associations are affected by treatment in a manner consistent with that of the response of fungal community composition in the bulk soil. 5. The above hypothesised changes in microbial community composition with treatment and soil depth would be partially mediated by changes in soil pH and water content.

Site Description and Sampling
Long-term climate manipulations were carried out in North Wales (53°03¢ N 3°28¢ W) on a peaty podzol. The mean annual precipitation is 1263 mm, and mean annual temperature 7.4°C. The vegetation is dominated by Calluna vulgaris, and the soil consists of a carbon-rich topsoil layer ( 8 cm deep, 87% organic matter) and an underlying layer of gley subsoil ( 4 cm deep, 37% organic matter). Climate manipulations started in 1999, and consisted of a summer drought treatment, a warming treatment as well as unmanipulated control plots. There are three replicate plots per treatment. The treatments are imposed using an automated retractable roof system, with drought plots having the roofs cover the plots during summer rainfall events ( 54% of summer rainfall is excluded) and the warming plots having the roofs cover the plots overnight to keep in the heat (0.2°C increase in mean annual temperature). A full description of the experimental set-up is provided in (Beier and others 2004).
Soil samples were collected in February 2003 and 2017 using a stainless steel auger. Samples were collected from both the organic topsoil and gleyed mineral subsoil for both years, and in 2017, the 2 cm interface between the two soil layers was analysed separately. The samples for each plot were bulked together for topsoil and subsoil in 2003, while three soil cores per plot from 2017 were analysed with the three depths separately. In total, there were 15 samples from 2003 and 78 from 2017, as some samples were not able to be included. Soil pH and EC were measured on frozen 2017 soil samples using a Corning 220 pH meter (VWR combination electrode for pH and EC 662-1805; Jenway 4510). Soil pH was measured from a 1:2.5 (w/v) soil-to-0.01 M CaCl 2 suspension after equilibration for 0.5 h.

DNA Extraction and Sequencing
DNA was extracted from 0.2 g frozen field moist soil using a PowersoilÒ DNA Isolation Kit (Mo Bio Laboratories Inc., Carlsbad, CA) according to the manufacturer's instructions. Amplicons were generated using a 2-step amplification approach, using Illumina Nextera tagged primers. Bacterial 16S V4 primers 515f GTGYCAGCMGCCGCGGTAA and 806r GGACTACNVGGGTWTCTAAT (Walters and Climate change alters microbial communities others 2016), and Fungal ITS primers GTGART-CATCGAATCTTTG and TCCTCCGCTTATTGA-TATGC (Ihrmark and others 2012) were each modified at 5' end with the addition of Illumina pre-adapter and Nextera sequencing primer sequences. Amplicons were generated using a highfidelity DNA polymerase (Q5 Taq, New England Biolabs). After an initial denaturation at 95°C for 2 min, PCR conditions were: denaturation at 95°C for 15 s; annealing at temperatures 55°C and 52°C for 16S and ITS reactions, respectively; annealing times were 30 s with extension at 72°C for 30 s; repeated for 25 cycles. A final extension of 10 min at 72°C was included.
PCR products were cleaned using a ZR-96 DNA Clean-up Kit (Zymo Research Inc., Irvone, CA) following manufacturer's instructions. MiSeq adapters and 8nt dual-indexing barcode sequences were added during a second step of PCR amplification. After an initial denaturation 95°C for 2 min, PCR conditions were: denaturation at 95°C for 15 s; annealing at temperatures 55°C; annealing times were 30 s with extension at 72°C for 30 s; repeated for 8 cycles with a final extension of 10 min at 72°C. Amplicon sizes were determined using an Agilent 2200 TapeStation system. Libraries were normalised using SequalPrep Normalization Plate Kit (Thermo Fisher Scientific), quantified using Qubit dsDNA HS kit (Thermo Fisher Scientific) and pooled together at equal concentrations. The pooled library was diluted to achieve 400 pM in a 40-ll volume after denaturation and neutralisation. Denaturation was achieved with 4 ll 2 M NaOH for 5 min followed by neutralisation with 4 ll 2 M HCl. The library was then diluted to its load concentration of 14 pM with HT1 buffer and 5% denatured PhiX control library. A final denaturation was performed by heating to 96°C for 2 min followed by cooling in crushed ice. Sequencing was performed on an Illumina MiSeq using V3 600 cycle reagents. The DNA sequences are available on the European Nucleotide Archive under primary accession reference PRJEB33721.

Molecular Bioinformatics
Illumina demultiplexed sequences for 16S and ITS were processed separately in R using DADA2 (Callahan and others 2016) to quality filter, merge, denoise and assign taxonomies as follows: Amplicons reads were trimmed to 270 and 220 bases, forward and reverse, respectively. Filtering settings were maximum number of Ns (maxN) = 0, maximum number of expected errors (maxEE) = c(3,5) and amplicon primer sequences removed using trimLeft = c(20,20). Sequences were dereplicated and the DADA2 core sequence variant inference algorithm applied. Forward and reverse reads were then merged using mergePairs function to produce amplicon sequence variants (ASVs). Sequence tables were constructed from the resultant ASVs and chimeric sequences were removed using removeBimeraDenovo default settings. ASVs were subject to taxonomic assignment using assignTaxonomy at default settings; training databases were GreenGenes v13.8 (DeSantis and others 2006; McDonald and others 2012) and Unite v7.2 (Kõ ljalg and others 2005) for 16S and ITS, respectively. Fungal taxa were assigned to trophic modes using FunGUILD (Nguyen and others 2016).

Root Fungal Colonisation
Soil cores for root assays were extracted in April 2015. Each plot had a single 8 cm diameter core taken to 8 cm depth. The cores were cut into 1 cm subsections, soaked in tap water and Calluna vulgaris roots removed by hand and washed to remove soil particles. Root length, diameter and number of tips were measured using WinRHIZO version 3.2 on a flatbed scanner. Proportional colonisation of ericoid mycorrhizae (ErM) and dark septate endophytes (DSE) was estimated using the magnified intersection technique (McGonigle and others 1990). Roots were bleached in 10% KOH for 20 h and then stained with a 5% vinegar-ink solution (Vierheilig and others 1998;Arndal and others 2013). Roots were cut to 1-2 cm in length and 2 mm passes made along each root length. At the end of each pass, all cells were examined for ErM colonisation or the presence of DSE hyphae. ErM colonisation was categorised as 0%, 0-1%, 1-10%, 10-50%, 50-90% and 90-100% colonisation based upon the classification system of Trouvelet and others (1986). DSE proportional colonisation was calculated as the number of colonised intervals divided by the total number of intervals. Root biomass, length and fungal colonisation data are available online at the NERC Environmental Data Centre (White and others 2019).

Statistics
All statistics were performed in R version 3.6.0 (R Core Team 2020). Bacterial taxa were rarefied to 25,000 reads 100 times using the vegan R package (Oksanen and others 2020) and the rounded average used for calculation of richness and diversity indices. The same procedure was used for fungal taxa with rarefaction to 10,000 reads. The cut-offs for rarefaction were identified based on evaluation of the read depths of the samples and removal of the samples with considerably lower read depths than the rest of the data (Salter and others 2014). In order to account for the structure of the data, we modelled the effect of climate treatment, depth and year of collection upon diversity with a Bayesian hierarchical model using the brms package (Bü rkner 2017). These models had plot identity nested in year as a grouping factor and normal priors on the population-level effects with mean 0 and standard deviation 5 for fungi and 10 for bacteria. Priors on the group-level effect were student T distributed with shape parameter 3, mean 0 and the above standard deviation. Model diagnostics-that is, R-hat, effective sample size and graphical posterior predictive checks-were assessed for all models and found to be acceptable. This same procedure was used for the analysis of changes in soil pH and EC with depth and treatment in the 2017 samples, with priors having a standard deviation of 1 for soil pH and 10 for soil EC. Two-dimensional NMDS ordination on the Bray-Curtis distances between taxa based on the rarefied data was used to characterise community composition. The response of community composition to treatment was analysed using PERMA-NOVA, the homogeneity of dispersion assumption was shown to hold at p > 0.1, and plot identity was included as a blocking factor. All figures were plotted using ggplot2 (Wickham 2009).
Co-occurrence networks were constructed using the SpiecEasi package in R (Kurtz and others 2015). Bacterial and fungal taxa that appeared in over 50% of the 2017 samples were used simultaneously to construct networks that contained intraand inter-kingdom co-occurrence relationships for each of the control, warming and drought treatments (Tipton and others 2018). The network was then plotted and its characteristics determined using the igraph package (Csardi and Nepusz 2006). Node degree and betweenness were used to identify the key taxa within the network.
Indicator taxa for the warming and drought treatments were identified through analysing the differential relative abundance of taxa by treatment using the DESeq2 method (Love and others 2014) and the apeglm shrinkage estimator (Zhu and others 2019). Year, climate treatment and soil depth were included as additive terms within the DESeq2 model. The presence of any interaction terms was tested using the likelihood-ratio test, and the absence of any taxa responding to interaction term compared to the reduced model taken as confirmation that no interaction occurred.
The impact of depth and treatment on DSE colonisation was modelled using a Bayesian regression model with a zero-inflated beta distribution using the brms package, and with plot identity as a random effect (Bü rkner 2017). The impact of depth and treatment upon ErM colonisation was modelled using the brms package as well, as a cumulative ordinal regression model assuming the latent variable to be normally distributed (Bü rkner and Vuorre 2019). Plot identity was modelled as a group-level effect, and an interaction between depth and treatment was assumed. The emmeans package was used to obtain contrasts between the treatments given soil depth (Lenth 2021).
The relative impact of soil pH, EC, soil moisture and temperature on the soil microbial community composition was established using multivariate hierarchical Bayesian regression models within brms. Soil moisture and temperature were taken from in situ sensors on the day of sampling. Soil temperature is measured at 5 cm depth with a T107 sensor from Campbell scientific. Soil moisture is measured as volumetric water content with CS616 sensors from Campbell scientific. The numeric predictors (pH, temperature, moisture) were first transformed by centring the mean at zero and dividing by twice the standard deviation so they were on a similar scale to any binary predictor. All models had plot identity included as a grouping factor. Models were compared using leave-one-out cross-validation to estimate pointwise out-of-sample prediction accuracy (Vehtari and others 2017). For the models predicting NMDS scores, the models were fit with both NMDS scores as response variables simultaneously so that the residual correlation between the two could be estimated. Data from 2017 only were used within these models due to the absence of pH and EC values from 2003.

Microbial Diversity
There were 8818 unique sequences from a total of 4,514,220 reads returned by the 16S primer, of which 8673 were matched to bacteria and 138 to archaea. There were 2539 unique ITS sequences from a total of 5,711,663 reads, of which 2377 matched to fungi. The bacterial data were rarefied to 25,000 reads, and the fungal to 10,000 reads. There were 28 samples of the 93 that failed to amplify enough fungal DNA for inclusion in the analysis, of which 12 samples had 0 reads and 14 samples had less than 200 reads. Of the failed Climate change alters microbial communities samples, 15 were from warming plots, compared to 5 and 8 from control and drought plots, respectively (Supplementary Table 1 While both bacterial and fungal richness showed trends with soil depth, sampling year and climate treatment, our models showed no clear effects of treatment, depth or year once plot identity was incorporated (Figure 1). Fungal richness tended to be lower in the subsoil samples regardless of year and treatment (50% quantile of subsoil predictions lower than all middle zone and topsoil predictions); however, there were no particular trends with treatment or year and all the estimates for every depth/year/treatment combination had their 95% quantile overlap. Bacterial richness within all depth/year/treatment combinations was similar, with the 50% quantiles overlapping in all cases. Interestingly, the group-level plot identity factor was estimated to be much greater in the bacterial richness model compared to the fungal richness model (45 compared to 4.6). The response of bacterial and fungal Shannon and inverse Simpson diversity indices to depth, year and treatment followed the same patterns as richness (Supplementary Figures 1 and 2). Overall, there was no clear response of microbial diversity to climate treatments, soil depth or duration of treatment.

Microbial Community Composition
The microbial community composition was impacted by depth and treatment in both sampling periods, with a clear separation by depth and some separation by treatment (Figure 2). The warming subsoil shows greater similarity to the topsoil than to the subsoil of the control and drought plots. There was a significant effect of both soil depth and climate treatment upon both bacterial and fungal composition in 2017, but no significant interaction  Figures 3 and 4). The majority of fungi were not able to be assigned to any trophic mode, likely related to 63% of our fungal taxa not being assigned to a genus, and there was limited evidence of change in trophic modes to treatment (Supplementary Figure 5).

Network Analysis
The nature of the co-occurrence patterns across the microbial data is dependent on whether the bacteria and fungi are considered together or separately. There were 161 bacterial taxa and 54 fungal Figure 1. Change in bacterial richness per 25,000 rarefied reads (a) and fungal richness per 10,000 rarefied reads (b) with soil depth (topsoil, subsoil), year (2003, 2017) and treatment (control, drought and warming). The Bayesian hierarchical model with plot identity as a group-level effect found that the 50% quantiles for every treatment/year/depth combination overlapped for bacteria (a), and although the subsoil fungal diversity tended to be lower than the topsoil diversity the 95% quantiles for every treatment/year/depth combination overlapped for fungi (b). taxa that appeared in over 50% of the samples, and were thus included in the networks (Figure 3,  Supplementary Figure 6). The network construction was run for bacteria and fungi together, bacteria only and fungi only. The majority of the abundant microbial taxa included within our network showed a distinct depth preference, but no effect of treatment (Supplementary Figure 7). There were many links within the joint network that would not have been found if bacteria and fungi were only considered separately (123 out of a total 428 links, with another 283 within bacteria and 22 within fungi). Inter-kingdom links include four links from a fungal ASV identified as Mortierella humilis to bacterial taxa and two from a fungal ASV identified as Hyaloscypha fuckelii to bacterial taxa, both fungi being probable saprotrophs and none of the bacterial taxa being identified to species level. The bacteria only network had 312 edges, and the fungi only network had 10 edges. The interactions between the different microbial communities in our sites act across the kingdom boundaries and ignoring across-kingdom interactions changes specific links and overall network stability (Supplementary Figure 6).

Indicator Taxa
Analysis of the different relevant abundance of bacterial and fungal taxa in drought and warming compared to control treatments revealed a small number of taxa that responded strongly to the treatments, with no taxa responding differently to the treatment in the different depths or years (Ta-ble 1, Supplementary Figure 8). There were more bacterial taxa that declined under the climate change treatments than increased, and over half of the taxa that declined did so under both drought and warming. However, in the fungal communities while the warming treatments caused a decline in more taxa than increased, there were far more fungal taxa that increased under the drought treatment than decreased.  The number of taxa (percentage of all taxa) that responded to warming, drought or both treatments compared to control at a p < 0.1 significance level.

Root Fungal Colonisation
The climate treatments impacted overall root biomass and the change in mycorrhizal colonisation of roots with depth as well as the overall microbial community composition. There was a decrease in root biomass and number of root tips with depth ( Supplementary Figure 9). Drought and warming had only a limited effect upon the number of root tips and fine root biomass once depth was accounted for (elpd of all models within standard errors of each other), but did decrease the overall root biomass by about 50% and about 33%, respectively (elpd difference 3.7 ± 2.2). The rate of ericoid mycorrhizal colonisation also tended to decrease with depth, with a decline in the number of 50-90% colonised sections and increase in the number of 1-10% colonised sections in both control and drought plots (95% quantile of depth effect -0.18 to 0.10 in control and -0.29 to 0.01 in drought). However, proportionally more roots were colonised by dark septate endophytes at the lower depths. These changes in colonisation with depth were not apparent in the warming plots (95% quantile included 0 for depth effect in warming plots in both models), and overall warming had lower ErM colonisation rates than control and drought (95% HPD -1.39 to -0.44 on probit scale, Figure 4), and higher levels of DSE colonisation at intermediate depths (95% HPD 0.18 to 0.68 on log-odd scale, Supplementary Fig-ure 10). Overall, compared to the control plots the warming plots had lower rates of fungal colonisation, while drought plots had higher rates of ErM colonisation in the topsoil.

Impact of Soil Chemical Properties on Microbial Communities
The change in pH and EC with depth was altered by the experimental treatment ( Figure 5). Under control conditions, pH increases with depth (subsoil was 0.12 units higher than topsoil, 95% interval: 0.03 to 0.22), but within the drought and warming plots there was no change in pH with depth (95% interval was -0.10 to 0.10 for drought and -0.16 to 0.03 for warming). EC was lowest in the subsoil in the control and drought plots (control subsoil was 21.2 units below topsoil, 95% interval: 8.6 to 33.0, drought was 31.6, 95% interval 14.1 to 48.1); however, in the warming plots there was less change in EC with depth (95% interval 1.4 to 34.9). However, these differences in pH and EC with depth and treatment were limited, with the addition of treatment to the model performing equivalently to the model with depth only (elpd difference of adding treatment for pH + 1.0 ± 3.0, for EC -0.7 ± 1.1). The subsoil of the warming plots had similar conditions to the topsoil throughout the site. The models suggest that the impact of treatment upon fungal richness was fully mediated by the changes in soil physicochemical environment, while bacterial richness was less impacted by the change in soil properties. No models did particularly well at predicting the bacterial richness, but the best model for investigating bacterial richness had solely climate treatment and soil depth as predictors (Bayes R 2 = 0.227 ± 0.071), with three other equivalent models by elpd (taking standard errors into account) adding subsets of pH, EC, moisture and temperature to the predictors. Bacterial richness was lowest in the intermediate depth and highest in the drought plots. Fungal richness was impacted by the physical properties of the site, with the best model including pH, EC, soil depth, temperature and moisture as predictors (Bayes R 2 = 0.651). The model with only climate treatment and soil depth as predictors was notably worse than all the other models (elpd difference of 3.4 compared to < 1 for all other models). Fungal richness decreased with depth and increased with EC and pH. There was also evidence for decreasing richness with moisture and increasing with temperature.
The microbial community composition was still impacted by the treatment after accounting for changes in soil temperature, moisture, pH and EC. The best model had treatment, pH and EC as predictors of bacterial and fungal composition (Figure 6). The relative impact of pH on bacterial composition was higher than that of EC, whereas the reverse was true for fungal composition. The treatments have resulted in changes in the soil physicochemical structure which have impacted the microbial community composition, but the impact of treatment is only partially mediated by the measured changes in soil chemistry.

Microbial Community Response to Warming and Drought
Our results show that fungal taxa are more responsive to drought and warming than bacterial taxa, with warming altering soil biological and chemical properties throughout the soil profile. We found that these microbial community responses were dependent upon the duration of treatment, in agreement with our hypothesis. Although we found limited evidence of treatment effects after and climate treatments in 2017. pH was higher in the subsoil compared to topsoil in control plots (0.12, 95% CI 0.03-0.22), but drought showed no difference (0.00, 95% CI -0.10-0.10) and pH in warming plots tended to decrease with depth (0.07, 95% CI -0.16-0.03). EC increased with depth in all three treatments (95% CI: control -8 to -33; drought -14 to -48; warming -1 to -35). Figure 6. The parameter estimates for the impact of soil depth, EC, pH and climate treatment on bacterial (top row) and fungal (bottom row) NMDS scores. The circles are the mean estimate, the thick bars the standard error and the thin bars the 90% CI. Data presented are from 2017 only. The impact of depth is represented as the impact of the subsoil (DepthsS) and topsoil (DepthT) relative to the transition zone, and the impact of treatment represented relative to control. 4 years, previous work on this site after a similar period of time found that drought treatment effects upon fungal communities were limited to the summer treatment period, and our sampling occurred outside this period in the winter season (Toberman and others 2008). Our results indicate that it takes longer than 4 years to develop legacy effects of climate change treatments, and based on comparison to previous studies of microbial growth and biomass at our site this legacy effect may not have been apparent after 13 years of treatment (Rousk and others 2013). Our observed impacts of climate change are at least partially moderated by changes in the soil chemical environment which are likely driven by alterations in soil hydrology and the plant community at the site in response to treatment (Domínguez and others 2015;Robinson and others 2016). The importance of the soil physicochemical properties in mediating climate change impacts upon the microbial community is in agreement with results from other climate change experiments (Deltedesco and others 2020). We have observed changes in the fungal community using both DNA sequencing of the soil and microscopic examination of the mycorrhizal colonisation of C. vulgaris roots. This indicates that changes in the bulk soil microbial community reflect changes in the functional capability of the soil microbiome to interact with plants which could alter biogeochemical cycling (Read and Perez-Moreno 2003;Tedersoo and Bahram 2019).
Previous studies on the microbial response to long-term climate change have largely focused on measuring the microbial biomass response, which is highly dependent on the soil type and climate conditions (Ren and others 2018;Zhou and others 2018). Drought and warming have been found to impact microbial composition more strongly and more persistently than microbial richness in both short and long-term experiments (Sayer and others 2017; Tó th and others 2017; de Vries and others 2018; Yu and others 2018; Birnbaum and others 2019). This is consistent with our finding that both bacterial and fungal richness show limited response to climate change treatments across soil depths and regardless of duration of treatment. Changes in composition in response to long-term drought have been shown in some cases to be driven by changes in the rarer taxa in a calcareous grassland (Sayer and others 2017). This supports our finding that the dominant taxa show limited response to treatment in contrast to their strong depth preference, as demonstrated by the qualitative description of the taxa in the network analysis by colouring according to either depth or treatment preference.

Soil Depth
A variety of observational studies have shown a distinct difference in microbial composition with depth, consistent with our results (Griffiths and others 2003;Serkebaeva and others 2013;Delgado-Baquerizo and others 2017;Seuradge and others 2017). The subsoil is generally less biologically active and shows less seasonal variation in its physical conditions and microbial communities (Griffiths and others 2003). Therefore, it is no surprise that many studies have found that the topsoil is more strongly linked to the aboveground land use and plant productivity (Delgado-Baquerizo and others 2017; Seuradge and others 2017). However, our results showed greater impacts of climate change deeper into the soil than we expected, impacting soil physicochemical properties, plant-soil interactions and microbial processes. This change in the stratification with depth of physicochemical properties and microbial communities could be linked to the distinct change in hydrological behaviour and roots over the course of the experiment (Robinson and others 2016). We found that the overall root biomass decreased under drought and warming by around a half and a third, respectively, in the top soil layers while fine root biomass did not, which agrees with previous results at this site that showed increased proportion of fine root biomass under drought (Robinson and others 2016). This shift in root allocation away from thicker roots to finer roots could potentially be related to a reduction in the structural strength of the soil or changes in the hydrology such that thicker roots are less important for soil penetration and water uptake while finer roots are maintained for nutrient absorption (Hodge and others 2009). Our soils are admittedly shallow and the changes with depth we have observed are over relatively short distances (Sowerby and others 2008). However, the results still hold relevance for a geographically extensive soil type and are likely to occur in many deeper soil types.
The change in stratification with depth of the soil physical, chemical and biological properties in response to warming may be related to aboveground changes in plant communities, particularly the increased cover of moss in the warming plots (Domínguez and others 2015). Moss acts as a layer of insulation on the soil surface, buffering changes in soil temperature and moisture (Turetsky and others 2012). Mosses have been found to reduce evapotranspiration, increase surface infiltration and influence the partitioning of heat fluxes (Beringer and others 2001; Blok and others 2011). The Climate change alters microbial communities thermal and hydrological influences of moss cover can lead to differences in belowground soil microbial communities, changing their biomass and activity (Gornall and others 2007;Benavent-Gonzá lez and others 2018), with relevance to global biogeochemical cycling (Porada and others 2014). The insulative properties of moss could be reducing the magnitude of daily and seasonal heat and moisture fluctuations at the soil surface, therefore reducing subtle temperature and hydrological differences between the topsoil and subsoil (Turetsky and others 2012). However, our results indicate that the observed changes are occurring in the subsoil, rather than in the topsoil, which could mean that the changes in water infiltration and heat fluxes are leading to increasing penetration of water and translocation of chemicals within the soil profile.

Evaluating the Concept of a Bacterial Versus a Fungal Response
The community composition response to drought and warming we have found shows some differences between bacterial and fungal taxa. However, we found more co-occurrence links between specific bacteria and fungi than we hypothesised, suggesting that caution should be used in interpreting the broad-scale results as a bacterial response versus a fungal response. Although the overall abundance and dominance of bacteria and fungi may change in response to climate change this obscures changes in fine-scale dynamics which could be creating specific drought and warming communities that are composed of a combination of both bacteria and fungi. We do find that more fungal taxa respond positively to drought and negatively to warming which indicates that the drought microbial community could consist of relatively more fungi. This is consistent with previous results showing higher tolerance of fungi to drought (Haugwitz and others 2014). Neglecting the inter-kingdom links could lead to erroneous assumptions about microbial community stability, as suggested by work in the human microbiome (Tipton and others 2018). It is possible that the correlations between bacteria and fungi are driven solely by changes in environmental conditions and do not represent a true microbial interaction, as these co-occurrence networks are prone to this kind of error (Carr and others 2019). However, this result does offer an intriguing avenue for future work in establishing whether specific inter-kingdom interactions such as those we have identified between fungal saprotrophs and bacterial taxa exist and can influence microbial community structure and activity in soils.

Fungal Root Colonisation
Our results from examination of the fungal root colonisation are in agreement with previous results that have indicated that ErM and DSE colonisation respond differently to experimental warming (Olsrud and others 2009; Binet and others 2017). However, the treatment responses we found differed from previous results, with Olsrud and others (2009) finding no effect of warming on ErM colonisation after 6 years of treatment in a subarctic birch forest. Other results from a Danish heathland with similar experimental set-up to our site on a sandy soil found that ErM colonisation was lower in the warming treatment at 5-10 cm depth compared to drought and control, and that DSE was lower at depth in the drought treatments (Arndal and others 2013). Overall, we found lower ErM colonisation in response to warming; however, we found higher ErM colonisation in response to drought which is different from Arndal and others' results. These contrasting results could reflect the difference in duration of treatment, differences in soil type and drainage behaviour or potentially the differences in the plant response to treatment. The implications of shifts in mycorrhizal type for plant community resilience and biogeochemical cycling in response to warming are unclear due to the lack of knowledge on the relative role of DSE versus ErM colonisation upon plant nutrient uptake and stress resilience (Newsham and others 2009). However, there are some suggestions that DSE colonisation may enable the uptake of organic nitrogen compounds and improve resilience to certain stressors (Newsham 2011;Hill and others 2019).

Implications for Ecosystem Function
Our analysis has revealed that microbial community shifts in response to climate change occur over decades, not just years, which could have continued, long-term effects upon ecosystem functions. While we have identified clear shifts in some functions-that is, mycorrhizal colonisation-it is unclear what impact other changes in microbial community composition could have upon ecosystem function and stress resilience. One potential ecosystem function of clear importance to both global carbon cycling and local ecosystem health is soil respiration, and microbial communities are known to be particularly important in determining soil respiration (Davidson and Janssens 2006). Our site has previously been shown to have experienced increases in respiration in response to drought and warming, which could be related to the microbial community composition and fungal root colonisation changes described here (Reinsch and others 2017). It is important to consider that changes in ecosystem function need not be driven by changes in overall microbial biomass or dominant taxa, as rare taxa can be disproportionately important to biogeochemical cycles (Jousset and others 2017). Therefore, the small number of specific taxa we have identified as showing differences between the climate treatments could be directly influencing ecosystem functions, as well as potentially driving shifts in microbial community network connectivity with associated impacts on ecosystem functioning (Wagg and others 2019). The changes we have observed are persisting over the winter season when treatments are not in effect, indicating the possibility of legacy effects of climate change on ecosystem functions mediated by changes in microbial communities.

CONCLUSIONS
We have found that bacterial and fungal community composition and diversity is altered by longterm drought and warming treatments. The impact of simulated climate change is greater after eighteen years of treatment compared to four years. The shift in bulk soil community composition is also reflected by a shift in fungal root colonisation of the dominant Calluna vulgaris plant species. These changes are at least partially driven by changes in the distribution of roots throughout the soil profile and changes in the soil chemical environment. In general, fungal taxa appear to be more responsive to the climate change treatments than bacterial taxa but the presence of cross-domain cooccurrence relationships and sensitive taxa in both domains cautions against interpreting our results as a fungal versus a bacterial response. These insights into soil microbial community response to drought and warming can inform our understanding of what drives differences in soil functional response to climate change.

AC KNOWLEDGEMENT S
We thank all those who have maintained the Clocaenog experimental site over the decades of operation, and Rachel Harvey and Timo Breure for their assistance in soil sampling. We acknowledge the interests of the Ecological Continuity Trust (ECT), whose national network of long-term experiments includes the experiment on which this research was conducted and who funded the microbial analysis. F.M.S. was supported by the award of a studentship grant by the NERC and BBSRC funded STARS (Soils Training and

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/. D A T A AV AI L ABI LI TY Data for this manuscript are publicly available under the European Nucleotide Archive (https://ww w.ebi.ac.uk/ena/browser/view/PRJEB33721) and the Environmental Information Data Centre (http s://doi.org/10.5285/3d468857-f5d0-4dc4-88f3-6be 6df19608b). response to simulated global change is phylogenetically conserved and linked with functional potential. ISME J 10:109-118.