Soil Bacterial and Fungal Richness and Network Exhibit Different Responses to Long-Term Throughfall Reduction in a Warm-Temperate Oak Forest

: Prolonged drought results in serious ecological consequences in forest ecosystems, particularly for soil microbial communities. However, much is unknown about soil microbial communities in their response to long-term consecutive droughts in warm-temperate forests. Here, we conducted a 7-year manipulated throughfall reduction experiment (TFR) to examine the responses of bacterial and fungal communities in terms of richness and networks. Our results show that long-term TFR reduced bacterial, but not fungal, richness, with rare bacterial taxa being more sensitive to TFR than dominant taxa. The bacterial network under the TFR treatment featured a simpler network structure and fewer competitive links compared to the control, implying weakened interactions among bacterial species. Bacterial genes involved in xenobiotic biodegradation and metabolism, and lignin-degrading enzymes were enhanced under TFR treatment, which may be attributed to TFR-induced increases in ﬁne root biomass and turnover. Our results indicate that soil bacterial communities are more responsive than fungi to long-term TFR in a warm-temperate oak forest, leading to potential consequences such as the degradation of recalcitrant organics in soil.


Introduction
Microorganisms play essential roles as decomposers and provide crucial ecosystem functions such as soil carbon sequestration [1,2]. There is a challenge in understanding the responses of these highly diverse and complex microbial communities to extreme droughts, which are projected to continually increase in frequency and intensity in the future [3]. Previous studies that have mostly conducted short-term drought experiments have demonstrated considerable effects on soil microbial communities across different ecosystems [4,5]. However, we still know little about the impacts of long-term consecutive droughts on soil microbial communities and the potential functions they mediate [6], particularly in warm-temperate forests.
Soil microorganisms include bacteria and fungi that have different roles in soil biogeochemical cycles. Due to the different cellular structures and physiological properties, bacterial and fungal communities are very likely differentially susceptible to drought [7]. Generally, bacteria are thought to be more sensitive to drought than fungi [8,9]. A recent study showed that drought had a greater influence on soil bacterial communities than fungal communities, leading to the degradation of recalcitrant carbon [10]. However, other studies have indicated that fungi are more sensitive to soil moisture [11,12]. This inconsistency is mainly due to different precipitation change settings [13] and different ecosystem types [14] or different microbial community measurement methods. Additionally, previous studies often focused on single characteristics (e.g., abundance or diversity) of microbial communities, rather than a large number of network interactions between different microbial taxa. This is an important knowledge gap, given that environmental change can recombine microbial network interactions [15]; furthermore, the characteristics of microbial networks can determine their response to climate change [16].
Soil microbes do not exist alone but often form complex interspecific networks that can affect microbial communities in their response to climate change [17]. Network analysis has proved to be a powerful way to explore the interactions between different taxa of complex microbial communities [18,19]. Despite its pitfalls [20], network analysis provides important new insights into the complex microbial community structure and interactions among taxa, and the responses of these to environmental change, which may be more difficult to understand using the routine diversity metrics used in microbial ecology [21]. For example, de Vries et al. (2018) revealed that short-term drought had more influence on soil bacterial than fungal networks [17]. However, few studies have been conducted to evaluate bacterial and fungal networks in response to prolonged drought under field conditions.
Oak forest is a critically important component of forest ecosystems, with the largest area and vegetation carbon storage in China, that greatly contributes to regulating soil carbon sequestration and regional climate change [22,23]. Here, we carried out a 7-year manipulation experiment designed with throughfall reduction (TFR) in a warm-temperate oak forest to explore the responses of soil microbial communities. Our previous results showed that a 4-year TFR experiment did not have significant impacts on microbial biomass, while little was known about soil microbial communities during the first four-year period [24]. Another study in this natural forest showed no significant change in soil bacterial diversity across four seasons from September of 2015 to June 2016 [25], indicating that significant changes in microbial communities are not necessarily observed in the short term in this area [26]. In this study, we consecutively sampled during the sixth and seventh years of the imposed drought and evaluated the impacts of prolonged drought on the abundance and diversity of bacterial and fungal communities. Recently developed Tax4Fun and network analysis was also performed to offer new insight into the impacts of long-term drought on soil microbial communities. We hypothesized that TFR would negatively affect soil bacterial richness more than fungal richness, as bacteria are generally more susceptible to drought than fungi [27,28]. Specifically, we predicted that bacterial taxa, such as Actinobacteria and Chloroflexi, known to be adapted to drought conditions [29,30], would increase in response to long-term TFR. By contrast, we predicted that bacterial genes related to the degradation of recalcitrant carbon would be enhanced under TFR treatment [31]. We also hypothesized that TFR treatment would weaken bacterial interactions, leading to a simpler network structure, fewer hubs and connectors, and lower competitive links compared to the fungal network [19].

