Ribosomal RNA Genes

The up-to-date reference 16S rRNA gene sequences were maintained as described earlier [3 (link)]. We attempted to select a sequence with the best quality for each species by using the following strategy. For cases in which multiple sequences were available for a type strain, the sequence extracted from its whole-genome assembly (WGA) was selected. As for PCR-derived sequences, the quality of sequencing was checked manually by secondary-structure-aware alignment using the EzEditor program [13 (link)]. Maximum-likelihood phylogenetic trees of each taxonomic group, such as phyla, classes, orders or families, were generated from manually aligned 16S rRNA gene sequences using RAxML software [14 (link)]. All 16S rRNA gene sequences were assigned taxonomically to the species level as a part of the complete taxonomic hierarchy which consisted of phylum, class, order, family, genus and species (subspecies if applicable).

The simulated reads used here were derived from the reference databases using the “Cross-validated classification performance” notebooks in our project repository. The reference databases were either Greengenes or UNITE (99% OTUs) that were cleaned according to taxonomic label to remove sequences with ambiguous or null labels. Reference sequences were trimmed to simulate amplification using standard PCR primers and slice out the first 250 bases downstream (3′) of the forward primer. The bacterial primers used were 27F/1492R [27 (link)] to simulate full-length 16S rRNA gene sequences, 515F/806R [28 (link)] to simulate 16S rRNA gene V4 domain sequences, and 27F/534R [29 (link)] to simulate 16S rRNA gene V1–3 domain sequences; the fungal primers used were BITSf/B58S3r [30 (link)] to simulate ITS1 internal transcribed spacer DNA sequences. The exact sequences were used for cross validation and were not altered to simulate any sequencing error; thus, our benchmarks simulate denoised sequence data [4 (link)] and isolate classifier performance from impacts from sequencing errors. Each database was stratified by taxonomy and 10-fold randomized cross-validation data sets were generated using scikit-learn’s library functions. Where a taxonomic label had less than 10 instances, taxonomies were amalgamated to make sufficiently large strata. If, as a result, a taxonomy in any test set was not present in the corresponding training set, the expected taxonomy label was truncated to the nearest common taxonomic rank observed in the training set (e.g., Lactobacillus casei would become Lactobacillus). The notebook detailing simulated read generation (for both cross-validated and novel taxon reads) prior to taxonomy classification is available at https://github.com/caporaso-lab/tax-credit-data/blob/0.1.0/ipynb/novel-taxa/dataset-generation.ipynb.
Classification performance was also slightly modified from a standard machine-learning scenario as the classifiers in this study are able to refuse classification if they are not confident above a taxonomic level for a given sample. This also accommodates the taxonomy truncation that we performed for this test. The methodology was consistent with that used below for novel taxon evaluations, so we defer its description to the next section.

Individuals were grouped as follows. Group 1 (five D. reticulatum fed a lab diet for 7 days); Group 2 (five D. reticulatum fed a lab diet for 14 days); Group 3 (five D. reticulatum fed a lab diet and infected on Day 7 with P. hermaphrodita with feces collected 7 days postinfection—14 days in total); Group 4 (three A. valentianus fed a lab diet for 7 days); Group 5 (three A. valentianus fed a lab diet for 14 days); Group 6 (three A. valentianus infected with P. hermaphrodita with feces collected 7 days postinfection—14 days in total). Feces were collected from each slug for DNA extraction.
DNA was extracted from feces using DNeasy PowerSoil Pro Kit (Qiagen) following the manufacturer's instructions. The presence of bacterial DNA was checked after extractions using PCR amplification of the hypervariable regions of the 16S rRNA gene. This was carried out using the primers 27f (5′‐AGAGTTTGATCMTGGCTCAG‐3′) and 1492r (5′‐TACGGYTACCTTGTTACGACTT‐3′) (Lane, 1991 ) with the following thermocycler conditions: 3 min at 95°C followed by 35 cycles of 15 s at 95°C, 30 s at 55°C, 1.5 min at 72°C, and a final step of 8 min at 72°C. Amplicons were visualized using agarose gel electrophoresis to confirm that PCRs had worked; in all cases, bands of the correct size were present, and no amplification of bacterial DNA could be seen in the extraction negative control or the PCR negative control.
DNA samples were sent for 16S rRNA metagenomic sequencing (Novogene). The V4 hypervariable region of the 16S rRNA gene was amplified using the primers 515F (5′‐GTGCCAGCMGCCGCGGTAA‐3′) and 806R (5′‐GGACTACHVGGGTWTCTAAT‐3′). All PCR reactions were carried out with Phusion® High‐Fidelity PCR Master Mix (New England Biolabs). Sequencing libraries were generated with NEBNext® Ultra^TM DNA Library Prep Kit for Illumina and quantified via Qubit and Q‐PCR. Libraries were sequenced on an Illumina NovaSeq. 6000 platform to generate 2 × 250 bp paired‐end reads.
Analysis of the raw reads occurred at Novogene using the following method. Paired‐end reads were merged using FLASH (V1.2.7) (Magoč and Salzberg, 2011 (link)). Quality filtering on the raw tags was performed under specific filtering conditions to obtain high‐quality clean tags according to the QIIME (V1.7.0) (Caporaso et al., 2010 (link)). The tags were compared with the reference database (SILVA database) using the UCHIME algorithm (Edgar et al., 2011 (link)) to detect chimera sequences. Detected chimera sequences were then removed to obtain Effective Tags. All Effective Tags were processed by UPARSE software (v7.0.1090) (Edgar, 2013 (link)). Sequences with ≥97% similarity were assigned to the same Operational Taxonomic Units (OTUs).
For each OTU, QIIME (Version 1.7.0) in the Mothur method was performed against the SSU rRNA database of SILVA Database for species annotation at each taxonomic rank (Threshold:0.8~1) (Quast et al., 2012 (link)). MUSCLE (Version 3.8.31) (Edgar, 2004 (link)) was used to obtain the phylogenetic relationship of all OTUs.
OTUs abundance information was normalized using a standard of sequence number corresponding to the sample with the least sequences. OTUs were analyzed for Alpha diversity (Wilcoxon test function) and Beta diversity (AMOVA—Analysis of Molecular Variance) to obtain richness and evenness information in samples. AMOVA was also used to compare the taxonomic compositions of infected and noninfected slugs in weighted PCoA. Analysis of Alpha and Beta diversity were all performed on the normalized data and calculated with QIIME (Version 1.7.0). Significant intragroup variation is detected via MetaStats based on their abundance.

When the bacteria were cultured in TSB medium to OD₆₀₀ 1.0, they were treated in PBS containing 0.4 mM H₂O₂ for 10 min and collected by centrifugation. Total RNA was extracted using RNAprep Pure Cell/Bacteria Kit (Tiangen Biotech, Beijing, China). Then, cDNA was synthesized with 1 μg RNA and PrimeScript™ IV 1st strand cDNA Synthesis Mix (Takara, Dalian, China). For the qRT-PCR, cDNA was mixed with specific primers (Table 2) and chamQ universal SYBR quantitative RT-PCR master mix (Vazyme, Nanjing, China). 16S rRNA gene was used as internal controls.