Effects of Ocean Acidification, Hypoxia, and Warming on the Gut Microbiota of the Thick Shell Mussel Mytilus coruscus Through 16S rRNA Gene Sequencing

Gut microbiota play a very important role in the health of the host, such as protecting from pathogens and maintaining homeostasis. However, environmental stressors, such as ocean acidification, hypoxia, and warming can affect microbial communities by causing alteration in their structure and relative abundance and by destroying their network. The study aimed to evaluate the combined effects of low pH, low dissolved oxygen (DO) levels, and warming on gut microbiota of the mussel Mytilus coruscus. Mussels were exposed to two pH levels (8.1, 7.7), two DO levels (6, 2 mg L−1), and two temperature levels (20, 30°C) for a total of eight treatments for 30 days. The experiment results showed that ocean acidification, hypoxia, and warming affected the community structure, species richness, and diversity of gut microbiota. The most abundant phyla noted were Proteobacteria, Bacteroidetes, and Firmicutes. Principal coordinate analysis (PCoA) revealed that ocean acidification, hypoxia, and warming change microbial community structure. Low pH, low DO, and increased temperature can cause shifting of microbial communities toward pathogen dominated microbial communities. Linear discriminant analysis effect size (LEfSe) showed that the significantly enriched biomarkers in each group are significantly different at the genus level. Phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) analysis revealed that the gut microbiome of the mussels is associated with many important functions, such as amino acid transport and metabolism, transcription, energy production and conservation, cell wall, membrane and envelope biogenesis, and other functions. This study highlights the complexity of interaction among pH, DO, and temperature in marine organisms and their effects on the gut microbiota and health of marine mussels.