Study Site and Experimental Design
The selected research site was at the Baotianman National Field Research Station of Forest Ecosystem (30 • 20 -33 • 36 N, 111 • 47 -112 • 04 E), Henan Province, in central China. The mean annual precipitation and air temperature are 894 mm and 15.1 • C, respectively (1400 m a.s.l.) [32]. The soil is Haplic Luvisol [33] and has a pH of 4.75 [24]. The typical vegetation is deciduous broadleaf tree species such as Quercus variabilis, Quercus aliena var. acuteserrata and Fagus engleriana, and some coniferous tree species such as Pinus armandii, Pinus tabulaeformis and Pinus massoniana.
The experimental platform was established in an oak secondary forest (60 years old), dominated by Q. aliena var. acuteserrata. In the spring of 2013, six 20 × 20 m plots were built, with three plots designed for manipulated throughfall reduction ("TFR") and three plots as controls ("control"). The TFR treatment was as described in Lu et al. (2017) [24]. Briefly, about 160 plastic sheets (0.5 m × 3 m), covering 50% of a plot area (20 m × 20 m), were suspended 1.5 m above the floor to exclude throughfall during the growing seasons (May-October) from 2013 until 2017. The 50% TFR imposed a~30% reduction in annual precipitation during the experiment and significantly reduced the soil moisture, but did not change the surface soil temperature [24]. In the spring of 2018, we adjusted the magnitude of the TFR from 50% to 70%. To prevent underground and surface water flow, we inserted a 0.7 m-deep plastic barrier around each TFR plot. We gathered the litter that fell on the plastic sheets every week and spread it back evenly within the plot to avoid potential changes in the input of litterfall. We set up buffer zones (~2 m wide) along the inside edge of each plot, and no measurements were taken inside the buffer zone.

Sampling and Measurements
In July 2018 and August 2019, we randomly collected five soil samples from each plot at 0-10 cm depths. After picking out fine roots (<2 mm), the five collected samples from each plot were completely mixed into one composite soil sample. One part of the composite samples was kept at 4 • C for chemical and microbial biomass analyses. One small part of the composite samples was stored at −20 • C for soil microbial analysis. The biomass of fine roots (FRB) was determined after oven drying at a constant weight.
A modified in-growth core method [34] was conducted to determine fine root production. Ten stainless steel cubes (20 × 20 × 20 cm) were created in each plot for two soil layers (0-20 and 20-40 cm) and were refilled with rootless native soil in May 2019. The fine roots in these in-growth blocks were collected at the end of October 2019. The turnover rate of fine roots is defined as the ratio of the fine root production to FRB [35].
Gravimetric soil moisture was calculated as the mass lost after drying a certain soil at 70 • C for 72 h. In addition, an Em50 data logger was installed with five 5TM combined soil temperature (ST) and moisture (SM) probes in each plot. The ST and SM were continuously measured at 30 min intervals at a depth of 0-5 cm. The soil pH was determined with a 1:2.5 ratio of air-dried soil to deionized water. Soil organic C (SOC) and total N (TN) were analyzed with an elemental analyzer. The total P (TP) was determined by alkaline fusion molybdenum-antimony colorimetry [36]. The soil microbial biomass C (MBC) and nitrogen (MBN) were analyzed by the chloroform fumigation extraction method [37]. The MBC and MBN were based on the differences in extractable C content between fumigated and unfumigated soils, with conversion coefficients of 0.45 [38] and 0.54 [37].

DNA Extraction and Illumina Sequencing
Microbial DNA from each soil sample was extracted using the E.Z.N.A. ® soil DNA Kit (Omega Bio-tek, Norcross, GA, USA). The extracted DNA concentration and purity were determined with a NanoDrop 2000 UV-vis spectrophotometer. We used the primer pairs 515F/907R to amplify the bacterial 16S rRNA gene [39]. The fungal ITS gene was amplified with the primers ITS1F and ITS2R [40]. The PCR amplification and Illumina sequencing were conducted by the Majorbio Company (Shanghai, China) using the Illumina MiSeq PE300 platform.

