Amynthas corticis is an epi-endogeic earthworm species with global distribution (Cunha et al., 2017). The earthworms were collected using an electroshocking method from grassland in Beijing (39°57′N, 116°17′E). Specifically, four soil probes were placed into the soil in a square plot with a resulting sampling area of 1 m². The probes were part of an electroshocking protocol that was used to stimulate earthworms emerging within the sampling area (Weyers et al., 2008). Electroshocking has been shown to have minimal effects on soil biota (Staddon et al., 2003). Species identification of the earthworm individuals was based on morphological characteristics such as body color, first dorsal pore, clitellum position, and the number of body segments (Zhao et al., 2015). For those individuals that were identified as A. corticis, only intact and mature adult earthworms with similar body size (about 3-g fresh weight) were retained for use in the experimental setup in order to circumvent the effects of minor factors (e.g., integrity and age) on results. Target individuals were brought to the laboratory, and non-targets were returned to the field. Then the target earthworms were surface-washed with sterile water and euthanatized by freezing at −40°C for 10 min (Zhang et al., 2010).

To clarify whether the functioning microbes in decomposition were primarily from the indigenous microbes of earthworm or microbes dwelling in the soil, two soils were collected from geographically distant sites: Beijing (39°57′N, 116°17′E), where the earthworms were originally collected from, and Xuchang (34°02′N, 113°81′E), which is a non-native environment for the collected earthworms. The soils collected from Beijing and Xuchang were denoted as native soil and non-native soil, respectively. Both soils were classified as fluvo-aquic (Zou et al., 2015; Zhao et al., 2019) but with different bacterial community compositions due to the geographical separation (Supplementary Table 1 and Supplementary Figure 1). Theoretically, if the microbes taking part in the decomposition are derived from soil, this will result in different decomposition dynamics and decomposer communities for the native soil and non-native soil treatments; however, if indigenous earthworm microbes are the primary decomposers, then similar decomposition dynamics and decomposer communities will be observed. At each sampling site, the soil was collected from the top 20 cm of the field in five 1 m × 1 m plots. In each plot, 20 soil cores (5-cm diameter) were randomly collected. All collected soil samples were thoroughly homogenized and immediately transported to the laboratory on ice. The samples were then stored at 4°C for less than 2 weeks before use.

For each microcosm, sieved (mesh size 2 mm) and homogenized soil (120-g dry weight equivalent) was placed in polypropylene bottles (12.5-cm height, 6-cm diameter). Three euthanatized earthworms were put in a nylon mesh bag (mesh size < 1 mm), and their weights were recorded. The mesh bags allowed for efficient sampling and nutrient exchange. Then the decomposition bag containing dead earthworms was buried in the middle layer of soil (4 cm from the bottom), mimicking a realistic situation. Soil moisture was adjusted to 40% of water-holding capacity by adding sterile deionized water. All microcosms were incubated at 20°C in the dark in a climate chamber. Based on a preliminary experiment, approximately 80% of earthworm mass is lost in native soil due to decomposition within 10 days. Therefore, to ensure sufficient samples for subsequent analysis, the incubation period was set as 8 days. To explore earthworm decomposition dynamics and bacterial community shifts, five replicates were destructively sampled on days 0, 1, 3, 5, and 8 for the following DNA extraction and DOM analyses. In total, 50 decomposition bags (2 soil types × 5 time points × 5 replicates = 50) and 50 soil samples (collected underneath the decomposing earthworm) were collected.

The mass loss was determined from the difference of dry weight between the initial and final mass of each sampling time. At each sampling time, decomposition bags were removed from the microcosms and were freeze-dried at −40°C for 24 h to a constant weight using a vacuum freeze dryer (Labconco, Kansas City, MO, United States). After drying, the undigested soil in earthworm gut was carefully removed from earthworm tissues using tweezers. Then the tissue samples were weighed and manually ground to fine powder with a mortar. The freeze-dried samples were also used for subsequent DOM analysis and DNA extraction. Because fresh earthworms were used for the decomposition experiment, the initial tissue dry mass (three euthanatized earthworms) in each microcosm was estimated based on the ratio of tissue dry mass to total wet mass (including tissue dry mass, dry mass of undigested soil in gut, and water). To determine the ratio, three replicates of fresh earthworms collected from the same batch were freeze-dried, and the tissue dry mass of these three earthworms was determined following the same procedure described above to calculate the proportion of tissue dry mass to total wet mass.

