Sample collection
Earthworms and surrounding soil samples were collected from the top 20 cm of soil based on the depth of the cultivated horizon from 28 provinces across China, including 35 arable lands and 16 forested lands. Samples were collected from June 15 to August 20, 2017. The latitude and longitude of sampling sites are supplied in Additional file 1: Table S1, and shown in Additional file 7: Fig. S1. Except for forested land samples from Chongqi (two replicates), each collection comprised three replicates of earthworm and soil samples, thus a total of 152 earthworm samples and 152 surrounding soil samples were obtained. For collection of earthworm samples, we firstly excavated multiple rectangular block of soils (length: 50 cm, width: 30 cm and high: 20 cm) at each location, and placed them onto the clean white cotton. The soil block was broken up and all earthworms carefully collected by hand. Because different soils have different earthworm species it is not possible to find a dominant species in all soils. We collected 3–30 individuals per species in a replicate, 60–240 individuals per site and 20–80 individuals per replicate. Using morphological analysis, we identified the dominant species based on abundance. Soil samples were collected using a five-point sampling method. The collected soils and earthworms were divided into two parts. One part was used to extract DNA, and the other was used to analyze chemical properties. The samples for the chemical analysis were preserved in sealed plastic bags and transported to the laboratory at 4 oC within 24 h. Samples for DNA extraction (earthworms and soils) were immediately stored in anhydrous alcohol and kept at -20 oC in the laboratory until use.
DNA extraction and high-throughput quantitative PCR (HT-qPCR)
Earthworm samples stored in anhydrous alcohol were carefully washed twice using sterile water to remove allogenic materials from the surface, rinsed for 2 min with 0.5% sodium hypochlorite solution to sterilize the surface, and rinsed three times with sterile water at 4 °C. Earthworms were immediately dissected using sterile forceps and surgical scissors in a bacteria-free operating environment to obtain earthworm gut contents. Gut contents (0.5 g) were collected into a 2-mL Eppendorf tube, to which 0.96 mL sterile 100 mM phosphate buffer pH 7.4 was added. The FastDNA Spin Kit for soil (MP Biomedicals, Illkirch, France) was used to isolate DNA from earthworm gut contents or from 0.5 g soils. DNA extraction followed the manufacturer’s instructions. The 100 µL DES (DNA elution solution) was used to elute the extracted DNA into a sterile 1.5-mL Eppendorf tube. Since gut contents often contain larger amounts of PCR inhibitors than soil, we used a universal DNA purification kit (HYK-002, Lulong Biotech, China) to purify the obtained DNA samples (earthworms and soils) before further analysis. The quality and concentration of DNA were checked using 1.0% agarose gel electrophoresis and the Quibt 3.0 Fluorometer (Invitrogen). Resulting DNA was stored at − 20 °C until further study.
Before HT-qPCR, DNA was diluted to 30 ng µL− 1 using sterile water. Although the extracted DNA was purified, to further confirm that inhibition was negligible, we conducted a spiking test with a dilution series of the template DNA from soils and earthworm gut contents to determine recovery ratios of PCR. In the spiking test, we added three known amounts of ARGs (tetC, tetG and sul2) into soil and gut DNA extracts. At the same time, we diluted a series of DNA templates (0, 10, 100 and 1000 fold) from soils and earthworm gut contents, including three different earthworm ecotypes. The spiked and diluted DNAs were detected using HT-qPCR. The abundance of added ARGs did not present significant difference between the spiked soils and gut DNA extracts (P > 0.05). No significant difference was observed in ratios of PCR from the dilution series of the template DNA between soils and gut contents (P > 0.05). These results indicated that the PCR inhibitors were negligible in extracted soils and gut contents.
A total of 285 primer sets (Additional file 2: Table S2) targeting resistance genes to major classes of antibiotics, 16S rRNA genes and 10 mobile genetic elements (MGEs) were used. The Wafergen SmartChip Real-time PCR system (Wafergen Inc., USA) was used to perform the HT-qPCR. Three technical repeats and a non-template negative control were included in the amplification of each primer set. The HT-qPCR amplification conditions have been previously described [28, 29]. We used SmartChip qPCR software (V 2.7.0.1) to analyze the HT-qPCR data. Data from amplifications with efficiencies outside the range of 0.9–1.1 were discarded. Only data with three positive replicates were considered as a successful amplification and used in the downstream analysis. The threshold cycle number (CT) for the detection limit was set at 31. We normalized the copy number of ARGs per 16S rRNA gene or bacterial cell to reflect the abundance of ARGs, avoiding biases arising from different bacterial abundance in each sample. Normalization was carried out using the bacterial 16S rRNA gene abundance in each sample. Calculations followed the equations below [28]:
Relative gene copy number = 10^((31 − CT(detection))/(10/3))
Relative abundance of ARGs = Relative ARG copy number/Relative 16S rRNA gene copy number
Normalized ARG copy number per bacterial cell = Relative abundance of ARGs × 4.1.
Based on the Ribosomal RNA Operon Copy Number Database, on average, a bacterial cell contains 4.1 copies of 16S rRNA gene, and this correction factor was used in this study [28, 29].
Real-time quantitative PCR (RT-qPCR)
To explore the relationship between bacterial abundance and ARG abundance, we used RT-qPCR with the Roche 480 system to determine the absolute bacterial 16S rRNA copy numbers of earthworm gut content and soil. Universal primers (5 ‘GGGTTGCGCTCGTTGC) and (5’ ATGGYTGTCGTCAGCTCGTG) were used to target the bacterial 16S rRNA gene region. The qPCR system and conditions were as previous studies [28, 30]. Each sample was replicated three times in all qPCR, including a non-template reaction as the negative control. The standard curve was obtained by detecting triplicate 10-fold serial dilutions of plasmid DNA with the target-gene (amplification efficiency 97 − 103%, r2 > 0.99).
Bacterial 16S rRNA gene amplification, sequencing and bioinformatic analysis
To examine the biogeographic characteristic of earthworm gut microbial communities and explore the contribution of bacterial communities to the change in ARGs, we characterized the bacterial communities of earthworm gut and soil using a set of primers (515F: GTGCCAGCMGCCGCGG − 907R: CCGTCAATTCMTTTRAGTTT) to amplify the 16S rRNA gene hypervariable V4-V5 region [28, 31]. Unique barcodes were included in each primer to distinguish different samples. The DNA was amplified and purified according to previous studies [32, 33]. We used a Qubit 3.0 fluorimeter to determine the concentration of purified products, and 24 products of equal concentration were pooled as a library. Finally, the Illumina Hiseq2500 platform (Novogene, Tianjin, China) was used for high-throughput sequencing of each library.
Quantitative Insights Into Microbial Ecology (QIIME v1.9.1) was used to analyze high-throughput sequencing data of earthworm gut and soil, following the online instructions [34]. Briefly, low-quality reads, ambiguous nucleotides and primer sequences were removed to obtain clean sequences. We assigned high quality sequences to operational taxonomic units (OTUs) according to 97% sequence similarity using UCLUST [35]. We discarded singleton OTUs to generate the final OTU table, and assigned taxonomy to OTUs based on the representative sequence using the RDP Classifier 2.2 [36]. The PyNAST aligner was used to align representative OTUs based on the GreenGenes 13.8 reference database [37]. We built a phylogenetic tree via the FastTree algorithm [38]. To minimize the error of potential sequencing, we discarded those OTUs with a total number of reads of < 100 in all samples before downstream analysis [39].
Sample analysis
Soils and earthworm gut contents were air dried in the dark and 2 mm sieved. We used soil/water (1/2.5) and gut content/water (1/5) suspensions to measure pH, and a C/N instrument (Vario MAX, C/N, Elementar) to determine total nitrogen and carbon. The clay content of soil was determined using a particle number analyzer (Malvern Instruments Ltd.). Soil types and properties were summarized in Additional file 3: Table S3.
Identification of earthworm species
After earthworm species were identified by morphological characteristics, we further used Cytochrome Oxidase I (COI) barcode gene sequencing as confirmation. DNA of earthworm tissues was isolated using the DNeasy Blood and Tissue Kit (QIAGEN, China (Shanghai) Co., Ltd) based on the manufacturer’s instructions. The isolated DNA was amplified using the primers LCO1490 (5’ GGTCAACAAATCATAAAGATATTGG) and HCO2198 (5’ TAAACTTCAGGGTGACCAAAAAATCA) to target the COI barcode region [5]. The amplification system and conditions were as previously designed [40]. PCR products were sequenced, and the sequences submitted to NCBI to identify the species of earthworms via the Basic Local Alignment Search Tool (BLAST). Information on earthworm species and ecological group at each location is listed in Additional file 4: Table S4.
The microcosm assay and long-term field experiments
Earthworms comprise abundant species (~ 8000) and the distribution of species varies greatly in different regions. According to the distinct feeding and burrowing habits of earthworms, they are usually classified into three ecological groups (epigeic, endogeic and anecic) [14]. Epigeic earthworms mainly ingest plant litter and live on the surface, and rarely burrow. Endogeic earthworms feed on a large number of soil food sources (e.g. mineral soil and humus), live below the surface, and commonly burrow. Anecic earthworms usually ingest organic residues, forage the surface, live in the mineral soil layer, and burrow. Earthworms of different ecotypes often have different feeding guilds, affecting the roles of earthworms in soil ecosystems. To reveal the change of ARGs among soils, earthworm gut contents and casts, we conducted a microcosm experiment using three different ecological earthworm species (epigeic: Eisenia andrei, endogeic: Amynthas robustus and anecic: Amynthas hupeiensis). The farmland soils (29o47’ N, 121o21’ E) from Ningbo were used in the microcosm experiment, and soil properties are summarized below. Three treatments were performed, including three kinds of different earthworm species representing different ecotypes. Each treatment had 15 earthworms added into 2.5 kg dried soil (75% of the field water-holding capacity) and was replicated four times. The microcosms were incubated in the greenhouse for 28 days. After 28 days, we collected soils, earthworms and casts samples, and earthworms were dissected to obtain gut contents. The obtained samples were used to detect the abundance and diversity of ARGs.
To explore the effect of earthworms on the ARGs in the soil ecosystem, we conducted the microcosm experiment with a full factorial design, including two soil conditions (control and amended with pig manure) and earthworm treatments (control and adding 15 earthworms; Amynthas sp.). Each treatment was replicated four times, and a total of 16 pots were used. Earthworm density was selected based on the following considerations. The typical earthworm density in the agricultural soil is 1–5 earthworms in 2.5 kg soil, but with the input of organic matter the number of earthworms will also increase. 15 earthworms in 2.5 kg soil can be also found in the soil, especially, which was applied organic fertilizer (e.g. pig manure). The soil (clay content = 7.35%, pH (CaCl2) = 4.75, CEC = 13.76 cmol kg− 1, total C 32.4 g kg− 1 and total N 3.77 g kg− 1) was collected from arable land (depth: 0–20 cm) at Ningbo (29o47’ N, 121o21’ E), Zhejiang province, China. Each pot contained 2.5 kg dried soil, and water content of the soil was adjusted to 75% of the field water-holding capacity. After incubating for two weeks, earthworms of similar size (~ 60 g fresh weight) and seeds of a common crop lettuce (Lactuca sativa) were added. The microcosms were incubated in the greenhouse for 75 days. We further collected soil samples from two long-term field experiments with or without earthworms, in the corn season. M. guillelmi was the dominant species of earthworm in natural fields. The two long-term field experiments were described as below: 1) Wheat-corn rotation was used. The field experiment lasted for one year. Each treatment was replicated three times. 2) Corn-vetch rotation was conducted lasting for two years. Each treatment was replicated five times. After harvest, the DNA was extracted from soil and earthworms as above, and the ARGs and bacterial communities of soil and earthworms were analyzed as described above.
Data processing and statistical analysis
The figure describing the sampling sites was produced via ArcMap 10.1. Data on ARGs and bacterial community were presented as mean ± SE. We set the significance level at 0.05 and indicated P values as follows: *P < 0.05, **P < 0.01 and ***P < 0.001. SPSS version 20.0 was selected to compare the significance of the data using the linear mixed model, paired-samples t-test and independent-sample t-test. We used the principal coordinate analysis (PCoA) to reveal patterns of earthworm ARG profiles via RStudio with labdsv 1.8-0 and ggplot 2. The differences of shared ARGs and dominant bacterial OTUs (with relative read abundance > 0.05) between soil and earthworms were compared using STAMP software. The violin plots were conducted using the ggplot2 packages of R with version of 3.4.4, and the Venn diagram was produced via the online version of Venny 2.1. We used the vegan 2.3-1 packages of R to perform a Nonmetric multidimensional scaling analysis depicting the difference between earthworm gut and soil bacterial communities [41, 42]. Partial redundancy analysis was conducted via the vegan 2.3-1 package of R to reveal the contribution of earthworm gut bacterial communities, soil properties (bacterial community and physical and chemical properties) and gut content properties (pH and total C) to the percentage of variation in ARG profiles. We used the vegan 2.3-1 package of R to conduct the PERMANOVA (Adonis test and Anosim test), the Mantel test and the Procrustes test and to calculate the Shannon index of bacterial communities. OriginPro 9.1 was used to produce the histogram and Scatter diagram, and to perform linear regression.