Sequencing Data Processing
The DNA sequencing reads were analyzed using fastp and merged using FLASH according to the following criteria: reads < 50 bp and ambiguous bases were discarded; the allowable overlap region mismatch ratio was 0.2; the allowable barcode mismatch number was 0; the allowable primer mismatch number was 2. Operational taxonomic units (OTUs) were built at 97% similarity using UPARSE (version 7.1). Singletons, and chimeric sequences identified by UCHIME were discarded [41]. The taxonomy of each OTU was analyzed using the RDP Classifier with a confidence level of 70%, and the SILVA database for bacteria or the UNITE dataset for fungi [42]. A total of 602,237 high-quality bacterial sequences (from 36,269 to 72,359) and 685,712 fungal sequences (from 39,916 to 74,370) were obtained. All the obtained data were rarefied to the minimum sequence of the sample before downstream analyses. The community composition and diversity were analyzed on a platform provided by the Majorbio Company (www.i-sanger.com). Tax4Fun was run to predict the functional pathways and critical enzymes related to carbon degradation in bacterial communities using functional annotation based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases [43].

Network Analysis
Bacterial and fungal networks were constructed for the control and TFR treatments, respectively. We used the OTUs that appeared in four out of six samples for network computation. The correlation between pairwise OTUs was measured with Spearman's Rho. A cutoff value for the similarity matrix was automatically generated according to the default settings. All the analyses were conducted using the MENAP (http://ieg2.ou.edu/ MENA/) [44]. The node topology roles were classified into four categories according to within-module connectivity (Zi) and among-module connectivity (Pi) [15]: (a) peripherals, (b) connectors, (c) module hubs and (d) network hubs. Module hubs (Zi > 2.5, Pi ≤ 0.62) and connectors (Zi ≤ 2.5, Pi > 0.62) were proposed as keystone species because of their crucial roles in the network structure [44]. The network topologies were shown by Cytoscape 3.7.1 [45].

Statistical Analyses
All the statistical analyses reported were performed in R and SPSS version 24.0 for Windows (SPSS, Chicago, IL, USA). We carried out one-way ANOVAs to evaluate the TFR effects on microbial attributes in 2018 and 2019. We also carried out two-way ANOVAs to test the effects of the TFR and year on the soil properties, microbial biomass and fine root biomass. The differences in bacterial and fungal community composition were visualized using ordination by nonmetric multidimensional scaling (NMDS) [46], and were tested using permutational multivariate analysis of variance (999 permutations) with "Adonis" in the vegan package [47]. The significance of differences in network characteristics between control and TFR were analyzed with the χ 2 test or Student's t-test. We carried out Pearson correlation to evaluate the relationships between soil microbial richness and abiotic or biotic factors.

Soil and Root Characteristics
During the prolonged drought treatment from 2018 to 2019, the TFR significantly decreased the soil moisture (p < 0.05), but the soil temperature showed no significant changes between the TFR and the control (Table S1, Figure S1). The TFR significantly increased the fine root biomass and turnover (p < 0.05) ( Table S1). The TFR had no significant effects on soil chemical properties, including the pH, SOC, TN and TP (Table S1).

Microbial Diversity and Composition
The TFR significantly decreased the bacterial richness (p < 0.05) during the prolonged drought treatment from 2018 to 2019 (Figure 1). The TFR significantly decreased the bacterial Shannon diversity index (p = 0.044), with only a slight decrease in the bacterial Pielou's evenness (p = 0.068) in 2018 ( Figure S2). For fungi, the TFR had no significant effects on the fungal richness, diversity or Pielou's evenness (Figure 1b, Figure S2). Pearson correlation analysis revealed that bacterial richness positively correlated with soil moisture (p < 0.05) but negatively correlated with fine root biomass (p < 0.05) ( Figure S12). We also found that fungal richness negatively correlated with fine root biomass (p < 0.05) ( Figure S12).

Functional Profiles of Bacteria
Using Tax4Fun as a predictive exploratory tool, the bacterial metabolic pathways were classified into six categories (level 1; Figure S11). The TFR treatment decreased the Genetic Information Processing pathways during the prolonged drought treatment from 2018 to 2019 (p < 0.05). The TFR decreased Cellular Processes in 2018 (p < 0.05) and slightly decreased them in 2019 (p = 0.07). However, the TFR increased the Metabolic Pathways in 2019 (p < 0.05) and slightly increased them in 2018 (p = 0.08). Specifically, we found that Xenobiotics Biodegradation and Metabolism (level 2) related to Metabolic Pathways (level 1) were enriched under long-term TFR treatment from 2018 to 2019 (p < 0.05) (Figure 3a).   The production of exoenzymes for carbon degradation from bacteria was predicted by using Tax4Fun (Figure 3b). The TFR treatment increased the relative abundance of monophenol monooxygenase (1.14.18.1) (p < 0.05) during the prolonged drought experi-ment from 2018 and 2019 and increased limonene-1,2-epoxide hydrolase (3.3.2.8) only in 2019 (p < 0.05).

