Earthworm-invaded boreal forest soils harbour distinct microbial communities

Earthworm invasion in North American forests has the potential to greatly impact soil microbial communities by altering soil physicochemical properties, including structure, pH, nutrient availability, and soil organic matter (SOM) dynamics. While most research on the topic has been carried out in northern temperate forests, little is known about the impact of invasive earthworms on soil microbial communities in hemiboreal and boreal forests, characterized by a slower decay of organic matter (OM). Earthworm activities can increase OM mineralization, altering nutrient cycling and biological activity in a biome where low carbon (C) and nitrogen (N) availability typically limits microbial and plant growth. Here, we characterized and compared microbial communities of earthworm-invaded and non-invaded soils in previously described sites across three major soil types found in the Canadian (hemi)boreal forest using a space-for-time approach. Microbial communities of forest floors and surface mineral soils were characterized using phospholipid fatty acid (PLFA) analysis and metabarcoding of the 16S rRNA gene for bacteria and archaea and of the internal-transcriber-spacer-2 (ITS2) region for fungi. In forest floors, the effects of earthworm invasion were minor. In mineral soil horizons, earthworm invasion was associated with higher fungal biomass and greater relative abundance of ectomycorrhizal fungi. Oligotrophic bacteria (Acidobacteriota and Chloroflexi) were less abundant in invaded mineral soils, where Gram(+) : Gram(−) ratios were also lower, while the opposite was observed for the copiotrophic Bacteroidota. Additionally, earthworm-invaded mineral soils harboured higher fungal and bacterial species diversity and richness. Considering the important role of soil microbial communities for ecosystem functioning, such earthworm-induced shifts in their community composition are likely to impact nutrient cycling, as well as vegetation development and forest productivity at a large scale, as the invasion progresses in these (hemi)boreal systems.

1 Introduction

Microorganisms play a key role in soils, regulating litter decomposition and nutrient and carbon (C) availability. By ingesting litter and mixing it with mineral soil, earthworms can provide microorganisms with better access to C sources and thus accelerate decomposition and nutrient cycling (Blouin et al., 2013; Curry and Schmidt, 2007; Edwards, 2004). In Canada and the northern United States, the last glaciation extirpated native earthworm populations, resulting in the development of forests devoid of earthworms. These northern forests are now being invaded by exotic earthworm species, which can substantially alter existing soil organic matter (SOM) dynamics (Frelich et al., 2019). Changes in microbial community composition and diversity resulting from earthworm presence have been observed under laboratory (Chang et al., 2016; Butenschoen et al., 2007; Gómez-Brandón et al., 2010; de Menezes et al., 2018) and field conditions, mostly in temperate forests (Price-Christenson et al., 2020; Dempsey et al., 2013; Groffman et al., 2015; McLean et al., 2006; McLean and Parkinson, 1997).

Comparable studies are minimal for boreal forests despite multiple reports of earthworm invasion occurring throughout North America (Cameron et al., 2007; Moore et al., 2009). In temperate forests, invasive earthworms can be found in more than 80 % of suitable habitats, as reported in an extensive survey of the Great Lakes region of the United States (Fisichelli et al., 2013; Frelich et al., 2019). In boreal forests, 50 % of the low-human-impact sites sampled in Alaska by Saltmarsh et al. (2016) were invaded, and Cameron and Bayne (2009) estimated that 49 % of the boreal forests of northeastern Alberta will be invaded by 2059. The impacts of invasive earthworms on soil morphology are similar in temperate and boreal forests and include the following: (1) thinning of the forest floors and (2) development of novel Ahu surface mineral horizons, characterized by higher pH, clay, and OC content (Bohlen et al., 2004; Lejoly et al., 2021; Lyttle et al., 2011). However, as N limitations for microbial activity are greater in boreal forests (Högberg et al., 2017), the influence of invasive earthworms on microbial communities could be much larger than in temperate forests.

Fungi, including free-living saprotrophs and root-associated symbionts such as mycorrhizal fungi (Lladó et al., 2017), make up a large fraction of biomass in boreal forest soils (Clemmensen et al., 2015). Most northern forest tree species form symbiotic associations with ectomycorrhizal (EcM) fungi, except for maple (Acer spp.), which associates with arbuscular mycorrhizal (AM) fungi (Brundrett et al., 1990), while Populus spp., often categorized as dual mycorrhizal, are mostly ectomycorrhizal in western Canada (Karst et al., 2021). Previous studies have reported a positive response of mycorrhizal fungi, including both AM (Dempsey et al., 2013; Drouin et al., 2016) and EcM (Jang et al., 2022) fungi, to earthworm invasion in temperate forests. To our knowledge, the response of EcM fungi to earthworm invasion has never been studied in the field for boreal forests (Addison, 2009; Cameron et al., 2012). Because EcM fungi can access C through symbiosis, their growth is mostly limited by N availability, while free-living saprotrophic fungi depend on both soil C and N availability (Högberg et al., 2021, 2017). Saprotrophic fungi are able to decompose a wide range of organic compounds, including lignin and cellulose, with the help of extracellular enzymes (Lebreton et al., 2021). Earthworms, by moving fresh litter down the soil profile (Blouin et al., 2013), could lessen C limitations for saprotrophic fungi and result in their increased relative abundance. However, Crowther et al. (2013) reported that increased grazing by soil fauna (isopods) reduced the abundance of the dominant saprotrophic fungi but increased fungal diversity as competition with other saprotrophic taxa decreased. With so few field studies, the specific impacts of earthworm invasion on soil fungal diversity and abundance are not well known.

