Microbiome Analysis Reveals Diversity and Function of Mollicutes Associated with the Eastern Oyster, Crassostrea virginica

Tissue samples were dissected from 29 oysters collected from an aquaculture farm in Ninigret Pond, Rhode Island, USA. The 16S rRNA gene-based amplicon sequencing was performed on multiple tissues of 19 oysters, and a microbiome-enriched metagenome was sequenced on the gut samples pooled from the remaining 10 oysters (see Data Set S1A in the supplemental material). Of the 19 oysters, amplicon samples from the pallial fluid and hemolymph were collected from 9 oysters, amplicon samples from the gut, mantle, gill, and inner shell were collected from all, and additional metagenome samples were prepared from the gut of 12 oysters (Materials and Methods). Individual oysters were tested for potential infection by common causative agents of oyster disease in the region, including P. marinus (causative agent of Dermo disease), H. nelsoni (causative agent of MSX disease), and H. costale (causative agent of SSO disease) (Materials and Methods). Overall, 14 of the 29 oysters were infected with P. marinus with a range of infection severities from 0.5 to 5 as measured by the Mackin index (30, 31), while little to no infection from H. nelsoni or H. costale was detected among all 29 samples (Data Set S1A). Due to the low number of individuals infected with H. nelsoni and H. costale, these infections were not considered a factor in subsequent analyses.

The alpha diversity of the oyster microbiome was compared among different tissue types using the unique amplicon sequence variant (ASV) counts, the Shannon index, and the Simpson index across all the profiled oyster samples (Fig. 1). The unique ASV counts revealed a significantly higher diversity in the inner shell samples than in the gut, gill, mantle, and pallial fluid samples (P value < 0.05, Wilcoxon test); the Shannon index revealed a higher diversity in the inner shell than the gut and mantle (P value < 0.05, Wilcoxon test); and the Simpson index showed statistically significant difference only between the inner shell and the mantle (P value < 0.05, Wilcoxon test). A high level of heterogeneity was observed among the different tissues, with the median number of unique ASVs per sample ranging from 188 in the gut to 838 in the inner shell (Data Set S1B). Similarly, when comparing the pooled ASVs from different oyster tissues, a total of 15,011 unique ASVs were identified among all the oyster samples, but only 496 ASVs (3.3%) were present in all six tissue types analyzed (Fig. S1). A high level of heterogeneity was also observed among individual oysters even when the same tissue type was examined, with the gut and mantle representing the lowest percentage (0.6%) and the pallial fluid representing the highest percentage (3.4%) of conserved ASVs (present in >80% of samples for a given tissue type) identified from each tissue type (Data Set S1B). Correlations between P. marinus infection and microbiome alpha diversity were examined within each tissue type. Significant differences were identified only in the inner shell samples, where the uninfected samples were significantly more diverse than the infected samples (P value < 0.05, Wilcoxon test) based on all three alpha diversity measures.

The beta diversity among different samples was measured based on phylogenetic isometric log-ratio (PhILR)-transformed ASV counts (Materials and Methods). Statistical analysis of the sample distances suggested that the gut samples were the most distinct of the profiled tissue types, as they were significantly different from all other tissue types (P value < 0.05, pairwise permutational multivariate analysis of variance [PERMANOVA]). The gill, mantle, and inner shell samples were also significantly different from one another, while the pallial fluid and hemolymph samples had no statistically significant difference from other tissue types with the exception of the gut (Data Set S1C). The comparison of P. marinus-infected against uninfected samples within each tissue type revealed no statistical significance in beta diversity (P value < 0.05, PERMANOVA).

Taxonomic assignment of the ASVs obtained from this study revealed several major bacterial classes that had a median relative abundance of greater than 10% within at least one tissue type, including Mollicutes, Chlamydiae, Spirochaetia, Fusobacteriia, and Gammaproteobacteria (Fig. 2). The relative abundance of these bacterial classes was examined among different tissue types using the ANCOM-BC (analysis of compositions of microbiomes with bias correction) approach (32). The gut samples had a significantly higher relative abundance of Mollicutes and Chlamydiae than any other tissue types examined, with an estimated log fold change ranging from 1.6 to 3.6 for the Mollicutes and from 2.0 to 4.6 for the Chlamydiae compared to the other tissues (Fig. 2A and B; Data Set S1D). Similarly, the mantle and gill samples had a significantly higher abundance of Spirochaetia than other tissues, with the mantle having 2- to 3-fold-higher abundance than other tissue types and 1.5-fold-higher abundance than the gill (Fig. 2C; Data Set S1D). The hemolymph had around a 2-fold-higher relative abundance of Mollicutes, Chlamydiae, and Fusobacteriia compared to the mantle (Fig. 2A, B, and D; Data Set S1D). The gill had a significantly lower abundance of Chlamydiae compared to all other tissue types, while the inner shell had a significantly higher abundance of Mollicutes than the mantle and gill and a lower abundance of Gammaproteobacteria than the pallial fluid samples (Fig. 2A and E).