Freeze-dried earthworm powder (0.02 g) was dissolved in ultrapure water with a ratio of 1:10. The suspensions were centrifuged at 8,000 g and 4°C for 10 min and then filtered through a 0.45-μm cellulose acetate membrane. The collected filtrate was diluted 20-fold with ultrapure water to an appropriate concentration (fluorescence intensity < 10,000) before excitation-emission matrix (EEM) spectroscopy analysis. The EEM analysis is a widely used tool to characterize the composition of complex DOM mixture (Chen et al., 2003), since different DOM components emit different fluorescence intensities with specific emission wavelengths when excited by specific excitation wavelengths. EEM fluorescence spectra were measured on an F-4500 fluorescence spectrophotometer (Hitachi, Tokyo, Japan) with a xenon excitation source. EEM scans were made at the excitation range from 200 to 400 nm and an emission range from 280 to 500 nm at 5-nm increments. The scanning rate was set at 12,000 nm min^–1. The spectrum of ultrapure water was recorded as the blank (Zhang et al., 2016). Raman scattering, which can result in distortion of the measured spectrum due to energy change, was removed by subtracting the fluorescence spectra of the blank. EEM fluorescence spectra were classified into five fractions to represent five major types of DOM (Chen et al., 2003): tyrosine-like proteins (region 1; excitation, 200–250 nm; emission, 280–330 nm), tryptophan-like proteins (regions 2; excitation, 200–250 nm; emission, 330–380 nm), fulvic acid-like organics (region 3; excitation, 200–250 nm; emission, 380–550 nm), microbial byproduct-like materials (region 4; excitation, 250–400 nm; emission, 280–380 nm), and humic acid-like organics (region 5; excitation, 250–400 nm; emission, 380–550 nm).

Total DNA was extracted from 0.2 g of earthworm powder samples and soil samples from each of the replicates, for each treatment and sampling time, using the QIAGEN DNeasy PowerSoil Kit (QIAGEN, Hilden, Germany). The quality of DNA was assessed using the ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, United States). Primer pairs 338F (5′-ACT CCT ACG GGA GGC AGC AG-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′), targeting the V3–V4 region of 16S RNA genes, were used for amplification (Tian et al., 2018). An additional 6-base barcode was added to 806R (Sampson et al., 2016). The barcode was unique to each sample and permitted the identification of individual samples within a single MiSeq sequencing run. We performed three PCR technical replicates per sample on each DNA extract. The 20-μl PCRs were conducted using 0.5 μl of each primer (10 mM), 10 μl of Premix Taq (1.25 U of DNA polymerase, 0.4 mM of dNTPs, 3 mM of Mg²⁺; Takara, Dalian, China), 1 μl of template DNA (10 ng), and 8 μl of nuclease-free PCR-grade water. The PCR was performed on a T100 Thermal Cycler (Bio-Rad, Hercules, CA, United States). The thermal cycling conditions were as follows: 94°C for 5 min, followed by 28 cycles of 94°C for 30 s, 55°C for 30 s, 72°C for 45 s, and a final extension at 72°C for 10 min. Negative controls, other ingredients not changed except the DNA templates replaced with sterilized water, were included for detecting potential contamination (Sun et al., 2019). PCR products were pooled and purified with the AxyPrep DNA Gel Extraction Kit (Axygen, Hangzhou, China). The barcode-tagged amplicons were mixed in equimolar concentrations for MiSeq library construction and then subjected to 250 base pair (bp) paired-end sequencing on the MiSeq platform (Illumina, San Diego, CA, United States) in a commercial laboratory (Majorbio, Beijing, China).