Similarly to fungi, bacteria and archaea decompose plant and soil organic macromolecules, but they are also involved in a wider range of ecosystem processes, including N fixation and transformations (Lladó et al., 2017). While vegetation composition is less important than for fungi (Baldrian, 2017), pH is a key driver of bacterial community composition (Lladó et al., 2017; Fierer and Jackson, 2006). Earthworm gut transit alters bacterial communities both structurally and functionally (Pedersen and Hendriksen, 1993; Wang et al., 2021; Medina-Sauza et al., 2019). Because of an overall increase in nutrient availability following earthworm invasion, bacteria would be favoured over fungi, the latter being characterized by slower growth and typically decomposing more complex organic compounds (McLean et al., 2006; Soares and Rousk, 2019). In line with this statement, decreased fungi : bacteria ratios have been reported following earthworm invasion in temperate forests (Dempsey et al., 2011), although not systematically (Chang et al., 2017), and remain to be investigated for boreal forests.

Soil microbial communities are commonly characterized using either DNA metabarcoding or phospholipid fatty acid (PLFA) analysis. These techniques are considered to be complementary as they show different sensitivities to land use and other environmental changes (Orwin et al., 2018). Both methods have been successfully used to study the impacts of endemic and exotic earthworms on soil microbial communities (Gómez-Brandón et al., 2010; Butenschoen et al., 2007; Chang et al., 2016; Price-Christenson et al., 2020; de Menezes et al., 2018; Dempsey et al., 2013). Although the DNA and PLFA methods used in combination may provide more robustness to the analysis and conclusions of a study, this approach has rarely been used in past studies assessing earthworm-induced changes to soil bacterial and fungal communities (Shao et al., 2019).

Here, we characterized soil microbial communities of the major (hemi)boreal forest types of North America, invaded and non-invaded by earthworms, using both metabarcoding and PLFA analysis. Besides the presence or absence of earthworms, our sampling design ensured that other soil forming factors were kept identical, therefore minimizing differences between invaded and non-invaded areas of the same site. Our objective was to identify shifts in bacterial, archaeal, and fungal communities associated with earthworm invasion separately for forest floors and mineral soils considering their intrinsically different physicochemical properties and distinct microbial communities (Prescott and Grayston, 2013). We combined two molecular methods (DNA and PLFA analyses) to detect differences in the relative abundances of both functional and taxonomical groups of microorganisms, and we used PLFA data for quantitative analysis of broad microbial groups. We hypothesized that earthworm invasion would shift microbial community composition in both mineral soils and forest floors. For fungal communities, we expected no changes in the ratio of ectomycorrhizal to saprotrophic fungi as earthworm invasion would positively impact them both. Despite differences in soil type and vegetation among sites, we expected shifts to be similar across sites, indicative of increased nutrient availability and pH, including lower fungi : bacteria ratios.

2 Materials and methods

2.1 Sample collection

Sites were selected from the most common soil types and vegetation covers of the Canadian boreal forest, including Brunisol (Soil Classification Working Group, 1998) under sugar maple (Acer saccharum Marshall – Valcartier), Podzol under black spruce (Picea mariana (Mill.) Britton, Sterns, and Poggenburg – Grands Jardins), and Luvisol under trembling aspen (Populus tremuloides Mich. – EMEND (Ecosystem Management Emulating Natural Disturbance) site; Table S1 in the Supplement). The sites used in the present study are a subset of those described in Lejoly et al. (2021) and correspond to sites with lower anthropogenic impacts. Although A. saccharum is not recognized as a boreal forest tree species, it was included in the study as it is often found in the temperate–boreal ecotone in hemiboreal forests. As the home range of A. saccharum is shifting northwards with climate change, sugar maple forests will eventually establish themselves in regions currently covered by boreal forest tree species, negatively impacting soil C storage (Boilard et al., 2023).

The presence of earthworms was determined in 2018 by hand-sorting of the litter and hot-mustard extraction from the mineral soil surface (Lawrence and Bowers, 2002) and from clear signs of earthworm activity such as the presence of surface casts or extensive bioturbation. Within each site, it was possible to identify earthworm-invaded and earthworm-free patches. At Valcartier and Grands Jardins, all three functional groups of earthworms (epigeic, endogeic, and anecic) were found, indicating a later stage of invasion, while only one epigeic species (D. octaedra) was observed at the EMEND site, suggesting a relatively recent invasion (Lejoly et al., 2021; Table 1). As demonstrated in our previous study, the sampling points were identified as differing only in terms of the presence of earthworms, keeping other soil forming factors (slope, aspect, vegetation, parent material) as similar as possible (Lejoly et al., 2021). Earthworm presence resulted in thinner forest floors and reworked surface mineral soils (Lejoly et al., 2021), matching previous findings for invaded North American temperate forests (Bohlen et al., 2004; Lyttle et al., 2011; Langmaid, 1964; Alban and Berry, 1994), for invaded Fennoscandian arctic forests (Wackett et al., 2018), and for the typical mull form of humus observed in temperate forests that developed in the presence of earthworms (Muys and Granval, 1997; Beyer et al., 1991). For Valcartier, the effects of earthworm invasion on the forest floor varied, ranging from minimal changes (no significant differences in terms of thickness and C stocks compared to the control but a change from mor to moder) to more drastic changes where the forest floor had almost disappeared.

At each site, three to six sampling points were selected for both earthworm-invaded and earthworm-free zones. The sampling points were randomly distributed within these zones, and sampling occurred once for each location over the months of June and July 2019. As recommended by Ferlian et al. (2018), we separately sampled the forest floors and the top layers (0–10 cm) of the mineral soils to independently test for the impact of earthworm invasion on these two different soil layers. All the samples were stored at −20 ^∘C for transportation and then at −80 ^∘C prior to freeze-drying.