The shotgun metagenomic sequencing of oyster gut microbiomes resulted in the collection of over 2.5 billion raw reads and 2.1 billion quality-filtered reads. Around 1.7 billion (80%) and 7.8 million (0.4%) of the quality-filtered reads were assigned to the oyster host and P. marinus, respectively. A total of 394,986,464 reads that were unmapped to the eastern oyster or P. marinus genomes were coassembled across the two rounds of metagenomic sequencing, resulting in 612,574 unique contigs. At least one nearly complete 16S rRNA gene was detected for each of the five major bacterial classes (along with the Bacteroidetes and Cyanobacteria), ranging from 1,396 to 1,562 bp in length. Six of the full-length 16S rRNA genes were mapped to the amplicon data with full coverage and 100% identity to the corresponding ASVs. These ASV-mapped full-length 16S rRNA genes were taxonomically assigned to the Mollicutes (named Mollicutes-1 and Mollicutes-2), Chlamydiae, Spirochaetia (named Spirochaetia-1 and Spirochaetia-2), and Bacteroidia (Data Set S2A).

The relative abundance of these 16S rRNA genes was probed across different oyster tissues and between P. marinus-infected and uninfected samples using ANCOM-BC and based on mappings to the ASV abundances (Fig. S2). The Mollicutes-1 showed a significant differential abundance (P value < 0.001, ANCOM-BC) between the P. marinus-infected and uninfected samples of the gut and the inner shell, where it was more abundant among the uninfected samples of the gut (log fold change of 4.3) but was slightly less abundant in the uninfected samples of the inner shell (log fold change of −1.2). A significantly higher relative abundance of the Mollicutes-2 (P value < 0.001, ANCOM-BC) and Bacteroidia (P value < 0.05, ANCOM-BC) was also observed in P. marinus-uninfected than the infected pallial fluid samples (log fold change of 5.2 and 2.3, respectively). In contrast, the Chlamydiae had significantly lower relative abundance in the uninfected gill, mantle, and inner shell samples compared to the P. marinus-infected samples (log fold changes of −0.3, −0.2, and −1.1, respectively; P value < 0.001, ANCOM-BC).

Binning of the coassembled contigs produced two MAGs with completeness greater than 80% and contamination less than 2%. Each MAG contained a full-length 16S rRNA gene, with the first from the class Mollicutes (Mollicutes-1) and the second from Chlamydiae. The Mollicutes MAG included 47 contigs with a total length of 0.62 Mb and a GC content of 28.7%. The estimated completeness and contamination were 97.4% and 1.8%, respectively. Similarly, the Chlamydiae MAG included 57 contigs with a total length of 1.08 Mb and a GC content of 41.0%, and it had a completeness of 86.2% and contamination of 0.4%. Besides the Mollicutes and Chlamydiae MAGs, a partial Spirochetes MAG was reconstructed, with a completeness of 43.2%, a contamination of 0.7%, and a full-length 16S rRNA gene (Spirochaetia-1) included in the MAG.

Phylogenomic reconstructions were performed based on conserved single-copy genes (CSCGs) identified between the two nearly complete oyster MAGs and corresponding reference genomes in the Mollicutes and Chlamydiae taxa (Materials and Methods). In total, 34 CSCGs and 179 CSCGs were used to build the Mollicutes and Chlamydiae phylogenies, respectively. Based on the phylogenomic reconstruction, the oyster Mollicutes MAG formed a basal branch to a group of marine Mycoplasma, with the nearest neighboring branches including Mycoplasma todarodis (isolated from a squid), Mycoplasma marinum (isolated from an octopus), and Mycoplasmatales bacterium DT_67 and DT_68, two MAGs obtained from deep-sea sinking particles (herein referred to as DT_67 and DT_68) (33) (Fig. 3A). Meanwhile, the Chlamydiae MAG was found within the order Parachlamydiales most closely related to Simkania negevensis, an obligate intracellular bacterium with a broad host range from amoebae to animals (34, 35) (Fig. S3).

