To investigate the response of the rhizofungal community to foliar herbivory, we used Chinese tallow [T. sebifera ‘Small’ (Euphorbiaceae)] and the armyworm S. litura (Lepidoptera: Noctuidae; a leaf-eating caterpillar) as our study system. Chinese tallow seeds were collected in November 2018 from a natural population near Wuhan, China (31°33′N, 114°07′E). Spodoptera litura larvae were obtained from Henan Jiyuan Baiyun Industry Co., Ltd., and reared on an artificial diet (Tu et al., 2010).

Triadica seeds were sown in an unheated greenhouse at Henan University, Henan, China (34°49′N, 114°18′E), in April 2019. Four weeks later, similar-sized plants (with four or five leaves on them) were transplanted into individual pots (height: 15 cm, upper diameter: 19 cm, lower diameter: 12 cm). These pots were filled with a commercial soil and topsoil (surface litter was removed and then the topsoil was passed through a 1 cm mesh screen) that came from an abandoned field that had been used to grow tallow trees at Henan University (1:1 v:v). To diminish the space effects, we randomly shifted the position of the plants every day. To protect against natural insect herbivores, we enclosed all the prepared tallow plants within a nylon mesh cage (6 m × 4 m × 2 m, 0.1 mm mesh size). Then, 50-day-old healthy and uniform plants were used for the following herbivory treatments.

To understand the effects of foliar insect herbivory on soil fungi dynamics, we established treatment groups consisting of three different herbivory intensities (no herbivory, NH; low herbivory, LH; and high herbivory, HH) and three different durations of elapsed time since herbivore removal (2 weeks since herbivore removal, 2 W; 4 weeks since herbivore removal, 4 W; and 6 weeks since herbivore removal, 6 W). The herbivory treatments were implemented as follows: LH: One S. litura larva was placed on one healthy plant and removed after 7 days of feeding (~10% leaf area consumed); HH: Four S. litura larvae were placed on one healthy plant and removed after 7 days of feeding (~60% leaf area consumed). The control group (NH) comprised healthy plants with no herbivores. Each treatment had three replicates. To investigate the effects of time since herbivore removal upon soil fungal communities, we timed our harvests for 2 W, 4 W, and 6 W post-herbivore removal.

Holding their shoots, the experimental plants were gently shaken to remove any loosely bound soil, and then the soil that remained tightly bound to the surface of the fine roots was collected by manually brushing the soils and designated the rhizosphere soil (Philippot et al., 2013; Fan et al., 2017). In total, 27 rhizosphere soil samples were collected. All of these samples were stored at −20°C, and their DNA was extracted within 2 weeks.

Using a Power Soil DNA kit (MO BIO, Carlsbad, CA, United States), total DNA from each soil sample was extracted and stored at −40°C. To quantify the DNA concentration, a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, United States) was used; the DNA was diluted to approximately 25 ng/μl using distilled water and stored at −20°C until analysis.

