Abstract

Host evolutionary history is a key factor shaping the earthworm cast microbiome, although its effect can be shadowed by the earthworm’s diet. To untangle dietary from taxon effects, we raised nine earthworm species on a uniform diet of cow manure and compared cast microbiome across species while controlling for diet. Our results showed that, under controlled laboratory conditions, earthworm microbiomes are species-specific, more diverse than that of the controlled diet, and mainly comprised of native bacteria (i.e. not acquired from the diet). Furthermore, diet has a medium to large convergence effect on microbiome composition since earthworms shared 16%–74% of their bacterial amplicon sequence variants (ASV). The interspecies core microbiome included 10 ASVs, while their intraspecies core microbiomes were larger and varied in ASV richness (24%–48%) and sequence abundance across earthworm species. This specificity in core microbiomes and variable degree of similarity in bacterial composition suggest that phylosymbiosis could determine earthworm microbiome assembly. However, lack of congruence between the earthworm phylogeny and the microbiome dendrogram suggests that a consistent diet fed over several generations may have weakened potential phylosymbiotic effects. Thus, cast microbiome assembly in earthworms seem to be the result of an interplay among host phylogeny and diet.

Introduction

Earthworms are important members of temperate terrestrial ecosystems, where they represent the largest animal biomass component of the soil (Edwards 2004). Earthworms show a wide range of feeding strategies from pure detritivores, which feed on dead organic matter to strictly geophagous that ingest soil and feed on soil organic matter (Edwards 2004). These different feeding strategies are expected to affect the composition and diversity of their gut microbiomes because of the diversity of microbial and nutrient inputs. The earthworm gut is an anaerobic environment in which ingested microorganisms may perish or flourish according to their physiological requirements (Drake and Horn 2007). In addition, modifications occurring during gut transit are critical to soil ecosystems since modified microbial communities return to soil as faeces (i.e. earthworm casts), modifying decomposition rates (Aira and Domínguez 2011).

Previous studies using technologies of microbial profiling with low taxonomic resolution (i.e. phospholipid fatty acids, ribosomal intergenic spacer analysis, or single-strand conformation polymorphism) and mimicking natural conditions, suggest that earthworm microbiome composition could be species-specific (Thakuria et al. 2009, Nechitaylo et al. 2010, Gómez-Brandón et al. 2011, 2012, Aira and Domínguez 2014). These results suggest that permanent symbiosis between host and microbiota could have an important role in earthworm microbiome assembly. However, recent studies indicate that the microbial communities of earthworm casts may reflect the composition of the ingested microbial communities (Zeibich et al. 2019 and references therein). This would suggest an environmental effect on earthworm microbiome assembly, where diet would be a prominent factor driving earthworm microbiome assembly. Moreover, recent studies point to the existence of a core microbiome in earthworms (Aira et al. 2015, Liu et al. 2018, 2019). However, comparative analyses of earthworm cast microbiomes across multiple species are largely obscured by differences in their diets, which confounds microbiome contributions from the earthworms (i.e. native) vs. diet. By feeding a uniform diet to multiple earthworm species in a controlled environment, one could effectively control the dietary input to the microbiome and test for host-specific microbial community differences.

Therefore, our general goal here is to assess how cast microbiome varies across earthworm species while controlling and accounting for diet. Our experimental design includes nine earthworm species with different ecological preferences, allowing us to effectively measure and partition the relative contributions of native earthworm microbiome across species while controlling for diet. We hypothesize that differences in microbiome composition and structure (i.e. beta-diversity—how bacterial taxa are distributed among earthworm species) will be mainly attributed to earthworm species effects, while high similarities in microbial composition would suggest microbial convergence resulting from feeding earthworms the same diet. We coupled next-generation 16S rRNA amplicon sequencing with bioinformatic and statistical methods to characterize the microbiomes of the earthworm casts and diet and to determine: (i) the composition and structure of the earthworm cast microbiomes, (ii) whether bacteria from the diet appear in the earthworm cast microbiome and whether this effect is earthworm species-specific, (iii) how large is the earthworm cast native bacterial community (i.e. bacteria not coming from the diet) and whether it varies across earthworm species, and (iv) the composition and relative proportion of the inter- and intraspecies core microbiomes.

Material and methods

Diet, earthworm species, and cast sampling

We fed all the earthworm species growing in the stock cultures and in Petri dishes the same diet consisting of cow (Bos taurus) manure collected from a farm near the University of Vigo (Galicia, NW Spain). We used this uniform diet (cow manure) to control for microbial variation associated with more natural diets and because at least eight of the chosen earthworm species feed on animal faeces in their usual natural environment, although some of them can also feed on decaying organic matter. Thus, although this is a controlled laboratory study, it mimics to a certain extent natural conditions in terms of earthworm diet and rearing conditions. Ten juvenile (< 10 mg) individuals of the nine earthworm species Dendrobaena hortensis, Eisenia andrei, E. fetida, Lumbricus friendi, L. rubellus, L. terrestris (Lumbricidae), Eudrilus eugeniae (Eudrilidae), Microscolex dubius (Acanthodrilidae), and Perionyx excavatus (Megascolecidae) were picked from ongoing stock cultures fed with cow manure. We chose these species for being representative of their genera and/or families. These stock cultures were reared under laboratory conditions for more than 1 year, feeding on cow manure ad libitum. For most of the studied earthworm species, this means that at least three generations occurred in the stock culture feeding on cow manure. By doing this, we ensure that any effect of previous diets on earthworm microbiome composition will be minimized. Juvenile earthworms were then placed individually in sterile plastic Petri dishes that were partly filled (75%) with vermicompost collected from their specific stock culture and fed ad libitum with cow manure. Petri dishes (N = 90) were randomly placed in an incubation chamber at 20°C and 90% relative humidity. To sample earthworm casts, earthworms were removed from the Petri dishes, washed three times with sterile distilled water, and placed on moistened sterile filter paper in clean, sterile Petri dishes (cast sampling dishes). All handling tasks were done under sterile conditions in a laminar flow cabinet. Sampling dishes were then placed in the incubation chamber for 24 h. Once the earthworms emptied their guts, fresh (newly egested) casts were collected from each sampling dish with a sterile spatula (flame sterilized between samplings) and stored in 1.5 ml Eppendorf tubes at −80°C. After that, earthworms were returned to their own breeding dishes. This process was also done under sterile conditions and repeated at least five times to yield 250 mg of fresh cast material per each earthworm specimen. Thus, in total, we had 90 earthworm casts samples (from individual earthworms) and 10 cow manure samples (N = 100).