To determine the similarity of the Mollicutes and Chlamydiae MAGs to genomes within their corresponding taxa, average nucleotide identity (ANI) and average amino acid identity (AAI) were calculated between the MAGs and all other species represented in their corresponding CSCG phylogenies (Data Set S2B and C). For the Mollicutes MAG, the highest pairwise ANI values were observed with M. marinum (66.8%), DT_67 (66.8%), DT_68 (66.4%), and M. todarodis (66.0%), all representing strains from marine sources, as well as Mycoplasma mobile (66.1%), a freshwater fish pathogen. Similarly, some of the highest pairwise AAI values were observed with M. marinum (51.9%) and M. todarodis (50.3%) with orthologs identified from 57.4% and 53.3%, respectively, of the protein-coding genes in the oyster Mollicutes MAG. While DT_67 and DT_68 had slightly higher AAI values of 54.4% and 53.6%, respectively, only a limited number of orthologs were identified, covering 22.5% and 22.2% of the protein-coding genes in the Mollicutes MAG, respectively (Data Set S2B). The oyster Chlamydiae MAG had an ANI between 62.8% and 64.3% to all reference genomes of the Chlamydiae phylum. The highest AAI value (47.9%) was observed with its nearest neighbor in the phylogenomic reconstruction, Simkania negevensis, with ortholog identifications for 53.6% of the protein-coding genes in the Chlamydiae MAG (Data Set S2C). Overall, both the Mollicutes and Chlamydiae MAGs had lower ANI and AAI values to the reference genomes than expected for strains of the same species (both generally ≥95%) (36, 37).

A genome-scale metabolic reconstruction was performed to illustrate the metabolic capacity encoded by the oyster Mollicutes MAG (Fig. 3B). Overall, the Mollicutes MAG encoded a highly reduced metabolism with few carbon utilization pathways and limited biosynthetic capability. Transport via the phosphotransferase system (PTS) was predicted for glucose, fructose, mannose, and N-acetylglucosamine (GlcNAc). A complete glycolysis pathway and all subunits of ATP synthase were present, while genes of the tricarboxylic acid (TCA) cycle were largely missing. The ability to convert pyruvate to fumarate was predicted based on the presence of a malate dehydrogenase and a fumarate hydratase. The conversion of acetyl coenzyme A (acetyl-CoA) to acetate was inferred by the identification of a phosphate acetyltransferase and an acetate kinase. An arginine deiminase (ADI) pathway was identified, potentially contributing to the production of ATP via the exchange of arginine and ornithine (38) and enabling the production of precursors for a de novo pyrimidine biosynthesis pathway. The Mollicutes MAG also encoded a putative chitinase with the catalytic domain homologous to an experimentally verified chitinase in Pyrococcus furiosus (PDB: 2DSK) and a chitin-binding domain at the C terminus. It was noted that a mechanism of converting pyruvate to acetyl-CoA was missing in the MAG. However, putative subunits of a pyruvate dehydrogenase complex were identified from other contigs in the metagenomic coassembly, with top BLAST hits to members of the Mycoplasma genus. Therefore, the potential conversion of pyruvate to acetyl-CoA was speculated in the metabolic reconstruction (Fig. 3B).

A comparative genomic analysis was performed between the oyster Mollicutes MAG and genomes from closely related Mycoplasma isolates, including two marine species, M. marinum and M. todarodis, as well as one freshwater species, M. mobile (Data Set S2D). Genes in the oyster Mollicutes MAG can be largely classified into four categories based on their conservation with other species from the comparative genomic analysis (Fig. S4). Over 40% of the genes in the Mollicutes MAG were conserved among all reference genomes, encoding functions in central metabolism (e.g., ATP synthase, amino acid metabolism, nucleotide metabolism, etc.) and information storage and processing (e.g., replication, translation, transcription). Around 10% of genes were conserved between the Mollicutes MAG and the two marine isolates but were missing from M. mobile. These included genes belonging to a de novo pyrimidine biosynthesis pathway, multiple carbohydrate utilization genes, and one gene that encodes an alpha-amylase domain-containing protein. Only 18 genes (less than 3%) were uniquely conserved between the Mollicutes MAG and M. mobile, including an arginine deiminase and genes associated with DNA repair and the repair of enzymes under oxidative stress. Finally, around 36% of the genes were unique to the Mollicutes MAG. These included genes involved in pyruvate metabolism (e.g., malate dehydrogenase and fumarate hydratase) and the ADI pathway (e.g., carbamate kinase, ornithine carbamoyltransferase, and the arginine/ornithine antiporter). Overall, the oyster Mollicutes MAG maintained conserved functions with known Mycoplasma species while carrying unique metabolic capability in its genome.

