Sample collection was performed in April 2021 in the area near Dadeng Island (the natural habitat of the green mussel P. viridis), Fujian Province, China (Figure 1). Thirty mussel individuals (10 from each of the three sampling sites, 5–7 cm in shell length) were randomly collected and gently washed with 0.1% phosphate-buffered saline (0.1% phosphate buffered saline, PBS). The mussels were immediately stored in coolers (4°C) and transferred to the laboratory within 2 hours. In the laboratory, the mussels were dissected under sterile conditions. The adductor muscle of each green mussel individual was incised with a sterile scalpel. The hemolymph (H) was then taken from the posterior adductor muscle using a sterile 1 ml syringe, transferred to sterile tubes, and immediately frozen in liquid nitrogen following the method of Musella et al. (2020). The byssus (B), mantle (M), gill (G), foot (F), digestive gland (D), and gonad (testis, T; ovary, O) were dissected from each green mussel individual using a sterile scalpel, rinsed with sterile PBS, transferred to sterile tubes, and immediately frozen in liquid nitrogen. All the samples were stored at −80°C.

The samples for the outer shell surface (U) and substratum (S) of P. viridis were obtained by swab sampling. Swab samples were taken by washing the outer shell surface and the substratum (nylon rope) with 0.1% PBS and then gently rubbing the surfaces using cotton swabs. The obtained swab samples were stored in sterile tubes at 4°C and then transferred to a −80°C refrigerator. A total of 1.5 L of seawater (W) was also collected from each of the three sampling sites (Figure 1). Seawater samples were stored in coolers (4°C) and transferred to the laboratory along with the mussel and swab samples. In the laboratory, the seawater samples were filtered with MF-Millipore membranes (0.22 μm, Biosharp, Anhui, China). The membranes were then immediately frozen in liquid nitrogen and stored at −80°C.

Total microbial DNA was extracted from the tissue and environmental samples using a HiPure Stool DNA Kit (Magen, Guangzhou, China) according to the manufacturer's protocol. Specifically, 200 μl of hemolymph, 20–30 mg of byssus, mantle, gill, foot, digestive gland, testis, or ovary, one filtered seawater sample, and two to three cotton swabs of the outer shell surface or substratum were used for total microbial DNA extraction. DNA concentrations were measured using a Qubit^® dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA).

The library construction and next-generation sequencing (NGS) were conducted by GENEWIZ (Suzhou, China). The sequencing library was constructed using a MetaVX Library Preparation Kit (GENEWIZ, Suzhou, China). Briefly, 20–30 ng of DNA was used to generate amplicons that covered the V3–V4 hypervariable regions of the bacterial 16s rRNA gene. The forward primer sequence was 5′-CCTACGGRRBGCASCAGKVRVGAAT-3′, and the reverse primer sequence was 5′-GGACTACNVGGGTWTCTAATCC-3′. The concentration was detected using a microplate reader (Tecan, Infinite 200 Pro), and the fragment sizes (expected value ~600 bp) were measured through 1.5% agarose gel electrophoresis.

The PCR amplification procedure was as follows. A 25 μl of PCR mixture was prepared containing 2.5 μl of TransStart buffer, 2 μl of dNTPs (2.5 mM each), 1 μl of forward primer, 1 μl of reverse primer, 0.5 μl of TransStart Taq DNA polymerase, and 17.5 μl of gDNA. The first round PCR was performed by the following program: 3 min of denaturation at 94°C followed by 14–16 cycles of denaturation at 94°C for 10 s, annealing at 57°C for 90 s, elongation at 72°C for 15 s, and a final extension at 72°C for 5 min. In addition, a splice with Index primers was added to the end of the PCR product of 16S rDNA by PCR to facilitate NGS sequencing. The 50 μl of PCR mixture was prepared containing 2.5 μl of TransStart buffer, 2 μl of dNTPs (2.5 mM each), 3 μl of INDEX primer N, 3 μl of Index primer S, 4 μl of 1x cocktail, 0.5 μl of TransStart Taq DNA polymerase, and 10 μl of ddH₂O. The second round of the PCR procedure was performed using the following program: 3 min of denaturation at 94°C followed by 10–12 cycles of denaturation at 94°C for 10 s, annealing at 60°C for 30 s, elongation at 72°C for 15 s, and a final extension at 72°C for 5 min at 4°C.

The libraries were quantified to a final concentration of 10 nM. Sequencing was performed on a MiSeq platform (Illumina, San Diego, CA, USA) using a 300 paired-end protocol according to the instrument manual. The MiSeq Control Software (MCS) was used to obtain raw sequencing data.