The resulting sequence reads were first sorted based on orientation, and the barcodes were extracted by using extract.barcodes.py in QIIME (Caporaso et al., 2010). Then the sequences were demultiplexed using the q2-demux plugin in QIIME 1.91 (Caporaso et al., 2010). PCR primers were trimmed using Cutadapt allowing for one mismatch (Martin, 2011). The trimmed reads were merged using the FLASH 2 (minimum overlap, 10; max overlap, 65; max mismatch density, 0.25), and quality filtering was performed with the fastq_filter command (Magoc and Salzberg, 2011). Paired reads were discarded if the average number of base error in the read (expected errors) was greater than 2. PCR chimeras were filtered out using USEARCH (Edgar, 2010). Then, identical sequences were dereplicated, and the singletons with possible erroneous sequences were removed using USEARCH (Edgar, 2010). The remaining high-quality sequences were clustered into operational taxonomic units (OTUs) using UCLUST at a 97% sequence similarity cutoff (Shi et al., 2019). The OTUs were assigned taxonomy in QIIME using the SILVA ribosomal RNA gene database project (Quast et al., 2012). All non-bacterial sequences were removed.

We estimated the temporal dynamics of mass loss and the bacterial alpha diversity of decomposing earthworms using linear regression analysis. The slopes of linear fits were compared by bootstrapping followed by t-test (Efron and Tibshirani, 1986; Zhou et al., 2008). The differences for mass loss, bacterial richness, Simpson index, and the amount of DOM at specific sampling time in native soil and non-native soil were compared using independent-samples t-test. To understand how the DOM composition changed during decomposition, principal coordinate analysis (PCoA) was performed based on the Bray–Curtis distance of the DOM fluorescence spectra. Also, permutational multivariate analysis of variance (PERMANOVA) was conducted to test the significant distinctions in the DOM constituent across different sampling times.

Principal coordinate analysis based on the Bray–Curtis matrix was performed to explore changes in the bacterial community composition, and the shift patterns were confirmed by PERMANOVA. To examine how the identity of bacterial taxa changed during earthworm decomposition, the shared genera between different sampling times and the unique genera at a specific sampling time were visualized via a network graph using Gephi v0.9.2 (Bastian et al., 2009). To further examine how the abundance of specific bacterial genera changed during earthworm decomposition, we used linear regression to test for the significant trends over time. Fast expectation-maximization microbial source tracking (FEAST), a source tracking method, was used to quantify the contribution of the likely source environments (source) to the microbial composition on decaying earthworm (Shenhav et al., 2019). Mantel test, based on Pearson’s correlations between the two distance matrices (e.g., Bray–Curtis distance of DOM composition), was performed to test whether the decomposition process depends on shifts in the microbial community (Ge et al., 2016).

Potential functions for bacterial taxa were evaluated with the program “functional annotation of prokaryotic taxa” (FAPROTAX). The FAPROTAX defines the metabolic functions in terms of prokaryotic taxa (e.g., genera or species) (Louca et al., 2016). Currently, 90 functional groups covering the process of carbon, nitrogen, and hydrogen cycling (e.g., methylotrophy, nitrate respiration, and hydrocarbon degradation) are included in the database. To identify the specific functional group that regulated the decomposition process, a correlation matrix was constructed by calculating Pearson’s correlation coefficients between the abundance changes for a bacterial functional group and the changes in DOM fluorescence intensity. The changes of bacterial functional group abundance and DOM were calculated by subtracting their values at day 0 from the values at each sampling time. When the abundance of functional group changed in concert with the changes in DOM amount, positive correlations would be detected. In contrast, negative correlations would be observed when the abundance of functional group and the amount of DOM changed in opposite ways. The Benjamini–Hochberg method was used to adjust the original p-values of correlation analysis to control the false discovery rate caused by multiple testing. PCoA, PERMANOVA, and Mantel test were conducted using the “vegan” package in R platform².