Table 1 Earthworm species found at each site after hand-sorting of the litter and hot-mustard extraction for the mineral soil. Data from Lejoly et al. (2021).

2.2 DNA extraction and sequencing

Freeze-dried soil samples were ground on a TissueLyser II (Qiagen, Hilden, Germany). Approximately 250 mg of ground sample was used for DNA extraction using the DNeasy PowerSoil DNA Isolation Kit (Qiagen, Hilden, Germany). Because of its stable nature, the extracted DNA also included non-viable microorganisms (Frostegård et al., 2011). For bacteria and archaea, primers 515F-Y (5^′-GTGYCAGCMGCCGCGGTAA-3^′) and 926R (5^′-CCGYCAATTYMTTTRAGTTT-3^′; Parada et al., 2016) were used for polymerase chain reaction (PCR) amplification of the V4–V5 regions of the 16S rRNA gene, while for fungi, primers ITS9F (5^′-GAACGCAGCRAAIIGYGA-3^′) and ITS4R (5^′-TCCTCCGCTTATTGATATGC-3^′; White et al., 1990) were used for PCR amplification of the internal-transcriber-spacer-2 (ITS2) region (Integrated DNA Technologies, Coralville, IA, USA). Amplification was carried out following the Platinum SuperFi Green PCR Master Mix protocol (Invitrogen, Carlsbad, CA, USA), with 30 PCR cycles and an annealing temperature of 60 ^∘C for ITS2 and 55 ^∘C for 16S. PCR products were purified using Sera-Mag Select magnetic beads (Cytiva, Marlborough, MA, USA). In preparation for sequencing, a second index PCR was performed with oligonucleotides containing the Illumina overhang adapters. PCR products were again purified as described above. All purified PCR products from both amplified regions were then pooled using 3–6 µL of each sample to build the amplicon library. DNA concentration and average size were determined using the Qubit dsDNA HS Assay Kit on a Qubit fluorometer (Life Technologies, Carlsbad, CA, USA) prior to running an HS DNA Bioanalyzer (Agilent Technologies, Wilmington, DE, USA) to normalize the pooled library to 4 nM before sequencing. The pooled libraries were sequenced on an Illumina MiSeq platform using a MiSeq Reagent v3 600-cycles Kit (Illumina Inc., San Diego, CA, USA) at the University of Alberta Molecular Biological Sciences Unit.

2.3 Bioinformatic analyses

All bioinformatic analyses were performed in “Quantitative insights into microbial ecology 2” (QIIME2 version 2021.4 software; Bolyen et al., 2019). Raw sequence data were demultiplexed and quality filtered using the q2-demux plugin followed by de-noising with DADA2 (Callahan et al., 2016). After trimming off the first 20 base pairs (bp) to remove primers, sequence reads were truncated where the average quality score dropped below 34 (at 285 and 215 bp for forward and reverse 16S reads and at 260 and 225 bp for forward and reverse ITS reads) and were de-replicated with paired-end settings to generate amplicon sequence variant (ASV) tables containing read counts. Using ASVs over operational taxonomic units results in finer-resolution (single nucleotide) differences between sequences (Callahan et al., 2017). However, this approach tends to inflate diversity measurements, with these being 1.3 and 2.1 times higher for fungal and bacterial richness, respectively (Glassman and Martiny, 2018). The mean amplicon size was 316 ± 37 bp, ranging from 240 to 433 bp, and 372 ± 2 bp, ranging from 350 to 426 bp, for ITS and 16S ASVs, respectively. After training the naïve Bayes classifiers (q2-feature-classifier), taxonomy was assigned to ASVs using classify-sklearn (Bokulich et al., 2018) with SILVA 138 SSURef NR99 full-length (99 %) and UNITE 8.3 (97 %) databases for 16S and ITS2, respectively (Quast et al., 2013; Abarenkov et al., 2010). Only prokaryotic (bacterial and archaeal) or fungal ASVs were kept, and the ASV tables generated by the QIIME2 were imported into R (version 4.0.5; R Core Team, 2022) for further analysis. Based on rarefaction curves, a total of 4437 and 2832 sequences for prokaryotic and fungal ASVs, respectively, were randomly selected for analysis of alpha diversity.