Following functional inferences obtained from the metabolic reconstruction of the oyster Mollicutes MAG, emphasis was placed on the composition and diversity of Mollicutes across different oyster microbiomes obtained from this and other studies (Fig. S5). A total of three additional data sets were identified from public databases, including the community profiling of gut samples from the adult eastern oyster (39), the profiling of adult eastern oyster biodeposits (40), and the profiling of eastern oyster larvae homogenate (25). Raw data from these prior studies were downloaded and reanalyzed with a standardized pipeline to minimize potential inconsistencies caused by variations in the computational procedures (Materials and Methods). According to the comparative analysis, the relative abundance of Mollicutes was over 5 and 8 log fold higher in the adult oyster gut than in adult biodeposits and larvae homogenate, respectively (ANCOM-BC, P value < 0.05). While a significantly higher abundance of Mollicutes was observed in the biodeposits than the larvae homogenate, the log fold change (1.7) was much lower than what was observed between the gut and larva samples. Interestingly, another bacterial class, the Spirochaetia, was also significantly enriched in adult oyster gut compared to the biodeposits (log fold change up to 2.9). In contrast, the Gammaproteobacteria were significantly more abundant in the larvae homogenate than the adult gut (log fold change between 2.7 and 3.0) and the biodeposits (log fold change of 3.9).

To further elucidate the taxonomic distribution of the Mollicutes ASVs recovered from oyster tissues and the surrounding water, a phylogenetic analysis was performed by placing the Mollicutes ASVs from this study and a previous study of the adult eastern oyster gut microbiome (39) into a reference tree of 16S rRNA gene sequences from the SILVA database, derived from environmental clones and organisms in pure culture (Materials and Methods). Branches in the tree were collapsed with labels indicating the dominant environmental source of the SILVA references. A bubble plot was included on the right of the phylogeny to indicate the number of unique Mollicutes ASVs identified from different tissue or water samples (Fig. 4).

The Mollicutes ASVs identified from oyster tissues were largely divided into four major clades of the phylogenetic tree (Fig. 4 and Data Set S2E): (i) Sydney rock oyster and water (clade I), including ASVs mapped to a clade of uncultured references taxonomically assigned to the family Mycoplasmataceae; (ii) crustaceans and soil (clade II), including ASVs mapped to a clade containing reference sequences from “Candidatus Bacilloplasma,” “Candidatus Lumbricincol,” and unidentified members of the family Mycoplasmataceae; (iii) water, sediment, and feces (clade III), including ASVs mapped to a clade of references from Izemoplasmatales, Acholeplasmatales, and RF39; and (iv) marine invertebrates (clade IV), including ASVs mapped to a clade of Mycoplasma. The highest number of unique Mollicutes ASVs across all eastern oyster tissues was observed in the clade IV (marine invertebrates), as shown both in this study and by Pierce and Ward (39). This clade also included the metagenome-assembled 16S rRNA genes, Mollicutes-1 and Mollicutes-2, and their corresponding ASVs (Data Set S2E). In contrast, the Mollicutes ASVs identified from surrounding water samples were mainly placed into clade III (water, sediment, and feces) of the phylogeny. Interestingly, the gut, compared to any other tissues, included the highest number of unique Mollicutes ASVs in clade IV but the lowest number of unique ASVs in clade III. Overall, a high level of phylogenetic diversity was revealed within the eastern oyster-associated Mollicutes.

