Earthworm invasion in North American forests has the potential to greatly impact soil microbiomes by altering soil physicochemical properties. 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 boreal forest using phospholipid fatty acid (PLFA) analysis and metabarcoding of the 16S rRNA gene (bacteria and archaea) and ITS2 region (fungi).

Samples of mineral soils, forest floors, and earthworm casts were collected in the summer of 2019 in multiple boreal forest sites across Canada.

Phospholipid fatty acids were extracted from freeze-dried samples following a modified Bligh and Dyer extraction method. PLFA identification was performed using the Sherlock Microbial Identification System version 4.5 following standard nomenclature. All unsaturated PLFAs were in cis configuration.

DNA was extracted from 250 mg of ground and freeze-dried soil sample using the DNeasy PowerSoil DNA Isolation Kit. For bacteria and archaea, primers 515F-Y (5’-GTGYCAGCMGCCGCGGTAA-3’) and 926R (5’- CCGYCAATTYMTTTRAGTTT-3’) were used for amplification of the V4-V5 regions of the 16S rRNA gene with 30 PCR cycles and an annealing temperature of 55 ºC. For fungi, primers ITS9F (5’-GAACGCAGCRAAIIGYGA-3’) and ITS4R (5’- TCCTCCGCTTATTGATATGC-3’) were used for amplification of the internal transcriber spacer 2 (ITS2) region with 30 PCR cycles and an annealing temperature of 60 ºC. 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. Sequence data were archived at NCBI-SRA (BioProject PRJNA850095).

Raw sequence data were demultiplexed and quality filtered using the q2‐demux plugin followed by denoising with DADA2 in QIIME2. After trimming off the primers, 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; at 260 and 225 bp for forward and reverse ITS reads), and dereplicated with paired end setting to generate amplicon sequence variants (ASVs) tables containing read counts. After training the naïve Bayes classifiers (q2‐feature‐classifier), taxonomy was assigned to ASVs using classify-sklearn with SILVA 138 SSURef NR99 full-length (99 %) and UNITE 8.3 (97 %) databases for 16S and ITS2, respectively. Fungal functional guilds were determined with FUNGuild database. Fungal ASVs associated with a unique guild ranked (highly) probable were grouped into ectomycorrhizal, pathogenic, and saprotrophic fungi. Additional ASVs were manually selected based on the literature and added when relevant.