Adult M. coruscus were acquired from Gouqi Island, China (30°72’N, 122°77’E) and transported to the laboratory on the same day. The mussels were cleaned and placed in a 10 L polycarbonate tank filled with seawater filtered through a 1.2 μm pore size acetate-fiber membrane (FSW). The water temperature was maintained at 21°C to reflect the collection site’s average winter temperature. The mussels received a daily mixed diet comprising Isochrysis zhanjiangensis and Platymonas helgolandica.

Mussel spawning was induced using a previously described method (Yang et al., 2008) with modifications. In brief, after overnight storage on ice, mussels were placed in 10 L polycarbonate tanks at 21°C until spawning began, when they were transferred immediately to 2 L beakers containing FSW. Egg and sperm FSW suspensions were collected and mixed gently to achieve fertilization, followed by incubation for 20 min. Excess sperm was washed off with FSW using a 20 μm mesh size nylon plankton net, and the fertilized eggs were incubated in a 2 L beaker with FSW at 18°C. Two days later, the swimming larvae were fed daily with 5×10⁴ cells/mL of I. zhanjiangensis and P. helgolandica and cultured at 21°C.

The Institutional Animal Care and Use Committee (IACUC) of Shanghai Ocean University approved the animal handling procedures and experimentation.

Pediveliger larvae (35 days post fertilization, larva at this stage can metamorphose into a juvenile) were randomly and evenly distributed in three randomized triplicate groups. The groups were subjected to three temperature treatments terminating at the same time: the larvae were placed directly in seawater at 25°C for 4 hours group Mc25.S); 4 days at 25°C (group Mc25.L), by raising the water one degree a day; and 21°C as control (group Mc21). Four degrees of warming represents the differential temperature rise predicted for 2100. The global mean temperature will have increased by 3.2–4.6 °C (Rogelj et al., 2016), and the increase in seawater temperature over a short time during marine heat waves (Field et al., 2014).

During the trial, the larvae were fed microalgae as usual. At the end of the experiment, 50 µg of larva were combined into a sample (every treatment group has three replicates), frozen, and stored at −80 °C for microbiome analysis.

According to the manufacturer’s instructions, a FastDNATM Spin Kit for Soil (MP Biomedicals, Santa Anna, CA, USA) was used to extract genomic DNA from the larval samples. The purity and concentration of the DNA were measured on a Nanodrop One instrument (Thermo Fisher Scientific, Waltham, MA, USA). Sterile, filtered (0.2 μm) water was added to standardize the DNA concentration to 5 ng/μl for all samples.

The universal bacterial primers, 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’) with a 12 bp barcode (synthesized by Invitrogen, Carlsbad, CA, USA) were used to amplify the bacterial 16S rRNA gene V3-V4 region. The 50 µl PCR reactions comprised 25 μl of 2× Premix Taq (Takara Biotechnology, Dalian, China), one μl of each primer (10 M), three μl of DNA (20 ng/μl) template, and 20 μl of distilled water. The reactions were run on a BioRadS1000 instrument (Bio-Rad, Hercules, CA, USA) with the following thermocycling conditions: initialization for 5 min at 94°C; 30 cycles of 30 s of denaturation at 94°C, 30 s of annealing at 52°C, and 30 s extension at 72°C; followed by a 10 min final elongation at 72°C.

Gel electrophoresis in 1% agarose was used to detect the PCR amplicons. GeneTools (version 4.03.05.0, SynGene, New Delhi, India) was used to assess the density of the PCR bands. The amplicons showing bright leading bands were mixed at equidensity ratios. An EZNA Gel Extraction Kit (Omega Bio-tek, Norcross, GA, USA) was used to purify mixed PCR products. A NEBNext^® Ultra™ DNA Library Prep Kit for Illumina^® (New England Biolabs, Ipswich, MA, USA) was used to generate the sequencing libraries with index codes according to the manufacturer’s protocol. A Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2100 system (Agilent Technologies, Waldbron, Germany) were used to assess the quality of the libraries. The Illumina Hiseq 2500 platform (Illumina Inc., San Diego, CA, USA) was used to sequence the libraries, generating 250 bp paired-end reads. The sequencing was carried out by Guangdong Magigene Biotechnology Co. Ltd (Guangzhou, China).