After 8 days of incubation, a decrease in the total mass of approximately 60% was detected (Figure 1A, 63% in native soil; Figure 1B, 55% in non-native soil). There were no significant differences in mass loss between native soil and non-native soil (p = 0.33). At each sampling time, the mass loss was also unaffected by soil origin (Supplementary Table 2). Decomposition rates, which were estimated as the slope of a linear regression of mass loss versus time, showed similar temporal trends in native soil (Figure 1A) and non-native soil (Figure 1B). In native soil, the decomposition rate before day 3 was 15.77 and then decreased to 3.25 (Figure 1A). In non-native soil, the decomposition rate before day 3 was 13.06 and then decreased to 2.80 (Figure 1B). A significant increase in decomposition rate was detected before day 3 (p < 0.01 for both soils); then the decomposition rate significantly decreased after day 3 (Figure 1A, p = 0.03; Figure 1B, p = 0.08). However, there was no significant difference in decomposition rate for soil either before day 3 (p = 0.19) or after day 3 (p = 0.40).

Throughout the study period, the total amount of earthworm DOM showed a clear temporal pattern in both native and non-native soils, with a distinct increase before day 3 (p < 0.01) and a subsequent decrease after day 3 (Figures 2A,B, p < 0.01). More specifically, in native soil, the fluorescent intensities of all five DOM types increased during the first 3 days (p < 0.01), with a peak of ∼2.2 × 10⁵ arbitrary units on day 3 and then a decrease to ∼1.3 × 10⁵ arbitrary units by day 8 (Figure 2D, p < 0.01). In non-native soil, the fluorescent intensities increased to a peak of ∼2.0 × 10⁵ arbitrary units during the first 3 days (p < 0.01) and then decreased to ∼1.1 × 10⁵ arbitrary units over time (Figure 2E, p < 0.01). Both in native soil and non-native soil, the extracted DOM was dominated by tyrosine-like proteins (region 1, about 45%), tryptophan-like proteins (region 2, about 35%), and microbial byproduct-like materials (region 4, about 18%, Figure 2 and Supplementary Table 3). The amount of total DOM and specific DOM types was not affected by soil origin for most sampling times, except day 3 (Supplementary Table 4, p > 0.05 for all except the tyrosine-like proteins, tryptophan-like proteins, and microbial byproduct-like materials at day 3).

The shifts in DOM composition over time were also visualized in the PCoA based on the Bray–Curtis distance, and the significance of differences was confirmed by PERMANOVA. The decomposition time had a significant effect on DOM composition in both soils (p < 0.01), as shown by the progressive changes of DOM composition over time in native soil (Figure 2G) and non-native soil (Figure 2H). Similar to the results of DOM amount, DOM composition was not affected by soil origin for most sampling times, except day 3 (Figure 2F and Supplementary Table 5, p > 0.05).

Bacterial community richness and Simpson index showed higher values before day 3 and then decreased thereafter (Figure 3). In native soil, the richness increased initially during the first 3 days (p = 0.02) and then decreased (p < 0.01, Figure 3A). In non-native soil, the richness did not change initially (p = 0.91) and then decreased during later stages of decomposition (p < 0.01, Figure 3B). The Simpson index is another characteristic of diversity, with a value range between 0 and 1. A value closer to zero indicates that one or a few species dominate, while a value closer to 1 indicates that the number of each species in the community is relatively equal and that more species are present in the sample. The Simpson index showed no significant change with time before day 3 in either native (p = 0.71, Figure 3C) or non-native (p = 0.19, Figure 3D) soils, whereas this index declined significantly after day 3 (p < 0.01 for both soils). The temporal changes of bacterial richness were marginally different in both soils either before day 3 (p = 0.06) or after day 3 (p = 0.07, Figure 3E). The Simpson index displayed a significant difference between the soils only after day 3 (p = 0.03, Figure 3F). For each sampling time, the bacterial richness was significantly affected by soil origin except on day 1 (Supplementary Table 6, p < 0.05 for all except day 1). The Simpson index was significantly affected by soil origin only on days 5 and 8 (Supplementary Table 6).