Bacterial and Fungal Networks
The TFR decreased the proportion of negative correlations in the bacterial network (p < 0.05), but the opposite occurred in the fungal networks (p < 0.05) (Figure 4, Table S2). The TFR decreased the average degree of nodes (the average links per node in the network) in the bacterial networks (p < 0.05), while the opposite occurred in the fungal networks (p < 0.05) (Figure 4, Table S2). We found no network hubs, and most of the nodes were peripherals in the bacterial and fungal networks ( Figure 5). For bacteria, four module hubs and five connectors were observed under the control, and the numbers decreased to two and three under the TFR treatment, respectively ( Figure 5). However, fewer module hubs and connectors were observed in the fungal communities in our study ( Figure 5).

Effects of TFR on Bacterial and Fungal Richness
Molecular analyses offer an insight into the soil microbial communities in their response to long-term TFR. In the study, bacterial richness was reduced after the 6-year prolonged consecutive TFR (Figure 1, Figure S2). Similar results have also been reported in tropical forest and subtropical forest ecosystems, with manipulated experimental drought [10,48]. Many studies have suggested that altered precipitation can influence soil bacteria via other physicochemical properties such as soil pH, SOC, soil nitrogen or soil phosphorus content [4,49,50]. In our study, most of the measured physicochemical properties had weak relationships with bacterial diversity (Figure S12), as they remained relatively stable under TFR treatment (Table S1). Additionally, we found negative correlations between bacterial richness and fine root biomass ( Figure S12). Similarly, in an invasion and drought interaction experiment, Fahey et al. (2020) found that root biomass had negative effects on microbial richness [46]. This relationship suggested that the higher fine root biomass associated with the drought may reduce bacterial richness [51][52][53], while it may have been an accompanying phenomenon with the long-term TFR (Table S1, Figure 1) in our study. Furthermore, soil moisture, the most crucial attribute, significantly correlated with bacterial richness (Figure S12), which demonstrates its crucial role in influencing microbial growth and activities [54,55].
Despite increasing evidence of soil microbes responding to environmental changes, there is still considerable uncertainty in the response of microbial taxa to drought [4,56]. In the study, bacteria were more influenced by TFR than fungi across different taxonomic levels (phylum, class, order and family). Earlier studies suggested that Proteobacteria and Bacteriodetes were sensitive to drought, while Firmicutes and Actinobacteria were resistant to drought conditions [57]. Nevertheless, we found that long-term TFR mainly affected relatively rare bacterial taxa. This is also consistent with Fahey et al. (2020), who found that experimental drought only affected rare bacterial and fungal taxa (i.e., <1% relative abundance) [46]. Our results suggest that rare bacterial taxa were more sensitive to TFR treatment than dominant taxa, probably resulting in the reduction of soil bacterial richness under the prolonged drought treatment in the present study [58]. Given that microbes may show quite different responses to drought based on adaptation to specific environmental conditions [8,59], more research in warm-temperate forests is needed.
Generally, the activity of a soil bacterial community is positively correlated with soil moisture and tends to decrease as microbes die or enter dormancy under drought [60,61]. In the present study, TFR decreased bacterial genes associated with Genetic Information Processing and Cellular Processes ( Figure S11), which could be explained by the decline in bacterial richness under TFR treatment, as reduced species richness allows for fewer metabolic activities [62]. However, we found that TFR increased the relative abundance of Xenobiotic Biodegradation and Metabolism, and exoenzymes for aromatic and lignin degradation (Figure 3), likely because the TFR-induced increase in fine root biomass and turnover (Table S1) stimulated bacteria to utilize this root material as a source of carbon, nitrogen or energy [4,63]. Previous studies also indicated that bacterial genes associated with the degradation of complex plant polysaccharides were enhanced under drought conditions [31,64]. Our results suggest that bacteria possibly switched to degrading more recalcitrant carbon within plant materials under long-term consecutive TFR.