INTRODUCTION
The phenomenon in which pH reduction occurs in seawater due to atmospheric CO 2 diffusion into the surface waters of oceans is called ocean acidification (Orr et al., 2005). It is estimated that due to anthropogenic activities ocean pH will decrease ∼0.3-0.5 units by the end of this century (Pachauri et al., 2014). It has been reported that ocean acidification has adverse effects on the marine calcified organisms with consequences for species diversity, trophic interaction, and other ecological processes (Pörtner et al., 2014). There is evidence that during the last century, average temperature has increased ∼0.7 • C, and it will further increase 1.8-4 • C by the end of the 21st century (IPCC, 2007). The temperature fluctuations can affect growth, metabolism, and development of marine bivalve mollusks, and in the case of thermal stress conditions, it can also cause increased mortality and decreased growth in the bivalve mussels (Dickinson et al., 2012). Dissolved oxygen (DO) concentration of <2 mg L −1 is considered a hypoxic condition in the aquatic ecosystems (Wu, 2002). The East China Sea has one of the largest coastal hypoxic areas in the world which covers an area of more than 12,000 km 2 (Chen C. et al., 2007). Hypoxia can cause reduced ingestion of food and also reduced growth (Welker et al., 2013). Hypoxia can also adversely affect the immune system of bivalves (Monari et al., 2007), leading to an increased risk of pathogen infections in bivalves (Wang et al., 2012(Wang et al., , 2014. The digestive physiology and energetic budget of bivalves can be affected by dissolved oxygen concentrations . The previous studies found that the digestive function and antioxidant response of Mytilus coruscus were affected by hypoxia, acidification, and warming of the ocean (Khan et al., 2020(Khan et al., , 2021. It is worth noting that in recent years, more and more scholars have paid attention to the potential combined effects of ocean hypoxia, acidification, and warming on marine life (Sampaio et al., 2021).
It has been proved that microbiota associated with host play a very important and essential role in animal health. Microbiota in the host provide many functions to the host, such as protection from diseases and nutrient processing (Sweet and Bulling, 2017). The gut microbiome is considered an organ contributing to the regulation of host metabolism (Rastelli et al., 2019). The balanced and healthy gut microbiota of animals play an important role to prevent infections and support the host immune system (Kau et al., 2011;Man et al., 2017). The studies have found that gut microbiota have functions in physiology, reproduction, development, energy balance, behavior, immunity, and life history of the host (Semova et al., 2012). Gut microbiota also play an important role in maintaining the homeostasis of the host, and changes in its composition can bring metabolic shifts that may result in host phenotype alterations (Turnbaugh et al., 2006;Visconti et al., 2019). High microbial abundance and diversity have been found in marine invertebrate animals (Olson and Kellogg, 2010). The alterations in the host-microbial communities due to stressful conditions are associated with an increased risk to diseases and poor health conditions because opportunistic pathogens and non-resident microbial species find their way to increase populations within the host (Vezzulli et al., 2013;Lokmer and Wegner, 2015). Mass mortality incidents of pacific oyster (Crassostrea gigas) have been reported in farming areas of Europe (Samain and McCombie, 2008), this loss was attributed to the complex interactions among environmental variables, microbial pathogens, and oysters (Pernet et al., 2014). The marine bivalve mollusks have high microbiome diversity and abundance (Olson and Kellogg, 2010;Pierce and Ward, 2018). The microbiota of invertebrates can be altered partially by fluctuations in the external environment (Tanaka et al., 2004;Beleneva and Zhukova, 2009). The stable host-microbiota associations in the digestive tract of a variety of animals like crustaceans (Freese and Schink, 2011), insects (Robinson et al., 2010), and oysters (Wegner et al., 2013) have been reported.
Various abiotic environmental factors, such as pH, DO, and temperature can affect the mussel microbiota communities (Rubiolo et al., 2019). The structure of the gut microbiome is highly changeable and can be altered as a result of environmental factors (Goodrich et al., 2016;Xie et al., 2016;Bestion et al., 2017). The changes in these environmental factors can directly affect mussel microbiome communities by altering their structure and diversity which will finally affect the ability of the mussel to convert food into energy and other physiological functions (Rubiolo et al., 2019). Low pH value has been reported to increase the growth of Vibrionaceae, and consequently, increased growth of vibrio genus might have adverse effects on mussel growth (Meron et al., 2011). Decreased pH affects bacterial communication, and this process is also called quorum sensing. It is a key process that affects the microbial interactions in the ecosystem, depending on autoinducing peptides (AIPs) in Gram-positive bacteria and acyl homoserine lactones (AHLs) in Gram-negative bacteria (Miller and Bassler, 2001;Manefield and Whiteley, 2007). Both have important ecological and biological functions and regulate symbiosis, competence, virulence, secondary metabolite production, extracellular enzymes, biofilm formation, and bioluminescence in the marine environment (De Kievit, 2009;Chong et al., 2012;Mangwani et al., 2012). Changes in environmental factors such as temperature, would affect microbial community composition either directly or indirectly (Mathai et al., 2020). The increased temperature is a critical factor in regulating the bacteria virulence (Lee et al., 2001;Kimes et al., 2012) and innate immunity of bivalves (Cheng et al., 2004). In marine biotic and abiotic habitats, temperature is an important factor to shape microbiota communities (Fuhrman et al., 2008). Due to global warming, alteration of mutualistic microbial communities to pathogen-dominated communities and increase in infectious diseases in the hosts have already been reported (Ritchie, 2006;Altizer et al., 2013). Few studies have been carried out to reveal the bacterial diversity of mollusks digestive gland (Rubiolo et al., 2018). Previous studies reported that high water temperature altered M. coruscus gut microbiota by favoring the accumulation of opportunistic pathogens in adult mussels (Li et al., 2018). The homeostasis and health controlling mechanisms of gut microbiota in aquatic animals under hypoxic conditions are still unclear (Fan et al., 2020). Especially the combined effects of acidification, warming, and hypoxia on the gut microbiota of the aquatic animals are still not well-known .
Mytilus coruscus is an ecologically and economically important bivalve species in the coastal waters of Yellow Sea and East China Sea (Dong et al., 2017). The habitats of the economically important bivalve mollusks face the conditions of hypoxia (Diaz and Rosenberg, 2008), thermal stress (Baumann and Doherty, 2013), and ocean acidification (Wallace et al., 2014), which are dangerous for the survival of the bivalves species in recent years (Beck et al., 2011). As these multiple stressors practically exist in the habitat of M. coruscus in the summer season, their effects on the mussels may be synergistic, antagonistic, or additive, it is urgent to understand the complex effects of these multiple stressors (Stevens and Gobler, 2018), and the interactive effects of these stressors need to be clarified Sui et al., 2017). Therefore, in this experiment, we elevated the combined effects of ocean acidification, hypoxia, and warming on the gut microbiota of the thick shell mussel M. coruscus, which can help us to understand the physiological responses of the mussels to these stressful conditions and to make plans for better health and protection of marine animals.

Mussel Collection
Mussels were collected from the Gouqi Island located in the Zhejiang province of China. Healthy mussels (shell length: 8.43 ± 0.6 cm; wet weight: 79.58 ± 4.64 g) were selected and cleaned by removing all kinds of attachments adhered to the mussel shell surfaces. After that, mussels were placed in the aquatic tanks (28L), equipped with filtering and circulating water devices. The standard natural seawater conditions of the collection site of mussels were established in all experimental tanks, i.e., pH 8.1 ± 0.1, DO 6 ± 0.1 mg L −1 , temperature 20 ± 0.2 • C. Before starting the experiment, the mussels were acclimatized for 1 week under the above conditions. Mussels were fed with microalgae Chlorella spp. (2.5 × 10 4 cells ml −1 ) twice a day.

Experimental Setup
In the current experiment, we tested the effects of three major environmental factors (pH, DO, and temperature) on the mussel gut microbiota. Two pH levels (8.1 and 7.7), two DO levels (6 and 2 mg L −1 ), and two temperature levels (20 and 30 • C) were set for the experiment. The pH level 8.1 represents the ongoing mean pH level at the sampling site while pH level 7.7 represents the forecast average pH at the end of 2,100 (Orr et al., 2005;IPCC, 2007), and extreme of present natural variability at the Shengsi island of East China Sea (Shang et al., 2019). At the sampling site of mussels, DO concentrations of <2-3 mg L −1 were reported by Chen C. et al. (2007) and Chen J. H. et al. (2007). DO level ≤ 2 mg L −1 is considered hypoxic level for aquatic organisms (Sui et al., 2017), while 6 mg L −1 is considered as normal DO level. Thus, in our experiment, DO 2 mg L −1 was considered as hypoxic treatment for mussels. In the present experiment, 20 • C denotes the mean temperature experienced by the mussels, while 30 • C represents the peak temperature level experienced by the mussels at the site of collection (Wu et al., 2016). All factors were used in a full factorial design. A total of 8 treatments were set in the experiment T1 (DO 6 mg L −1 × pH 8.1 × 20 • C), T2 (DO 6 mg L −1 × pH 8.1 × 30 • C), T3 (DO 6 mg L −1 × pH 7.7 × 20 • C), T4 (DO 6 mg L −1 ×pH 7.7 × 30 • C), T5 (DO 2 mg L −1 × pH 7.7 × 20 • C), T6 (DO 2 mg L −1 × pH 7.7 × 30 • C), T7 (DO 2 mg L −1 × pH 8.1 × 20 • C), and T8 (DO 2 mg L −1 × pH 8.1 × 30 • C). The experiment was conducted in triplicates with a density of 20 mussels per tank (28 L). A recirculatory system with filters was used to prevent metabolic wastes accumulation in tanks by the mussels. Before starting the experiment, the mussels were acclimatized gradually to the experimental conditions. DO (from 6 to 2 mg L −1 ) and pH (from 8.1 to 7.7) were gradually reduced within 2 days, and the temperature was increased step by step from 20 to 30 • C in 5 days. The entire experiment lasted for 30 days.

Experimental System
For maintenance of required temperature and DO levels in the experimental tanks, we used electronic thermostats and DO regulators (Loligo systems Apps, Tjele, Denmark), respectively. Low DO level was maintained by injecting N 2 directly into tanks using an oxygen regulator linked to a computer. Required experimental low pH level was maintained by introducing pure CO 2 directly using WTW pH 3310 connected to pCO 2 /pH feedback STAT systems (DAQ-M) and SenTix 41 pH electrodes (Loligo Systems Inc.), operated by Cap CTRL software (Loligo Systems Inc.). A Salinometer (S/Mill-E, Atago) was used to measure salinity. Salinity was kept constant at 25 psu. For pH measurement, a portable pH meter (pH-201, MSITECH (Asia-Pacific) Pte. Ltd., Singapore) gauged with the NBS scale was used. The titration method was adopted to determine total alkalinity (A T ) values for all treatments. To calculate pCO 2 , calcite saturation state (Ωca) and aragonite saturation state (Ωar) from A T , pH NBS , temperature, and salinity, CO 2 sys was used with constants, K1 and K2 (Mehrbach et al., 1973) modified by Millero (2010). The seawater chemical parameters of all treatment groups in the experiment are summarized in Table 1.

Sampling of Tissues and Preparation
Prior to sampling, mussels were stopped feeding for 12 h. Three mussels per replicate were randomly collected for dissection at the sampling day (day 30). Guts were carefully removed and then cleaned with phosphate buffer (50 mM, pH 7.4), dried using tissue paper, placed in tubes, and kept on ice. For further analysis, the tissues of three mussels from each replicate were pooled representing one sample of the relative treatment and immediately preserved at −80 • C. Thus, for each treatment, three samples were collected.

DNA Extraction and Amplification
A DNA extraction kit was used for total genomic DNA extraction by following the protocols of the manufacturer. With NanoDrop and agarose gel, the DNA quality and quantity were determined. Dilution of the extracted DNA was done to the concentration of 1 ng µL −1 and for further processing preserved at −20 • C. The diluted DNA was used as a template for PCR amplification of bacterial 16S ribosomal ribonucleic acid (rRNA) genes with the barcoded primers and Takara Ex Taq (Takara). For bacterial diversity analysis, V3-V4 variable regions of 16S rRNA genes were amplified with universal primers 343F (5 ′ -TACGGRAGGCAGCAG-3 ′ ) and 798R (5 ′ -AGGGTATCTAATCCT-3 ′ ) (Nossa et al., 2010).

Library Construction
By using gel electrophoresis, amplicon quality was visualized, purified with AMPure XP beads (Agencourt), and amplified for another round of PCR. After purification with the AMPure XP beads again, a Qubit dsDNA assay kit was used for the final amplicon quantification. For subsequent sequencing, the same amounts of purified amplicon were pooled.

Data Analysis
A FASTQ file format was used for raw sequencing data. Using Trimmomatic software, paired-end reads were pre-treated (Bolger et al., 2014) to recognize and break off unclear bases (N). Using the sliding window trimming approach, low-standard sequencing with average quality score <20 was also removed.
With the help of Flash software (Reyon et al., 2012) pairedend reads were grouped. The parameters of assembly were: 10 base pairs of minimal overlapping, 200 base pairs of maximum overlapping, and 20% of maximum non-similarity rate. Further denoising of sequences was performed by abandoning reads with ambiguous or homologous sequences or base pairs below 200 bp reads and reads with 75% of bases above Q20 were retained. Reads with chimera were then detected and removed. QIIME software (version 1.8.0) was used to achieve these two steps (Caporaso et al., 2010).
To produce operational taxonomic units (OTUs), clean reads were subjected to the removal and clustering of primer sequences, for this purpose Vsearch software with a 97% similarity cutoff was used (Edgar et al., 2011). With the help of the QIIME package, the representative read of each OTU was selected. Using an RDP classifier (confidence threshold of 70%), all representative reads were annotated and blasted against the Silva database version 123 . All representative reads were annotated and blasted against the Unite database using blast (Lobo, 2008).
The richness index of Chao1 and the diversity indices of Shannon and Simpson were calculated based on OTUs, and the differences of microbial diversity indices were determined by the Kruskal-Wallis test. A principal coordinate analysis (PCoA) was performed to reflect the differences in microbial construction between samples. A linear discriminant analysis of effect size (LEfSe) was performed to assess the difference in the relative abundance of bacteria, and the logarithmic transformation linear discriminant analysis (LDA) score threshold was 3.5. PICRUSt function prediction of the mussel gut microbes is based on the cluster of orthologous groups of proteins (COG) database.

Unique Biomarkers Detected in the Mussel Guts
The LEfSe was conducted to determine taxonomic differences in relative abundance of bacteria and unique biomarkers in all experimental treatments. Mucispirillum is the bacterial genus level biomarker in T1, distinguished this group from all the other groups. Acinetobacter is genus level biomarker in T3 that is unique to this group. Endozoicomonas, Kistimonas genera distinguished T4 from all other groups. Pseudomonas is the unique genus of T5. Polaribacter and Candidatus − Rhabdochlamydia are the unique genera of T6. Pseudorhodobacter and Mycoplasma are the special genera of T7. Similarly, Atopostipes and Flammeovirga are the special genera of T8 (Figure 3).

DISCUSSION
In order to better forecast and understand the physiological responses and other ecological mechanisms of marine organisms in a continuously changing ocean environment due to the combined stress of ocean acidification, hypoxia, and warming, we need multi-stressor exposure experiments. The environmental stressors can act independently or exert synergistic, antagonistic effects on organism performances (Todgham and Stillman, 2013). Very less is known about the variations in microbial compositions in bivalve mollusks (Lokmer and Wegner, 2015). This study is the first report to provide the evidence of the combined effects of ocean acidification, hypoxia, and warming on the gut microbiota of mussels. This study highlights the importance of multi-stressor exposure experiments to reveal the combined effects of three global stressors that practically exist in the oceans at the same time. The intestinal microbiota populations arrange an ecological network via interspecific interactions and rely on the ecological network to regulate its stability and dynamic balance (Coyte et al., 2015). Strong relationships and more sub-modules are the key network properties for effectiveness and strength of the network (Barabasi and Oltvai, 2004;Kitano, 2004). Bacterial communication is mediated by quorum sensing, which involves the production, release, and detection of autoinducers (small hormone-like molecules). Using these chemical communications, bacteria can detect their environment and respond by changing in population size and/or community structure, which in some instances may be harmful for the host organisms (Waters and Bassler, 2005). The resident microbes by competing with pathogen colonization may contribute to the host protection and therefore play a role in the host immune system (Kamada et al., 2012;Chu and Mazmanian, 2013;Desriac et al., 2014). In the current study, the results also explain such phenomena. Evidence from Goodrich et al. (2016) and Xie et al. (2016) explained that gut microbiota are strongly pliable and changeable under varying environmental conditions, whereas the core gut microbiota remain essentially stable (Nicholson et al., 2012). Bivalve microbiota are highly diverse in nature and are easily affected by environmental conditions, such as water temperature, pH, salinity, DO, nutrients, and infections (Green and Barnes, 2010;Lokmer and Wegner, 2015). In the current study, different pH, DO, and temperature levels affected the relative abundance and community composition of gut microbiota of the M. coruscus among the treatments. The stability of the microbiota communities can be detected in organisms as a result of any ecological or environmental disturbance, which is very important for doing community functions . In the present experiment results, we did not find any perpetuity of microbial communities among experimental treatments as a result of different pH, DO, and temperature levels. The differences in the relative abundance and community composition of microbiome populations among treatments indicate the unstable state of microbiome populations in the results; it also indicates that environmental conditions can influence the relative abundance and community composition of gut microbiota. The high-level interaction disturbances among the host, environment, and microbiota can destabilize the balance of the community composition and lead to poor states which have drastic consequences for the host. In order to achieve the microbial stability, the host has the capacity for balancing the community composition Brown et al., 2013). In the results, the unstable community structure may indicate that the host attempts to stabilize the microbiome community structure and composition and arrange a stable ecological network. However, the current research cannot determine whether the change of bacterial community is due to the physiological changes of mussels caused by the environment or the direct response of bacteria to the environment.
The results showed that dominant microbiota in the gut of the M. coruscus are Proteobacteria, Becteroidetes, and Firmicutes at the phylum level. In addition, the results are similar to Yang et al. (2020), who investigated the effects of ocean acidification and microplastics on the gut microbiota of this species. Similarly, FIGURE 2 | Representing microbial diversity indices of Chao 1 (A), Shannon (B), and Simpson (C) and principal coordinate analysis of gut microbiome (D). [T1 (DO 6 mg L −1 × pH 8.1 × 20 • C), T2 (DO 6 mg L −1 × pH 8.1 × 30 • C), T3 (DO 6 mg L −1 × pH 7.7 × 20 • C), T4 (DO 6 mg L −1 × pH 7.7 × 30 • C), T5 (DO 2 mg L −1 × pH 7.7 × 20 • C), T6 (DO 2 mg L −1 × pH 7.7 × 30 • C), T7 (DO 2 mg L −1 × pH 8.1 × 20 • C), T8 (DO 2 mg L −1 × pH 8.1 × 30 • C)]. Song et al. (2018) reported that Bacteroidetes and Firmicutes are the two dominant bacteria phyla that colonize the digestive tract of aquatic animals. Firmicutes flourish within mussels because they have capacity to produce spores that show resistance to dehydration and also severe environmental conditions (Sathe et al., 2017). Unstable microbiota community structure in the intestine and microecological disorders can be reflected by intestinal Proteobacteria. Proteobacteria symbiotic bacteria are present in the intestine of healthy animals. Under certain environmental conditions, in the intestine these Proteobacteria may become intestinal microorganisms that are responsible for causing inflammation (Shin et al., 2015). Previous studies reported that microbial communities in Rapana venosa and Achatina fulica are dominated by Tenericutes, Bacteroidetes, and Proteobacteria (Pawar et al., 2014;Yang et al., 2019).
The differences in the dominant microflora in the digestive tract of different organisms in these studies and the current study may be due to differences in their diets and different environmental conditions. Yang et al. (2020) reported that pH variations have no significant effects on the gut microbiota of the M. coruscus. But Zeng et al. (2020) reported that the gut microbiota of moth larva have been significantly affected by dietary pH. In the current experiment results, low pH has no significant effects in most treatments groups. Some authors indicated that high temperature increases heterogeneity of the microbiota communities among individuals (Erwin et al., 2012;Lokmer and Wegner, 2015). Li et al. (2019a) reported that increased temperature affected the microbial diversity in the hemolymph of mussels. The microbiota with lower diversity in unhealthy oysters were found and the FIGURE 4 | Representing gut microbiota presumptive functions and their relative abundances of the mussel Mytilus coruscus in different experimental treatments. [T1 (DO 6 mg L −1 × pH 8.1 × 20 • C), T2 (DO 6 mg L −1 × pH 8.1 × 30 • C), T3 (DO 6 mg L −1 × pH 7.7 × 20 • C), T4 (DO 6 mg L −1 × pH 7.7 × 30 • C), T5 (DO 2 mg L −1 × pH 7.7 × 20 • C), T6 (DO 2 mg L −1 × pH 7.7 × 30 • C), T7 (DO 2 mg L −1 × pH 8.1 × 20 • C), and T8 (DO 2 mg L −1 × pH 8.1 × 30 • C)]. diversity was largely affected by heat stress (Lokmer and Wegner, 2015). Hypoxic conditions can disturb the homeostasis of the gut microbiota by quenching important communities and destroying their network structure (Fan et al., 2020). Results of these experiments coincide with the current experimental findings. In current results, we also found that increased temperature, low pH, and low DO have affected the microbiome diversity, relative abundance and community structure, and composition of gut microbiota populations. Probably mussels need to adjust their physiological activities when facing unfavorable environmental conditions and thus the gut microbiota changes.
According to Li et al. (2019a), PCoA analysis revealed that temperature is an important factor affecting the microbial community, 21 • C treatments groups were separated from 27 • C treatments groups, suggesting that disturbance of microbial community may be a sign of an unhealthy status. Similar results were also reported by Lokmer and Wegner (2015). In the current experiment, the PCoA analysis revealed that all treatments having different pH, DO, and temperature influenced microbial populations and shaped them by separating various groups from each other.
Microbial diversity significantly increased in the gut microbiota of mussels (Mytilus galloprovincialis) by elevated temperature as revealed by Chao1 and Shannon and Simpson indices (Li et al., 2019b). We also observed significantly increased gut microbiota as a result of elevated temperature in some treatments groups as indicated by Chao1 and Shannon and Simpson indices. Li et al. (2018) reported that significant reduction in microbial diversity was noted in live mussels (M. coruscus) gut microbiota as a result of heat stress. In our results, we also found decreased microbial diversity in some groups due to high temperature as revealed by Chao 1 and Shannon and Simpson indices. Similarly, some groups showed decreased microbial diversity due to hypoxia and decreased pH in terms of Chao 1 and Shannon and Simpson indices. High microbial species diversity in healthy animals enhances their ability to tolerate unfavorable environmental conditions (Lozupone et al., 2012). Low microbial diversity highly favors unhealthy status in many organisms (Chang et al., 2008;Green and Barnes, 2010).
The LEfSe analysis in the current study showed that the most unique taxa of all experimental treatments have pathogenic potential. This analysis indicates that acidification, hypoxia, and warming have caused establishment and dominancy of pathogenic taxa in the gut of mussels, as stated by Lokmer and Wegner (2015), in C.gigas environmental stressors such as high temperature favored the shift of bacterial communities toward pathogen-dominated microbiota communities and also enhanced the secondary opportunistic pathogens colonization. Moreover, we found no unique phylum, family, class, order, and genus levels taxa in T2 group compared with other experimental groups. It may be due to that this treatment has favored the environmental conditions for species taxon or sub taxa of phylum, family, class, and order. The deviations from the normal microbiome may indicate pathogenic infections or changes in metabolic processes. This can help us to identify and predict the effects of global change on mussel health (Shade and Handelsman, 2012).
Microbial community function prediction has revealed that intestinal microbiota are associated with carbohydrates, purines, pyrimidines, and amino acid metabolism in the host (Human Microbiome Project Consortium, 2012). The results from PICRUSt analysis predicted that microbiota of the thick shell mussel M. coruscus is mainly associated with amino acid transport and metabolism, transcription, energy production and conversion, cell wall, membrane and envelop biogenesis, translation, ribosomal structure and biogenesis, and other functions. Therefore, further multi-stressor exposure studies are necessary to reveal the potential mechanisms involved in these pathways.

CONCLUSION
The present study is the first report on the combined effects of three major global environmental problems, i.e., ocean acidification, hypoxia, and warming on the gut microbiome of the mussel M. coruscus. The current study demonstrates that ocean acidification, hypoxia, and warming have effects on the relative abundance, community composition, species richness, and diversity of the gut microbiota communities of mussels. Our results revealed that low pH, low dissolved oxygen (DO), and high temperature favor the proliferation of opportunistic and pathogenic bacteria, and predicted many possible functions of gut microbiota which need further investigations to explore their mechanisms. Further studies are needed for the development of strategic plans to save our marine organisms from the adverse effects of these ongoing global environmental problems.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: NCBI [accession: PRJNA744459].

AUTHOR CONTRIBUTIONS
YW, MH, and XZ contributed to the conception and design of the study. FU, YS, XC, and HK performed the experiments. FU and YS performed the statistical analysis. FU, AZ, JF, and WL wrote the first draft of the manuscript. JP, XZ, MH, and YW wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.