Zhang et al. Microbiome (2021) 9:196 https://doi.org/10.1186/s40168-021-01150-6 RESEARCH Open Access Gammaproteobacteria, a core taxon in the guts of soil fauna, are potential responders to environmental concentrations of soil pollutants Qi Zhang1, Zhenyan Zhang1, Tao Lu1, Yitian Yu1, Josep Penuelas2,3, Yong-Guan Zhu4,5 and Haifeng Qian1* Abstract Background: The ubiquitous gut microbiotas acquired from the environment contribute to host health. The gut microbiotas of soil invertebrates are gradually assembled from the microecological region of the soil ecosystem which they inhabit, but little is known about their characteristics when the hosts are under environmental stress. The rapid development of high-throughput DNA sequencing in the last decade has provided unprecedented insights and opportunities to characterize the gut microbiotas of soil invertebrates. Here, we characterized the core, transient, and rare bacterial taxa in the guts of soil invertebrates using the core index (CI) and developed a new theory of global microbial diversity of soil ecological microregions. Results: We found that the Gammaproteobacteria could respond indiscriminately to the exposure to environmental concentrations of soil pollutants and were closely associated with the physiology and function of the host. Meanwhile, machine-learning models based on metadata calculated that Gammaproteobacteria were the core bacteria with the highest colonization potential in the gut, and further identified that they were the best indicator taxon of the response to environmental concentrations of soil pollution. Gammaproteobacteria also closely correlated with the abundance of antibiotic resistance genes. Conclusions: Our results determined that Gammaproteobacteria were an indicator taxon in the guts of the soil invertebrates that responded to environmental concentrations of soil pollutants, thus providing an effective theoretical basis for subsequent assessments of soil ecological risk. The results of the physiological and biochemical analyses of the host and the microbial-community functions, and the antibiotic resistance of Gammaproteobacteria, provide new insights for evaluating global soil ecological health. Keywords: Gut microbiota, Pesticide, Indicator taxa, Machine learning, Soil invertebrate, Antibiotic resistance genes * Correspondence: hfqian@zjut.edu.cn 1 College of Environment, Zhejiang University of Technology, Hangzhou 310032, P. R. China Full list of author information is available at the end of the article © The Author(s). 2021 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 http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Zhang et al. Microbiome (2021) 9:196 Page 2 of 17 Background [17], predict the occurrence of fusarium wilt in soil Gut microbiotas help to maintain the health and physio- based on soil microbiotas [18], and predict microplastic logical homeostasis of metazoans, in which epithelial type and its disposal methods based on characteristic en- membranes are the most important interfaces between vironmental microbiotas [19]. The emergence of re- an organism and its environment [1]. In particular, sym- search interest in the gut microbiotas of soil fauna, the biotic gut microbes promote the digestion of complex substantial development of sequencing technology, the substrates and reduce intestinal colonization by patho- advocacy of publicly available data, and machine learn- gens through training of the host immune system and by ing, together, provide an effective route to obtain and affecting body maturation [2]. The complexity of mam- merge the considerable public sequencing data related malian gut microbiotas hinders the analysis of host- to the gut microbiotas of soil fauna, for the identification microbiota interactions, so invertebrates harboring sim- of not only the core gut microbiotas of soil fauna but ple gut microbiotas are an important experimental also of indicator species that respond to soil pollutants. model for determining the relative contributions of indi- Oxytetracycline (OTC) is the most commonly used vidual microbes to the physiological processes of the antibiotic in animal husbandry and several studies indi- host [3]. cated that OTC is a common residual antibiotic in soil Soil fauna account for 23% of known animals and are to which manure has been applied [20, 21]. The strobi- essential components of soil ecosystems, involved in the lurin fungicide azoxystrobin (AZ) is the highest selling decomposition of litter, the transfer of energy, and the fungicide, with $1.165 billion in global sales in 2016 formation of soil microstructure [4, 5]. Recent studies (http://cn.agropages.com/, accessed on October 19, have found that the gut microbiotas of soil invertebrates 2017). Both the antibiotic OTC and the fungicide AZ are important for microbial research and have been as- are broad-spectrum antimicrobials, and are the most sociated with the enrichment and dissemination of genes representative of soil antimicrobials. Therefore, we used conferring antibiotic resistance [6], the colonization and AZ and OTC as typical treatments to investigate the transfer of pathogens [7], and the potential for the cyc- changes of the gut microbiota in the model soil inverte- ling of soil elements [8], thereby expanding the eco- brate Folsomia candida (Collembola) and analyzed func- logical and evolutionary potential of the hosts [9]. The tional genes using eukaryotic transcriptome sequencing, gut tract of soil invertebrates is dominated by transient which could identify the interactions between the gut microbiota [5], depending on time, environmental condi- microbiota and host physiological functions. We then tions, nutrient availability, stage of growth and health of defined the core gut microbiota of F. candida using a the host, and even circadian rhythms [10]. Therefore, large-scale analysis of data for the microbial community characterizing the core gut microbiotas of soil inverte- and identified taxa that could be used as potential indi- brates is challenging. cators of environmental concentrations of soil pollutants Soil pollution is one of the most serious hidden dan- using random forest machine learning. We further inves- gers to global agricultural food security [11]. Soil inver- tigated the antibiotic resistance genes (ARGs) in F. can- tebrate faunae, as some of the first organisms exposed to dida gut using HT-qPCR and identified the taxa related soil pollutants, are often used as indicator species in to the ARGs. toxicological experiments [12, 13]. The residual environ- mental concentrations of soil pollutants are usually low Materials and methods and rarely greatly affect the physiology of soil inverte- Details are provided in Supplementary Information and brate fauna in the short term. The microbiotas within Figure S1. the guts of soil invertebrates have more sensitive charac- teristics than their hosts or even soil microbiotas [14, Laboratory experiments 15]. The process to identify core microbiota is inevitably Experimental design affected by different experimental methods, PCR The test species, Folsomia candida (“Berlin strain,” ori- primers, sequencing depth, and other technical opera- ginally obtained from Aarhus University, Denmark), was tions [16]. Thus, core microbiota that can stably exist reared in our laboratory for 2 years following within the gut of a species have the potential to become Organization for Economic Co-operation and Develop- an indicator organism for general study. The develop- ment (OECD) guideline 232 and was placed in a suitable ment of a single test, however, is labor intensive and breeding environment, maintained in Petri dishes con- would not be universal, so identifying a general indicator taining a mixture of charcoal and plaster of Paris (1:8 w/ taxon is difficult. w). These Petri dishes were kept in a thermostatic box A method combining metadata analysis and machine (Safe Co., Ningbo, China) at 75% relative humidity and a learning has recently been used to predict the develop- temperature of 20 ± 2 °C with a 16:8-h light: dark photo- ment of cancer based on human intestinal microbiotas period (800 Lux). Ultrapure water was added once a Zhang et al. Microbiome (2021) 9:196 Page 3 of 17 week to ensure substrate moisture and the F. candida adults and larvae in all microcosms after 28 days of ex- were fed twice a week with dry yeast powder. posure using the imageJ 1.52, and 13 adults were col- To obtain F. candida, all of the same age, we trans- lected for the locomotory test. Two adults were then ferred 60-70 active adults to a new substrate to lay eggs isolated for determining ROS content and CYP450 en- for 2 days, and then removed these adults. As the juve- zymatic activity respectively, following the manufac- niles hatched, they were immediately transferred to a turer’s instructions for the nematode ELISA kit based on new substrate for culture. Before the exposure experi- the double-antibody sandwich method (Jiangsu Enzyme ment, 7-9-day-old collembolans were transferred to the Industrial Co., Ltd., Yancheng, China) [15]. test soil (2.27, 16.47, and 81.25% of clay, silt, and sand, In the other experimental setup, sixty 14-16-day-old respectively; pH 5.16; water-holding capacity, 46.8%; pre-incubated collembolans were added to control, AZ, total N content, 3.8 g/kg) from a vegetable field (29° 49′ OTC, and AO experimental microcosms with eight rep- N, 121° 20′ E; Zhejiang, China) to adapt the test collem- licates. These experimental microcosms were conducted bolans to the new culture environment. in same-sized beaker with 65 g moist soil. After expos- When the mortality of these collembolans in the new ure for 28 days, all collembolans were collected for gut soil culture environment was less than 1% during the microbiome, resistome, and eukaryotic transcriptome pre-incubation process, we used these collembolans for analysis. Fifty adults per replicate were washed three subsequent exposure experiments. To simulate the nat- times with ultrapure water and dissected under a stereo ural field environment, we did not feed in the pre- microscope using very precise tweezers to obtain the incubation and exposure experiments and added distilled collembolan gut samples [6]. water twice weekly to maintain soil moisture. The gut samples were transferred to 2-mL round- We firstly conducted a pre-exposure experiment, to bottomed centrifuge tubes containing 20 μL of protein- measure reproduction rate and mortality, for selecting ase K and 180 μL of a lysis buffer solution for the extrac- the most suitable concentration for the formal experi- tion of DNA [15]. A total of 2 g of soil per sample, ments, and the specific details are described in the Sup- without collembolans, was collected for analyzing the plementary Information. The formal laboratory soil microbial community. Fifty adults each replicate experiments were divided into two parts, one for deter- were provided for extraction of RNA. These collected mining the death, reproduction, locomotion, reactive collembolans were firstly washed three times with ultra- oxygen species (ROS) concentration, and cytochrome pure water and then immediately frozen in liquid nitro- P450 (CYP450) enzymatic activity of F. candida, and the gen for RNA extraction. other for transcriptomic measurements (gene expres- sion), the gut microbiota (bacterial and fungal communi- Effects of soil pollution on F. candida gut microbiota ties), and the gut resistome of F. candida. Before starting A total of 50 guts of adult F. candida, per replicate, were the experiments, we added ultrapure water to restore used to isolate high-quality DNA using a DNeasy® Blood the microorganisms in the soil for a week at room & Tissue Kit (QIAGEN, GER). The V4 hypervariable re- temperature, and the moisture content was maintained gion of the 16S rRNA gene was amplified using universal at 60% of the maximum (60.23%), as detailed in OECD primers (forward primer 515F 5′-GTGCCAGCMGCC guideline 232. The laboratory experiments were divided GCGG-3′ and reverse primer 806R 5′-GGACTACNVG into two parts, one for determining the death, GGTWTCTAA-3′) [14], and region 1 of the internal reproduction, locomotion, ROS concentration, and transcribed spacer (ITS) gene was amplified using the CYP450 enzymatic activity of F. candida, and the other forward primer ITS1F (5′-CTTGGTCATTTAGAGGAA for transcriptomic measurements (gene expression), the GTAA-3′) and the reverse primer ITS2 (5′-GCTGCG gut microbiota (bacterial and fungal communities), and TTCTTCATCGATGC-3′) [21]. The 50-μL reactions (25 gut resistome of F. candida. In one part of the experi- μL of TaKaRa ExTaq DNA polymerase, 1 μL of DNA ment, soil microcosms were established for exposing F. (range 10-15 μg/mL), 1 μL of universal forward, 1 μL re- candida to the pollutants: control (no pollutants), AZ verse primer, and 22 μL of PCR-grade water) were amp- (0.3 mg AZ/kg dry soil), OTC (10 mg/kg dry soil), and lified following reaction conditions previously described AO (combined 0.3 mg AZ/kg dry soil and 10 mg OTC/ [14]. The PCR products were then purified, pooled, and kg dry soil). Twenty 14-16-day-old pre-incubated col- sequenced using the Illumina MiseqPE300 platform lembolans were introduced into these experimental mi- (Meiji, Shanghai, China). crocosms, each treatment had four replicates and were conducted in sterile glass beakers (inner diameter 5.3 Effects of soil pollution on F. candida gut resistome cm, 6.5 cm high) containing 30 g of moist soil at 18 °C We used a total of 384 primer sets (Data set S2) (includ- with a diurnal light cycle (16:8 h light: dark). We added ing 320 antibiotic resistance genes (ARGs), 57 mobile sterile water twice a week and recorded the numbers of genetic elements (MGEs), and the 16S rRNA gene) to Zhang et al. Microbiome (2021) 9:196 Page 4 of 17 investigate the composition and abundance of ARGs in sample sequence was normalized to 2000 reads, and a the collembolan gut using the SmartChip Real-time PCR total of 415 samples from the 17 independent studies, system (Wafergen, USA). The PCR reaction mixture of including the 17 soil pollutants (e.g., the fungicide azox- each well compose with collembolan DNA template, ystrobin; the insecticide cypermethrin; the herbicide gly- primers, sterile water, and LightCycler 480 SYBR Green phosate; the antibiotics tetracycline, sulfamethoxazole, I Master mix. The reaction system (95 °C (10 min) and and oxytetracycline; the antibiotic substitute Macleaya 40 cycles of 95 °C (0.5 min) and 60 °C (0.5 min)) was cordata extract; the heavy metals arsenic, silver nitrate, used for HT qPCR. SmartChip qPCR software was used silver nanoparticles, and nano-copper oxide; and the to analyze the raw data and a threshold cycle (CT) of 31 emerging pollutants micro-, nano-, and tire-tread plas- was used to detect ARGs. Only when 3 technical repli- tics), five kinds of soil invertebrates (collembolans, cates and 3 biological replicates were amplified at the enchytraeus, earthworms, mites, and ants), and three ex- same time, did we consider an ARG to have been de- posure methods (oral exposure, soil microcosm, and tected. The relative abundance of ARGs was calculated field experiment). All metadata were divided into “Con- using the equation below [22]. trol” and “Pollution” based on the sample information Copy number of gene = 10 ((31-CT (measurement))/(10/3) for each experiment (Data set S1). Effects of soil pollution on F. candida gene expression at Construction and validation of the predictive model the transcriptomic level Machine-learning algorithms were used to predict the A total of 50 collected adults per replicate were immedi- heterogeneous microbial data for identifying taxa in the ately snap frozen in liquid nitrogen and stored at −80 °C guts of the soil invertebrates that were most strongly as- to ensure the RNA integrity. The RNA was isolated from sociated with the stress of the soil pollutants. We uti- each replicate pool using an RNA extraction kit (HiPure lized three machine-learning algorithms, random forest Universal RNA Midi Kit, Magen, Guangzhou, China). (RF) [23], logistic regression (LR) [24], and support- We checked the concentration, purity, and integrity of vector machine (SVM) [25]. The receiver operating char- the RNA and isolated the mRNA. The enriched mRNA acteristic curve (ROC) and the area under the curve in- was then randomly broken into fragments of ~300 bp dicated that the RF algorithm performed well on our and reverse transcribed to produce stable double- merged data. We therefore used the RF algorithm to stranded cDNA using ReverTra Ace qPCR RT Kit build the predictive model. Classification models based (TOYOBO, Osaka, Japan), which was sequenced using on each taxonomic level could distinguish between the the Illumina Novaseq 6000 platform (Meiji, Shanghai, bacterial communities in the guts of the soil inverte- China). brates in the control and pollutant treatments using the randomForest package in R with default parameters Data collection and description and processing of the 16S (version 4.6-14). The results indicated that the average rRNA gene metadata accuracy rate was similar across all taxa, so we selected Based on a previous search principle [19], we searched the RF model at the genus level with the lowest esti- the Web of Science Core Collection and Science Direct mated rate of out-of-bag (OOB) errors (17.58%) [26]. for the terms “gut microbiota of soil fauna or ‘species’ name” “gut microbial community of soil fauna or ‘spe- Statistical methods cies’ name” and “gut microbiota of soil invertebrate or The means ± standard errors (SEs) of each treatment ‘species’ name.” A total of 33 studies were collected from were calculated. A two-tailed Welch’s t test was used to these databases, but only 20 independent experiments identify significant differences between groups. A princi- were publicly available with incomplete 16S rRNA gene pal coordinate analysis (PCoA) based on unweighted sequences (Data set S1), including three 16S rRNA gene UniFrac distances for the guts of the soil invertebrates hypervariable regions, V4, V4-V5, and V3-V4, which and for the bacterial and fungal communities in the sur- were amplified using the respective primer pairs 515F/ rounding soil was performed using the Majorbio Cloud 806R (forward primer: 5′-GTGCCAGCMGCCGCGG Platform (www.majorbio.com), and the output was visu- TAA-3′, reverse primer: 5′-GGACTACNVGGGTWTC alized using OriginPro 9.1. The Adonis function (9999 TAA-3′), 515F/907R (forward primer: 5′-GTGCCAGC permutations) was used in a PERMANOVA to identify MGCCGCGGTAA-3′, reverse primer: 5′-CCGTCAAT differences among the treatments using the vegan 2.4-3 TCMTTTRAGTTT-3′), and 334F/806R (forward pri- package in R version 3.6.1. Function prediction analysis mer: 5′-CCTACGGGAGGCAGCAG-3′, reverse primer: of fungi using FUNGuild was performed on Majorbio 5′-GGACTACHVGGGTWTCTAAT-3′). We used Cloud Platform (www.majorbio.com). Vsearch (version 2.7.1) and QIIME2 (version 2021.2) to Heatmaps were generated using TBtools (Toolbox Bi- filter, classify, and merge these raw data sets. Each ologists v0.655), and histograms, line and box charts, Zhang et al. Microbiome (2021) 9:196 Page 5 of 17 and linear regressions were produced using GraphPad data. The SEM was produced using Amos Graphics v22 Prism 8.00. The weighted gene co-expression network (IBM Corp., Armonk, NY, USA). Meta-analysis and sen- analysis (WGCNA) was performed using the Majorbio sitivity analysis were performed using the STATA statis- Cloud Platform (www.majorbio.com). The co- tical software package version 15.0 (Stata Corp., College occurrence network analysis of laboratory experiments, Station, TX, USA). In addition, to characterize positive based on the relative abundance of all bacterial and fun- and negative co-occurrences separately, we used the co- gal classes, using pairwise Spearman’s rank correlations hesion among taxa to reveal the interactions, similarity, (r) in the psych package in R (r > 0.6 or r < −0.6, P < and differences between both positive and negative spe- 0.05), was performed using Gephi v0.9.2. The shared cies interactions in the niches of microbial taxa, using network based on the frequency of each bacterial class the previous equation to calculate the positive and nega- among all independent studies was produced using tive cohesion values [27]. Gephi v0.9.2. Structural equation models (SEMs) were built to calculate the direct and indirect effects among Results the gut bacteria (Shannon index), treated groups, bacter- Effects of AZ, OTC, and AO on the physiology and iaShannon/fungiShannon (B/F) index, cytochrome P450 (en- biochemistry of F. candida. zymatic activity), the HAA index, Gammaproteobacteria We designed a locomotory map (radii of 5, 10, 15, 20, (relative abundance), and transcriptome (PC1 of the and 30 cm) to characterize the motility of F. candida TPM value using Bray-Curtis distances). The signifi- after exposure to pollutants (Fig. 1A). Thirteen individ- cance of each path-coefficient was analyzed by calculat- uals of the same size were treated with 0.3 mg AZ/kg ing its critical ratio (P < 0.05). The goodness-of-fit index dry soil, 10 mg OTC/kg dry soil, or AO for 28 days and (GFI) and the Bentler comparative fit index (CFI) indi- were then put at the center of the map. The number of cating the goodness-of-fit of the models to the original individuals in areas A1-A5 after 2, 4, 6, 8, and 10 min Fig. 1 Effects of azoxystrobin (AZ), oxytetracycline (OTC), and their combination (AO) on the physiology and biochemistry of Folsomia candida. Locomotory test map (A). The number of the individuals in different areas of the map after 2-10 min of exposure to AZ, OTC, and AO (B). The high area activity (HAA) index (C), with different K values (D). E The numbers of adults and larvae, ROS concentrations, and the enzymatic activity of CYP450 after exposure to AZ, OTC, and AO. * (P < 0.05) and ** (P < 0.01) indicate significant differences between the bacterial and fungal communities (two-tailed Welch’s t test) Zhang et al. Microbiome (2021) 9:196 Page 6 of 17 were recorded (Fig. 1B). The AZ and OTC concentra- concentrations of pollutants would enrich potential soil tions, which did not cause mortality, were chosen ac- pathogens, inducing ecological risks. cording to the environmental residual levels and pre- experiment results (Figure S2). F. candida motility was Effects of AZ, OTC, and AO on gut microbiotas of F. calculated using Eq. 1 (Supplementary Information), and candida an index of high area activity (HAA) was calculated for Analysis of the 16S rRNA and ITS gene sequences from each time point. F. candida motility for the treatments 20 adult F. candida guts showed that richness was sig- was in the order control > AZ > TOC > AO based on nificantly higher in the fungal community than the bac- the slopes of the fitted lines (Fig. 1C, D). terial community (P < 0.001, df = 30, two-tailed Welch’s The numbers of adults and larvae did not differ signifi- t test) (Figure S5), indicating that the potential cantly between the treated groups and the control. The colonization in F. candida guts was much lower for concentration of reactive oxygen species (ROS) was fungi than bacteria. The addition of pollutants did not higher in OTC-, AZ-, and AO-treated groups than the affect the diversity or abundance of the fungal commu- control (Fig. 1E: P < 0.01, P < 0.01, P < 0.05, df = 6, n = nity (Shannon and Chao indices: P > 0.05, df = 6, two- 4 per group, two-tailed Welch’s t test), and the enzym- tailed Welch’s t test) (Figure S6A, B), but unweighted atic activity of cytochrome P450 (CYP450) was signifi- UniFrac distances indicated that the diversity and abun- cantly higher only in the AO group (P < 0.05, df = 6, n = dance of the fungal community differed between the AO 4 per group, two-tailed Welch’s t test). group and the control (P <0.001, df = 26, two-tailed Welch’s t test) (Figure S6C), perhaps due to the inter- action between fungi and bacteria. We established a co- Effects of AZ, OTC, and AO on the microbiota of the occurrence network of the bacterial and fungal species surrounding soil. to test this hypothesis and found that the index of stabil- A total of 0.5 g of soil surrounding the F. candida speci- ity (negative: positive cohesion) for all microbial commu- mens was carefully separated from the microcosms, both nities was highest in the AO group, representing the with and without treatment, under a microscope to most unstable community network (Figure S7). avoid interference by the larvae for determining the ef- Exposure to the pollutants significantly affected the di- fects of AZ, OTC, and AO on the diversity and compos- versity and structure of the bacterial community, espe- ition of the soil microbiota. The soil bacterial and fungal cially in the AZ and AO groups (Fig. 2A, B: AZShannon communities were assessed using deep amplicon se- and AZChao: P < 0.01, df = 6, two-tailed Welch’s t test; quencing of the 16S rRNA and internal transcribed spa- AOShannon: P < 0.01, df = 6, two-tailed Welch’s t test; cer (ITS) genes, respectively. Exposure to AZ, OTC, and PERMANOVA for treatments: R2 = 0.6398, P = 0.001, AO for 28 days did not significantly affect the alpha- Adonis analysis). Totals of 2 (1 up and 1 down), 10 (8 (Shannon and Chao indices) or beta-diversity of the soil up and 2 down), and 6 (4 up and 2 down) OTUs differed bacterial community (P > 0.05, df = 6, n = 4 per group, significantly between the OTC, AZ, and AO groups as two-tailed Welch’s t test; PERMANOVA for treatments: compared to the control, respectively (Fig. 2C: P < 0.05, R2 = 0.2764, P = 0.182, Adonis analysis) (Figure S3A, B). df = 6, n = 4 per group, MetagenomeSeq difference ana- The number and classification of operational taxonomic lysis). The associated phylogenetic relationships are units (OTUs) did not differ significantly between the shown in Fig. 2D. The OTUs shared between the treat- treated groups and the control (P > 0.05, df = 6, n = 4 ment groups were classified into Proteobacteria per group, metagenomeSeq difference analysis). (phylum), with Alphaproteobacteria and Gammaproteo- Fungal diversity was significantly higher in the AZ and bacteria (classes) identified as common potential indica- AO groups, although the fungi were richest in the AO tor taxa after exposure to the pollutants (Fig. 2E), similar group (Figure S3C), and the treatments affected the to previous studies [22, 29]. Proteobacteria may be a po- composition of the fungal community (principal coord- tential microbial signature of dysbiosis in soil inverte- inate analysis (PCoA) of unweighted UniFrac distance) brates, just as in human guts [30]. Interestingly, the [28], PERMANOVA for treatments: R2 = 0.2392, P = relative abundance of Gammaproteobacteria gradually 0.024, Adonis analysis) (Figure S3D). The number and increased with the intensity of pollutant stress (OTC < classification of OTUs differed significantly in each AZ < AO) (Fig. 2F: P > 0.05, P < 0.01, and P < 0.05, df = group (Figure S4A-C), and the OTUs with highest abun- 6, n = 4 per group, two-tailed Welch’s t test). The co- dances belonged to plant or animal pathogens (function occurrence network indicated that the Gammaproteo- prediction analysis of fungi using FUNGuild, a tool for bacteria, the most common taxon across all treatments, parsing fungal OTUs) (Figure S4D). This finding indi- was significantly correlated with the abundances of 178 cated that soil microbiotas would respond differently to other bacteria and 16 fungi, indicating that this taxon different soil pollutants and that low environmental occupied an important central position for correlating Zhang et al. Microbiome (2021) 9:196 Page 7 of 17 Fig. 2 Effects of azoxystrobin (AZ), oxytetracycline (OTC), and their combination (AO) on the gut microbiota of Folsomia candida. Alpha diversity (Shannon and Chao indices) (A). PCoA of OTU file data using unweighted UniFrac distances) for different groups of gut bacteria (B). Significantly different OTUs between the control and treated groups (C). Evolutionary relationships of the OTUs (D). Classification of significantly different OTUs at the phylum and class levels and significantly different OTUs shared between the treated groups (E). The relative abundance of significantly different OTUs shared in all treated groups (F). Co-occurrence network analysis identifing the relationships among the bacteria, fungi, and Gammaproteobacteria (G). The size and color of the nodes represent the relative abundance and the associated type of OTU data (black, bacteria; green, fungi; and blue, Gammaproteobacteria). The colored lines indicate the links among the bacteria, fungi, and Gammaproteobacteria. * (P < 0.05) and ** (P < 0.01) indicate significant differences between the bacterial and fungal communities (two-tailed Welch’s t test), NS indicates “not significant” the fungal and bacterial communities (Fig. 2G, Spear- microorganisms (Fig. 3A). A PCoA of gene expression man’s > 0.6, P < 0.05). The associated topological value (TPM) using the Bray-Curtis distances indicated data also showed that Gammaproteobacteria played the differences among the groups (Fig. 3B: PERMANOVA most important role in the co-occurrence network for treatments: R2 = 0.6327, P = 0.001), and Bray-Curtis through interactions with other bacteria and fungi (Fig- dissimilarity indicated that the OTC and AO groups ure S8). were significantly separated from the control (P < 0.001 and P < 0.01, respectively, two-tailed Welch’s t test). We Relationships between Gammaproteobacteria and host classified a total of 24,436 genes (21,359 known, 3077 physiology, biochemistry, and function unknown) into five unique module eigengenes based on We collected 50 F. candida after 28 days of exposure to differences in the strength of the interaction between AZ, OTC, and AO to extract RNA for transcriptome se- genes, conforming to a scale-free distribution, to deter- quencing based on the F. candida genome in NCBI mine whether gut Gammaproteobacteria were associated (GCF_002217175.1) [31] after transferring the gut con- with host functions (R2 = 0.82, soft power (β) = 6, mini- tents to sterile water to avoid contamination from gut mum module size = 30, weighted gene co-expression Zhang et al. Microbiome (2021) 9:196 Page 8 of 17 Fig. 3 (See legend on next page.) Zhang et al. Microbiome (2021) 9:196 Page 9 of 17 (See figure on previous page.) Fig. 3 Correlations among the transcriptome, gut microbiome, and gut Gammaproteobacteria of Folsomia candida (A). PCoA based on the Folsomia candida gene expression value (TPM) data using Bray-Curtis distances showing the different clusters in the control and the treated groups (OTC, AZ, and AO) (B). Heatmap of the level of expression (the color key indicates the TPM value) (C), and associated KEGG functional pathway (the color scale and red text indicate the primary and various secondary pathways, respectively) (D). Ordinary least squares (OLS) linear regression between Gammaproteobacteria (relative abundance) and physiological and biochemical indicators of F. candida (number of juveniles, HAA index, enzymatic activity of CYP450, and ROS concentration) (E). Structural equation model (SEM) of the relationships among the gut bacteria (Shannon index), treated groups, bacteriaShannon/fungiShannon (B/F) index, cytochrome P450 (enzymatic activity), the HAA index, Gammaproteobacteria (relative abundance), transcriptome (PC1 of the TPM value using Bray-Curtis distances), and the goodness-of-fit index (GFI) and the Bentler comparative fit index (CFI) indicating the goodness-of-fit of the models to the original data. Dashed lines indicates the “not significant correlation” (F). The direct, indirect, and total standardized effects of Gammaproteobacteria on the indicators (G) network (WGCNA) analysis) (Figure S9A-C). The MEb- invertebrates. Previous studies reported that the gut rown (gene membership of brown) co-expression mod- microbiotas of soil invertebrates were mainly composed ule was selected for further analysis, because it was most of the phyla Proteobacteria, Firmicutes, Actinobacteria, relevant to the control phenotype, to distinguish be- and Bacteroidetes [32], which are anaerobic and faculta- tween the control and treated groups (module signifi- tive anaerobic bacteria due to the special anaerobic en- cance = 0.451, module membership vs gene significant: vironment of the gut. R2 = −0.697, P < 0.012, WGCNA analysis) (Figure S9D- We unified the determination of the core gut micro- F). Seventy-nine MEbrown genes were significantly biotas of soil invertebrates through collection of data negatively correlated with the relative abundance of from 33 independent experiments by searching the Web Gammaproteobacteria (Spearman’s > 0.6, P < 0.05, of Science Core Collection and Science Direct, 20 of Spearman analysis), and gene expression was lower in which were publicly available and contained incomplete the treated groups than the control (Fig. 3C). These 16S rRNA gene sequences (Data set S1). We then reana- genes were annotated to functions of immunity, diges- lyzed these sequences and homogenized them to 2000 tion, development, detoxification, vitamin carbohydrate reads for each sample. The merged OTU table contained metabolism, protein assembly and synthesis, and cell 415 gut samples from soil invertebrates in 17 independ- growth and apoptosis based on the KEGG pathway data- ent experiments: nematodes, springtails, earthworms, base (Fig. 3D). Locomotive activity (HAA index) and mites, and ants (Data set S1). The gut microbiotas of soil ROS concentration for the host were also significantly invertebrates are contributed by host factors such as diet correlated with the relative abundance of Gammaproteo- [33], habits [34], region [6], and soil type [35], so refining bacteria (Fig. 3E: R2 = −0.4568, P = 0.0041; R2 = 0.3395, the core gut microbiotas of soil invertebrates at a higher P = 0.0395, ordinary least squares linear regression ana- level of classification, the class level, was necessary. First, lysis). These findings indicated that the relative abun- the results of the taxa-detection rate for each experiment dance of gut Gammaproteobacteria was negatively identified seven shared classes: Actinobacteria, Alpha- correlated with the normal physiological functions of the proteobacteria, Ktedonobacteria, Acidimicrobiia, Acido- host, similar to the index of oxidative stress (ROS level), bacteriae, Planctomycetes, and Gammaproteobacteria indicating that the host was stressed. Structural equation (Fig. 4A). We then analyzed the frequency and relative modeling (SEM) found that the gut Gammaproteobac- abundance of each class across all samples (Fig. 4B, C). teria were connected with the expression of host func- The two indicators, frequency and relative abundance, tional genes by interfering with the gut bacterial had different rank orders. A core index (CI) was estab- community and were also directly associated with host lished for calculating the symbiotic potential of core locomotion and level of oxidative stress (Fig. 3F, G: χ2/df microbiotas from all taxa based on a metadata analysis = 0.961, P = 0.450; GFI = 0.913, CFI = 1.00, RAMSEA = using Eq. 2 (Supplementary Information), and we used 0.000, SEM analysis). the normalized core index (NorCI) to determine the core gut microbiota of the soil invertebrates (Eq. 3, Sup- Defining the core microbiotas of the guts of soil plementary Information) (Fig. 4D). According to the pre- invertebrates based on a metadata analysis vious classification method, based on the NorCI Global databases of soil invertebrates have recently been thresholds of 0.83 and 0.3 (Supplementary Information) developed, but are concentrated mainly on species iden- and the frequency of all samples, we further divided the tification (Global Biodiversity Information Facility bacterial classes into core (Gammaproteobacteria, (GBIF), www.gbif.org). The recent rapid development of Alphaproteobacteria, and Planctomycetes; 50.68%), tran- molecular sequencing technology has provided an effect- sient (Bacilli, Ktedonobacteria, Actinobacteria, and Acid- ive way to combine the diversity of global species of obacteria; 35.63%), and rare (lower-CI taxa; 13.69) gut symbiotic microbiota and the gut microbiotas of soil microbiota (Fig. 4E). Interestingly, CI and NorCI for the Zhang et al. Microbiome (2021) 9:196 Page 10 of 17 Fig. 4 Defining the core microbiomes of the guts of soil invertebtrates based on a metadata analysis. Network analysis showing the unique and shared classes in the 17 independent experiments. The node size indicates the number of connections, the node color indicates the following: black nodes, different independent experiments; blue nodes indicate bacterial classes not shared between all experiments; red nodes indicate bacterial classes shared by all independent experiments (A). The frequencies (B) and average abundances (C) of the 20 most abundant classes in all samples, respectively. Core index (CI) for core, transient, and rare bacteria (at the class level) in the guts of the soil invertebrates (f, class frequency; n, number of independent experiments with class data; s, sequence length of the class; S, total sequence length from all samples) and NorCI, the normalized core index using all classes detected to normalize CI (D). Relative abundances of the core (50.68%), transient (35.63%), and rare (13.69%) classes in all samples (E). All data from the independent studies were merged into a single data set for all analyses invertebrates and the surrounding soil (Figure S10A) experiments. The samples were from earthworms, showed that the core taxa in soil were very similar to the Enchytraeus, collembolas, ants, and mites exposed to 15 transient taxa in the invertebrate gut. A PCoA based on kinds of soil pollutants (e.g., the fungicide azoxystrobin; the OTU data using Bray-Curtis distances, however, the insecticide cypermethrin; the herbicide glyphosate; identified a pronounced separation between the soil in- the antibiotics tetracycline, sulfamethoxazole, and oxy- vertebrates and the bacterial communities of the sur- tetracycline; the antibiotic substitute Macleaya cordata rounding soil (Figure S10B). Surprisingly, CI was also extract; the heavy metals arsenic, silver nitrate, silver highest for Gammaproteobacteria in the surrounding nanoparticles, and nano-copper oxide; and the emerging soil, but the PCoA based on Gammaproteobacteria OTU pollutants micro-, nano-, and tire-tread plastics), ma- data demonstrated that the gut Gammaproteobacterial nure, and lime. community differed from the community in the sur- To determine whether gut bacterial taxa could be used rounding soil. This finding may have been due to the as biomarkers for responses to soil pollutants, we built special anaerobic environment of the gut (Figure S10C). three machine-learning models, random forest (RF), support-vector machines (SVM), and logistic regression Screening the most important gut bacterial taxa that (LR), in which the accuracy rate (area under the curve) respond to pollutant pressure using a machine-learning indicated that RF was the best model for predicting the method and associated community function classification of samples (Fig. 5A). We therefore calcu- We screened data for a total of 415 gut samples (193 lated the accuracy of the classification based on the bac- controls and 222 pollution groups) from 17 independent terial data at the levels of phylum, class, order, family, Zhang et al. Microbiome (2021) 9:196 Page 11 of 17 Fig. 5 Screening of gut microbiota that responded to pollutant pressure using a machine-learning method. The machine-learning modules of the random forest (RF), support-vector machine (SVM), and logistic regression (LR) analyses were built using the OTU data for the bacterial community (A). The associated AUC and ROC curves indicated that the RF model was the most accurate model. The RF module constructed using all samples at the phylum, class, order, family, genus, and OTU levels. Dashed lines indicates the 95% confidence interval (B). Prediction using the test data in the RF module. Black indicates the control, and blue indicates the treated groups in the comparison between raw (O) and predictive (P) information (C). Each test used a tenfold cross-validation method to verify the accuracy of the model predictions, and the 18 most abundant bacterial genera were identified by applying the RF classification of the relative abundances of the control and treated samples (D). RF classification of the relative abundance of the control and treated samples as based on the data for all genera for calculating their mean decreases in accurracy, combined with the frequency of each genus for screening the most important indicator genus, Paraburkholderia, belonging to Gammaproteobacteria (E). The associated relative abundance of Gammaproteobacteria in the control and the treated groups (F). P values were determined using two-tailed Welch’s t tests Zhang et al. Microbiome (2021) 9:196 Page 12 of 17 genus, species, and OTUs and selected the most accur- accumulating within the gut, could act as indicators of ate classification of the genus data to establish the RF soil disturbance by pollution. A PCoA based on the model (Fig. 5B). The predication of this model based on OTU data from the rarefied sequencing reads for each the test data set was 100% accurate for the pollution gut sample using the Bray-Curtis distances presented samples and 84.2% accurate for the control samples (Fig. significantly separated clusters between control and pol- 5C). We next conducted a tenfold cross-validation with lutant treatments (Fig. 6A: PERMANOVA for treat- five repeats to evaluate the importance of genera as po- ments: R2 = 0.1256, P < 0.001, Adonis analysis). Bacterial tential indicators (Fig. 5D), which combined with the diversity (Shannon index) was higher in the pollutant frequency data for each genus, and confirmed the bac- treatments (Fig. 6B: P < 0.001, df = 415, two-tailed terial biomarker taxa (Fig. 5E) belonging to Gammapro- Welch’s t test), indicating that pollutant residues in the teobacteria and their enrichment in contaminated guts soil interfered with the structure and diversity of the gut (Fig. 5F). Meanwhile, taking the heterogeneity of differ- microbial community. The topological properties ana- ent experimental pollutants, we used a meta-analysis lysis exhibited an insignificant change between control and a sensitivity analysis based on a random model and and treatment group (Figure S12). The negative:positive found an upward trend of Gammaproteobacteria in the cohesion of each bacterial community and the contaminated group (Figure S11). This finding sup- Gammaproteobacteria-correlated community network in ported the hypothesis that Gammaproteobacteria, by the control and pollutant treatments were calculated Fig. 6 Effects for soil pollutant in stability of microbial networks of soil invertebrate gut. Principal coordinate analysis (PCoA) of the 17 independent experiments with 415 gut samples (193 control and 222 treated), including 17 kinds of soil pollutants (Data set S1) using Bray-Curtis dissimilarity performed on OTU data from the control and treated microbial communities (A). R2 and P values were determined using Adonis analysis. Shannon index for the gut mirobial communities in the control and treated groups (B). P values were determined using two-tailed Welch’s t tests. Negative:positive cohesion (calculated using the formula in (28)) in the whole (C) and Gammaproteobacteria-associated communities (D) networks with and without soil pollution, indicating that the treated groups formed a more stable community network than the control as shown in the schematic diagram (E) Zhang et al. Microbiome (2021) 9:196 Page 13 of 17 using the equation described by Hernandez et al. [27], to community was significantly affected by soil pollution, determine the function of the Gammaproteobacteria but the fungal OTUs were not shared among the pollu- community in the gut. The higher negative:positive cohe- tant treatments. Interestingly, the gut bacterial commu- sion in the bacterial communities of the contaminated nity had a stress response to the pollutants similar to the guts demonstrated that soil pollutants stimulated the gut host physiological response and also shared significantly to develop a more stable community network (Fig. 6C, different OTUs in each treatment group, suggesting that E), especially the Gammaproteobacteria, which greatly the gut and soil bacteria possessed common response promoted the stability of the gut bacterial community characteristics to soil pollution, indicating that the gut (Fig. 6D). bacteria may not have been disturbed by specific targets of the soil pollutants. The diversity of colonization by Discussion the fungal community in the gut of F. candida was The rapid development of next-generation sequencing strongly affected by dietary structure, but the potential technology and bioinformatics has enabled large-scale, for colonization was much lower for the fungi than the cross-influence analyses of complexity, which has greatly bacteria. enlarged our understanding of the interaction between The gut bacteria of soil invertebrates may thus be suit- microbial communities and their ecological niches [36– able potential indicators of soil pollution. We classified 38]. Interestingly, microbiotas directly and closely associ- the OTUs shared at the phylum and class levels as Pro- ated with personal, public, and planetary health [39] teobacteria and Gammaproteobacteria (Fig. 2E), which have gradually become indicators for predicting human were significantly enriched in the F. candida gut in the and ecological health [17, 18, 37, 38, 40]. Our previous pollutant treatments. Interestingly, the change in relative studies have found that the gut microbiota of soil inver- abundance of Gammaproteobacteria was consistent with tebrates was more sensitive than soil microbiotas to soil the HAA index after exposure to OTC, AZ, and AO, pollution [14, 15, 34]. Identifying common potential in- supporting Gammaproteobacteria as a potential indica- dicator taxa in invertebrate guts may therefore be benefi- tor of the response to soil pollution. We constructed the cial to the assessment of soil ecological risk. We used a co-occurrence network, including both bacteria and microcosm experiment for prediction, and used meta- fungi, to test the contributions of Gammaproteobacteria data and machine learning for verification, to identify in- to the entire gut microbiota and found that Gammapro- dicator taxa from the guts of soil invertebrates in teobacteria correlated with most other bacteria and fungi response to environmental concentrations of soil pollu- to maintain the interoperability and intercommunica- tion and their relationships with the associated gut tions throughout the microbial community (Fig. 2G). microbiome and resistome and with characteristics of F. candida is sensitive to soil pollutants, which will in- host physiology and biochemistry. evitably affect the level of molecular functions. We used We exposed the model soil invertebrate F. candida to a WGCNA analysis of the transcriptomic data to identify common soil pollutants (a fungicide (AZ) and an anti- the functional characteristics of the host under the pres- biotic (OTC)) in our microcosm experiment, separately sure of soil pollution for determining the MEbrown co- targeting fungal and bacterial communities. The survival expressed gene module, significantly correlated with the and fecundity of F. candida were not significantly af- control phenotypes, which differed from the pollution fected, but the significant increase in ROS confirmed phenotypes (Figure S8). Interestingly, the expression of that the host was negatively affected by these pollutants genes significantly correlated with Gammaproteobacteria (Fig. 1D). The enzymatic activity of CYP450 was also was downregulated in the pollutant treatments. These higher in the pollutant treatments, indicating a biochem- genes were assigned to immunological, digestive, meta- ical stress response and detoxification mechanism by F. bolic, and other functional pathways essential to the candida [15]. This result indicated that soil pollutants at host. We used SEM to identify the relationships between environmentally relevant residual concentrations could the host function and Gammaproteobacteria abundance, not be indicated via ordinary physiological phenomena. which again demonstrated that Gammaproteobacteria The HAA index also indicated that the effect of AZ and were negatively linked to host function (transcriptome OTC was higher when combined, than individually (Fig. data) by affecting the diversity of the gut bacterial com- 1C), unlike the effect on the F. candida gut microbiota. munity (Fig. 3G, H). The relative abundance of Gamma- A specific taxon rather than the entire bacterial commu- proteobacteria was significantly and directly correlated nity was thus likely an indicator of the response to soil with the levels of oxidative stress and locomotive ability pollution. of the host. These results strongly supported Gamma- We could not differentiate the diversity and structure proteobacteria as a potential indicator taxon in response of the soil bacterial community between the control and to soil pollution, although limited to the environmental pollutant treatments. The structure of the fungal concentrations in this study. Zhang et al. Microbiome (2021) 9:196 Page 14 of 17 Proteobacteria have been considered as gut microbiota in Orchesella cincta and Folsomia candida, Enterobacte- associated with dysbiosis in humans, blooms of which riaceae, Pseudomonadaceae, and Sphingomonadaceae in have been correlated with obesity, diabetes, and some Caenorhabditis elegans, all of which are classified to the immunological disorders [30]. Proteobacteria are faculta- Gamma- and Alpha-proteobacteria [49], and these taxa tive anaerobic bacteria, whose unique oxygen consump- also strictly conformed to our definition. The definition tion affects the gut environment, making it conducive to of the core taxa of the same genetic relationship was colonization by many types of anaerobic bacteria and de- conserved at the family level, but the cross-species core stroys the stability of the original gut microbial commu- taxa definition in this study can only be conserved at the nity. The phylum Proteobacteria is generally considered class level. The metadata analysis also indicated that gut to be comprised of 116 validated bacterial families and core taxa have similar core taxa classification character- has the largest phylogenetic composition, with highly di- istics, but the specific species composition and structure verse morphological and physiological functions for are significantly different, indicating that the assembly of maintaining a competitive advantage in adapting to the gut microbiotas of soil invertebrates depended on complex and diverse ecological niches (e.g., soil [36], the host and deterministic processes, perhaps associated plants [41], freshwater [42], seawater [43], and the at- with the structure and function of the gut. These results mosphere [44]). In particular, one group of Proteobac- provide a basis for future research on the microbial di- teria, the Gammaproteobacteria, often aggressively versity of soil ecosystems. occupies the ecological niche of the symbiotic relation- We divided a total of 415 samples, from 17 independ- ship between plants and insects [45, 46]. We therefore ent experiments, into two parts, control and pollutant hypothesized that Gammaproteobacteria would occupy treatments. The diversity and structure of the bacterial an important core position in the gut microbiotas of soil community differed significantly between the controls invertebrates. We analyzed 16S rRNA gene sequences and pollutant treatments, and the bacterial communities from the guts of 17 soil invertebrates available in public of the contaminated guts were more stable (Fig. 6), indi- databases but found no consensus among previous re- cating that higher bacterial-community diversity can ports for defining the core microbiota. A current strict maintain the stability of community function. Previous definition for an indicator species is that all samples studies have reported that soil invertebrates will only ac- must include a specific taxon (at the OTU level) [47], cumulate a small amount of pollutants in their guts [14], but this definition is not suitable for many distantly re- leaving the gut community under sub-stressed condi- lated host species [10, 48]. Various independent studies tions, which could reduce competition between taxa and have included different test species, environments, and increase bacterial diversity. The Intermediate Disturb- methods, especially those based on metadata analyses, ance Hypothesis [50] states that appropriate interference which may be less appropriate, because the number of prevents competitive exclusion in communities [51], in- shared OTUs may be much lower than at other levels of dicating that coerced systems are more stable than unco- classification. erced systems. The gut bacterial communities in our We thus took into account both the percent occur- study had a positive defensive strategy under pollution rence in samples and the relative abundance of each stress, and that Gammaproteobacteria played a key role taxon [16] and used class as the core taxonomic level. in maintaining the stability of the communities. Due to We also constructed a shared network based on each in- the complex anaerobic environment of the gut and its dependent experiment, identifying the classes shared high species heterogeneity, finding an effective method among all studies (Fig. 3A), to avoid interference of the for determining the action of Gammaproteobacteria on variation of different independent experiments. Based on the host, such as the commonly used verification the above factors, we developed an equation using CI to methods, isolation of gut microbiota, and artificial con- represent the roles of each taxon and identified the core struction of functional microbiota [52, 53], is difficult. taxa based on the metadata analysis (Eq. 2, Supplemen- The construction of sterile gut environments in our tary Information). Gammaproteobacteria had the highest study has especially identified large obstacles to the veri- CI, supporting their role as an important core taxon in fication of the specific identity of Gammaproteobacteria. the guts of soil invertebrates. Other core gut taxa The combination of metadata analysis and machine (Alphaproteobacteria and Planctomycetes) with a total learning is a good alternative to laboratory-validated relative abundance 50.68% were identified using 0.83 as methods by building models using large-scale computa- the CI’ threshold. A CI’ threshold of 0.3 was used to tional analysis. The RF model we used in this study has identify the transient taxa (35.63%): Bacilli, Ktedonobac- been widely used in the prediction of different commu- teria, Actinobacteria, and Acidobacteria (Fig. 2E). Previ- nities, such as distinguishing between the gut bacterial ous studies also defined the core taxa of common soil communities of mice that feed on normal diets and invertebrates. They included Rickettsia and Pseudomonas those with high salt content [54], distinguishing between Zhang et al. Microbiome (2021) 9:196 Page 15 of 17 soil fungal and bacterial communities with and without analysis; HAA: Index of high area activity; ROS: Reactive oxygen species; fusarium wilt [55] and distinguishing between the root CYP450: Cytochrome P450; OTUs: Operational taxonomic units; WGCNA: Weighted gene co-expression network; KEGG: Kyoto Encyclopedia microbial communities of indica and japonica rice [18]. of Genes and Genomes; PERMANOVA: Permutational analysis of variance; Interestingly, we also found that the healthy and con- RIN: RNA integrity number; TPM: Transcripts per million; SEM: Structural taminated guts were accurately distinguished using an equation modeling; GBIF: Global biodiversity information facility; NorCI: Normalized core index; OOB: Out-of-bag RF model, which can be built at any taxonomic level with similar average accuracy rates, so we selected the RF model at the genus level with the lowest estimated Supplementary Information The online version contains supplementary material available at https://doi. rate of out-of-bag (OOB) errors (17.58%), indicating that org/10.1186/s40168-021-01150-6. Gammaproteobacteria was the core taxon that may re- spond to soil pollutants and was enriched in the con- Additional files 1: Extended Methods. Test soil, species, and taminated guts of the soil invertebrates (Fig. 5E, F). pollutants; Laboratory experimental design; DNA extraction, DNA amplification, library preparation, sequencing and bioinformatic analysis; Most interestingly, previous meta-analyses of gut RNA isolation, transcript sequencing, library preparation and bioinformatic microbiota showed that Gammaproteobacteria (Escheri- analysis; Data collection and description and processing of the 16S rRNA chia coli) may be the main ARBs in the human gut [56]. metadata; Construction and validation Prediction model; Equation for the high area activity (HAA) index and core index (CI). We detected that Gammaproteobacteria was the only Additional files 2: Figure S1. Graphic representation of the experiment bacterial class that was significantly and strongly corre- and analysis method. Figure S2. Number of adults in the control and all lated with the relative abundance of ARGs (Figure S13), treatments. Figure S3. Alpha diversity (Shannon and Chao indexes) and indicating that Gammaproteobacteria may affect the gut beta diversity of the bacterial and fungal communities in the soil surrounding Folsomia candida. Figure S4. Number of significantly resistome of F. candida. The soil invertebrate gut is an different fungal OTUs in all treatments compared to the control group. important reservoir for ARGs, transferring them across Figure S5. Richness (Chao index) of the bacterial and fungal the soil food web, thus allowing its use as one of the in- communities in the gut of Folsomia candida. Figure S6. Alpha diversity (Shannon and Chao indexes) and beta diversity of the fungal community dicators for risk assessment of soil ecology and health in the gut of Folsomia candida. Figure S7. Stability of the networks of [57]. Therefore, Gammaproteobacteria not only respond interaction between bacteria and fungi in all groups. Figure S8. The to the environmental concentration but can also reflect degrees and closeness centrality of bacterial classes (relative abundance) in bacteria-bacteria and bacteria-fungi co-occurrence networks from all la- the dynamic changes among ARGs, which may allow boratory samples. Figure S9. WGCNA analysis of Folsomia candida. Fig- evaluation of soil ecological health risks from two per- ure S10. Core index (CI) and normalized CI (NorCI) of the bacterial spectives. It is therefore a reasonable and flexible strat- communities in the soil surrounding the soil invertebrates across eight in- dependent experiments. Figure S11. Meta-analysis and sensitivity ana- egy to construct Gammaproteobacteria general-purpose lysis of Gammaproteobacteria relative abundance in soil invertebrate primers, to provide an indication of soil ecological and guts. Figure S12. The topological properties (degrees, closeness central- health risks. ity, and betweenness centrality) of bacterial classes in the co-occurrence networks of control and pollution group from all metadata samples. Fig- ure S13. The identification of antibiotic resistance bacteria (ARB) in soil Conclusion invertebrate gut. Figure S14. Annotation of all mapped genes using the GO, KEGG, COG, NR, Swiss-Prot, and Pfam databases. We combined our microcosm experiment, metadata Additional files 3: Data set S1. Uploading sequence data ID. Data set analysis, and a machine-learning method to identify the S2. Antibiotic resistance gene primers. core taxa in the guts of soil invertebrates and provided a quantitative method for identifying core taxa in micro- Acknowledgements bial communities based on metadata analysis that is suit- Not applicable. able for different habitats and species. Interestingly, we also determined that Gammaproteobacteria were a po- Authors’ contributions QZ designed and performed all practical aspects of the study with guidance tential indicator taxon in the guts of the soil inverte- from HFQ. QZ wrote the first draft of the manuscript and HFQ, YGZ, and JP brates that responded to environmental concentrations contributed substantially to revisions. YTY and ZYZ searched all related of soil pollutants, thus providing an effective theoretical literature and collected the publicly accessible databases. QZ, ZYZ, and TL were responsible for the machine learning model construction and related basis for subsequent assessments of soil ecological risk. data analysis. QZ performed the visualization of all data and the artistic Additionally, the results from the physiological and bio- design of all figures. HFQ and TL acquired funding for this project. All chemical analyses of the host, and the microbial- authors read and approved the final manuscript. community functions and antibiotic resistance of Gam- Funding maproteobacteria, also provide new insights for evaluat- Funding was provided by the National Natural Science Foundation of China ing global soil ecological health. (21976161, 21777144, and 41907210). Abbreviations Availability of data and materials ARG: Antibiotic resistance gene; CI: Core index; OECD: Organization for All meta-analysis sequencing data used here were either obtained from pub- Economic Cooperation and Development; ITS: Internal transcribed spacer; licly accessible databases or directly from the authors of the studies, and the RF: Random forest; LR: Logistic regression; SVM: Support-vector machine; Uploading sequence data IDs are shown in Data set S1. The laboratory ex- ROC: Receiver operating characteristic curve; PCoA: Principal coordinate periment sequencing data has been deposited in the National Center for Zhang et al. Microbiome (2021) 9:196 Page 16 of 17 Biotechnology Information Sequence Read Archive (SRA) database (accession 13. Ye X, Xiong K, Liu J. Comparative toxicity and bioaccumulation of number: PRJNA743602). fenvalerate and esfenvalerate to earthworm Eisenia fetida. J Hazard Mater. 2016;310:82–8. https://doi.org/10.1016/j.jhazmat.2016.02.010. Declarations 14. Zhang Q, Zhu D, Ding J, Zheng F, Zhou S, Lu T, et al. The fungicide azoxystrobin perturbs the gut microbiota community and enriches Ethics approval and consent to participate antibiotic resistance genes in Enchytraeus crypticus. Environ Int. 2019;131: Not applicable. 104965. https://doi.org/10.1016/j.envint.2019.104965. 15. Zhang Q, Zhang Z, Zhou S, Jin M, Lu T, Cui L, et al. Macleaya cordata Consent for publication extract, an antibiotic alternative, does not contribute to antibiotic resistance Not applicable. gene dissemination. J Hazard Mater. 2021;412:125272. https://doi.org/10.101 6/j.jhazmat.2021.125272. Competing interests 16. Astudillo-García C, Bell JJ, Webster NS, Glasl B, Jompa J, Montoya JM, et al. The authors declare that they have no competing interests. Evaluating the core microbiota in complex communities: a systematic investigation. Environ Microbiol. 2017;19(4):1450–62. https://doi.org/1 Author details 0.1111/1462-2920.13647. 1 College of Environment, Zhejiang University of Technology, Hangzhou 17. Cammarota G, Ianiro G, Ahern A, Carbone C, Temko A, Claesson MJ, et al. 310032, P. R. China. 2Global Ecology Unit CREAF-CSIC-UAB, CSIC, Bellaterra, Gut microbiome, big data and machine learning to promote precision 08193 Barcelona, Catalonia, Spain. 3CREAF, Campus Universitat Autònoma de medicine for cancer. Nat Rev Gastro Hepat. 2020;17(10):635–48. https://doi. Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain. 4Key org/10.1038/s41575-020-0327-3. Laboratory of Urban Environment and Health, Institute of Urban 18. Yuan J, Wen T, Zhang H, Zhao M, Penton CR, Thomashow LS, et al. Environment, Chinese Academy of Sciences, Xiamen 361021, P. R. China. Predicting disease occurrence with high accuracy based on soil 5 State Key Lab of Urban and Regional Ecology, Research Center for macroecological patterns of Fusarium wilt. ISME J. 2020;14(12):2936–50. Ecoenvironmental Sciences, Chinese Academy of Sciences, Beijing 100085, P. https://doi.org/10.1038/s41396-020-0720-5. R. China. 19. Wright R, Langille M, Walker T. Food or just a free ride? A meta-analysis reveals the global diversity of the plastisphere. ISME J. 2020;15(3):789–806. Received: 11 July 2021 Accepted: 12 August 2021 https://doi.org/10.1038/s41396-020-00814-9. 20. Alavi N, Babaei AA, Shirmardi M. Assessment of oxytetracycline and tetracycline antibiotics in manure samples in different cities of Khuzestan Province, Iran. Environ Sci Pollut R. 2015;22(22):17948–54. https://doi.org/1 References 0.1007/s11356-015-5002-9. 1. Lee W-J, Hase K. Gut microbiota–generated metabolites in animal health and disease. Nat Chem Biol. 2014;10(6):416–24. https://doi.org/10.1038/ 21. Dahshan H, Abd-Elall AMM, Megahed AM. Veterinary antibiotic resistance, nchembio.1535. residues, and ecological risks in environmental samples obtained from 2. Huttenhower C, Gevers D, Knight R, Abubucker S, Badger JH, Chinwalla AT, poultry farms. Egypt Environ Monit Assess. 2015;187(2):2–10. https://doi. et al. Structure, function and diversity of the healthy human microbiome. org/10.1007/s10661-014-4218-3. Nature. 2012;486(7402):207–14. https://doi.org/10.1038/nature11234. 22. Zhu D, An X-L, Chen Q-L, Yang X-R, Christie P, Ke X, et al. Antibiotics disturb 3. Aksoy E, Telleria EL, Echodu R, Wu Y, Okedi LM, Weiss BL, et al. Analysis of the microbiome and increase the incidence of resistance genes in the gut multiple tsetse fly populations in Uganda reveals limited diversity and of a common soil collembolan. Environ Sci Tech. 2018;52(5):3081–90. species-specific gut microbiota. Appl Environ Microb. 2014;80(14):4301–12. https://doi.org/10.1021/acs.est.7b04292. https://doi.org/10.1128/AEM.00079-14. 23. Liaw A, Wiener M. Classification and regression by RandomForest. Forest. 4. De Deyn GB, Raaijmakers CE, Zoomer HR, Berg MP, de Ruiter PC, Verhoef 2001;2(3):18–22. HA, et al. Soil invertebrate fauna enhances grassland succession and 24. Yarnold P, Grimm L. Reading and understanding multivariate statistics; 1995. diversity. Nature. 2003;422(6933):711–3. https://doi.org/10.1038/nature01548. 25. Cortes C, Vapnik V. Support-vector networks. Mach Learn. 1995;20(3):273–97. 5. Berg M, Stenuit B, Ho J, Wang A, Parke C, Knight M, et al. Assembly of the https://doi.org/10.1007/BF00994018. Caenorhabditis elegans gut microbiota from diverse soil microbial 26. Kurogi S, Shigematsu R, Ono K. Properties of direct multi-step ahead environments. ISME J. 2016;10(8):1998–2009. https://doi.org/10.1038/ismej.2 prediction of chaotic time series and out-of-bag estimate for model 015.253. selection. Neural Inf Process. 2014;8835:421–8. https://doi.org/10.1007/ 6. Zhu D, Delgado-Baquerizo M, Su J-Q, Ding J, Li H, Gillings MR, et al. 978-3-319-12640-1_51. Deciphering potential roles of earthworms in mitigation of antibiotic 27. Hernandez DJ, David AS, Menges ES, Searcy CA, Afkhami ME. Environmental resistance in the soils from diverse ecosystems. Environ Sci Technol. 2021; stress destabilizes microbial networks. ISME J. 2021;5(6):1722–34. https://doi. 55(11):7445–55. https://doi.org/10.1021/acs.est.1c00811. org/10.1038/s41396-020-00882-x. 7. Zhu Y-G, Zhao Y, Zhu D, Gillings M, Penuelas J, Ok YS, et al. Soil biota, 28. Lozupone C, Knight R. UniFrac: a new phylogenetic method for comparing antimicrobial resistance and planetary health. Environ Int. 2019;131:105059. microbial communities. Appl Environ Microb. 2005;71(12):8228–35. https:// https://doi.org/10.1016/j.envint.2019.105059. doi.org/10.1128/AEM.71.12.8228-8235.2005. 8. Xiang Q, Zhu D, Chen Q-L, O’Connor P, Yang X-R, Qiao M, et al. Adsorbed 29. Zhu D, Zheng F, Chen Q-L, Yang X-R, Christie P, Ke X, et al. Exposure of a sulfamethoxazole exacerbates the effects of polystyrene (∼2 μm) on gut soil collembolan to Ag nanoparticles and AgNO3 disturbs its associated microbiota and the antibiotic resistome of a soil collembolan. Environ Sci microbiota and lowers the incidence of antibiotic resistance genes in the Technol. 2019;53(21):12823–34. https://doi.org/10.1021/acs.est.9b04795. gut. Environ Sci Tech. 2018;52(21):12748–56. https://doi.org/10.1021/acs.est. 9. Feldhaar H. Bacterial symbionts as mediators of ecologically important traits 8b02825. of insect hosts. Ecol Entomol. 2011;36(5):533–43. https://doi.org/10.1111/j.13 30. Shin N-R, Whon TW, Bae J-W. Proteobacteria: microbial signature of 65-2311.2011.01318.x. dysbiosis in gut microbiota. Trends Biotechnol. 2015;33(9):496–503. https:// 10. Thomas T, Moitinho-Silva L, Lurgi M, Björk JR, Easson C, Astudillo-García C, doi.org/10.1016/j.tibtech.2015.06.011. et al. Diversity, structure and convergent evolution of the global sponge 31. Faddeeva A, Studer RA, Kraaijeveld K, Sie D, Ylstra B, Mariën J, et al. microbiome. Nat Commun. 2016;7(1):11870. https://doi.org/10.1038/ Collembolan transcriptomes highlight molecular evolution of hexapods and ncomms11870. provide clues on the adaptation to terrestrial life. PLoS One. 2015;10(6): 11. Kopittke PM, Menzies NW, Wang P, McKenna BA, Lombi E. Soil and the e0130600. https://doi.org/10.1371/journal.pone.0130600. intensification of agriculture for global food security. Environ Int. 2019;132: 32. Sun M, Chao H, Zheng X, Deng S, Ye M, Hu F. Ecological role of earthworm 105078. https://doi.org/10.1016/j.envint.2019.105078. intestinal bacteria in terrestrial environments: a review. Sci Total Environ. 12. Xing Y, Meng X, Wang L, Zhang J, Wu Z, Gong X, et al. Effects of 2020;740:140008. https://doi.org/10.1016/j.scitotenv.2020.140008. benzotriazole on copper accumulation and toxicity in earthworm (Eisenia 33. Xiang Q, Zhu D, Chen Q-L, Delgado-Baquerizo M, Su JQ, Qiao M, et al. fetida). J Hazard Mater. 2018;351:330–6. https://doi.org/10.1016/j.jhazmat.201 Effects of diet on gut microbiota of soil collembolans. Sci Total Environ. 8.03.019. 2019;676:197–205. https://doi.org/10.1016/j.scitotenv.2019.04.104. Zhang et al. Microbiome (2021) 9:196 Page 17 of 17 34. Zhang Q, Zhu D, Ding J, Zhou S, Sun L, Qian H. Species-specific response of 53. Morella NM, Weng FC-H, Joubert PM, Metcalf CJE, Lindow S, Koskella B. the soil collembolan gut microbiome and resistome to soil oxytetracycline Successive passaging of a plant-associated microbiome reveals robust pollution. Sci Total Environ. 2019;668:1183–90. https://doi.org/10.1016/j. habitat and host genotype-dependent selection. P Natl Acad Sci. 2020; scitotenv.2019.03.091. 117(2):1148–59. https://doi.org/10.1073/pnas.1908600116. 35. Zhu D, Chen Q-L, Li H, Yang X-R, Christie P, Ke X, et al. Land use influences 54. Wilck N, Matus MG, Kearney SM, Olesen SW, Forslund K, Bartolomaeus H, antibiotic resistance in the microbiome of soil collembolans Orchesellides et al. Salt-responsive gut commensal modulates TH17 axis and disease. sinensis. Environ Sci Tech. 2018;52(24):14088–98. https://doi.org/10.1021/acs. Nature. 2017;551(7682):585–9. https://doi.org/10.1038/nature24628. est.8b05116. 55. Zhang J, Liu YX, Zhang N, Hu B, Jin T, Xu H, et al. NRT1.1B is associated with 36. Lauber CL, Hamady M, Knight R, Fierer N. Pyrosequencing-based assessment root microbiota composition and nitrogen use in field-grown rice. Nat of soil pH as a predictor of soil bacterial community structure at the biotechnol. 2019;37(6):676–84. https://doi.org/10.1038/s41587-019-0104-4. continental scale. Appl Environ Microb. 2009;75(15):5111–20. https://doi. 56. Li X, Stokholm J, Brejnrod A, Vestergaard GA, Russel J, Trivedi U, et al. The org/10.1128/AEM.00335-09. infant gut resistome associates with E. coli, environmental exposures, gut 37. Hermans SM, Buckley HL, Case BS, Curran-Cournane F, Taylor M, Lear G. microbiome maturity, and asthma-associated bacterial composition. Cell Using soil bacterial communities to predict physico-chemical variables and Host Microbe. 2021;29(6):975–87. https://doi.org/10.1016/j.chom.2021.03.017. soil quality. Microbiome. 2020;8(1):79. https://doi.org/10.1186/s40168-020- 57. Zhu D, Xiang Q, Yang X-R, Ke X, O’Connor P, Zhu Y-G. Trophic transfer of 00858-1. antibiotic resistance genes in a soil detritus food chain. Environ Sci Technol. 38. Lurgi M, Thomas T, Wemheuer B, Webster NS, Montoya JM. Modularity and 2019;53(13):7770–81. https://doi.org/10.1021/acs.est.9b00214. predicted functions of the global sponge-microbiome network. Nat Commun. 2019;10(1):992. https://doi.org/10.1038/s41467-019-08925-4. 39. Prescott SL, Wegienka G, Logan AC, Katz DL. Dysbiotic drift and Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in biopsychosocial medicine: how the microbiome links personal, public and published maps and institutional affiliations. planetary health. Biopsychosocial Med. 2018;12(1):7. https://doi.org/10.1186/ s13030-018-0126-z. 40. Astudillo-García C, Hermans SM, Stevenson B, Buckley HL, Lear G. Microbial assemblages and bioindicators as proxies for ecosystem health status: potential and limitations. Appl Microbiol Biot. 2019;103(16):6407–21. https:// doi.org/10.1007/s00253-019-09963-0. 41. Redford AJ, Fierer N. Bacterial succession on the leaf surface: a novel system for studying successional dynamics. Microb Ecol. 2009;58(1):189–98. https:// doi.org/10.1007/s00248-009-9495-y. 42. Zhang Q, Zhang Z, Lu T, Peijnenburg WJGM, Gillings M, Yang X, et al. Cyanobacterial blooms contribute to the diversity of antibiotic-resistance genes in aquatic ecosystems. Comms Bio. 2020;3(1):737. https://doi.org/10.1 038/s42003-020-01468-1. 43. Hanno Teeling BMF. Dörte Becher. Substrate-controlled succession of marine bacterioplankton populations induced by a phytoplankton bloom. Science. 2012;336(6081):608–11. https://doi.org/10.1126/science.1218344. 44. Whon TW, Kim MS, Roh SW, Shin NR, Lee HW, Bae JW. Metagenomic characterization of airborne viral DNA diversity in the near-surface atmosphere. J Virol. 2012;86(15):8221–31. https://doi.org/10.1128/JVI.002 93-12. 45. Gottel NR, Castro HF, Kerley M, Yang Z, Pelletier DA, Podar M, et al. Populus deltoides roots harbor distinct microbial communities within the endosphere and rhizosphere across contrasting soil types. Appl Environ Microb. 2011;77(17):5934–44. https://doi.org/10.1128/AEM.05255-11. 46. Gatehouse LN, Sutherland P, Forgie SA, Kaji R, Christeller JT. Molecular and histological characterization of primary (Betaproteobacteria) and secondary (Gammaproteobacteria) endosymbionts of three mealybug species. Appl Environ Microb. 2012;78(4):1187–97. https://doi.org/10.112 8/AEM.06340-11. 47. Segata N, Baldini F, Pompon J, Garrett WS, Truong DT, Dabiré RK, et al. The reproductive tracts of two malaria vectors are populated by a core microbiome and by gender- and swarm-enriched microbial biomarkers. Sci Rep. 2016;6(1):24207. https://doi.org/10.1038/srep24207. 48. Schmitt S, Tsai P, Bell J, Fromont J, Ilan M, Lindquist N, et al. Assessing the complex sponge microbiota: core, variable and species-specific bacterial communities in marine sponges. ISME J. 2012;6(3):564–76. https://doi.org/1 0.1038/ismej.2011.116. 49. Bahrndorff S, de Jonge N, Hansen JK, Lauritzen JMS, Spanggaard LH, Sørensen MH, et al. Diversity and metabolic potential of the microbiota associated with a soil arthropod. Sci Rep. 2018;8(1):2491. https://doi.org/10.1 038/s41598-018-20967-0. 50. Joseph H. Connell. Diversity in tropical rain forests and Coral Reefs. Science. 1978;199(4335):1302–10. https://doi.org/10.1126/science.199.4335.1302. 51. Marull J, Tello E, Fullana N, Murray I, Jover G, Font C, et al. Long-term bio- cultural heritage: exploring the intermediate disturbance hypothesis in agro-ecological landscapes (Mallorca, c. 1850–2012). Biodivers Conserv. 2015;24(13):3217–51. https://doi.org/10.1007/s10531-015-0955-z. 52. Motta EVS, Raymann K, Moran NA. Glyphosate perturbs the gut microbiota of honey bees. P Natl Acad Sci. 2018;115(41):10305–10. https://doi.org/10.1 073/pnas.1803880115.