Twenty-nine live oysters, 1 to 3 years in age, were sampled from an aquaculture farm along with 1 liter of surface water in Ninigret Pond, Rhode Island (USA), on 25 August 2017. In this farm, eastern oysters are farmed using a rack and bag culture system, with racks suspended approximately 1.5 m from the bottom of the pond. Environmental conditions at the site were measured as part of another study (15). Oysters were transported to the laboratory in coolers on ice. Upon arrival at the laboratory, the outer shells of individual oysters were scrubbed to remove visible sediments, and the length, width, and mass of each oyster were recorded (see Data Set S1A in the supplemental material). The oysters were sprayed with 70% ethanol followed by the immediate dissection to obtain the pallial fluid, hemolymph, mantle, gill, gut, and inner shell samples. Hemolymph and pallial fluid were collected for 9 of the 29 oysters (fluid volumes recorded in Data Set S1A). The pallial fluid was collected prior to shucking by creating a small notch in the shell and using a sterile syringe and needle through the notch. Following shucking, hemolymph was obtained with a sterile syringe from the adductor muscle sinus (50). Pallial fluid and hemolymph samples were centrifuged at 16,100 × g for 8 min, the supernatant was removed, and the pellet was resuspended in 500 μl of RNAlater (Invitrogen, Carlsbad, CA). For all of the oysters, samples of the gut (consisting of the digestive gland, intestine, and stomach tissue), gill, mantle, and inner shell (taken from the inner surface of the top and bottom shells using sterile cotton swabs) were stored individually in 1 ml of RNAlater (Invitrogen, Carlsbad, CA). A mixed sample of gill, mantle, and rectum tissues was also preserved for each oyster in 100% ethanol for protozoan pathogen quantification (described below). The water sample was prefiltered through a 153-μm mesh filter, followed by filtration through a 5-μm-pore-size (diameter of 47 mm) filter (Sterlitech, Kent, WA) and then a 0.2-μm-pore-size Sterivex filter (Millipore, Burlington, MA). All samples were stored in a −80°C freezer before further processing. Detailed tracking of all samples and their application in the 16S rRNA community profiling and metagenomic sequencing is provided in Data Set S1A.

For each of the 29 oysters, DNA was extracted from the 100% ethanol-preserved gill, mantle, and rectum samples (described above) using a previously described protocol with modifications in order to assess the abundance of the eastern oyster protozoan pathogens Perkinsus marinus (causative agent of Dermo disease), Haplosporidium nelsoni (causative agent of MSX disease), and Haplosporidium costale (causative agent of SSO disease) (51). First, a total of 5 mg sample was pooled from all three tissue types and was sterilely placed into a 1.5-ml microcentrifuge tube with 160 μl of preheated urea buffer and incubated at 60°C in a dry bath for 1 h. Samples were vortexed intermittently throughout the incubation. After incubation, the tubes were vortexed again and placed into a 95°C dry bath for 15 min. Samples underwent centrifugation for 5 min at 15,000 × g, and 100 μl of the supernatant was removed and mixed with 1 μl of 100× Tris-EDTA (TE) buffer, 50 μl ammonium acetate, and 400 μl of 95% ethanol. After mixing, the samples were left at −20°C overnight to precipitate. The following day, samples underwent centrifugation for 20 min at 15,000 × g to pellet the precipitate and were washed three times with 70% ice-cold ethanol to remove any remaining impurities. DNA was eluted in 100 μl of buffer (10 mM Tris-HCl, 0.1 mM EDTA, pH 8.0). Quantity and quality of the DNA were assessed with the NanoDrop 2000c (Thermo Fisher Scientific, Waltham, MA), and the samples were diluted to 100 ng in 3 μl to ensure consistency between quantitative PCRs (qPCRs).

Subsequently, the DNA was used as a template in a multiplex, real-time PCR testing for P. marinus, H. nelsoni, and H. costale. P. marinus primers and dually labeled probes from reference 31 were used in conjunction with MSX/SSO primers, MSX probe, and SSO probe (primer and probe sequences in Table 1). Each reaction with a 20-μl reaction mixture was carried out using 300 nM P. marinus primers, 450 nM Haplosporidium primers, and 75 nM specific probe for each pathogen with 3.55 μl water, 10 μl iQ Multiplex Powermix (Bio-Rad, Hercules, CA), and 3 μl (100 ng) template DNA. The thermal cycler protocol was as follows: 95°C for 60 s, followed by 40 cycles of 95°C for 15 s, then 60°C for 30 s. The qPCR copy number was used to determine the Mackin rating of each oyster following protocols adapted from reference 31. Copy numbers of the small-subunit (SSU) rRNA gene for H. costale and H. nelsoni were converted to a semiquantitative scale (0 for no infection, 1 for initial infection, 2 for moderate infection, and 3 for severe infection) adapted from reference 52.