For fungal ASVs, functional guilds were determined with the FUNGuild database using the FUNGuildR R package (Nguyen et al., 2016; http://github.com/brendanf/FUNGuildR, last access: 24 April 2023). Taxa associated with a unique guild ranked as probable or highly probable were selected and grouped into three categories: ectomycorrhiza, pathogens (including parasites), and saprotrophs. Taxa associated with more than one guild were manually checked and added to the above categories, when suitable, after carefully reviewing the literature, as recommended by Tedersoo et al. (2022). As general ITS primers are not suitable for AM fungi, they were not included in functional guild analysis (Stockinger et al., 2010). The proportion of each fungal guild was calculated as the frequency of reads assigned to a particular guild divided by the total reads counted across all three guilds. Sequence data were archived at the National Center for Biotechnology Information – Sequence Read Archive (NCBI-SRA) (BioProject PRJNA850095).

2.4 Phospholipid fatty acid (PLFA) analysis

Polar lipids were extracted from freeze-dried samples (3 g of mineral soil and 0.25 to 0.5 g of forest floor) following a modified Bligh and Dyer extraction method (Quideau et al., 2016). Prior to extraction, the 1,2-dinonadecanoyl-sn-glycero-3-phosphocholine (19:0 $/$ 19:0) nonadecanoate surrogate standard (Avanti^® Polar Lipids Inc., Alabaster, AL, USA) was added to determine the final recovery. Extracted phospholipids were purified on solid-phase extraction columns (Agilent Technologies, Wilmington, DE, USA) and methylated in a mild alkaline environment to obtain fatty acid methyl esters (FAMEs). FAMEs were analyzed on an Agilent 6890 series capillary gas chromatograph (GC) (Agilent Technologies, Wilmington, DE, USA) equipped with a 25 m Ultra 2 (5 %-phenyl)-methylpolysiloxane column, a flame ionization detector (Hewlett Packard, Santa Clara, CA, USA), and He as the carrier gas. PLFA concentrations were determined against the instrument standard methyl decanoate Me10:0 (Aldrich, St. Louis, MO, USA), added prior to GC identification and quantification (Quideau et al., 2016). Results were calculated as a concentration per gram of soil (nmol g⁻¹) using the surrogate standard. PLFA identification was performed using the Sherlock Microbial Identification System version 4.5 software (MIDI, Inc., Newrak, NJ, USA) following standard nomenclature (Maxfield and Evershed, 2014). All unsaturated PLFAs were in cis configuration. PLFAs with <14 and >20 C chain lengths are considered of non-microbial origin and were therefore removed before statistical analysis.

Because interpreting specific PLFAs as indicators of specific microbial groups must be done with caution (Frostegård et al., 2011), we only used individual PLFAs widely recognized in the literature as markers of broad microbial groups: 15:0 as general bacteria; i14:0, i15:0, a15:0, i16:0, a16:0, i17:0, and a17:0 for Gram-positive(+) bacteria; 16:1ω7, 18:1ω7, cy17:0, and cy19:0 for Gram-negative(−) bacteria; 18:2ω6 as a fungal marker; and 20:4ω6, 20:5ω3, and 20:1ω9 for eukaryotes, including fungi, as well as micro- and mesofauna such as protists, nematodes, and microarthropods (Stromberger et al., 2012; Zelles, 1999; Saetre, 1998; Orwin et al., 2018).

2.5 Microbial ratio calculations

We calculated PLFA ratios corresponding to broad taxonomic groups: the fungi : bacteria ratio as the sum of fungal PLFAs divided by the sum of Gram(+), Gram(−), and general bacterial PLFAs and the Gram(+) : Gram(−) ratio as the sum of Gram(+) PLFAs divided by the sum of Gram(−) bacteria, with the latter ratio being negatively correlated with carbon availability (Fanin et al., 2019). In addition, we examined potential changes in PLFAs that have been reported to vary in response to environmental stress, specifically changes in the membrane composition of microorganisms (Watzinger, 2015). The cyclo ratio was calculated as the sum of cy19:0ω9 and cy19:0ω7 divided by 18:1ω7, with an increase in cyclo-PLFA indicating a response to pH, osmotic, and/or thermal stresses (Yang et al., 2015; Guillot et al., 2000; Mykytczuk et al., 2010). The 10-methyl (10Me) ratio was obtained by dividing 10Me16:0 by 16:0 considering the fact that an increase in 10-methyl branching is linked to an increase in membrane fluidity, which is associated with lower temperatures (Poger et al., 2014).

For bacterial ASVs, the Proteobacteria : Acidobacteriota ratio was calculated as the ratio of the number of reads associated with those two phyla, with higher values typically being associated with a higher nutrient status (Orwin et al., 2018). For fungal ASVs, the ectomycorrhizal : saprotrophic ratio was also calculated to assess the relative response of these two fungal guilds to earthworm invasion; for instance, this ratio will decrease if earthworm invasion negatively affects EcM fungi and/or positively affects saprotrophic fungi.

2.6 Statistical analyses

All statistical analyses were performed in R version 4.0.5 (R Core Team, 2022). The threshold for significance was set at α = 0.1 to account for the higher probability of a type-two error associated with the low sample size. In addition to the unfortunate loss of three samples, we recognize that regional studies such as ours necessarily have low replication. Samples (forest floor and mineral soil) were divided into two categories: invaded and non-invaded, corresponding to the factor “invasion”. Because of a significant three-way interaction, all statistical analyses were conducted separately for forest floors and mineral soils.

Relative abundances of fungal and bacteria taxa at different taxonomical levels were calculated as the number of reads of the target category divided by the total number of reads. Two indices were calculated for alpha diversity based on rarefied ASV tables: species richness and effective number of species, corresponding to a transformed Shannon diversity index (Jost, 2006). For fungi, those indices were calculated for each functional guild, as well as for all fungal ASVs, while for bacteria, they were only based on all ASVs.

To determine whether earthworm invasion affected microbial communities functionally and taxonomically, two-way (factors: invasion and site) analyses of variance (ANOVAs) were performed separately for forest floor and mineral soil samples on fungal and bacterial taxa, diversity indices, and microbial – bacteria and PLFA – ratios presented in the previous section. The ANOVAs were run on bacterial phyla and families present in at least 60 % of the samples and accounting for >0.5 % of average relative abundance, as well as on dominant (>1 % of average relative abundance) fungal phyla and fungal guilds (as determined with FUNGuild). When necessary, data were first transformed with Tukey's ladder of powers using the transformTukey function from the rcompanion package to ensure normality of residuals (Package “rcompanion”, Mangiafico, 2023). The homogeneity of variance was checked with the Bartlett test. To account for the unbalanced design for both invasion and site, a type-II sum of squares was used for the ANOVAs using the car R package (Fox and Weisberg, 2019).

To determine whether earthworm invasion altered soil microbial community composition, permutational analyses of variance (PERMANOVAs) were performed on all three datasets (PLFA, 16S, and ITS) after a Hellinger transformation using Bray–Curtis distance matrices for the forest floors and mineral soils separately with the adonis function from the vegan package (Oksanen et al., 2022). Additionally, we ran indicator species analyses on Hellinger-transformed data for the forest floors and mineral soils separately to identify specific taxa associated with invaded or control samples using the multipatt function from the indicspecies R package (De Caceres and Legendre, 2009).

3 Results

3.1 Earthworm invasion effects on fungi

A total of 1 735 408 reads were obtained for fungal ASVs, ranging from 70 to 136 741 reads with an average of 34 028 reads per sample. After quality control and trimming, we obtained a total of 10 761 ASVs. On average, 80 % of unique ASVs belonged to at least one fungal guild. The following numbers of ASVs were identified for each guild: ectomycorrhizal (EcM; 208), saprotrophic (450), and pathogenic (126). Before rarefaction, 8 out of 51 samples were removed from the calculations of diversity indices because of low read counts: 3 invaded forest floors, 3 invaded mineral soils, and 2 control mineral soils. The relative abundances of fungal phyla, classes, orders, and families are reported in Figs. S1 to S4 in the Supplement.

In the forest floors, neither the fungal community composition (Table 2 and Fig. S5) nor the relative abundance of the dominant fungal phyla and guilds (Fig. 1 and Tables S2 and S3) differed with earthworm invasion status. However, the indicator species analysis revealed that the EcM genus Lactarius was associated with control forest floors (Table 3). Species richness and diversity were not affected by earthworm invasion in forest floors (Table 4). In the mineral soils, fungal community composition differed by earthworm invasion status, and the interaction with the site was marginally significant (Table 2 and Fig. S5). The post hoc test revealed that the differences between invaded and non-invaded mineral soils were significant for the Luvisol but not for the Podzol nor the Brunisol (Table 2). The relative abundance of Ascomycota, the most abundant fungal phylum, was significantly lower in earthworm-invaded mineral soils (Table S2). The relative abundance of dominant fungal guilds was also affected by earthworm invasion, with lower relative abundance of saprotrophic fungi (40 % compared to 60 % in control soils) and higher relative abundance of EcM fungi (58 % compared to 35 % in control soils; Fig. 1 and Table S3). As a result, the fungal communities shifted from saprotroph- to EcM-dominated, and the ectomycorrhizal : saprotrophic ratio increased from 1.25 to 2.42 (Table 5). However, diversity and species richness remained higher for saprotrophic fungi (16 and 30, respectively) than for EcM fungi (7 and 17, respectively).

Figure 1 Average (±1 SE) relative abundances (as number of reads) of fungal guilds. Different letters indicate significant differences (p value < 0.1) between control and earthworm-invaded forest floors or mineral soils and were obtained with the post hoc honestly significant difference (HSD) Tukey test after two-way analysis of variance (n = 7–18). Results of individual sites can be visualized in Fig. S6.

Table 2 Permutational analysis of variance performed on Hellinger-transformed data for fungal and bacterial amplicon sequence variants and phospholipid fatty acids (PLFAs) using Bray–Curtis distance matrix for earthworm-invaded (EW) and non-invaded (CONT) forest floors (LFH) and mineral soils (MIN). Df stands for degrees of freedom, and Sum sq stands for sum of squares. Post hoc pairwise comparisons between sites are also displayed when the interaction factor was non-significant. When the interaction was significant, the pairwise comparisons between invaded and control soils within each site are displayed instead.

The indicator species analysis identified the EcM genera Amphinema and Tomentella as indicators of invasion in mineral soils (Table 3). Fungal species richness and diversity were higher in invaded mineral soils for EcM (+104 % and 53 %, respectively), saprotrophic (+71 % and 61 %, respectively), and pathogenic (+150 % and 145 %, respectively) fungi, while only fungal species richness was higher when considering all fungi regardless of their functional guild (Tables 4 and S3).

Table 3 Indicator species for ectomycorrhizal and saprotrophic fungi, as well as bacteria, associated with control or earthworm-invaded forest floors and mineral soils. Indicator species were obtained using the multipatt function in the indicspecies R package on Hellinger-transformed data. The amplicon sequence variants (ASVs) were grouped at the finest identified taxonomic level between order and genus; taxa only identified to the class level were excluded. The statistic is a combination of specificity (A = 1 when only found in the one group) and fidelity (B = 1 when found in all samples of the one group). Relevant indicators were selected with B>0.5 and relative abundance >0.05 % in the taxonomical group for which they were selected as an indicator. No indicator species were found for pathogenic fungi.

Table 4 Averages (± 1 SE) of fungal and bacterial species richness and diversity (as the effective number of species) indices (n = 7–15). Different letters represent significant differences between control and earthworm invaded soils (p value < 0.1) and are presented separately for the forest floors and mineral soils.

Table 5 Averages (±1 SE) of selected fungal (n = 7–18), bacterial (n = 7–18), and phospholipid fatty acid (PLFA; n = 8–19) ratios. Different letters represent significant differences between control and earthworm-invaded soils (p value < 0.1) and are presented separately for the forest floors and mineral soils. The 10Me ratio was obtained by dividing the PLFA 10Me16:0 by 16 : 0; the cyclo ratio was calculated by dividing the sum of cy19:0ω9 and cy19:0ω7 by 18:1ω7.

3.2 Earthworm invasion effects on bacteria and archaea

For bacterial ASVs, a total of 2 302 500 reads was obtained, with an average of 45 147 per sample, ranging from 35 530 to 54 272 reads. After quality control and trimming, we obtained a total of 17 086 ASVs. While the term “bacterial communities” is used here, archaea were not excluded. The relative abundances of bacterial phyla, classes, orders, and families can be found in Figs. S7 to S10.

In the forest floors, bacterial community composition did not differ by earthworm invasion status (Table 2 and Fig. S11). However, specific taxonomic groups were affected, with the relative abundance of Actinobacteriota being significantly lower in invaded soils (24 %) than in the control soils (32 %; Fig. 2 and Table S4). Verrucomicrobiota increased significantly in earthworm-invaded forest floors, reaching 3 % of relative abundance, mainly due to the increase in the Chthoniobacteraceae family. Within the Bacteroidota phylum, the Chitinophagaceae family was also positively affected by earthworm invasion, with the genus Parafilimonas identified as an indicator of invasion in the forest floor (Table 3).

In the mineral soils, earthworm invasion altered the composition of bacterial communities according to the permutational analysis of variance (PERMANOVA) results (Table 2 and Fig. S10). Notably, the relative abundance of Acidobacteriota significantly decreased from 27 % to 21 % following earthworm invasion, resulting in a higher Proteobacteria : Acidobacteriota ratio (Fig. 2 and Table 5). While Actinobacteriota were not affected by earthworm invasion at the phylum level, the Gaiellaceae family responded positively. The Chitinophagaceae family was also positively affected by earthworm invasion in the mineral soil, mirroring the significant increase in Bacteroidota. Although their relative abundance was much lower (<10 %), it is still worth noting that the Gemmatimonadota were positively affected by earthworm invasion in the mineral soil. Within the Chloroflexi, the uncultured bacteria 1921-2, 1921-3, and B12WMSP1 were indicators of non-invaded mineral soils. The Hyphomicrobiaceae and Nitrosomonadaceae families were positively affected by earthworm invasion. Similarly, the relative abundance of Pirellulaceae was significantly higher in earthworm-invaded mineral soils, of which Pirellula was selected as an indicator together with two other Planctomycetota, namely Fimbriiglobus and the uncultured bacterium OM190 (Table 3). Lastly, these shifts were associated with higher species richness and diversity in invaded mineral soils (Table 4).

Figure 2 Heatmap of changes in relative abundance (as number of reads) of bacterial phyla and families. Only significant differences (p value < 0.1), determined by two-way analyses of variance (ANOVAs) followed by post hoc HSD Tukey tests, are displayed (n = 7–18). Orange and red hues represent a significantly higher relative abundance, and blue hues represent a significantly lower relative abundance in earthworm-invaded soils compared to controls. Results of the two-way ANOVAs associated with the significant differences can be found in Table S3, and individual sites can be visualized in Fig. S12.

3.3 Earthworm invasion effects on microbial and faunal PLFAs

Results of the permutational analysis of variance (PERMANOVA) showed that invasion did not significantly affect the overall PLFA composition of the forest floors and mineral soils (Table 2 and Fig. S13). While the total PLFA biomass was not affected by earthworm invasion, the fungal PLFA 18:2ω6 was significantly higher in earthworm-invaded mineral soils, with a concentration almost twice as high as in the control soils of the same soil layer (Table 6 and Fig. 3). Consequently, the fungi : bacteria ratios were slightly but significantly higher in earthworm-invaded mineral soils (0.11) compared to the controls (0.09) as the concentrations of bacterial PLFAs were not affected by earthworm invasion. Additionally, the concentration of eukaryotic PLFAs (20:4ω6, 20:5ω3, and 20:1ω9) was higher in earthworm-invaded mineral soils. Also in the mineral soils, the Gram(+) : Gram(−) ratio was significantly lower after earthworm invasion, as were the 10Me and cyclo ratios, indicative of environmental stress (Tables 5 and S5). In earthworm-invaded forest floors, the observed changes are limited to higher concentrations of Gram(+) bacteria PLFA i14:0 and a17:0, and of the eukaryotic PLFA 20:4ω6, while PLFA ratios are not affected by earthworm invasion. For all PLFA concentrations and ratios, the interaction between invasion and site was non-significant (Table S5).

Figure 3 Average (±1 SE) concentrations (nmol per gram of soil) of major phospholipid fatty acid groups. Different letters indicate significant differences (p value < 0.1) between control and earthworm-invaded forest floors or mineral soils and were obtained with the post hoc HSD Tukey test after two-way analysis of variance (n = 8–19). Results for individual sites can be visualized in Figs. S14–S15.

Table 6 Average (±1 SE) concentrations (nmol per gram of soil) of fungal phospholipid fatty acids (PLFAs) associated with specific microbial (general bacteria, Gram(+), Gram(−), and fungi) and eukaryotic, including micro- and mesofauna, groups (n = 8–19). Different letters represent significant differences between control and earthworm-invaded soils (p value < 0.1) and are presented separately for the forest floors and mineral soils.

4 Discussion

4.1 Earthworm invasion and shifts in microbial community composition

The effects of earthworm invasion on soil microbial communities were similar across sites spanning the Canadian hemiboreal and boreal forests, regardless of soil type. As expected, the three sites, corresponding to three distinct soil types, harboured different microbial communities, but the interaction between invasion and site was almost always non-significant. Despite potentially variable amplitudes, earthworm invasion was associated with similar shifts in microbial community composition, as illustrated by permutational analyses of variance, as well as in relative and absolute abundances of specific microbial groups obtained from amplicon sequencing and PLFA analysis. An exception to this was found for the fungal community composition, where differences were only significant for the Luvisol.

Our results corroborate previously reported changes in bacterial communities associated with earthworm activity. De Menezes et al. (2018) reported an increase in Verrucomicrobiota, and Gong et al. (2018) reported an increase in Gemmatimonadota after incubation with earthworms in pasture and arable soils, respectively. Chitinophagaceae, involved in the degradation of chitin and other complex sugars such as cellulose and hemicelluloses (Rosenberg, 2014), were also previously positively correlated with earthworm presence (de Menezes et al., 2018). The observed increase in Proteobacteria : Acidobateriota ratios, indicative of higher nutrient status (Smit et al., 2001), has been previously associated with earthworm activity in arable soils (Gong et al., 2018).

While we initially hypothesized that the fungi : bacteria ratios would decrease following invasion, we observed the opposite and measured an overall increase in fungal biomass in mineral soils, estimated by the fungal PLFA 18:2ω6. This should be added to the list of contradictory results previously documented. Dempsey et al. (2011) reported a decrease in fungi : bacteria ratios in a temperate hardwood forest on an area basis (m²) but not per gram of soil, which is how we calculated our ratio. Similarly, in an incubation with arable soil, fungal biomass decreased in the presence of earthworms (Butenschoen et al., 2007). Chang et al. (2017) observed an increase in fungi : bacteria ratios but only in the presence of multiple earthworm species in deciduous temperate forests, consistent with our findings considering that most of our sites were also invaded by multiple earthworm species. Similarly, in a Chinese subtropical forest, the presence of the exotic earthworm Pontoscolex corethrurus, native to South America, also increased fungi : bacteria ratios (Shao et al., 2019).

The analysis of fungal guilds revealed that the positive impact of earthworm invasion on fungal biomass was limited to EcM fungi present in the mineral horizons of the studied boreal forest soils. Previous field studies observed a similar response for EcM and AM fungi following earthworm invasion in temperate forests (Dempsey et al., 2013; Drouin et al., 2016; Jang et al., 2022). In comparison, in a mesocosm experiment, Cameron et al. (2012) reported that the presence or density of invasive earthworms did not affect EcM fungal community composition or colonization of white spruce roots. However, only a subset of EcM fungi able to withstand disturbance was likely to be present in the pots used in this mesocosm study, and the seedlings were grown in conditions not conducive to growth of fine roots, a required process for changes in EcM fungal colonization to occur. Thus, the conditions of Cameron et al. (2012) may have underestimated the effects of earthworms. The observed decrease in the relative abundance of saprotrophic fungi goes against our initial hypothesis. Saprotrophic fungi and bacteria both decompose OM compounds, but bacteria are typically more competitive in degrading simple C compounds (Rousk and Bååth, 2011). An increase in those simple C compounds in the presence of earthworms would explain the observed decrease in the relative abundance of saprotrophic fungi, giving bacteria a competitive advantage. However, as the fungal biomass (estimated by the fungal PLFA 18:2ω6) increased, this might not correspond to an absolute decrease in saprotrophs.

Similarly to fungi : bacteria ratios, there are also mixed responses to earthworm presence in Gram(+) : Gram(−) ratios. In line with our findings, Butenschoen et al. (2007) found that endogeic earthworms favoured Gram(−) over Gram(+) bacteria in arable soils. On the other hand, Dempsey et al. (2013) saw a decrease in Gram(−) and an increase in Gram(+) bacteria in the surface mineral soils after earthworm invasion in a northern hardwood forest. Other studies found variable responses of the relative abundance of Gram(−) and Gram(+) bacteria among earthworm species (Chang et al., 2016, 2017). Taken together, the absence of consensus among studies about the impact of earthworms on the fungi : bacteria and Gram(+) : Gram(−) ratios shows that the relationships are complex and not well understood. There is, however, evidence that pH and C availability, among other factors, can alter the relative abundance of those broad microbial groups (Rousk et al., 2010; Fanin et al., 2019; Zhou et al., 2017).

In addition to soil type and dominant canopy, the three studied sites also differed from one another in terms of earthworm invasion. The Luvisol was only invaded by epigeic earthworms, in high abundance, suggesting that the introduction of invasive earthworms is quite recent, while the Podzol and the Brunisol were invaded by all three functional groups of earthworms, suggesting a later stage of invasion (Lejoly et al., 2021). Additionally, while the Luvisol and the Podzol sites are dominated by EcM-associating tree species, the Brunisol site is dominated by AM-associating sugar maple trees. Because of this design, we cannot determine which factor – soil type, dominant mycorrhizae, or earthworm assemblages – was predominant. However, it is clear that the direction of the changes observed was mostly the same for all three sites, although the amplitude of change might vary.

4.2 Microbial community composition and potential environmental changes

The current invasion of North American forests by exotic earthworms can have cascading effects on ecosystem functioning (Frelich et al., 2019). Although microbial communities in forest floors were relatively insensitive, this soil layer often disappears with earthworm invasion (Lejoly et al., 2021). In contrast, the more pronounced changes occurred in the mineral soil, whose relative relevance to soil nutrient cycling will increase because of the loss of forest floor. Given that soil microbial communities play a central role in nutrient cycling and SOM dynamics (Simpson et al., 2007), the shifts observed in the mineral soil are likely to have far-reaching consequences for nutrient availability and ecosystem functioning.

Although less pronounced than in the mineral soils, taxonomic shifts in the forest floors match the lower TOC and TN content previously observed at these sites (Table S1; Lejoly et al., 2021). Notably, the EcM fungal genus Lactarius, identified as an indicator of the absence of earthworms in the forest floor, is often classified as nitrophilic (Lilleskov et al., 2002, 2011). Rodriguez-Ramos et al. (2021) also found the same taxon to be negatively affected by wildfire and salvage logging, which both alter the integrity of the forest floor, as does earthworm invasion. Verrucomicrobiota, including the Chthoniobacteraceae family, increased in invaded forest floors and are often considered to be oligotrophic (Orwin et al., 2018; Hu et al., 2022; Sun et al., 2017; Ramirez et al., 2012). In comparison, the relative abundance of Actinobacteriota, classified as copiotrophs (Ramirez et al., 2012; Orwin et al., 2018), decreased in invaded forest floors. This is in line with the common view that, by feeding on more palatable substrates, earthworms lead to the accumulation of more recalcitrant C compounds in forest floors (Curry and Schmidt, 2007).

Numerous studies have previously shown that earthworms increase C availability through the incorporation of fresh litter into the soil, priming microbial communities with labile C and thus enhancing the decomposition of more recalcitrant OM (Bohlen et al., 2002; Medina-Sauza et al., 2019; De Graaff et al., 2010; Fontaine et al., 2004). Our results corroborate these observations, with lower Gram(+) : Gram(−) ratios in earthworm-invaded mineral soils, indicating higher C availability (Fanin et al., 2019). Additionally, the lower relative abundance of the oligotrophic Acidobacteriota and Chloroflexi and the higher relative abundance of the copiotrophic Bacteroidota further illustrate that the microbial community is more adapted to higher C availability in earthworm-invaded mineral soils (Fierer et al., 2007; Zeng et al., 2022). Bacteroidota, as well as Chitinophagaceae, are also often positively correlated with high N additions in boreal forests (Högberg et al., 2014), which aligns with the numerically higher TN content of invaded mineral soils (Table S1). Additionally, the EcM fungal genera Aphimena and Tomentella, indicators of earthworm invasion in the mineral soil, have been positively correlated with soil fertility and pH in boreal forests (Sterkenburg et al., 2015; Haas et al., 2018). While we did not measure nutrient availability, the differences in microbial community composition of earthworm-invaded and control mineral soils might indicate a shift in C and N cycling associated with earthworm invasion.

In the mineral soils, the response of the PLFA ratios suggests that earthworm invasion affected environmental conditions for microbial activity. The decrease observed for the cyclo ratio (i.e. the sum of PLFAs cy19:0ω9 and cy19:0ω7 divided by 18:1ω7) has been previously associated with a decrease in environmental stress of osmotic or acidic origin (Yang et al., 2015; Guillot et al., 2000; Mykytczuk et al., 2010), which is consistent with the higher pH found in earthworm-invaded mineral soils (Table S1; Lejoly et al., 2021). The increase in pH can also explain the higher bacterial diversity and richness (Fierer and Jackson, 2006) and the lower Gram(+) : Gram(−) ratios (Frostegård et al., 1993) associated with earthworm-invaded mineral soils. In addition, the lower 10Me ratios observed in invaded mineral soils could indicate a decrease in membrane fluidity, linked to higher temperatures (Poger et al., 2014; Kieft et al., 1994; Zhang and Rock, 2008). The thinning of the forest floor layer observed in earthworm-invaded soils (Lejoly et al., 2021) could increase soil temperature during the growing season, similarly to the mechanical removal of forest floor in managed forests (Tan et al., 2005).

While some of our observations align with previous findings for earthworm invasion in temperate forests, there are some interesting divergences, including the relative overall increase in fungi compared to bacteria and the decrease in Gram(+) : Gram(−) ratios in invaded mineral soils. In boreal forests, fungi : bacteria ratios are typically higher than in temperate forests as fungi have slower growth and lower nutrient requirements compared to bacteria (Soares and Rousk, 2019); fungal communities therefore play a major role in C sequestration (Francisco et al., 2016; Clemmensen et al., 2013; Chen et al., 2020). As EcM fungi can slow down C cycling (Averill and Hawkes, 2016), their increase in invaded mineral soils could positively affect soil C storage in boreal forests, which contain more than twice as much soil organic carbon (SOC) as temperate forests per unit area (Lal, 2005). However, because AM fungi were excluded from this study, we cannot conclusively state that fungal communities became EcM-dominated in all invaded soils. Considering the inherent differences between the two forest ecosystems, extrapolating the findings of temperate forests to boreal forests is therefore not recommended. Our results suggest that hemiboreal and boreal forests are affected by invasive earthworms differently from temperate forests and confirm the importance of implementing additional studies in the boreal biome, including enzyme analyses and/or RNA techniques to target the activity of microbial communities.

5 Conclusions

Using a combined PLFA and metabarcoding approach, our study presented the first evidence of microbial community composition shifts associated with earthworm invasion in hemiboreal and boreal forest soils. Both fungal and bacterial communities were affected, but changes differed between forest floors and mineral soils, with more pronounced shifts in the latter. This result is important because as the thickness and C stocks of forest floors drastically decrease with earthworm invasion (Lejoly et al., 2021), the relative importance of forest floor microbial communities for soil nutrient cycling will likely decrease. In contrast, the microbial communities of earthworm-invaded mineral soils, whose C stocks are not affected by invasion (Lejoly et al., 2021), were quite distinct from those of pre-invasion mineral soils. Earthworm-invaded mineral soils harboured greater fungal biomass, associated with higher relative abundance of EcM fungi and higher microbial diversity. Shifts in bacterial communities included lower Gram(+) : Gram(−) ratios, higher relative abundance of copiotrophic phyla, and lower relative abundance of oligotrophic phyla, all indicating higher nutrient and C availability in earthworm-invaded soils. The similar shifts observed across three intrinsically different soil types of North American hemiboreal and boreal forests constitute strong evidence that earthworm invasion has universally applicable effects on soil microbial community composition in these ecosystems.