The Trimmomatic V0.33 (Bolger et al., 2014) was used to filter the raw paired-end reads to generate high-quality, clean reads. FLASH V1.2.11 (Magoč and Salzberg, 2011) was used to merge the clean reads dependent on the overlap between the paired-end reads. Mothur V1.35.1 (Kozich et al., 2013) used the unique barcodes and primers to assign the sequences to the samples, then, the barcodes and primers were removed to obtain effective clean tags. USEARCH V10 (Edgar, 2010) was used for sequence analysis. Sequences sharing greater than or equal to 97% similarity were assigned to the same operational taxonomic unit (OTU). Each OTU was represented by the most frequently occurring sequence, which was screened and annotated. We used the SILVA (Quast et al., 2013) database to annotate the taxonomic information for each representative sequence (using the default confidence threshold ≥ 0.5). After removing contaminating OTUs, we determined the total number of OTU sequences and the OTU types, respectively. A normalized (subsampled) OTU table was obtained based on the sample with the fewest sequences. Using the OTU table, we calculated the annotation ratio for each classification level to determine the sequence composition of each sample at each classification level. The OTU table was also used to construct a Venn graph, histogram, and heat map displayed using R software (v2.15.3).

The complexity of species diversity in the samples was represented by the Chao1 richness index and Simpson’s and Shannon diversity indices, calculated using QIIME software (v1.9.1) (Bolyen et al., 2019) and displayed using the R software (v2.15.3). The vegan package of R software was used to carry out non-metric multidimensional scaling (NMDS) based on the weighted and unweighted distance matrix. Linear discriminant analysis Effect Size (LEfSe) was used to find the biomarker of each group. The obtained raw sequences were submitted to the NCBI database and received the accession number PRJNA1029461.

GraphPad Prism 9.0 (GraphPad Software, Boston, MA, USA) and one-way ANOVA tests were used to analyze the larval microbiota’s relative abundance and taxonomical composition at the different levels and Alpha Diversity.

The Animal Ethics Committee at Shanghai Ocean University approved all the protocols for mussel acclimation and experimentation in the present study (registration number SHOU-DW-2020-032).

The 16S rDNA gene sequencing yielded 1,233,013 sequences clustered into 6,830 distinct OTUs ( Supplementary Table S1 ). The larvae hosted 29 unique phyla, 72 classes, and 198 orders of bacteria. The most abundant phylum in the Mc21 group was Proteobacteria (73.1%), followed by Actinobacteria (21.2%), Bacteroidetes (4.5%), and many other unclassified bacteria (1.2%). The most abundant phyla in the Mc25.S group were Proteobacteria (71.0%), Actinobacteria (23.0%), Bacteroidetes (4.6%), and many other unclassified bacteria (1.4%). The most abundant phyla in the Mc25.L group were Proteobacteria (35.7%), Actinobacteria (3.8%), and Bacteroidetes (59.2%), and many other unclassified bacteria (1.2%) ( Supplementary Table S2 , Supplementary Figure S1 ).

We identified 18 different families with an abundance > 1%, and 235 families with an abundance < 1% were classified as “others” ( Figure 1A ). Most reads (16.60–31.53% of the OTUs) were obtained from the Mc21 and Mc25.S groups and were assigned to Burkholderiaceae, Nocardiaceae, and Xanthomonadaceae. In contrast, most reads (5.00–54.55% of total OTUs) obtained from the Mc25.L group were assigned to Burkholderiaceae, Flavobacteriaceae, and Rhodobacteraceae ( Figure 1A ; Supplementary Table S3 ). The analysis of the relative abundance showed that at the family level, the microbial community structure between the Mc21 and Mc25.S groups was similar but different from the Mc25.L group. The relative abundance of Burkholderiaceae, Nitrincolaceae, Nocardiaceae, and Xanthomonadaceae decreased significantly, and that of Flavobacteriaceae and Rhodobacteraceae increased significantly in the Mc25.L group compared with those in the Mc21 and Mc25.S (P < 0.05) ( Figures 1B, C ; Supplementary Table S3 ). Notably, the abundance of Pseudoalteromonadaceae decreased significantly irrespective of whether the larvae were in warmer seawater for 4 hours or 4 days. (P < 0.05) ( Figures 1B, C ; Supplementary Table S3 ).