DNA extractions were performed with the Qiagen Allprep PowerFecal DNA/RNA kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol with slight modifications, including the addition of 2 μl of proteinase K (Thermo Fisher Scientific, Waltham, MA) and 13.5 μl of 20% SDS (Thermo Fisher Scientific, Waltham, MA) to the lysis buffer. Mechanical lysis was performed via two rounds of homogenization with a TissueLyser (Qiagen, Hilden, Germany) for 2.5 min each at 30 Hz, and the samples were incubated for 5 min at room temperature in between. Samples included approximately 30 mg each of gut, gill, and mantle tissue as well as all the collected hemolymph, pallial fluid, and inner shell samples from each oyster (Data Set S1A). Additionally, one preparation of the ZymoBIOMICS Microbial Community Standard (Zymo, Irvine, CA) was extracted with the same protocol and used as a reference for quality control in the amplicon sequencing. The DNA from the 0.2- and 5-μm-pore-size filters from the water samples was extracted with the Qiagen QIAamp PowerFecal DNA kit following the manufacturer’s protocol with the same modifications as above. However, the TissueLyser was used only for a total of 1 min (30 s at 30 Hz, incubated for 5 min at room temperatures, and then 30 s at 30 Hz). A PVC pipe cutter was used to open the Sterivex filter, and a sterile scalpel and forceps were used to remove half of the filter sample for DNA extraction (53).

16S rRNA gene profiling was performed on the oyster tissue samples, the ZymoBIOMICS Microbial Community Standard sample, and the water samples collected at the time of sampling. Template DNA was amplified with V4 primers 515F/806R with Illumina adapter overhangs (Table 1) and 2× Phusion HF master mix (Thermo Fisher Scientific, Waltham, MA) (54, 55). The PCR was configured in a touchdown protocol and contained the following program: 3-min initial denaturation at 94°C, followed by 10 cycles of a 45-s denaturation at 94°C, 60-s annealing at 60°C with every cycle decreasing 1°C, and 30 s of elongation at 72°C, followed by 25 cycles of a 45-s denaturation at 94°C, 60-s annealing at 50°C, and 30 s of elongation at 72°C; 10 min of final extension at 72°C; and an indefinite hold at 4°C. Sequencing libraries were prepared at the Rhode Island Genomics and Sequencing Center (RIGSC) using a standardized protocol as follows: 50 ng of AMPure XP-cleaned PCR product was used as a template in a second PCR (5 cycles) with 2× Phusion HF master mix in order to attach full indices and adapters with the Illumina Nextera index kit (Illumina, San Diego, CA). The resultant PCR products were cleaned with AMPure XP, profiled with an Agilent Bioanalyzer DNA1000 chip (Agilent Technologies, Santa Clara, CA), quantified using a Qubit fluorometer (Invitrogen, Carlsbad, CA), and normalized for pooling. Then, the pooled library was quantified with qPCR with a Roche LightCycler 480 using the KAPA Biosystems Illumina kit (KAPA Biosystems, Wilmington, MA). The final pooled library underwent two rounds of 2 × 250-bp sequencing on the same Illumina MiSeq instrument (Illumina, San Diego, CA).

In order to enhance the metagenomic binning of coassembled contigs, a microbial cell enrichment procedure was performed on the gut samples through the adaptation of a procedure previously applied to marine sponges (56). The gut samples of 10 oysters were pooled and homogenized with 10 ml of ice-cold, sterile artificial seawater (28 ppt) in an autoclaved mortar and pestle and vortexed for 10 min. The homogenized sample was centrifuged at 100 × g at 4°C for 30 min. The supernatant was extracted and centrifuged again at 10,967 × g at 4°C for 30 min. After removing the supernatant, the resultant pellet was washed three times with ice-cold, sterile artificial seawater. The pellet from the final round of centrifugation was diluted with ice-cold, sterile artificial seawater and pushed through a 3-μm-pore-size (diameter of 47 mm) filter (Millipore, Burlington, MA) with a syringe. The flowthrough was subsequently pushed through a 0.2-μm-pore-size (diameter of 22 mm) filter (Sterlitech, Kent, WA). DNA extraction was then performed on the 0.2-μm filter using the Qiagen QIAamp PowerFecal DNA kit following manufacturer protocols. Mechanical cell lysis was performed via two rounds of homogenization with a TissueLyser (Qiagen, Hilden, Germany) for 1 min each at 30 Hz, and the sample was incubated for 5 min at room temperature in between.

Two rounds of metagenomic sequencing were performed (Data Set S1A). In the first round, 12 gut samples were individually barcoded and prepared for shotgun metagenomics using a read length of 2 × 150 bp across three lanes on an Illumina HiSeq 4000. In the second round, one microbiome-enriched oyster gut sample was prepared using the microbiome enrichment protocol described above and barcoded for shotgun metagenomics using a read length of 2 × 150 bp on an Illumina NextSeq 500 (mid-output). Preparation of metagenomic libraries was performed at the RIGSC using standardized protocols including DNA fragmentation with an S220 ultrasonicator (Covaris, Woburn, MA), quantification with a Qubit fluorometer (Invitrogen, Carlsbad, CA), and library preparation (including end repair, A-tailing, adapter ligation, and size selection) on the Wafergen Apollo 324 with PrepX reagents and consumables (TaKaRa Bio, Kusatsu, Japan), followed by quality control using an Agilent Bioanalyzer high-sensitivity chip (Agilent Technologies, Santa Clara, CA) and quantification of the final library with qPCR with a Roche LightCycler 480 using the KAPA Biosystems Illumina kit (KAPA Biosystems, Wilmington, MA). Sequencing of the metagenomic libraries was performed at Duke University’s Sequencing and Genomic Technologies Core.