Effects of TFR on the Interactions in Bacterial and Fungal Networks
Documenting the interactions between coexisting microbial taxa may help to ascertain their responses to drought [17]. In the study, we conducted a network analysis to evaluate the interactions between microbial species under long-term TFR treatment. Our results confirmed our prediction that bacterial interactions would be weakened by TFR treatment, generating simpler networks. First, the decreased complexity of the bacterial networks under TFR was supported by the decreased average degree (Figure 4, Table S2) and the longer harmonic geodesic distance [44]. These network changes indicated that relationships between different bacterial species could be weakened by the TFR treatment in this study. Second, fewer negative connections were observed in the TFR treatment than in the control (Figure 4, Table S2), indicating that the extent of competition among species was reduced [19]. Third, there were lower numbers of module hubs and connectors in the TFR compared to the control ( Figure 5). The lower connectivity and negative correlations may suggest that some bacterial species filled specific fragmented niche spaces where no direct competitors existed [65], and may also suggest that many soil bacteria were in an inactive or dormant state under TFR treatment [66]. Fungi, however, developed higher connectivity and contained more negative correlations under TFR (Figure 4, Table S2), and their interactions might have been strengthened by the TFR treatment.
It has been shown that increasing complexity can improve community stability [19,67]. It has also been proved that competition can promote the stability of microbial communi-ties [17,68,69]. Based on the above reasons, the simple network structure (lower average degree and fewer module hubs and connectors) and few competitive connections indicated that soil bacterial networks may be less stable under TFR treatment.

Conclusions
In this study, long-term experimental throughfall reduction (TFR) in a warm-temperate forest caused a reduction in bacterial richness but had little effect on fungal communities. The bacterial network became less connected and contained fewer negative connections under TFR treatment, and thus, bacterial interactions could be weakened by TFR. The higher abundance of bacterial genes associated with xenobiotic biodegradation and metabolism, and lignin-degrading enzymes under TFR suggests that bacterium-mediated lignin degradation may be strengthened by long-term TFR. Our results suggest that soil bacteria are more susceptible to long-term TFR than fungal communities in warm-temperate forests, with potential consequences for degradation of recalcitrant organics in soil.

Supplementary Materials:
The following are available online at https://www.mdpi.com/1999-4 907/12/2/165/s1, Figure S1: (a) Seasonal variation of daily (black bars) and monthly (grey bars) precipitation under ambient environment; (b) daily averages of automatically (30 min) measured volumetric soil moisture content and (c) soil temperature at 5 cm under control and throughfall reduction (TFR) treatment.; Table S1: Effects (P values) of troughfall reduction (TFR) and years (Y), and their interactions on soil properties, microbial biomass and fine root biomass. p < 0.05 are highlighted in bold; Figure S2: α-Diversity estimates of (a and b) bacterial and (c and d) fungalcommunities under different treatments in 2018 and 2019. OTU diversity and evenness estimates were represented by Shannon-Wiener and Pielou's evenness index, respectively; Control, ambient precipitation; TFR, throughfall reduction. Different letters represent significant difference at p < 0.05; Figure S3: Nonmetric multidimensional scaling (NMDS) ordinations of (a) bacterial and (b) fungal community dissimilarities under different treatments; Figure S4: Relative abundance of (a) bacterial and (b) fungal phyla under different treatments in 2018 and 2019. The groups with relative abundances higher than 1% are shown, while those with less than 1% relative abundance are integrated into "other"; Figure S5: Relative abundance of (a) bacterial and (b) fungal classes under different treatments in 2018 and 2019. The groups with relative abundances higher than 1% are shown, while those with less than 1% relative abundance are integrated into "other"; Figure S6: Relative abundance of (a) bacterial and (b) fungal orders under different treatments in 2018 and 2019. The groups with relative abundances higher than 1% are shown, while those with less than 1% relative abundance are integrated into "other"; Figure S7: Relative abundance of (a) bacterial and (b) fungal families under different treatments in 2018 and 2019. The groups with relative abundances higher than 1% are shown, while those with less than 1% relative abundance are integrated into "other"; Figure S8: Relative abundance of bacterial classes that showed a significant (p < 0.05) response to the TFR in 2018 and 2019; Figure S9: Relative abundance of (a) bacterial and (b) fungal orders that showed a significant (p < 0.05) response to the TFR in 2018 and 2019; Figure S10: Relative abundance of (a) bacterial and (b) fungal families that showed a significant (p < 0.05) response to the TFR in 2018 and 2019; Figure S11: Relative abundance of the predicted bacterial gene of metagenome related to KEGG pathways at level 1 under different treatments in 2018 and 2019; Control, ambient precipitation; TFR, throughfall reduction; * p < 0.05; Figure S12: Pearson correlations between the soil microbial diversity and soil properties and fine root biomass. SM: soil moisture; FRB: fine root biomass; SOC: soil organic carbon; TN: total nitrogen; TP: total phosphorus. Significant correlations are reported as: *, p < 0.05; **, p < 0.01; Table S2: Major properties of microbial networks under control and throughfall reduction (TFR) treatment, and associated random networks.