Raw reads were processed with the software Qiime (1.9.1) (Caporaso et al., 2010) to splice forward and reverse reads according to their overlap. Data with low quality were then removed from the spliced sequences using Cutadapt (1.9.1) (Martin, 2011), and chimeras were removed with Vsearch (1.9.6) (Rognes et al., 2016).

After quality filtering, an average of 84,710 reads per sample were clustered into operational taxonomical units (OTUs) using a 97% similarity threshold criterion. After normalization of the entire dataset, all 68,306 remaining OTU sequences were searched against the Silva 138 16SrRNA database (http://www.arb-silva.de/). Finally, the species taxonomic analysis of the representative OTU sequences was completed by using the Bayes algorithm of the Ribosomal Database Program classifier, and the community composition of each sample was examined at several taxonomic levels.

Alpha diversity indices (the Chao1 and Shannon indexes) were estimated through sequence random sampling according to the effective sequence number and calculated using Qiime (1.9.1). Rank-abundance was constructed in R (v 3.3.1) based on the OTU analysis results. Beta diversity was estimated by unweighted UniFrac analysis and principal coordinates analysis (PCoA). Unweighted UniFrac distances were calculated using Qiime (1.9.1), and the significance of data separation in the PCoA was based on the distances in the Bray-Curtis matrix. Finally, the functional profile analysis of the microbial community in each sample was performed using the software PICRUSt (v1.0.0) against the Kyoto Encyclopedia of Genes and Genomes (KEGG) categories.

Values were expressed as the mean ± standard deviation. The differences in alpha diversity among samples were statistically analyzed through Kruskal–Wallis tests, and a p < 0.05 was defined as indicating statistical significance.

Microbial communities were sequenced from 33 samples [27 tissue samples, including three each from hemolymph, byssus, mantle, gill, foot, digestive gland, ovary, testis, and outer shell surface samples, and six environmental samples (three substratum and three seawater samples)]. The NGS was conducted on the V3–V4 hypervariable region of the 16S rRNA gene to obtain the microbial composition and structure of each sample. The sequencing runs generated a total of 2,795,441 paired-end reads (46,806–132,794 reads per sample; Supplementary Table S2). A total of 2,254,088 clean reads passed quality filtering (40,762–100,570 per sample; Supplementary Table S2) and were clustered into 1,951,958 OTUs. According to the heatmap of the 30 most abundant OTUs (Figure 2), the samples could be divided into four clusters. The first cluster comprised samples from seawater; the second was composed of foot, outer shell surface, and substratum; the third comprised byssus, digestive gland, gill, mantle, ovary, and testis; and the fourth was the hemolymph, with specific OTUs that were rare in other samples.

Amplicon sequencing revealed the overall microbial composition of P. viridis tissues (Figure 3). At the phylum level, Proteobacteria (mean relative abundance (r.a.) ± SD, 54.40 ± 20.94%), Bacteroidetes (23.94 ± 14.94%), and Spirochaetae (6.66 ± 19.20%) were dominant in mussel tissues, followed by Firmicutes that accounted for 4.30% (Figure 3A). In addition, phylogenetic trees of the top 30 genera showed that most of the microorganisms belonged to Proteobacteria and Bacteroidetes (Supplementary Figure S2). Overall, the bacterial community composition at the phylum level was similar among different mussel tissues, except for hemolymph that was dominated by Spirochaetae. At the class level, Gammaproteobacteria (43.72 ± 26.90%), Bacteroidia (23.79 ± 14.88%), and Alphaproteobacteria (10.69 ± 12.10%) were dominant in mussel tissues (Figure 3B). At the family level, the most representative families were Comamonadaceae (16.54 ± 16.64%), Flavobacteriaceae (9.75 ± 16.51%), Enterobacteriaceae (8.09 ± 9.65%), and Spirochaetaceae (6.66 ± 19.20%) (Figure 3C and Supplementary Table S3). At the genus level, the dominant genera were Sphaerotilus (9.78 ± 13.81%), the unclassified genus of Spirochaetaceae (6.66 ± 19.20%), and Tolumonas (4.71% ± 8.29%) (Figure 3D).

Regarding the microbial community composition of environmental samples, Proteobacteria (58.15 ± 7.88%) and Bacteroidetes (33.53 ± 7.78%) were the dominant phyla (Figure 3A). At the class level, Alphaproteobacteria (34.54 ± 8.10%), Bacteroidia (33.39 ± 7.77%), and Gammaproteobacteria (23.60 ± 6.03%) were dominant in environmental samples (Figure 3B). The most representative families were Flavobacteraceae (24.27 ± 9.65%) and Rhodobacteraceae (19.08 ± 5.63%) (Figure 3C). The dominant genera were the unclassified genera of Rhodobacteraceae (8.01 ± 5.44%) and Tenacibaculum (6.94 ± 11.35%) (Figure 3D).

The bacterial community composition at the family level differed among the mussel tissues (Figure 3C and Supplementary Table S3). Spirochaetaceae, Vibrionaceae, and Mycoplasmataceae were the dominant bacteria in hemolymph, while the representative families in the mantle, digestive gland, gill, and gonad were Comamonadaceae, Enterobacteriaceae, Aeromonadaceae, and Weeksellaceae. Comamonadaceae and Methylophilaceae were dominant in mussel byssus. It was interesting to note that the foot and outer shell surface had similar dominant bacteria (Flavobacteriaceae and Rhodobacteraceae) with the substratum. These results implied that the microbiota of the tissues exposed or partly exposed to the environment (the outer shell surface and foot) was distinct from that of the tissues inside the shells. In particular, the foot (the organ responsible for exploring the substratum), the outer shell surface (the part exposed to the surrounding environment), the substratum, and seawater all shared the same dominant taxa (i.e., Flavobacteriaceae and Rhodobacteraceae), indicating the close connection of microbiota of foot and shells with that of the surrounding environment.

Alpha diversity analysis results showed that species richness in mussel tissues was significantly different from the diversity in the environmental samples (Kruskal–Wallis test, p < 0.05). According to the Chao1 index which evaluates the total number of species in a sample, the total number of microbial species was much higher in the substratum samples than in other samples, with the species number in the seawater sample being the lowest (Figure 4A). Among mussel tissues, the outer shell surface had the highest microbial diversity. In contrast, the microbial diversity of the byssus was the lowest. The Shannon diversity index (Figure 4B) showed that the substratum had the highest microbial diversity, while hemolymph had the lowest. The same trend was seen in the rank-abundance curve, where the abundance and evenness of microbial species in the substratum were the highest, and the hemolymph was the lowest (Figure 4C).

For the beta diversity analysis, a Bray-Curtis PCoA was performed to explore the microbial composition in different mussel tissues and environmental samples. The mussel tissues and environmental samples could be divided into four categories (Figures 4D–F). The microbial communities in the samples of the foot, outer shell surface, and substratum were similar. The seawater sample comprised a separate category. Byssus, gill, digestive gland, mantle, ovary, and testis were classified as one category, while hemolymph clustered apart from all these tissues owing to its unique microbial composition. Interestingly, the microbial diversity in the foot and outer shell surface was closer to that of the substratum rather than to seawater, suggesting that the mussel foot and shell have microbial connections with the environment, especially the substratum.

The functional profiles of the microbial communities for tissues and environmental samples were predicted using PICRUSt based on phylogenetic information. A total of 15,039 KOs (KEGG Orthology) were obtained through PICRUSt function prediction analysis, and 296 pathways were annotated. Pathways with high abundance included transporters, general function prediction only, and ABC transporters pathways, while pathways with low abundance included the biosynthesis of type II polyketide products and the biosynthesis of 12-, 14-, and 16-membered macrolides. The top 100 pathways ranked by abundance were utilized to cluster all 33 samples (Figure 5). The samples could be divided into two groups. In group I, seawater was clustered with the gill, byssus, digestive gland, mantle, and ovary. In group II, the foot, outer shell surface, substratum, testis, and hemolymph were clustered, with the hemolymph isolated from the other samples. As shown in Figure 5, most pathways exhibited a higher relative abundance in group I than in group II. For specific samples, the hemolymph microbiome was characterized by the pathways involved in osmotic regulation (taurine and hypotaurine metabolism), polycyclic aromatic hydrocarbon degradation, purine metabolism, and other environmental information processing. The outer shell surface and foot samples were clustered with substratum, and they were enriched in the degradation of aromatic compounds such as limonene and pinene degradation, RNA degradation, and glycosyltransferases. The byssus and gill microbiomes were characterized by the pathways involved in benzoate degradation, geraniol degradation, and fatty acid metabolism. The results showed a specific profile of seawater characterized by the enrichment in pathways involved in photosynthesis protein and nitrogen cycle (including glycine, serine, and glutamate metabolism and arginine and proline metabolism).