Demultiplexed read pairs from the two MiSeq runs were analyzed independently using the Quantitative Insights Into Microbial Ecology 2 (QIIME2) software version 2018.4 (57). First, the raw reads were imported using the qiime tools import function, and the read quality was inspected for the selection of trimming parameters using the qiime demux summarize function. Using the DADA2 plugin, samples were subjected to quality filtering, merging, denoising, and chimera detection with the qiime dada2 denoise-paired function trimming the forward and reverse reads at 20 bp at the 5′ end, the forward reads at 230 bp at the 3′ end, and the reverse reads at 200 bp at the 3′ end (58). This step led to the establishment of a count table that maps the occurrence of ASVs in each sample. The ASV count tables from the two different MiSeq runs were merged using the qiime feature-table merge function in QIIME2 using the sum overlap method specified with the argument –p-overlap-method sum. All subsequent steps were performed on the combined MiSeq run data. Following quality control and pooling of the two sequencing runs, samples with less than 10,000 total read pairs were dropped from further analysis due to low sequencing depth (Data Set S1A).

Taxonomic assignment was performed with QIIME2’s sklearn classifier by mapping to the SILVA database release 132 (59). Mitochondrial and chloroplast sequences were removed from the ASV count table based on the SILVA taxonomic assignments. The resulting ASV count table and taxonomy data were exported and analyzed in R version 4.0.2. The taxonomic assignments were used to validate the identification of all 8 bacterial members (3 Gram negative and 5 Gram positive) in the ZymoBIOMICS Microbial Community Standard sample, suggesting adequate cell lysis during DNA extraction. The Shannon and Simpson indices were computed based on the ASV count table with the diversity function in the vegan package version 2.5-6 (60). Statistical significance in the ASV counts, Shannon index, and Simpson index comparisons was evaluated using the pairwise.wilcox.test function in the core stats library in R with a Holm P value adjustment. Unless specified, the comparisons of alpha diversity, beta diversity, and taxonomic relative abundance were performed on all the oyster samples that underwent amplicon sequencing.

To account for the compositional nature of the ASV count table (61, 62), a phylogenetic isometric log-ratio transformation was performed with the PhILR package version 1.16.0 (63). The ASV phylogeny used in the PhILR analysis was constructed following multiple sequence alignment with MAFFT (64), alignment masking with the qiime alignment mask function, phylogenetic reconstruction with FastTree (65), and midpoint rooting in QIIME2. The PhILR-transformed ASV counts were used to compute between-sample distances using the Euclidean distance measure. Pairwise PERMANOVAs were conducted on the between-sample distances with the R package pairwiseAdonis version 0.4 to assess the statistical significance in comparisons of microbiome composition (66). Differential abundance of specific taxa was evaluated among the different tissue types and between P. marinus-infected and uninfected samples using the ancombc function with default parameters in the analysis of compositions of microbiomes with bias correction (ANCOM-BC) package, version 1.0.2, and the heatmap visualizations were constructed using the geom_tile function in ggplot2 to demonstrate the log fold changes and statistical significance (32). For the comparative analysis to prior studies, raw reads from three published eastern oyster microbiome data sets were retrieved from the NCBI Sequence Read Archive (SRA) linked to the BioProject accession numbers PRJNA518081 (25), PRJNA504404 (40), and PRJNA386685 (39) and subjected to the analytical protocols as specified above.

The data from both rounds of metagenomic sequencing underwent trimming and adapter removal with Trimmomatic version 0.38 (leading and trailing bases below a quality score of 3 were removed, a 4-bp sliding window was used with an average quality score of 15, and a minimum read length of 130 bp was required), and the last 10 bp of all reads were removed with Cutadapt version 1.9.1 (67, 68). Due to differences in the chemistry associated with the Illumina NextSeq and HiSeq platforms, the reads from the two rounds of metagenomic sequencing underwent slightly different quality control measures. The reads from the sample that was sequenced on the Illumina NextSeq 500 underwent additional steps in Trimmomatic and Cutadapt in order to remove poly(G). In Trimmomatic, a sequence composed of 50 guanine residues was added to the adapter sequences, and a sequence of 10 guanine residues was used as an adapter sequence in Cutadapt (allowing for an error rate of 20%) to enable the removal of poly(G) sequences at the 3′ end of reads.