The succession of the bacterial community on the decaying earthworms was also clearly detected, as indicated by PCoA performed on the entire dataset (Figures 4A–C). The communities of the first 3 days were clustered together and well separated from the bacterial communities of days 5 and 8 in both native soil and non-native soil (Figures 4B,C). Decomposition time had a significant effect on bacterial community composition in the native soil treatment (Figure 4B, p < 0.01) as well as in the non-native soil treatment (Figure 4C, p < 0.01). The bacterial communities displayed significant differences between native soil and non-native soil on most sampling days (Supplementary Table 7, p < 0.05 for all except day 8).

To examine how the identity of bacterial taxa changed during decomposition, the shared genera between different sampling times, and the unique genera at specific sampling times, were visualized (Figure 5). In the native soil (Figure 5A), the decaying earthworms supported more unique genera during the first 3 days, 30, 44, and 56 genera. While fewer unique genera, one (Proteobacterium) and two (Oscillochloris and Terrisporobacter), were present on days 5 and 8, respectively. The unique genera in the non-native soil (Figure 5B) showed a similar pattern to that in native soil with the number of unique genera before day 3 being 23, 54, and 29 genera, and then dropping to eight (e.g., Brevibacillus and Romboutsia) and seven (e.g., Streptosporangium and Allocatelliglobosispora) on days 5 and 8, respectively. As decomposition proceeded, the number of genera shared by the two adjacent days also decreased. The number dropped to four and two, from 29 and 22 genera in native soil (Figure 5A), and decreased to three and one, from 14 and 37 genera in non-native soil (Figure 5B).

To further examine how the abundance of specific bacterial genera changed during earthworm decomposition, we used linear regression to test for significant trends over time. Overall, 28% of 744 genera in the native soil treatment and 22% of 729 genera in the non-native soil treatment exhibited significant differential abundance as decomposition proceeded (Figures 6A,B and Supplementary Table 8). Specifically, the relative abundances of more genera on decaying earthworm decreased as decomposition proceeded, 194 in native soil and 149 in non-native soil. In contrast, only a few genera increased in abundance during decomposition, 11 in native soil and nine in non-native soil (representative genera shown in Figures 6C,D). These results suggested the abundance of more taxa on decaying earthworms decreased with decomposition time.

Source tracking was used to quantify the contribution of the likely source environments (source) to the microbial community on decaying earthworms (Figures 6E,F). The “soil” referred to the soil-resident microbes derived from incubated soil, which was defined as the soil bacterial community at day 0. The “earthworm” defined the indigenous bacteria in the earthworm, which was denoted as the bacterial community present on earthworm at day 0. Despite some similarities (shared OTUs), there were still considerable differences between soil and earthworm source (Supplementary Figure 3). In the native soil treatment, the earthworm-derived microbes accounted for 79.9% on day 1 and then decreased to 23.4% by day 8 (Figure 6E). In the non-native soil treatment, the earthworm-derived microbes accounted for 84.3% on day 1 and then decreased to 45.4% by day 8 (Figure 6E). The proportion of earthworm-derived microbes in the native soil treatment tended to be lower than in the non-native soil treatment (Figure 6E, p < 0.05 on day 8). Among the different compartments analyzed in this study, the contribution of soil-derived microbes was lower than that of earthworm-derived microbes (Figures 6E,F, p < 0.05 for all pairs). In the native soil treatment, the soil-derived microbes accounted for 4.8% on day 1 and then decreased to 0.8% (Figure 6F). In the non-native soil treatment, the soil-derived microbes accounted for 1.6% on day 1 and then decreased to 0.5% by day 8 (Figure 6F). In contrast to earthworm-derived microbes, the proportion of soil-derived microbes in the native soil treatment tended to be greater than that in the non-native soil treatment (Figure 6F, p < 0.05 on days 1 and 3). This result suggests that earthworm-derived microbes contributed primarily to the decomposition of earthworm tissues.