The nuclear ribosomal internal transcribed spacer region (ITS rRNA gene) was amplified by polymerase chain reaction (PCR), using the fungal primer set ITS5-F (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and ITS2-R (5′-GCTGCGTTCTTCATCGATGC-3′; Bokulich and Mills, 2013; Shi et al., 2020b). For the PCR, the samples were initially denatured at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, primer annealing at 50°C for 30 s, and extension at 72°C for 30 s; a final extension step of 5 min at 72°C was added, to ensure complete amplification of the target region. Each sample was amplified in triplicate, and the ensuing PCR products from each sample were pooled and purified. Next, the PCR products from each sample were combined in equimolar ratios in one tube and sequenced using an Illumina NovaSeq 6,000 platform at NovoGene, Beijing, China.¹

Using QIIME 2 software (version 2019.7; Bolyen et al., 2019), the raw data sequences were processed and analyzed by following the workflow at https://qiime2.org. Briefly, to obtain the sub-operational taxonomic unit (OTU) table, quality control of the raw sequencing data was performed using the Deblur tool (Amir et al., 2017). To remove low-quality regions from the sequences, each sequence was truncated at 120 bp according to a sequence quality plot. De novo chimera filtering by the VSEARCH tool (Rognes et al., 2016) was used to identify and filter out chimeras. Although the high sequence variability of the ITS marker gene can lead to a more accurate taxonomic identification than is possible using SSU rRNA, it also makes sequence alignment a lot more unreliable. To access reliable phylogenetic information from the ITS region, ghost-tree (Fouquier et al., 2016) was used to create a phylogenetic tree. For this, a pre-built version (ghost_tree_90_qiime_ver7_dynamic_10.10.2017) was used and the command “qiime vsearch cluster-features-closed-reference” was applied, to filter the feature table so as to retain those IDs that matched the ghost tree, and nearly 65% of high quality sequences (1,005,160/1,546,429) were retained. The downstream data analysis was based on the retained sequences. Subsequent taxonomy assignment was done using the Sklearn-based taxonomy classifier with the dynamic Unite database from 10 October 2017.² To compare fungal diversity metrics between the six treatments, 7,017 high-quality sequences were randomly selected for each sample.

The assembly of rhizosphere soil fungal communities under different foliar herbivory treatments was estimated using iCAMP, which is based on the turnover of each bin across samples (Ning et al., 2020). The five assembly processes—drift (DR), homogenizing dispersal (HD), dispersal limitation (DL), heterogeneous selection (HeS), and homogeneous selection (HoS)—were quantified in this study. Briefly, running the iCAMP procedure entails three major steps: the first is phylogenetic binning (dividing the observed data into different groups, namely “bins”); the second is the bin-based null model analysis to partition deterministic and stochastic contributions to HoS, HeS, HD, DL, and DR; and the third is to robustly integrate the results of different bins to assess the relative importance of each process (Ning et al., 2020). The minimal size of all bins in this study was 12, with 49 fungal bins obtained. A confidence index (percentage of null values less extreme than the observed value, i.e., a non-parametric one-side confidence level) was applied for the null model significance test.

Following the approach described by Weiss, co-occurrence networks were constructed using the SparCC package (Weiss et al., 2016). To improve the reliability of these networks, the OTU table was first filtered. Here we constructed six networks corresponding to the six treatment groups: NH, LH, HH, 2 W, 4 W, and 6 W. First, any singletons were removed, and then we retained only those OTUs present in more than 1/3 of all samples per treatment. In this way, 228 OTUs in NH, 203 OTUs in LH, 169 OTUs in HH, 187 OTUs in 2 W, 184 OTUs in 4 W, and 188 OTUs in 6 W were ultimately obtained and used to build the networks [following Weiss et al. (2016)].

Robustness testing is a powerful method to measure a network’s stability and has been used to investigate microbial community stability in multiple studies (Fan et al., 2017; Shi et al., 2019a; Wu et al., 2021). To do that in our study, each network natural connectivity was estimated by “attacking” the nodes (May, 1973) or edges (Jordan, 2009) in the SparCC network. In addition, following Banerjee et al. (2018), in our study network hubs, module hubs, and connectors were defined as keystone species (Shi et al., 2020b). To identify these hubs and connectors, the z-score (within-module degree) and c-score (participation coefficient) of each node within each network were calculated according to the method described by Guimera and Amaral (2005). On the basis of threshold values of the within-module degree (z-score) and participation coefficient (c-score) of nodes, those nodes with a z-score > 2.5 and c-score > 0.6 were designated network hubs, while those with a z-score > 2.5 and c-score < 0.6 were defined as module hubs; nodes with a z-score < 2.5 and c-score > 0.6 were deemed connectors; and nodes with a z-score < 2.5 and c-score < 0.6 were denoted as peripherals. The respective roles of these types of node in community networks has been described at length by Shi et al. (Shi et al., 2020b).

Using high-throughput sequencing, we obtained between 7,071 and 44,737 high-quality fungal sequences per sample (mean value: 27,921), of which >99% were classified into 1,569 distinct OTUs. Over 70% of the total OTU sequences belonged to the Acomycotas (Figure 1). The number of observed OTUs was not influenced by the stage of herbivore removal, but it did decrease with increasing herbivory intensity (Figure 2). The soil fungal communities were similar under the differing herbivory intensities and the time since herbivore removal (Table 1). These findings suggested that the rhizosphere soil fungal community was not affected by how long the plants were free of herbivore pressure or the herbivore intensity, but the number of observed species was influenced by herbivory intensity.

Not only the fungal community as a whole but also the individual species members exhibited variability in sensitivity to time since herbivore removal and herbivore intensity (Table 2). The relative abundance of many species of Agaricomycetes and Sordariomycetes was significantly increased the longer the plants remained free of herbivory. Among the three herbivory intensities, the HH treatment increased the relative abundance of several fungal species and decreased one species in the Agaricales. Furthermore, when compared with LH, the relative abundance of several species in the Ascomycota was decreased by the HH treatment.

The iCAMP analysis revealed that drift was the most important process among the five responsible for rhizosphere fungal assembly, with an average relative importance of 79.4%–96.3% (Figure 3; Supplementary Table S1). However, foliar herbivory evidently altered the relative importance of the five assembly processes. The LH and HH treatments both increased the relative importance of drift and decreased the ratio of homogeneous selection, and HH also increased the ratio of dispersal limitation. Interestingly, the ratio of homogeneous selection was enhanced under the 6 W treatment; this suggested that how processes contribute to fungal assembly depends on the time elapsed since herbivore removal. In short, both insect herbivory and times since herbivore removal could alter the relevance of the five processes for rhizosphere soil fungal assembly.

Using iCAMP enabled us to better detect and infer the assembly mechanisms, by providing information on the relative importance of different ecological processes in individual lineages (i.e., bins). Here, the observed 1,569 OTUs were divisible into 49 phylogenetic bins (Supplementary Table S2), yet we found that the relative importance of a given assembly process was independent of the relative abundance in bins (Figure 4). Our results also revealed that the relative importance of each bin not only differed within the same treatment but also among the treatments (Figure 5; Supplementary Table S3).

Robustness analysis is a powerful way to test a network’s stability. Accordingly, robustness of fungal networks to the treatments was also tested here, by altering the amplitude of natural connectivity via the deletion of nodes and edges (Figure 6). Nodes and edges were discarded in declining order of nodes’ betweenness. Doing this, we found that the natural connectivity of fungal networks presented sharp slopes in HH and 2 W treatment samples, indicating they had weak stability (Figure 6).

To identify the keystone species, z-scores and c-scores were calculated for the nodes in the network for each treatment. There were 4, 2, and 1 keystone species in the NH, LH, and HH treatment groups, respectively, and 7, 0, and 2 keystone species in the 2 W, 4 W, and 6 W treatment groups, respectively (Supplementary Table S4). This pattern suggested foliar insect herbivory as well as time since herbivore removal could slightly reduce the number of keystone species present in the soil fungal community of tallow plants.

In contrast to our hypothesis, herbivory intensity and elapsed time since herbivore removal did not significantly influence the fungal community dissimilarity. This suggests the rhizosphere fungal community may not sensitive to herbivory from leaf-eating insects. However, the fungi are sensitive to partial defoliation is evident by the shift upwards or downwards seen in certain fungal species under the experimentally imposed herbivory conditions. In particular, the number of observed rhizofungal species was found to be significantly lower in the HH than the NH and LH treatments. A previous study showed that aboveground herbivory could affect soil pH (Custer et al., 2020). In addition, myriad studies have reported that the microbial community is strongly related to soil pH (Chu et al., 2010; Fierer, 2017; Xue et al., 2018; Gong et al., 2019). Therefore, the lower number of observed species in HH treatment may be related to altered soil pH, indirectly caused by foliar insect herbivory. We thus suspect that insect herbivory aboveground would enhance root exudation levels, and this could lead to a changed rhizosphere soil pH that influences fungal diversity (Rudrappa et al., 2008; Hoysted et al., 2018).

Although the drift process drove most of the variation in fungal community composition across the samples from all treatments, slight changes in homogenous selection were discernible in some samples. For example, LH and HH treatments both reduced the relative contribution of the HoS process, while the DL process increased under HH. HoS refers to selection occurring under homogeneous abiotic and biotic conditions in both space and time (Dini-Andreote et al., 2015; Ning et al., 2020). So far as we know, deterministic processes (including HoS) which are mainly driven by biotic and abiotic factors are crucial for microbial community structure, therefore, the disturbed HoS process probably result from altered rhizosphere environments by foliar herbivory (Fan et al., 2017; Zhalnina et al., 2018; Shi et al., 2019b). This suggests herbivory has the potential for reducing the rhizosphere selection ability and limiting the community-level dispersal. This might be harmful for plant health, because the growth of plants is strongly related to their recruited microbes, and diminishing their dispersal ability would likely inhibit their recruitment ability (Bakker et al., 2018). Moreover, we found that the relative contribution of HoS increased at 6 W after herbivore removal. Probably, rhizosphere fungi possess self-regulatory ability after herbivory, which might be helpful for plant growth.

The iCAMP approach let us gain insight into the contribution of each fungal group by bin. Our results suggested that the contribution of each phylum or class in each bin was independent of their relative abundance (Figure 4), and the bin groupings contributed differentially to the five assembly scenarios in each treatment. Although other research has focused on the relative contribution of these processes to variation in soil microbial communities, they often neglect the critical role of each phylum or class. But, the role of these species is very important, because it is beneficial to supplement our knowledge of how to manage key microbial assembly processes by mediating the activity of the highly-contributing microbes (Ning et al., 2020).

Microbial community stability is crucial to ensure ecosystem functioning (Wu et al., 2021). Our study found that aboveground herbivores were quite capable of disturbing rhizofungal community stability. In particular, a high herbivory level significantly reduced the network stability, and this may portend a certain loss of rhizosphere ecological functioning. Further, the stability of the rhizofungal community was lowest at 2 W post-herbivory, and it tended to be stable at 6 W after herbivore removal. This result implied that the plant rhizosphere might exercise a self-recovery function. Numerous studies indicate microbial community stability can be affected by one or more local environmental stresses. For instance, Wu et al. (2021) found that enhanced warming could significantly reduce the microbial stability in the active layer, and Fan et al. found lower network stability in the microbial network of the rhizosphere than in bulk soil (Fan et al., 2017); however, few studies have considered the possibility of temporal or lag effects of biotic or abiotic disturbances to plants. This is the first study we know of that tested rhizofungal community stability under different herbivory intensities and at differing times after herbivore removal from the host plants. Given that plant roots have ability in filtering rhizosphere microbial communities and root exudation can mediate plant-microbe interactions (Fan et al., 2017; Zhalnina et al., 2018), while the root traits can be altered by aboveground herbivory (Karssemeijer et al., 2020), we suspect that foliar herbivory impacts on rhizofungal community network may mediated by the cascading effects on plants (Friman et al., 2021).

Collectively, our findings reveal that leaf herbivory from an aboveground insect herbivore can affect rhizofungal diversity, with potential weakening of fungal community stability under high herbivory intensity and soon after herbivore removal. The reduced fungal diversity and stability and disrupted assembly processes due to foliar herbivory may impose detrimental effects on plant health in various systems. However, these detrimental effects could be diminished with the extension of herbivore removal, it means the rhizofungal community may possess a self-recovery ability. Our results indicate that aboveground insect herbivores not only influence the rhizosphere microbial community diversity and stability, but might also trigger cascading effects on rhizosphere ecosystem function and cause diminished vitality in this active region belowground. Foliar herbivory should therefore be taken into consideration in future field and greenhouse studies that manipulate soil microbes to improve plant health. Terrestrial ecosystem includes a vast variety of herbivore guilds, plant species and soil types, whether foliar herbivores trigger similar fungi responses in other systems remain unclear, and needs to be further investigated in future studies.