The quality-filtered reads were then mapped with BBMap version 37.36 to a combined reference database containing the C. virginica genome (RefSeq assembly accession GCF_002022765.2), the P. marinus genome (RefSeq assembly accession GCF_000006405.1), all complete bacterial genomes downloaded on 18 April 2018, and a collection of manually curated draft genomes of bacteria isolated from oysters (strain names and NCBI RefSeq accession numbers in Data Set S2F). Reads mapped to a bacterial genome or unmapped to any genomes in the reference database were used for the metagenomic coassembly and binning. Paired reads were considered mapped to a bacterial genome if at least one read from the pair was mapped. Metagenomic coassembly was performed using MEGAHIT version 1.1.1 with default parameters (69). Binning of the coassembled contigs was performed with default parameters in MetaBat2 version 2.12.1 (70) using read mappings performed with BBMap, with the resulting SAM files converted to BAM format and sorted with SAMtools release 1.5 (71). Completeness and contamination of the MAGs reconstructed from the oyster microbiome were assessed using CheckM version 1.0.5. The ssu_finder function in CheckM was used to identify 16S rRNA genes in the metagenomic coassembly (72). Candidate 16S rRNA genes greater than 1,000 bp in length were retained.

All ASVs taxonomically assigned to the class Mollicutes from this and a prior study (39), two full-length Mollicutes 16S rRNA genes assembled from the metagenomes collected in this study, and 11 full-length 16S rRNA genes obtained from the Firmicutes outgroup (Data Set S2E) were placed into a reference tree with SILVA’s ACT server, using the SINA aligner (73). The minimum identity was set at 85%, and the number of neighbors per query was set at 5. The resultant phylogeny was collapsed based on the environment where the reference SILVA sequences were derived in iTOL (74).

Phylogenomic reconstructions were performed on the Mollicutes MAG and the Chlamydiae MAG with reference genomes from the two corresponding taxa (Data Set S2B and C). To ensure consistency in gene calling, all genomes were first analyzed using Prodigal version 2.6.3 (75) with genetic code 4 for the Mollicutes and standard genetic code for the Chlamydiae. Conserved single-copy genes (CSCGs) were identified through the analysis of bidirectional best BLAST hits as described in a previous study (76). The identified CSCGs were individually aligned with MUSCLE version 3.8.31, the alignments were concatenated, and a phylogeny was reconstructed with RAxML version 8.2.10 using the JTT substitution model and the GAMMA model of rate heterogeneity. Pairwise average nucleotide identity (ANI) values were computed with OrthANI version 1.40 (77). The average amino acid identity (AAI) was computed by the mean protein identity values of all bidirectional best BLAST hits identified based on the following thresholds: E value less than or equal to 1e−3, sequence identity greater than or equal to 30%, and coverage of 70% or higher to both sequences in the alignment.

A metabolic reconstruction of the Mollicutes MAG was developed based on an initial mapping of the proteome to the KEGG database through the KAAS server (78), and the annotation of transporters was based on homology searches to the Transporter Classification Database (79). The draft reconstruction underwent manual curations through comparisons to three Mycoplasma isolates closely related to the Mollicutes MAG (i.e., M. mobile, M. todarodis, and M. marinum) and using reference annotations from an existing metabolic model of Mycoplasma pneumoniae (80). The metabolic reconstruction was incorporated into a genome-scale model using PSAMM (81). Metabolic pathway gaps were identified using the PSAMM gapcheck function, which in turn were used to guide the manual curation of gene functions from the MAG and other contigs identified from the metagenomic coassembly. For example, the pyruvate dehydrogenase complex was initially not identified in the Mollicutes MAG, but its subunits were found on other contigs in the coassembly with top BLAST hits (in the NCBI nonredundant protein database) to members of the Mycoplasma genus. In this case, the pyruvate dehydrogenase reaction was included in the metabolic model to highlight the potential presence of this function. Some gap reactions, such as the transport of acetate, were included in the metabolic network to represent the potential export of acetate as a metabolic product.

Raw sequencing data are available from the NCBI Sequence Read Archive under the BioProject accession PRJNA658576.