Shifts in the DOM composition were significantly correlated to the changes in bacterial community composition for earthworms in native (Figure 7B, R = 0.54, p = 0.001) and non-native (Figure 7C, R = 0.66, p = 0.001) soils. The effects of the bacterial community were still maintained when we merged the two datasets from different soils (Figure 7A, R = 0.53, p = 0.001). In addition, the bacterial community compositions were significantly correlated to mass loss in native (Supplementary Figure 4B, R = 0.15, p = 0.04) and non-native (Supplementary Figure 4C, R = 0.25, p = 0.002) soils. The effects of bacterial community on mass loss were still maintained when we merged the two datasets from the different soils (Supplementary Figure 4A, R = 0.25, p = 0.002).

To understand how bacterial functional groups changed during earthworm decomposition, and thus link functional group succession to DOM shift, bacterial taxa were assigned to different functional groups using FAPROTAX database. Then Pearson’s correlation analysis was conducted to explore the relationships between the abundance change of each functional group and the amount change of each DOM type over time (Figure 8 and Supplementary Table 9).

Among the 90 functional groups in the FAPROTAX database, 42 and 40 groups were present in at least one of the samples in the native soil and non-native soil treatments, respectively. After false discovery rate correction, 35 functional groups in native soil and 34 functional groups in non-native soil, whose relative abundance changes were significantly correlated to the amount change of a specific DOM type at a specific sampling time (adjusted p < 0.05), were defined as those intimately associated with temporal DOM shift during earthworm decomposition. Among these functional groups, 29 were shared in both native soil and non-native soils. The direction of correlation between the abundance change of specific functional group and the amount change of each DOM type was consistent, within a specific sampling time in each soil treatment. Moreover, the direction of correlations between changes in specific functional group abundance and DOM amount was consistent at different sampling times except for the fermentation group.

We further examined functional group succession during decomposition. Generally, relative to the later time points, a larger number of functional groups were identified to be associated with the change of a specific DOM type in both native soil and non-native soil during the early decomposition phase (Figure 8). This result is consistent with the changes of bacterial richness during earthworm decomposition, further indicating that more bacterial taxa were involved in metabolizing fresh organic matter. Also, based on the phase in which the functional group imposes an effect, the functional groups could be classified into five categories: those associated with DOM changes (1) during early decomposition phase (Figure 8 and Supplementary Figure 5A, 18 and 15 in native and non-native soils, respectively), (2) during middle decomposition phase (Figure 8 and Supplementary Figure 5B, 11 and 10), (3) during late decomposition phase (Figure 8 and Supplementary Figure 5C, one and three), (4) during the entire decomposition process (Figure 8 and Supplementary Figure 5D, two and three), and (5) those irregularly associated with DOM changes during different decomposition phases (Figure 8 and Supplementary Figure 5E, three and three). For example, during the early decomposition phases, the functional groups involved in the nitrogen cycling and methylotrophy were predominant in the native soil treatment, while functional groups related to sulfur cycling, ureolysis, and chemoheterotrophy were predominant in the non-native soil treatment (Supplementary Figure 5A). The functional groups associated with different decomposition phases are shown in Supplementary Figure 5. Particularly, for some shared functional groups in the two soils, there was a time lag in the non-native soil treatment (e.g., methylotrophy and nitrite respiration).

We also examined how the functional groups may mediate the accumulation and consumption of a specific DOM type based on the direction of the correlations. For the 29 functional groups whose abundance changes significantly correlated with the amount change of a specific DOM type in both soils, 16 showed the same direction of correlation in native and non-native soils (nine positive and seven negative), while 13 functional groups showed the opposite directions between the two soils (Figure 8 and Supplementary Figure 5). The positively correlated groups in both soils included those associated with nitrogen cycling, methylotrophy, and dark hydrogen oxidation; and the negatively correlated groups in both soils included those associated with sulfur cycling, chemoheterotrophy, and animal parasites. We also found that six functional groups were significantly and positively correlated with DOM changes in the native soil only, while five functional groups were significantly correlated with DOM changes only in the non-native soil, and most of those correlations (four of five functional groups) were negative.