High-throughput sequencing and analysis of 16S rRNA amplicons

DNA was extracted using the MO-BIO PowerSoil® kit following the manufacturer’s protocols as in previous studies (Aira et al. 2015, Domínguez et al. 2021). All laboratory procedures were performed under a laminar flow hood to prevent contamination of the samples with microorganisms from the surrounding environment. We amplified and sequenced a fragment of the 16S rRNA gene covering the V4 region with the dual-index sequencing strategy described by Kozich et al. (2013). Sequencing was done on an Illumina MiSeq at the Center for Microbial System, University of Michigan. A total of seven samples from L. friendi and two from L. rubellus did not amplify and were not included in the analysis.

DADA2 (version 1.9) was used to infer the amplicon sequence variants (ASVs) present in each sample (Callahan et al. 2016a). Exact sequence variants provide a more accurate and reproducible description of amplicon-sequenced communities than is possible with OTUs defined at a constant level (97% or other) of sequence similarity (Callahan et al. 2017). Bioinformatics processing largely followed the DADA2 pipeline tutorial (https://benjjneb.github.io/dada2/tutorial.html), but truncating forward and reverse reads at 240 nt and 200 nt, respectively, with no ambiguous bases allowed, and each read required to have less than two expected errors based on their quality scores for 16S data. ASVs were independently inferred from the forward and reverse reads of each sample using the run-specific error rates, and then read pairs were merged. Chimeras were identified in each sample, and ASVs were removed if identified as chimeric in a sufficient fraction of the samples in which they were present. Taxonomic assignment was performed against the Silva database (version 132) using the implementation of the RDP naive Bayesian classifier available in the DADA2 R package (min boot 80; Wang et al. 2007). Sequence data have been uploaded to the GenBank SRA database under accession PRJNA749153.

Statistical analysis

We analyzed and plotted all the data in R version 3.5.1 using the phyloseq (McMurdie and Holmes 2013) and ggplot2 packages (Wickham 2009). Before any statistical analysis, we applied a prevalence filtering on our data keeping only ASVs present in at least two samples and also removing singletons and doubletons (Callahan et al. 2016b) present in the filtered data set. By doing this, we exclude taxa that probably are artefacts of the data collection process (Callahan et al. 2016b). We used the full data set for α-diversity calculations because α-diversity estimates are highly dependent on the number of singletons. All other analyses were done using the filtered data set. Rarefaction curves indicated that sampling depth was optimal for all samples for both the full data set (7154 ASVs and 1530 249 sequences; Figure S1, Supporting Information) and the filtered data set (2952 ASVs and 1466 133 sequences; Figure S2, Supporting Information). We normalized ASV counts using the variance-stabilizing transformation to render the data homoscedastic for downstream analyses (we chose the variance-stabilizing transformation over the rlog transformation due to the relatively equal read-counts across samples; Love et al. 2014). We used raw ASV counts when analysing differential ASV abundances with negative binomial models (Love et al. 2014).

An approximate maximum-likelihood phylogenetic tree was inferred using FastTree 2.1 (Price et al. 2010). We tested whether microbiomes from diet and earthworm cast microbiomes (i.e. joining data from all earthworm species) have the same structure, and whether microbiomes from each earthworm species have the same structure across species. To do this, we built dendrograms (Ward method) with distance matrices (weighted and unweighted Unifrac, Bray–Curtis and Jaccard) and ran UniFrac weighted and unweighted tests (10 000 iterations) as implemented in mothur (1.41.1; Schloss et al. 2009). We calculated corrected P-values with the Benjamini–Hochberg FDR in R. These tests are best suited for hypothesis testing as shown by Schloss (2008).

Taxonomic α-diversity was calculated as the number of observed ASVs and estimated as Chao1 richness, and diversity was estimated with the Shannon index. Phylogenetic diversity was calculated as Faith’s phylogenetic diversity (Faith 1992). We used mothur commands phylo.diversity (Faith index) and summary.single (observed richness, Chao1 richness, and Shannon index). We tested for differences between diet and earthworm microbiomes on α-diversity with linear model analysis. We checked the no,rmality of residuals and homogeneity of variance across groups. Tukey’s test was used for post hoc comparisons, and Benjamini–Hochberg FDR was used as a multiple test correction method using the “multcomp” package in R (Hothorn et al 2008).

We studied differential abundance of ASV and bacterial phyla between diet and earthworm microbiomes using negative binomial models as implemented in the package DESeq2 (Love et al. 2014). Differential abundances of ASVs and other bacterial taxa were determined according to Wald tests (DESeq2 default test) and obtained P-values adjusted by false discovery rate (P-adj < .05). Because multiple pairwise Wald tests were conducted for each diet and earthworm microbiome and among microbiomes of earthworm species, we further adjusted “raw” P-values using the Benjamini–Hochberg FDR method.

We defined native bacterial ASV as those ASVs present in earthworm microbiome samples after removing ASVs present in the diet microbiome. We estimated the number of native ASVs of each individual per earthworm species by randomization (Monte Carlo, 10 000 simulations) as the total number of ASVs minus the number of ASVs from one randomly selected diet microbiome. The statistical significance of the different proportions of native and diet ASVs between earthworm species was estimated using Monte Carlo (10 000 simulations) analysis (Manly 1997). For each of the 10 000 simulations, we generated the t-statistic of a generalized linear model run with earthworm species as a fixed factor. We used a quasibinomial error distribution and log link function to reduce the deviance in the model. For the null hypothesis, we followed the same procedure, but with data randomly allocated to earthworm species. We compared the estimated distribution of the t-test statistic in the original data with the randomized t-test distribution (null hypothesis). We calculated approximated P-values as the probability of randomized t-test values exceeding the distribution of t-test statistics based on the estimated data using 10 000 randomizations.

To check the existence of a core microbiome, we removed any ASVs of diet from earthworm samples because we presumed that those ASVs were diet dependent (i.e. what would not be in earthworms fed a different diet). Then, we first checked for the existence of an earthworm core microbiome. To do this, we used a prevalence criterion in which an ASV should be present in 100, 98, 95, 90, 85, or 80% of the samples from all earthworm species. We also looked for core microbiomes for each earthworm species using the most restrictive prevalence criteria (100% of samples for each earthworm species).

To test the relationship between the evolutionary history of earthworm and the amount of ASVs shared among them, we ran a Mantel test using as dissimilarity matrices the DNA distance (dist.alignment from seqinr library, Charif and Lobry 2007) of earthworm species and the number of shared ASVs (function mantel from library vegan, 999 permutations; Oksanen et al. 2020). To build the phylogenetic tree of earthworm species, we extracted total genomic DNA from tissues of earthworm species Eu. eugeniae, M. dubius, and P. excavatus using the DNAeasy Tissue kit (Qiagen). We amplified regions of the nuclear 28S rRNA and 18S rRNA and mitochondrial 16S rRNA, 12S rRNA, NADH dehydrogenase (ND1), cytochrome oxidase subunit I and II (COII and COII), and tRNA Asn, Asp, Val, Leu, Ala, Ser, and Leu genes (Domínguez et al. 2015). Nucleotide sequences from tRNAs were combined and all sequences from each gene region were aligned in MAFFT v6 (Katoh et al. 2005, Katoh 2008) under the global (G-INS-i) algorithm with representative sequences from the other six earthworm species obtained from Domínguez et al. (2015). The general time reversible model of evolution, with proportion of invariable sites and gamma distribution, was selected for each data partition (GTR + G + I). Maximum likelihood analysis of the eight concatenated partitions was performed in RAxML v8 (Stamatakis 2014) with 10 000 searches. Clade support was assessed using the nonparametric bootstrap procedure (Felsenstein 1985) with 10 000 bootstrap replicates.

Results

Composition, structure, and diversity of earthworm cast microbiomes

Diet and earthworm microbiomes showed significant differences in their bacterial community structure (unweighted UniFrac distances), forming two main clusters (Fig 1A, diet vs. earthworm microbiome, P < .0001). Moreover, the structure of earthworm microbiomes was species-specific in most cases, with bacteria from each earthworm species forming monospecific clusters (Fig. 1A). Some dendrograms, however, did not show monospecific clusters like M. dubius and L. rubellus in weighted UniFrac, or L. rubellus and L. friendi in Bray–Crutis and Jaccard dendrograms (Figure S3, Supporting Information). Nonetheless, all pairwise comparisons between clusters were significantly different with UniFrac unweighted and weighted tests, P < .05 FDR corrected. We also found significant pairwise comparisons between earthworm species clusters for the other distance measures used (weighted UniFrac, Bray–Curtis and Jaccard). It is important to note that earthworm microbiome assembly (Fig. 1A) was not mirrored by earthworm evolutionary relationships (Fig. 1B), with only two earthworm species (E. andrei and E. fetida) showing topological congruence. Earthworm microbiome showed increased abundances for most observed bacterial ASVs, with up to 83% of the bacterial ASVs showing higher abundances compared to diet microbiome (Fig. 2). This effect was also earthworm species-specific, with different ASVs increasing in different earthworm species (Figure S4 and Table S1, Supporting Information). In general, this increase in ASV abundance was large, with a 67% of the ASVs showing a significant log 2-fold change of 5–10, and a 23% of the ASVs (mean over all earthworm species) showing log 2-fold changes >10 (Table S1, Supporting Information). Most of the main bacterial phyla (Proteobacteria, Bacteroidetes, Planctomycetes, Actinobacteria, and Verrucomicrobia) had increased abundances in microbiomes from all earthworm species relative to diet microbiome (Fig. 1A; Figure S5, Supporting Information). Several bacterial phyla (Epsilonbacteraeota, Fibrobacteres, Kiritimatiellaeota, Lentisphaerae, and Spirochaetes) showed decreased abundances in earthworm microbiomes relative to diet microbiome (Fig. 1A; Figure S5, Supporting Information).

Structure of bacterial communities of diet microbiome (B. Taurus faeces) and earthworm microbiomes and evolutionary relationships of earthworm species. (A) Differences in structure between diet microbiome and microbiomes of earthworm species D. hortensis, E. andrei, E. fetida, Eu. eugeniae, Lumbricus friend, L. rubellus, L. terrestris, M. dubius, and P. excavatus. The dendrogram represents the dissimilarity of bacterial communities at the ASV level (variance-stabilized matrix of ASV counts, unweighted UniFrac distances, Ward method). Bars represent the relative abundance of dominant bacterial phyla. Low abundance bacterial phyla (< 1%) were grouped together. (B) Evolutionary relationships between the nine earthworm species. We have manually added the B. taurus tip as outgroup for illustrative purposes since its faeces were the earthworms diet.
Figure 1.

Structure of bacterial communities of diet microbiome (B. Taurus faeces) and earthworm microbiomes and evolutionary relationships of earthworm species. (A) Differences in structure between diet microbiome and microbiomes of earthworm species D. hortensis, E. andrei, E. fetida, Eu. eugeniae, Lumbricus friend, L. rubellus, L. terrestris, M. dubius, and P. excavatus. The dendrogram represents the dissimilarity of bacterial communities at the ASV level (variance-stabilized matrix of ASV counts, unweighted UniFrac distances, Ward method). Bars represent the relative abundance of dominant bacterial phyla. Low abundance bacterial phyla (< 1%) were grouped together. (B) Evolutionary relationships between the nine earthworm species. We have manually added the B. taurus tip as outgroup for illustrative purposes since its faeces were the earthworms diet.

Differences in abundance of bacterial ASVs between diet and earthworm microbiome. Relationship between differential ASV abundance (log2 fold change) and normalized ASV abundance. Negative and positive values indicate taxa under-represented and over-represented in earthworm microbiome relative to diet microbiome, respectively. Red circles indicate taxa showing significant (P-adj < .05) differential abundance between earthworm and diet microbiome, and red triangles indicate taxa with log2 fold change falling out of the limits of the y-axis. Inserted plots represent exemplar bacterial phyla and ASVs over- and under-represented, with asterisks showing significant differential abundances (log normalized counts) between D. hortensis (Dh), E. andrei (Ea), E. fetida (Ef), Eu. eugeniae (Eu), Lumbricus friend (Lf), L. rubellus (Lr), L. terrestris (Lt), M. dubius (Md), and P. excavatus (Pe) microbiomes and diet microbiome (dm) (P-adj < .05).
Figure 2.

Differences in abundance of bacterial ASVs between diet and earthworm microbiome. Relationship between differential ASV abundance (log2 fold change) and normalized ASV abundance. Negative and positive values indicate taxa under-represented and over-represented in earthworm microbiome relative to diet microbiome, respectively. Red circles indicate taxa showing significant (P-adj < .05) differential abundance between earthworm and diet microbiome, and red triangles indicate taxa with log2 fold change falling out of the limits of the y-axis. Inserted plots represent exemplar bacterial phyla and ASVs over- and under-represented, with asterisks showing significant differential abundances (log normalized counts) between D. hortensis (Dh), E. andrei (Ea), E. fetida (Ef), Eu. eugeniae (Eu), Lumbricus friend (Lf), L. rubellus (Lr), L. terrestris (Lt), M. dubius (Md), and P. excavatus (Pe) microbiomes and diet microbiome (dm) (P-adj < .05).

These marked changes in composition and structure between microbiomes of diet and earthworm species did not result in differences in α-diversity richness (P = .6) or diversity (P = .2), although earthworm microbiomes were up to a 21% more phylogenetically diverse (P = .003, Fig. 3A). Across earthworm species, we found that all earthworm microbiomes had higher phylogenetic diversity compared to diet with three species (E. andrei, E. fetida, and P. excavatus) having statistically significantly higher levels of diversity (Fig. 3B) and richness (Figure S6, Supporting Information).

Phylogenetic α-diversity of bacterial communities (A) of diet (cow manure) and earthworm microbiome, and (B) in diet and each earthworm species microbiome: D. hortensis, E. andrei, E. fetida, Eu. eugeniae, L. friendi, L. rubellus, L. terrestris, M. dubius, and P. excavatus. Different letters over box plots indicate significant differences between treatments (Tukey HSD test, FDR corrected). (C) Chord diagram represents the number of shared ASVs between microbiomes of earthworm species pairs using all ASVs of filtered data set. Samples were collapsed at earthworm species level.
Figure 3.

Phylogenetic α-diversity of bacterial communities (A) of diet (cow manure) and earthworm microbiome, and (B) in diet and each earthworm species microbiome: D. hortensis, E. andrei, E. fetida, Eu. eugeniae, L. friendi, L. rubellus, L. terrestris, M. dubius, and P. excavatus. Different letters over box plots indicate significant differences between treatments (Tukey HSD test, FDR corrected). (C) Chord diagram represents the number of shared ASVs between microbiomes of earthworm species pairs using all ASVs of filtered data set. Samples were collapsed at earthworm species level.

Native bacteria dominate the earthworm cast microbiome

Most of the ASVs found in earthworm microbiomes (92%–96% for all earthworm species; Table S2, Supporting Information) were not present in the diet (so they can be considered native) and varied significantly among earthworm species (Monte Carlo analysis P = .002; Table S2, Supporting Information). These native ASVs comprised the majority of sequences found in each earthworm species casts (85%–98% for all earthworm species with a mean of 94%). Moreover, we found that microbiomes from different earthworm species shared a variable number of ASVs. Thus, 658, 472, 281, 177, 98, 56, 52, and 28 ASVs were shared by two, three, four, five, six, seven, eight, and nine earthworm species, respectively (effect of number of earthworm species Χ27 = 154.66, P < .00001; Table S2, Supporting Information). More importantly, we found that in pairwise comparisons earthworm species pairs shared 39% of the ASVs (range 16%–74%; Fig. 3C), which comprised a mean of 58% of the sequences (range 29%–86%). This result did not vary when using the full data set (i.e. not filtered) and after removing the ASVs coming from diet microbiome in both full and filtered datasets (Figure S7, Supporting Information). Additionally, we found that evolutionarily close species (i.e. E. andrei and E. fetida) did not share more ASVs than more distantly related species (i.e. D. hortensis and E. andrei or P. excavatus and E. andrei; Fig. 3C; Figure S7, Supporting Information; Mantel test, r = −0.34, P = .98). In addition, only a small fraction of shared ASVs differed between earthworm species in pairwise comparisons (mean 12%, range 4%–28%; Table S3, Supporting Information). Again, these results were similar when we used the full data set and when removed the ASVs coming from diet in both full and filtered dataset (Figure S7, Supporting Information).

The earthworm core cast microbiome

We did not find any ASVs present in all the earthworm sample replicates (interspecies microbiome)—the most restrictive definition of a core microbiome. However, using less restrictive definitions of the core microbiome (ASV prevalence of 98, 95, 90, 85, and 80% across cast samples), we found that the earthworm core microbiome comprised 3, 4, 5, 9, and 10 ASVs that accounted for 8, 9,10, 12, and 13% of the sequences, respectively (Fig. 4). We also found distinct core microbiomes for each earthworm species. Although these core microbiomes comprised a small fraction of ASVs (mean 74, range 46–134 ASVs, which corresponds to a mean of 7.8 and a range of 5.6%–14.8% of ASVs), they represented a large proportion of all the sequences (mean 52%, range 43%–63% sequences; Table S2, Supporting Information). Most of the core ASVs were earthworm species-specific (235 of 360), i.e. they only appeared in one earthworm species core microbiome. From these ASVs, the most frequent ASVs belonged to the genera Flavobacterium, Pir4 lineage, Pirellula, Pseudomonas, Luteolibacter, and Chthoniobacter, and the families Chitinophagaceae, Burkholderiaceae, and Verrucomicrobiaceae. We also found that some core ASVs appeared in several earthworm species: 6, 11, and 26 ASVs were shared by five, four, and three earthworm species, respectively (Table S2, Supporting Information).

An earthworm core microbiome. (A) Percentage of core and noncore AVS in earthworm core microbiome. (B) Mean relative abundance (%) of ASVs (genus or most inclusive taxonomy found) from the core microbiome of casts from the nine earthworm species (D. hortensis, E. andrei, E. fetida, Eu. eugeniae, L. friendi, L. rubellus, L. terrestris, M. dubius, and P. excavatus). The top x-axis is for ASV1 and the down x-axis is for the other ASVs. We defined the earthworm core microbiome as those formed by ASVs present in 98, 95, 90, 85, and 80% of samples (N = 83).
Figure 4.

An earthworm core microbiome. (A) Percentage of core and noncore AVS in earthworm core microbiome. (B) Mean relative abundance (%) of ASVs (genus or most inclusive taxonomy found) from the core microbiome of casts from the nine earthworm species (D. hortensis, E. andrei, E. fetida, Eu. eugeniae, L. friendi, L. rubellus, L. terrestris, M. dubius, and P. excavatus). The top x-axis is for ASV1 and the down x-axis is for the other ASVs. We defined the earthworm core microbiome as those formed by ASVs present in 98, 95, 90, 85, and 80% of samples (N = 83).

Discussion

In this study, we aimed to untangle the microbial differences associated to earthworm species with different evolutionary histories. To remove variation associated to diets and habitat preferences, we reared the earthworms in the laboratory and fed them a uniform diet of cow manure. We then collected fresh casts and analyzed their bacterial composition using 16S rRNA amplicon next-generation sequencing. Our main results showed that under controlled conditions the earthworm cast microbiome: (i) were mainly composed by native bacteria, (ii) tended to be more diverse than the diet microbiome, and (iii) were significantly different in structure among host species. We also detected inter- and intraspecies core microbiomes of variable ASV richness and sequence abundance.

Composition, structure, and diversity of earthworm microbiomes

Microbiomes from the nine earthworm species were mostly dominated by ASVs that did not come from the diet, which resulted in marked differences between diet and earthworm microbiomes in composition, structure, and diversity. Bacteroidetes, e.g. was more abundant in all earthworm microbiomes, while Firmicutes significantly varied among earthworm species (Figure S5, Supporting Information). Moreover, our results confirmed that the high abundances of Actinobacteria and Proteobacteria previously found for E. andrei can be extended to the other earthworm species, as well as the increases in other bacterial phyla like Bacteroidetes, Chloroflexy, Planctomycetes, and Verrucomicrobia (Aira et al. 2015). At the ASV level, we found slight and unique differences in abundance between diet and earthworm microbiomes in each earthworm species. At the phylum level, the composition of bacterial communities of earthworm microbiomes from all earthworm species were relatively similar to those of epigeic earthworm species, such as E. andrei, Eu. eugeniae, andL. rubellus (Knapp et al. 2009, Aira et al. 2015, 2016, Wüst et al. 2009, Schulz et al. 2015) and to that of the anecic earthworm L. terrestris (Wüst et al. 2011), despite of differences in diets—some even used earthworms reared in their natural habitats. Earthworms are ecologically classified into three general ecotypes based on their feeding habits (Bouché 1977, Lee 1985); epigeic and anecic earthworms feed preferentially on dead organic matter, whereas endogeic earthworm feed on mineral soil. Consistent with previous work, our analyses found that microbiomes from the four epigeic earthworm species above had different proportions of bacteria belonging to phyla Proteobacteria, Firmicutes, Actinobacteria, Planctomycetes, Chloroflexi, Gemmatimonadetes, Bacteroidetes, and Cytophagales, among others. Bacterial communities of casts from L. rubellus and L. terretris were dominated by bacterial phyla Proteobacteria (> 30%), Actinobacteria, Acidobacteria, Bacteroidetes, and Firmicutes (17%–30%; Knapp et al. 2009, Wüs et al. 2009, Wüst et al. 2011), as we have observed here despite of marked differences in diet and rearing condition among these studies.

Differences between diet or environment and animal microbiomes have been also found in earthworms and nematodes (Gómez-Brandón et al. 2011, 2012, Aira and Domínguez 2014, Aira et al. 2015, Berg et al. 2016, Dirksen et al. 2016). Altogether, these results suggest that the symbiont consortia involved in earthworm digestive processes could be species-specific, and that, consequently, earthworm microbiomes should reflect only partially the composition of the ingested microbial communities. One could expect a nonsegregated clustering of earthworm microbiomes if diet were to have a strong effect on earthworm microbiome composition. This could occur whether bacteria from diet massively colonize earthworm microbiome or by a convergence in earthworm microbiome composition due to feeding of the same diet, although bacterial functional redundancy could mask this effect. However, we found that host species was the main determinant of earthworm microbiome composition, as previously found using other more or less natural experimental conditions (e.g. Thakuria et al. 2009, Gómez-Brandón et al. 2012, Aira and Domínguez 2014). In fact, bacteria from diet were hardly found in earthworm microbiomes, as it occurred in earthworms fed with sewage sludge (Domínguez et al. 2021). Similarly, native bacteria mostly dominated the gut microbiome of geophagous earthworm Aporrectodea caliginosa fed with soil (up to 89%; Aira et al. 2022). This result suggests that most of the bacteria found in earthworm microbiomes are therefore “native,” and contrasts with previous findings indicating that more than 95% of the gut bacteria come mainly from the ingested diet (Zeibich et al. 2019, and references therein). However, since our experimental design specifically controlled for the confounding effect of diet, more field studies should be performed to validate these results in natural environmental conditions.

Native bacteria dominate the earthworm cast microbiome

We found that most of the bacteria present in earthworm microbiomes were native (i.e. they did not come from the diet ingested) and accounted for most of the sequences (mean 94% for ASV for all earthworm species). Animal species-specific bacteria, i.e. bacterial lineages found only in the animal and not in the feeding environment, have been described in sponges and nematodes (Reveillaud et al. 2014, Berg et al. 2016, Dirksen et al. 2016) and also recently in the earthworm species A. caliginosa and E. andrei (Domínguez et al. 2021, Aira et al. 2022). The proportion of native bacteria was similar among all earthworm species. Moreover, an important part of the native earthworm cast bacterial community is species-specific, with a different proportion of ASVs shared among species. The proportion of ASVs shared seemed to be more related to the earthworm ecology than to the earthworm phylogeny. Thus, epigeic earthworm species (i.e. E. andrei, E. fetida, D. hortensis, or P. excavatus) shared more ASVs among them than with anecic earthworm species L. terrestris and vice versa. Similar results, i.e. a large proportion of shared bacterial taxa among different species, have also been reported for Drosophila species (Chandler et al. 2011). Contrary to our results, Zeibich et al. (2019) found that most phylotypes present in the gut content of the earthworm L. terrestris feeding on a mixture of roots, grass, and leaves were also present in soil samples, and that those phylotypes comprised the great majority of the sequences. These data also contrast with a recent study describing the bacterial communities of gut and casts in the geophagous earthworm A. caliginosa (Aira et al. 2022). In the aforementioned study, we found that most of the bacteria were unique to each microbiome soil (82%), gut (89%), and casts (75%), with a small overlap bacterial community composition between soil and gut and soil and casts (Aira et al. 2022). Vermicompost is mainly comprised of earthworm casts in a variable degree of ageing (Domínguez et al. 2010), so any accidental ingestion would result in the incorporation of the same bacterial pool already present in the casts. Moreover, bacteria from vermicompost are horizontally and not randomly acquired during cocoon formation in earthworms (Aira et al. 2018). These bacteria will colonize the earthworm gut during embryonic development (Davidson and Stahl 2008), which complicates discerning the origin of these bacteria—i.e. adult ingestion or embryonic acquisition. In addition, a recent study (Domínguez et al. 2021) has shown that earthworm casts share a relatively small amount of bacterial ASVs with vermicompost (15%), which supports the rationale that earthworms do not or hardly ingest vermicompost.

The earthworm core microbiomes

We found an earthworm core microbiome for the nine earthworm species, which was formed by up to 10 bacterial ASVs that accounted for 13% of the sequences. It is important to note that only three bacterial ASVs comprised the most restrictive definition of the interspecies core microbiome, i.e. ASV sample prevalence of 98%. Earthworm species-specific core microbiomes were larger and, more importantly, comprised roughly 50% of all sequences of each earthworm microbiome. These core ASVs are probably transferred vertically, as occurred with nephridial bacterial symbionts (Davidson et al. 2013). Alternatively, earthworms may horizontally acquire these core ASVs during cocoon formation, which is a selective process (Aira et al. 2018). Most core bacteria were earthworm species-specific, suggesting that diet has hardly any influence in their composition. Future studies should test if these core ASVs are maintained across different diets. If not, it would emphasize the role of diet as a factor structuring earthworm microbiomes, as seen in mammals and humans (Muegge et al. 2011, Nishida and Ochman 2018). Moreover, different ASVs from the genera Flavobacterium, Pirellula, Pir4 lineage, Pseudomonas, Chthoniobacter, and Luteolibacter repeatedly appeared in several earthworm core microbiomes. This indicates a high specificity for these bacteria and pointed again to vertical transmission.

Our results for E. andrei partially agree with a previous description of its core microbiome (Aira et al. 2015), where Proteobacteria (43%), Actinobacteria (40%), and Firmicutes (14%), with minor contributions from Bacteroidetes (3%) comprised the core microbiome. We did not find any Firmicutes in E. andrei core microbiome, and proportions of Proteobacteria (35%), Bacteroidetes (22%), and Actinobacteria (2.5%) differed. Furthermore, in our current study Planctomycetes made a large contribution (12%) to the core microbiome. It is important to note that the core microbiome of E. andrei and E. fetida, two ecologically close but biologically and phylogenetically different species (Domínguez et al. 2005, 2015, Pérez-Losada et al. 2005), were different but shared 33 ASVs, which comprise the 28% (E. andrei) and 44% (E. fetida) of the core ASVs microbiome. Similarly, the core microbiomes of L. rubellus and L. terrestris, two species from same genus (James et al. 2010), only shared the 17% of the ASVs. There is a report of a core microbiome for cocoons of E. andrei and E. fetida, which comprised 88 and 23 bacterial ASVs (Aira et al. 2018). In this case, Proteobacteria (50% and 52%) and Bacteroidetes (24% and 35%) comprised most of the core microbiome, with minor contributions of Actinobacteria (8% and 8%) and Verrucomicrobia (8% and 4%, for E. andrei and E. fetida, respectively). These values are similar to those found here for gut microbiomes. In our case, the two core microbiomes only shared 7% (E. andrei) and 26% (E. fetida) of the core ASVs. The differences found at genus or ASV levels arose because the bacterial communities of earthworm cocoons were mainly dominated by vertically transmitted symbionts and nonrandomly selected environmental bacteria (grape marc and vermicompost from grape marc, Aira et al. 2018). In the case of L. rubellus, the core microbiome differs from that described for the same species by Pass et al. (2015). Thus, in our case Proteobacteria and Bacteroidetes (87%) comprised most of the core instead of Proteobacteria and Actinobacteria (80%). However, Pass et al. (2015) used whole earthworms (not only their guts or casts) that were feeding on soil; in our case, however, the core only included the gut and earthworms were fed cow manure. All this evidence underlines the distinctive roles of host evolutionary histories and diet in shaping the composition of core microbiomes, and how core microbiome could change during earthworm ontogeny as nutritional or metabolic requirements change during growth, reproduction, and senescence.

Our data clearly show that the earthworm microbiome is highly diverse and is dominated by native bacteria, whose composition is highly earthworm species-specific. This, together with the large core microbiomes found for each earthworm species, suggest that phylosymbiosis (i.e. congruence between host-associated microbial communities and the phylogeny of host species; Brooks et al. 2016) processes could be prevalent in earthworm microbiome assembly. However, the lack of congruence between the earthworm phylogeny and the microbiome dendrogram suggests that having the same diet over several generations has weakened the effect of phylosymbiosis on earthworm microbiome assembly. This rationale is also supported by the variable degree of bacterial taxa shared among earthworm species, which again was not related to earthworm evolutionary history. Thus, earthworm microbiome assembly seems to be the result of an interplay among host phylogeny and diet.

The conclusions of our study are bound by our experimental design. Our main goal here is to disentangle native earthworm microbiome from diet and environment, hence we reared the animals in a controlled environment and fed them the same diet. Most of the earthworm species studied here usually feed on dead-organic matter and even manure if available, but they usually live in soil. Thus, although our experimental design mimics conditions in which earthworm are cultured (but at small scale), it may not resemble their living conditions in the wild. Consequently, results should be interpreted with caution when extrapolated to natural earthworm populations.

Authors’ contributions

J.D. and M.A. conceived the idea, performed the research, and wrote the main manuscript text. M.A. analyzed data and prepared the figures. All authors edited and reviewed the manuscript.

ACKNOWLEDGEMENTS

The authors would like to thank Hugo Martínez and Alberto Da Silva for help with sample collection and DNA extraction.

Funding

This study was partly cofunded by the Spanish Ministerio de Ciencia e Innovación (PID2021-124265OB-100) and the European Union’s Horizon 2020 research and innovation programme (LABPLAS_101003954).

Conflict of interest statement. None declared.

Data availability statement

DNA sequence data have been uploaded to the GenBank SRA database under accession PRJNA749153.

References

Aira
M
,
Bybee
S
,
Pérez-Losada
M
et al.
Feeding on microbiomes: effects of detritivory on the taxonomic and phylogenetic bacterial composition of animal manures
.
FEMS Microbiol Ecol
.
2015
;
91
:
1
10
.

Aira
M
,
Domínguez
J.
Changes in nutrient pools, microbial biomass and microbial activity in soils after transit through the gut of three endogeic earthworm species of the genus Postandrilus qui and Bouché, 1998
.
J Soils Sediments
.
2014
;
14
:
1335
40
.

Aira
M
,
Domínguez
J.
Earthworm effects without earthworms: inoculation of raw organic matter with worm-worked substrates alters microbial community functioning
.
PLoS ONE
.
2011
;
6
:
e16354
.

Aira
M
,
Olcina
J
,
Pérez-Losada
M
et al.
Characterization of the bacterial communities of casts from Eisenia andrei fed with different substrates
.
Appl Soil Ecol
.
2016
;
98
:
103
11
.

Aira
M
,
Pérez-Losada
M
,
Crandall
KA
et al.
Composition, structure and diversity of soil bacterial communities before, during and after transit through the gut of the earthworm Aporrectodea caliginosa
.
Microorganisms
.
2022
;
10
:
1025
.

Aira
M
,
Pérez-Losada
M
,
Domínguez
J.
Microbiome dynamics during cast ageing in the earthworm Aporrectodea caliginosa
.
Appl Soil Ecol
.
2019
;
139
:
56
63
.

Aira
M.
,
Pérez-Losada
M
,
Domínguez
J.
Diversity, structure and sources in bacterial communities of earthworm cocoons
.
Sci Rep
.
2018
;
8
:
6632
.

Berg
M
,
Stenuit
B
,
Ho
J
et al.
Assembly of the Caenorhabditis elegans gut microbiota from diverse soil microbial environments
.
ISME J
.
2016
;
10
:
1998
2009
.

Bouché
MB.
Strategies lombriciennes
. In:
Lohm
U
,
Persson
T
(eds),
Soil Organisms as Components of Ecosystems
.
Stockholm
:
Ecological Bulletin
,
1977
,
122
32
.

Brooks
AW
,
Kohl
KD
,
Brucker
RM
et al.
Phylosymbiosis: relationships and functional effects of microbial communities across host evolutionary history
.
PLoS Biol
.
2016
;
14
:
e2000225
.

Callahan
B
,
McMurdie
P
,
Holmes
SP.
Exact sequence variants should replace operational taxonomic units in marker-gene data analysis
.
ISME J
.
2017
;
11
:
2639
43
.

Callahan
BJ
,
McMurdie
PJ
,
Rosen
MJ
et al.
DADA2: High-resolution sample inference from illumina amplicon data
.
Nat Methods
.
2016a
;
13
:
581
3
.

Callahan
BJ
,
Sankaran
K
,
Fukuyama
JA
et al.
Bioconductor workflow for microbiome data analysis: from raw reads to community analyses
.
F1000Research
.
2016b
;
5
:
1492
.

Chandler
JA
,
Morgan Lang
J
,
Bhatnagar
S
et al.
Bacterial communities of diverse Drosophila species: ecological context of a host–microbe model system
.
PLos Genet
.
2011
;
7
:
e1002272
.

Charif
D
,
Lobry
J.
SeqinR 1.0-2: a contributed package to the R project for statistical computing devoted to biological sequences retrieval and analysis
. In:
Bastolla
U
,
Porto
M
,
Roman
H
,
Vendruscolo
M
(
eds.
),
Structural Approaches to Sequence Evolution: Molecules, Networks, Populations, Series Biological and Medical Physics, Biomedical Engineering
.
Berlin
:
Springer Verlag
,
2007
,
207
32
.

Davidson
SK
,
Powell
R
,
James
S
.
A global survey of the bacteria within earthworm nephridia
.
Mol Phylogenet Evol
.
2013
;
67
:
188
-200
.

Davidson
SK
,
Stahl
DA
.
Selective recruitment of bacteria during embryogenesis of an earthworm
.
ISME J
.
2008
;
2
:
510
8
.

Dirksen
P
,
Marsh
SA
,
Braker
I
et al.
The native microbiome of the nematode Caenorhabditis elegans: gateway to a new host-microbiome model
.
BMC Biol
.
2016
;
14
:
38
.

Domínguez
J
,
Aira
M
,
Breinholt
JW
et al.
Underground evolution: new roots for the old tree of lumbricid earthworms
.
Mol Phylogenet Evol
.
2015
;
83
:
7
19
.

Domínguez
J
,
Aira
M
,
Crandall
K
et al.
Earthworms drastically change fungal and bacterial communities during vermicomposting of sewage sludge
.
Sci Rep
.
2021
;
11
:
15556
.

Domínguez
J
,
Aira
M
,
Gómez-Brandón
M.
Vermicomposting: earthworms enhance the work of microbes
. In:
Insam
H
,
Franke-Whittle
I
,
Goberna
M
(eds),
Microbes at Work: From Wastes to Resources
.
Berlin Heidelberg
:
Springer
,
2010
,
93
114
.

Domínguez
J
,
Ferreiro
A
,
Velando
A.
Are Eisenia fetida (Savigny, 1826) and Eisenia andrei Bouché, 1972 (Oligochaeta, lumbricidae) different biological species?
.
Pedobiologia
.
2005
;
49
:
81
7
.

Drake
HL
,
Horn
MA.
As the worm turns: the earthworm gut as a transient habitat for soil microbial biomes
.
Annu Rev Microbiol
.
2007
;
61
:
169
89
.

Edwards
CA
.
Earthworm Ecology
. 2nd edn.
Boca Raton
:
CRC Press
,
2004
.

Faith
DP.
Conservation evaluation and phylogenetic diversity
.
Biol Conserv
.
1992
;
61
:
1
10
.

Felsenstein
J.
Confidence limits on phylogenies: an approach using the bootstrap
.
Evolution
.
1985
;
39
:
783
91
.

Gómez-Brandón
M
,
Aira
M
,
Lores
M
et al.
Epigeic earthworms exert a bottleneck effect on microbial communities through gut associated processes
.
PLoS ONE
.
2011
;
6
:
e24786
.

Gómez-Brandón
M
,
Lores
M
,
Domínguez
J.
Species-specific effects of epigeic earthworms on microbial community structure during first stages of decomposition of organic matter
.
PLoS ONE
,
2012
;
7
:
e31895
.

Hothorn
T
,
Bretz
F
,
Westfall
P.
Simultaneous inference in general parametric models
.
Biomet J
.
2008
;
50
:
346
63
.

James
SW
,
Porco
D
,
Decaëns
T
et al.
DNA barcoding reveals cryptic diversity in Lumbricus terrestris L., 1758 (Clitellata): resurrection of L . herculeus (Savigny, 1826)
.
PLoS ONE
.
2010
;
5
:
e15629
.

Katoh
K
,
Kuma
K
,
Toh
H
et al.
MAFFT version 5: improvement in accuracy of multiple sequence alignment
.
Nucleic Acids Res
.
2005
;
33
:
511
8
.

Katoh
T.
Recent developments in the MAFFT multiple sequence alignment program
.
Briefings Bioinf
.
2008
;
9
:
286
98
.

Knapp
BA
,
Podmirseg
SM
,
Seeber
J
et al.
Diet-related composition of the gut microbiota of Lumbricus rubellus as revealed by a molecular fingerprinting technique and cloning
.
Soil Biol Biochem
.
2009
;
41
:
2299
307
.

Kozich
JJ
,
Westcott
SL
,
Baxter
NT
et al.
Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform
.
Appl Environ Microbiol
.
2013
;
79
:
5112
20
.

Lee
KE
.
Earthworms: Their Ecology and Relationships with Soils and Land Use
.
Sydney
:
Academic Press
,
1985
.

Liu
D
,
Lian
B
,
Wu
C
et al.
A comparative study of gut microbiota profiles of earthworms fed in three different substrates
.
Symbiosis
.
2018
;
74
:
21
9
.

Love
MI
,
Huber
W
,
Anders
S.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
.
2014
;
15
:
550
.

Manly
BJF
.
Randomization, Bootstrap and Monte Carlo Methods in Biology
. 2nd edn.
London
:
Chapman and Hall
,
1997
.

McMurdie
PJ
,
Holmes
SP.
phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data
.
PLoS ONE
.
2013
;
8
:
e61217
.

Muegge
BD
,
Kuczynski
J
,
Knights
D
et al.
Diet drives convergence in gut microbiome functions across mammalian phylogeny and within humans
.
Science
.
2011
;
332
:
970
4
.

Nechitaylo
TY
,
Yakimov
MM
,
Godinho
M
et al.
Effect of the earthworms Lumbricus terrestris and Aporrectodea caliginosa on bacterial diversity in soil
.
Microb Ecol
.
2010
;
59
:
574
87
.

Nishida
AH
,
Ochman
H.
Rates of gut microbiome divergence in mammals
.
Mol Ecol
.
2018
;
27
:
1884
97
.

Oksanen
J
,
Guillaume Blanchet
F
,
Friendly
M
et al.
vegan: community ecology package. R package version 2.5-7
.
2020
. .

Pass
DA
,
Morgan
AJ.
,
Read
DS
et al.
The effect of anthropogenic arsenic contamination on the earthworm microbiome
.
Environ Microbiol
.
2015
;
17
:
1884
96
.

Pérez-Losada
M
,
Eiroa
J
,
Mato
S
et al.
Phylogenetic species delimitation of the earthworms Eisenia fetida (Savigny, 1826) and Eisenia andrei Bouché, 1972 (Oligochaeta, lumbricidae) based on mitochondrial and nuclear DNA genes
.
Pedobiologia
.
2005
;
49
:
317
24
.

Price
MN
,
Dehal
PS
,
Arkin
AP.
FastTree 2 – approximately maximum-likelihood trees for large alignments
.
PLoS ONE
.
2010
;
5
:
e9490
.

Reveillaud
J
,
Maignien
L
,
Murat Eren
A
et al.
Host-specificity among abundant and rare taxa in the sponge microbiome
.
ISME J
.
2014
;
8
:
1198
209
.

Schlos
PD.
Evaluating different approaches that test whether microbial communities have the same structure
.
ISME J
.
2008
;
2
:
265
75
.

Schloss
PD
,
Westcott
SL
,
Ryabin
T
et al.
Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities
.
Appl Environ Microbiol
.
2009
;
75
:
7537
41
.

Schulz
K
,
Hunger
S
,
Brown
GG
et al.
Methanogenic food web in the gut contents of methane-emitting earthworm Eudrilus eugeniae from Brazil
.
ISME J
.
2015
;
9
:
1778
92
.

Stamatakis
A.
RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies
.
Bioinformatics
.
2014
;
30
:
1312
3
.

Thakuria
D
,
Schmidt
O
,
Finan
D
et al.
Gut wall bacteria of earthworms: a natural selection process
.
ISME J
.
2009
;
4
:
357
66
.

Wang
Q
,
Garrity
GM
,
Tiedje
JM
et al.
Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy
.
Appl Environ Microbiol
.
2007
;
73
:
5261
7
.

Wickham
H
.
ggplot2: Elegant Graphics For Data Analysis
.
Berlin
:
Springer
,
2009
.

Wüst
PK
,
Horn
MA
,
Drake
HL.
Clostridiaceae and enterobacteriaceae as active fermenters in earthworm gut content
.
ISME J
.
2011
;
5
:
92
106
.

Wüst
PK
,
Horn
MA
,
Henderson
G
et al.
Gut-associated denitrification and in vivo emission of nitrous oxide by the earthworm families Megascolecidae and Lumbricidae in New Zealand
.
Appl Environ Microbiol
.
2009
;
75
:
3430
6
.

Zeibich
L
,
Schmidt
O
,
Drake
HL.
Fermenters in the earthworm gut: do transients matter?
.
FEMS Microbiol Ecol
.
2019
;
95
:
fiy221
.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Supplementary data