At the genus level, 21 genera were identified with an abundance > 1%, and 443 genera with an abundance < 1% were classified as “others” ( Figure 2A ). In the Mc21 and Mc25.S groups, the dominant genera were Alteromonas, Brevundimonas, Delftia, Microbacterium, Neptuniibacter, Neptunomonas, Pseudoalteromonas, Rhodococcus, Stenotrophomonas, Tenacibaculum, and Thalassotalea (average reads per group 1.04–31.40% of the total OTUs) ( Figure 2A ; Supplementary Table S4 ). Delftia, Lewinella, Neptuniibacter, Pseudoalteromonas, Rhodococcus, Ruegeria, Stenotrophomonas, Tenacibaculum, Thalassotalea, and Vibrio were the dominant genera in the Mc25.L group (average reads per group 1.09–52.65% of the total OTUs) ( Figure 2A ; Supplementary Table S4 ). By comparing the abundance in the different experimental groups, the results showed the similarity of the microbiota between the Mc21 and Mc25.S groups, except for Pseudoalteromonas (P < 0.05) ( Figure 2C ; Supplementary Table S4 ). In the Mc25.L group, the larval microbiota had more differences. The abundance of Delftia, Neptunomonas, Pseudoalteromonadaceae, Rhodococcus, and Stenotrophomonas decreased significantly, and Tenacibaculum increased considerably compared with that in the Mc21 and Mc25 groups (P < 0.05) ( Figures 2B, C ; Supplementary Table S4 ). These results indicated that rapid warming in a short time does not impact the M. coruscus larval microbial community structure. Conversely, a gradual increase in temperature over four days can change the population levels of certain bacterial species and possibly alter the microbial community composition.

Based on the Venn diagrams, the core microbial communities, i.e., occurring in all treatment groups, included 145 families and 216 genera, representing 57.31% of the total families and 47.68% of the entire genera ( Figure 3 ). Eighteen families and 21 genera with an abundance > 1% belonged to these core families and genera bacterial communities ( Supplementary Tables S3 , S4 ). In contrast, unique bacterial communities had an abundance of < 1% in each treatment group ( Supplementary Tables S3 , S4 ).

The microbiota richness and evenness were investigated between the control and the temperature treatment groups. There were no significant differences among the treatment groups in the Shannon and Simpson diversity indices (P > 0.05, Figures 4A, B ). However, the Chao1 index of the Mc25.L group was significantly higher than that of the Mc25.S and Mc21 (P < 0.05, Figure 4C ), indicating that the gradual increase in temperature caused a significant increment in the microbiota species abundance.

The NMDS analysis was employed as a proxy for temporal stability. It showed that the samples were homogenous in the Mc25.L group, and the larvae microbiota in this group changed compared with the other two groups, which showed some overlap ( Figure 5 ).

Linear discriminant analysis effect size (LEfSe) analysis was performed to reveal the significant ranking of abundant modules. The linear discriminant analysis (LDA) score and the cladogram showed differences between Mc25.S, Mc25.L, and Mc21 ( Figure 6 ). At the phylum or order level, the biomarkers demonstrating significant differences between the three groups were Caulobacterales in Mc25.S; PB19, Cellvibrionales (genus Halioglobus), Rhodobacterales (genus Litoreibacter and Kiloniella), Caulobacterales (family Hyphomonadaceae), Oceanospirillales (family Alcanivoracaceae), Micavibrionales (family Micavibrionaceae) in Mc25.L; and Alteromonadales (family Pseudoalteromonadaceae and genus Colwellia), Firmicutes (phylum Firmicutes and order Lactobacillales) in Mc21 ( Figure 6 ; Supplementary Table S5 ).