Microbiome Dysbiosis in Mytilus chilensis Is Induced by Hypoxia, Leading to Molecular and Functional Consequences

1. Introduction

Dissolved oxygen (DO) levels are critical for the survival and health of marine organisms [1]. In coastal marine ecosystems, reference thresholds have been established to consider the biological and ecological effects of dissolved oxygen levels in the ocean [2,3,4,5,6,7,8]. Normoxia is defined as the condition in which DO levels range between 9.0 and 3.0 mg L^-1, whereas hypoxia occurs when DO levels fall between 2.0 and 0.1 mg L^-1 [2,3,4,5,6,7,8]. Hypoxia on the south-central coast of Chile is a pressing issue severely affecting the region’s bivalve mollusks, with potentially severe consequences for the marine ecosystem [9,10,11,12,13]. Long-term research reveals that hypoxic zones occupy substantial portions of the water column during upwelling seasons, with an intensifying trend linked to climate change and altered oceanographic processes [14]. These trends portend a grim outlook for marine organisms reliant on stable oxygen availability [14].

Hypoxia-induced stress can severely impact bivalves, ranging from physiological to molecular levels [15,16]. Physiologically, this stress can increase clearance and respiration rates and reduce food intake, potentially resulting in stunted growth in these organisms[15,17]. Additionally, hypoxia stress can negatively affect the immune system of bivalves, increasing their susceptibility to pathogen infections[18,19]. At the cellular level, hypoxia can trigger autophagy, increase oxidative stress, and reduce cell viability [20]. Molecularly, prolonged hypoxia can affect protein metabolism, inflammation-related genes, and programmed cell death [20,21]. In extreme cases, hypoxia stress can lead to mass mortality and stranding of bivalves [10].

The microbiota, or microbiome, refers to the community of microorganisms inhabiting a specific environment with distinct physical and chemical characteristics [22,23,24]. The microbiota associated with the host plays a crucial role in animal health by providing vital functions such as disease protection and nutrient processing [15,22]. Moreover, microbiota composition influences the host’s physiology, stress tolerance, and fitness [25]. The microbiome is considered an organ that regulates host metabolism and is essential for maintaining a healthy balance in the host immune system due to its relationship with specific diseases [26,27,28,29]. Recent studies have shown that a diverse and balanced microbiota can indicate better metabolic health [30]. Greater microbiota diversity is associated with improved lipid profiles, lower levels of pro-inflammatory cytokines, and higher levels of anti-inflammatory cytokines [30]. Furthermore, correlations between microbiota diversity, enzyme activity and genetic pathways related to metabolism and health have also been observed [30].

Intestinal microbiota significantly influences the host’s physiology, reproduction, development, energy balance, behavior, and life history [31]. The intestinal microbiota of bivalve mollusks plays an essential role in their health and nutrition [32]. The diversity and abundance of microorganisms in their digestive tract assist in food digestion, strengthen the immune system and may influence their growth and development [32,33,34]. However, stress-induced alterations in microbial communities, such as those caused by hypoxia, may increase disease risk and compromise bivalve health [35,36]. For example, the proliferation of opportunistic pathogens, such as those from the Vibrio and Arcobacter genera, could significantly increase the host’s susceptibility to diseases, contributing to increased mortality [37]. Furthermore, fluctuations in the external environment, such as abiotic factors, can alter the structure, species richness, and diversity of intestinal microbiota [15].

Mytilus chilensis is a bivalve species of ecological and economic importance in the coastal waters of the Los Lagos Region in Chile [38,39]. Climate change, manifested in declining oxygen levels in the water, induces systematic changes in bivalve mollusks and their bacterial symbionts [15,40,41]. Economically valuable bivalve species, including M. chilensis, are increasingly exposed to hypoxic conditions, threatening their viability and sustainability [42,43]. The aquaculture industry in southern Chile, heavily reliant on seed collection from the Reloncaví Fjord and grow-out operations around Chiloé Island, faces recurrent hypoxia episodes exacerbated by seasonal upwelling and anthropogenic eutrophication [44,45].

The Reloncaví system, comprising the Reloncaví Fjord and Reloncaví Sound, is particularly vulnerable to hypoxia due to the influx of suspended allochthonous organic matter from rivers, especially during late winter and early spring [45,46,47]. This period is marked by glacial meltwater contributions, which dominate over precipitation-driven runoff [45,46,47]. Major riverine inputs, including the Puelo, Petrohué, and Cochamó Rivers, deliver substantial organic material, fueling microbial decomposition and oxygen consumption [48,49]. Recent risk assessments identify the Reloncaví estuarine system as a hotspot for high inorganic nutrient concentrations, intense phytoplankton blooms, and elevated chlorophyll levels, particularly in late winter (August–September) [50,51,52]. These conditions, driven by eutrophication from intensive salmon aquaculture, promote the formation of low dissolved oxygen water (LDOW) zones, where stratification and particulate organic matter deposition exacerbate oxygen depletion [52]. The prolonged water residence times in the Reloncaví system, coupled with high biological oxygen demand, including phytoplankton respiration and bacterial remineralization of organic material, underscore the urgent need to investigate hypoxia’s impact on M. chilensis, particularly at the microbiome level [45,52,53,54,55].

Therefore, our objective is to understand the influence of hypoxia on the intestinal and gill microbiota of the native mussel M. chilensis. Specifically, a comparative evaluation of the bacterial communities in the intestine and gills of M. chilensis exposed to hypoxia was conducted using 16S rRNA sequencing with nanopore technology. Our study is the first to investigate the effects of hypoxia on M. chilensis from a hologenome concept. This knowledge could enhance our understanding of host-specific microbiomes and their role in supporting host ecology. Additionally, it can help elucidate the physiological responses of M. chilensis to hypoxia and infer potential health and disease changes that may arise from future stress factors.

2. Materials and Methods

2.1. Experimental Design (Mussel Acclimatization, Hypoxia Challenge, and Sampling for Microbiological Analysis)

Blue mussels (M. chilensis) utilized in this study were sourced from the experimental laboratory at the Marine Biological Station of the Universidad de Concepción, Chile. From an initial pool of 480 individuals, 36 mussels were selected and distributed into three experimental replicates (n=12 mussels per replicate). The experimental design incorporated a structured sampling regime: three mussels were sampled under hypoxic conditions at day 10 (n=9 total mussels across replicates), three mussels were sampled following reoxygenation (normoxic conditions) at day 20 (n=9), three mussels were sampled following reoxygenation at day 40 (n=9), and three mussels were sampled under hypoxic conditions at day 50 (n=9). In total, 18 mussels were sampled under hypoxia, and 18 were sampled under reoxygenation.

The experiment spanned 50 days and consisted of alternating hypoxic (dissolved oxygen concentrations of 2.0 mg/L) and normoxic (dissolved oxygen concentrations of 7.2 ± 0.2 mg/L) phases, designed to simulate natural tidal fluctuations in oxygen levels experienced by mussels in their cultivation environment. The experiment commenced with a 10-day hypoxic exposure, followed by a 10-day reoxygenation period. This cycle was repeated, resulting in two 10-day hypoxic exposures and two 10-day reoxygenation periods, culminating in a 50-day experiment. The duration of each hypoxic exposure was determined based on methodologies established in prior studies [56].

Prior to the hypoxia challenge, mussels were acclimatized for 38 days in filtered seawater (12.5 ± 0.94 °C) under continuous flow, aeration, and feeding. Following acclimatization, dissolved oxygen concentrations within the recirculation system were monitored daily and adjusted as necessary using nitrogen gas injection to maintain the target hypoxic level of 2.0 mg/L.

Numerous studies have explored physiological adaptations to environmental stressors, such as cyclic oxygen fluctuations experienced by intertidal bivalves during hypoxic and reoxygenation cycles [20,21,57,58,59,60,61,62,63,64]. This study employed a comprehensive approach to investigate the systemic effects of hypoxia and reoxygenation on M. chilensis, with a specific focus on microbiological shifts within gill and digestive gland tissues. These tissues were selected due to their critical roles in bivalve physiology. Gills, as the primary interface between the organism and its environment, are central to respiration and filter feeding [65,66]. Given their constant exposure to ambient conditions, gills are particularly susceptible to physiological stress induced by hypoxia. Consequently, microbial changes within gill tissues were investigated as potential biomarkers of hypoxia-induced stress [67].

In addition to gills, the digestive gland was analyzed due to its multifunctional role in nutrient assimilation, metabolic regulation, and immune response [68,69,70,71]. Under hypoxic conditions, bivalves often exhibit valve closure, significantly reducing filtration and respiration rates [72,73]. Since hypoxia impacts filtration-dependent nutrient processing and metabolic activity, shifts in the digestive gland microbiota were examined to elucidate the broader physiological consequences of hypoxia-induced dysbiosis. By examining microbial dynamics in both gill and digestive gland tissues, this study aimed to identify shared perturbation patterns and assess the physiological implications of hypoxia exposure in M. chilensis.

To reduce inter-individual variability in the microbiota associated with gill and digestive gland tissues, samples from three mussels were pooled to create a single biological replicate, with each sequencing sample representing nine pooled individuals [74]. This pooling strategy, utilizing nine individuals per sequencing sample, was implemented to enhance the detection of consistent microbial patterns while optimizing sequencing resources. Gill and digestive gland samples were collected at four time points: short-term hypoxia (day 10), reoxygenation (day 20, normoxia), reoxygenation (day 40, normoxia), and long-term hypoxia (day 50). Specifically, Figures 1A,B, 3, 4 and 7 present microbiota composition from pooled samples at control time points (days 20 and 40, n=9 mussels each) and hypoxic phases (days 10 and 50, n=9 mussels each). Conversely, Figures 1C, 2, 5, and 6 display microbiota composition from combined normoxic samples at days 20 and 40 (n=18 total mussels) and combined hypoxic samples at days 10 and 50 (n=18 total mussels).

The samples were preserved in molecular-grade ethanol, transported at 4 °C, and stored at -80 °C until further processing for microbiological analysis.

2.2. DNA Isolation and 16S Amplification

Total bacterial DNA was isolated from homogenized gill and digestive gland tissues of mussels, using the phenol-chloroform extraction method.

Primarily, approximately 20–30 mg of gill and digestive gland tissues were sectioned and homogenized using ceramic beads. The homogenates were then mixed with lysis buffer and incubated at 37°C for 2 hours. Each sample was processed separately, with three biological replicates per condition. Briefly, the samples were thawed at room temperature, washed, minced, and vortexed with 1 mL of lysis buffer containing 10 mM Tris-HCl, 400 mM NaCl, 100 mM EDTA, 0.4% SDS, and 100 μg/mL Proteinase K (pH 8.0). The mixture was incubated at 37°C under constant agitation until complete lysis was achieved, followed by mechanical disruption using ceramic beads to ensure thorough tissue breakdown.

After incubation, 1 volume of phenol-chloroform was added to each sample, followed by centrifugation at 12,000 rpm for 5 minutes at room temperature. The aqueous phase was carefully collected, and an equal volume of chloroform was added, followed by an additional centrifugation at 12,000 rpm for 5 minutes. The resulting aqueous phase was mixed with molecular grade absolute ethanol and transferred to a DNeasy Blood & Tissue column (Qiagen, MD, USA) to continue purification following the manufacturer’s protocol.

The quality and purity of the extracted DNA were assessed using a Nanodrop One spectrophotometer (Thermo Scientific, MA, USA), and its integrity was verified through electrophoresis in a 1% agarose gel prepared in TAE buffer (Tris-Acetic Acid-EDTA). DNA concentration was further quantified by fluorescence using a Qubit 4 fluorometer (Thermo Scientific, MA, USA) with the dsDNA BR Assay Kit (Thermo Scientific, MA, USA).

For 16S rRNA gene amplification, the isolated DNA was diluted to a concentration of 50 ng/μL and used as a template in a 25 μL PCR reaction containing LongAmp Taq DNA polymerase (New England Biolabs, MA, USA) and universal 16S primers: 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 1492R (5′-GGTTACCTTGTTACGACTT-3′) [75]. The thermal cycling conditions included an initial denaturation step at 95°C for 1 minute, followed by 25 cycles of 95°C for 20 seconds, 56°C for 30 seconds, and 65°C for 2 minutes, with a final extension at 65°C for 5 minutes. The resulting 16S rRNA PCR amplicons were confirmed by electrophoresis in a 1.2% agarose gel prepared in TAE buffer.

2.3. Library Preparation and Nanopore Sequencing

Nanopore sequencing is an advanced technique for characterizing microbial communities by sequencing the 16S rRNA gene amplicon [76]. Following the PCR amplification of the 16S rRNA gene, the resulting amplicons were pooled according to the experimental groups and purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA) to remove primer dimers and nonspecific amplification products. The purified amplicons were then quantified using a Qubit 4 fluorometer (Thermo Scientific, Waltham, MA, USA) to ensure the appropriate library concentration for sequencing.

Library preparation was conducted using the 16S Barcoding Kit (SQK-16S024, Oxford Nanopore Technologies, Oxford, UK) following the manufacturer’s protocol. The amplicons were barcoded through a PCR reaction using LongAmp Taq polymerase (New England Biolabs, MA, USA) and purified according to the instructions provided by the supplier.

The quality and size distribution of the prepared libraries were evaluated using the 2200 TapeStation system (Agilent, CA, USA) with DNA ScreenTape (Agilent, CA, USA). This ensured that the amplicons were within the expected size range and minimized the presence of adapter dimers or residual primers. The final library concentration was determined using the Qubit 4 fluorometer with the High Sensitivity D5000 ScreenTape (Agilent, Santa Clara, CA, USA). To ensure accuracy and mitigate potential batch effects, a mock microbial community (ZymoBiomics Microbial Community Standard, Zymo Research, CA, USA) was included in the analysis as a quality control standard.

According to the manufacturer’s guidelines, libraries were pooled at equimolar concentrations for multiplexing and loaded onto a Spot-ON Flow Cell for sequencing using the MinION platform (Oxford Nanopore Technologies, UK). The sequencing efficiency and run quality were monitored in real-time using the MinKNOW software (Oxford Nanopore Technologies), enabling comprehensive and high-resolution microbial community profiling.

2.4. Data Processing and Taxonomic Assignment

A rigorous data processing pipeline was applied to the nanopore sequencing reads to ensure robust and reliable results. First, base-calling was performed using the Guppy software (version 6.3.2), which ensured high-quality nucleotide base identification. Subsequently, a strict quality filter (Q-score ≥ 7) was applied to remove low-quality reads, thereby ensuring the resulting dataset’s integrity.

The resulting FASTQ files were processed using Porechop to remove adapter sequences and to demultiplex the reads, assigning them to their respective samples [77]. The demultiplexed reads were then used as input for taxonomic classification through the Emu algorithm, specifically designed to annotate full-length 16S rRNA sequences generated from nanopore sequencing [78]. Emu utilizes an expectation-maximization approach, which enables precise and reliable taxonomic assignments, minimizing false positives and negatives [78].

To further refine taxonomic assignment, a customized 16S rRNA database was constructed by combining reference databases with sequences derived from previous studies. A minimum abundance threshold of 0.01 was established to exclude low-representative taxa or potential artifacts. Classified reads were grouped into operational taxonomic units (OTUs) at a 97% similarity threshold, which enabled the operational definition of the different bacterial species present in the samples.

2.5. Community Profiling and Statistical Testing

The resulting OTU table was analyzed using the Microbiome Analyst software to obtain a comprehensive overview of the microbial community structure and diversity. Initially, singleton OTUs—those present in only one sample—were removed to reduce noise and enhance the robustness of the analyses. Subsequently, a logarithmic transformation was applied to the abundance data to normalize the distribution and improve the interpretability of the results.

Principal Coordinates Analysis (PCoA) was performed based on a Bray-Curtis distance matrix to assess differences in the microbial community composition across sample groups. This analysis facilitated the visualization of sample relationships and the identification of clusters with similar microbial compositions. In addition, an Analysis of Similarities (ANOSIM) was conducted to determine whether statistically significant differences existed in community structure between the compared groups.

Furthermore, a rarefaction curve was constructed using the Vegan package in R to evaluate sampling coverage and microbial community richness [79]. This curve assessed whether the number of sequences obtained was sufficient to capture the total community diversity and if there were differences in richness between the sample groups.

2.6. Data Processing and Heat-Tree Visualization of Microbial Communities

To perform the heat-tree analysis, the R programming language (version 4.3.3) and the integrated development environment RStudio were employed. Several R libraries were utilized, including Metacoder for hierarchical taxonomic data analysis and visualization, Dplyr for efficient data manipulation, and Vegan for ecological diversity analysis [79,80]. These tools were selected for their robustness and widespread use in microbiome data analysis.

The microbiome dataset used in this study was derived from 16S rRNA sequencing data processed through the MicrobiomeAnalyst pipeline [81]. The dataset comprised two primary files: one containing taxonomic read abundance per sample and another with metadata detailing sample attributes, such as experimental treatments and tissue sources. Both files were imported into R, and an initial exploratory analysis was conducted to assess data integrity and structure.

To ensure data reliability and mitigate the influence of sequencing artifacts, a rigorous filtering process was implemented. Taxa with low read counts, specifically those with fewer than five reads, were removed to minimize the impact of sequencing errors. Furthermore, taxa exhibiting zero abundance across all samples were excluded. A prevalence threshold of 20% was applied, ensuring that only taxa present in at least 20% of the samples were retained. This filtering process was conducted separately for different tissue types, including gills and digestive glands, to account for potential tissue-specific variations in microbial communities.

The heat-tree visualization was generated using the Metacoder package, which enables the hierarchical representation of microbial communities [82]. In this visualization, node size corresponds to taxon abundance, while node color indicates statistical differences across experimental conditions. Prior to visualization, the dataset was transformed into a format compatible with the Metacoder package. Taxonomic hierarchies were structured according to lineage, and a taxmap object was created to organize and analyze taxonomic relationships.

2.7. Linear Discriminant Analysis Effect Size (LEfSe) and Correlation Network Analysis

The linear discriminant analysis effect size (LEfSe) method was employed to identify significant differences in bacterial species abundance between gill samples from challenged and controlled individuals. A significance threshold of FDR-adjusted p-value < 0.05 and a logarithmic, linear discriminant analysis (LDA) score ≥ 4.0 were set as cut-off values to identify differentially abundant taxa. The top 15 discriminative features were visualized using a dot plot highlighting the primary bacterial taxa driving the differences between groups.

To further explore microbial interactions, a network correlation analysis was performed to evaluate the co-occurrence patterns of microbial taxa in the samples. Networks were constructed using the Sparse Correlations for Compositional Data (SparCC) algorithm, which is particularly suited for microbiome data due to its ability to handle compositional structures. The correlation network was estimated using 100 bootstrap permutations, with a significance threshold of P < 0.05 and a minimum correlation coefficient of 0.3. These parameters were chosen to ensure the identified associations’ robustness and minimize the inclusion of spurious correlations. This approach helps identify potential interactions and dependencies among bacterial species within the microbiota.

2.8. Prediction of Metagenomic Functional Potential

PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) was employed to infer the functional potential of the microbiome in the gill and digestive gland tissues of M. chilensis under hypoxic conditions based on marker gene sequences [83]. The R package ggpicrust2 was used to facilitate the functional analysis and visualization [84]. Given that PICRUSt2 is optimized for short-read sequencing data, a preprocessing step was necessary to adapt our long-read amplicon sequences generated through nanopore technology. The software HyperEx extracted the V3-V4 regions from the full-length 16S rRNA sequences, generating a FASTA file suitable for downstream processing.

The resulting FASTA file and a BIOM file generated from the operational taxonomic unit (OTU) table were used as input for PICRUSt2 analysis. The functional pathways were constructed using the MetaCyc database to provide a comprehensive functional profile of the microbial communities [85]. The relative abundances of pathways were visualized using the R software. Statistical Analysis of Metagenomic Profiles (STAMP) was employed [86] to test for significant differences in pathway contributions between groups. A chi-square test corrected with the Benjamini-Hochberg false discovery rate (FDR) was applied to control for multiple comparisons, with a significance threshold set at FDR-adjusted p-value < 0.05.

For visualization purposes, the final plots included only data with a minimum relative abundance of 0.5% to focus on the most relevant functional changes. This approach enabled a detailed evaluation of pathway-level functional shifts associated with the microbiome’s response to hypoxic stress, providing insights into the potential functional adaptations within the microbial community.

2.9. Data Availability

The nanopore data for DNA analysis were deposited in the National Center for Biotechnology Information-Short Reads Archive (NCBI-SRA) under the BioProject accession number PRJNA1099139.

3. Results

3.1. Alpha and Beta Diversity Analysis of M. chilensis Microbiota under Normoxia and Hypoxia

Alpha and beta diversity analyses of the M. chilensis microbiota under normoxic and hypoxic conditions revealed significant differences in microbial composition and structure across both studied tissues (gills and digestive gland) (Figure 1).

Rarefaction curves (Figure 1A) showed an apparent decrease in microbial diversity under hypoxia compared to normoxia, particularly pronounced in the gills. The digestive gland showed lower diversity only during the first hypoxic sampling.

Principal Coordinate Analysis (PCoA) (Figure 1B) revealed a distinct clustering of normoxic and hypoxic samples for both tissues, indicating substantial shifts in microbial community composition.

Alpha diversity indices (Figure 1C) confirmed a significant decrease in microbial diversity in both tissues under hypoxia. This decline was more pronounced in the gills, suggesting a greater sensitivity of the microbial community present to low oxygen compared to the digestive gland.

Figure 1. Alpha and beta-diversity analysis for M. chilensis microbiota exposed to normoxia (blue) and hypoxia (conditions). The figure shows the rarefaction curves for all de samples (A), the principal coordinate analysis (PCoA) (B), and different alpha diversity index estimated for gills and digestive glands under the experimental conditions (C).

3.2. Taxonomic Shifts in the Microbiota of M. chilensis Under Normoxia and Hypoxia

Figure 2 presents a heat-tree analysis, where branch color reflects the log2 of the mean ratio between bacterial taxa’s relative abundance under normoxia (blue) and hypoxia (red). This approach visualizes taxonomic changes across hierarchical levels, highlighting variations induced by oxygen conditions.

In gills, an apparent decrease in the relative abundance of specific taxa under hypoxia was observed, while others exhibited an increase. This difference is represented by branch color: taxa with decreased abundance under hypoxia are blue, and those with increased abundance are red. Dominant taxa under normoxia included Roseobacteraceae, Oceanicola, Roseobacter, Thalassobius, Phaeobacter, Octadecabacter, Sulfitobacter, Leisingera, Ruegeria, Loktanella, Tateyamaria, Roseovarius, Litoreibacter, Pelagimonas, Pacifibacter, Shimia, Actibacterium, Saggitula, Rhodobacteraceae, Aliiroseovarius, Yoonia, Pseudophaeobacter, Planktotalea, Sedimentitalea, Amylibacter, Aestuariibius, and Neptunicoccus. Conversely, taxa representative of hypoxia included Campylobacterales, Poseidonibacter, Arcobacteraceae, Campylobacteraceae, Campylobacter, Aliarcobacter, Arcobacter, Malaciobacter, Halarcobacter, Thiovulaceae, Elizabethkingia, Weeksellaceae, Saccharicrinis, Paludibacter, Crocinitomix, Crocinitomicaceae, Labilibacter, Prolixibacteraceae, Marinifilum, Olleya, Antarticcibacterium, Firmicutes, and Bacillales.

The digestive gland displayed a similar response, with differential patterns in taxonomic abundance under hypoxia. Hypoxia significantly altered the microbiota composition, with several taxa exhibiting notable changes in relative abundance. Dominant taxa under normoxia included Octadecabacter, Sulfitobacter, Thalassobius, Phaeobacter, Roseobacter, Leisingera, Antarctobacter, Roseovarius, Tateyamaria, Shimia, Pacificibacter, Litoreibacter, Pelagimonas, Aliiroseovarius, Amylibacter, Sedimentitalea, Planktotalea, Rhodobacteraceae, Lacipirellulaceae, Polaribacter, and Aquimarina. In contrast, taxa enriched under hypoxia were Acidovorax, Delftia, Curvibacter, Burkholderia, Cupriavidus, Paraburkholderia, Neisseria, Massilia, Desulfovibrio, Solidesulfovibrio, Pseudodesulfovibrio, Acidithiobacillus, Acidiferrobacteraceae, Nevskiales, Steroidobacteraceae, Ectothiorhodospiraceae, Marinobacterium, Pseudoalteromonas, Shewanella, Vibrio, Francisella, Yersinia, Serratia, Pectobacteriaceae, Enterobacteriaceae, Ancylomarina, Marinifilaceae, Marinilabiliales, Bacteroidales, and Flavobacteriales.

3.3. Analysis of Bacterial Genus Relative Abundance in the Microbiota of M. chilensis under Normoxia and Hypoxia

Heatmap analysis revealed distinct shifts in the relative abundance of bacterial genera in the gill and digestive gland microbiotas of M. chilensis under normoxic and hypoxic conditions (Figure 3 and Figure 4). The microbial communities clustered into two distinct groups in both tissues, reflecting differential responses to oxygen availability. Figures 3A and 4A illustrate these patterns, with normoxic (blue) and hypoxic (red) conditions highlighting specific genera that were differentially distributed, suggesting a tissue-specific microbial adaptation to hypoxia.

Cluster 1 consisted of bacterial genera that showed a marked decrease in abundance under hypoxic conditions. In the gills (Figure 3B), these genera exhibited a significant reduction in relative abundance during hypoxia exposure, with a similar trend observed in the digestive gland (Figure 4B). Representative genera in the gill cluster included Salmonella, Thalassobius, Roseobacter, Cocleimonas, Nitratireductor, Planktomarina, Marinicella, Yoonia, Neptunicoccus, Tenacibaculum, Cellulophaga, Leucothrix, Tritonibacter, Jannaschia, Lacinutrix, Boseongicola, Sedimentitalea, Pseudahrensia, Pseudoruegeria, and Phaeobacter. In the digestive gland, prominent genera included Thiomicrorhabdus, Mobilisporobacter, Klebsiella, Helicobacter, Spongiibacter, Paraglaciecola, Citrobacter, Shewanella, Cellvibrio, Vibrio, Solidesulfovibrio, Rheinheimera, Serratia, Aquella, Glaciecola, Franconibacter, Latilactobacillus, Oceaniserpentilla, Lacticaseibacillus, and Labilibaculum.

In contrast, Cluster 2 comprised bacterial genera that increased in abundance under hypoxia. In the gills (Figure 3C), these genera showed a significant rise in relative abundance in response to hypoxic exposure, and a similar pattern was evident in the digestive gland (Figure 4C). The predominant gill genera under hypoxia included Poseidonibacter, Arcobacter, Phocaeicola, Marinifilum, Halarcobacter, Methylothermus, Azospirillum, Labilibaculum, Bacteroides, Saccharicrinis, Francisella, Salegentibacter, Draconibacterium, Paludibacter, Elizabethkingia, Polaribacter, Pedobacter, Malaciobacter, Roseimarinus, and Aliarcobacter. The digestive gland featured hypoxia-associated genera such as Sideroxydans, Planktotalea, Marinicella, Nitratireductor, Mariniblastus, Ahrensia, Litoreibacter, Profundibacter, Pseudahrensia, Seohaeicola, Jannaschia, Pelagimonas, Lentibacter, Olleya, Thalassobius, Octadecabacter, Sulfitobacter, Amylibacter, Fuerstiella, and Aliiroseovarius.

3.4. Linear Discriminant Analysis

Linear Discriminant Analysis (LDA) revealed distinct bacterial community compositions under normoxic and hypoxic conditions in both the gills and digestive gland of M. chilensis (Figure 5).

In the gills (Figure 5A), 31 bacterial species were significantly more abundant under normoxia, including Aquimarina macrocephali, Aquimarina muelleri, and Flavobacteriaceae bacterium. Conversely, 19 species exhibited increased abundance under hypoxia, such as Poseidonibacter parvus, Poseidonibacter lekithochrous, and Arcobacter nitrofigilis.

Similarly, the digestive gland (Figure 5B) displayed 14 bacterial species enriched under normoxia, including Polaribacter sp. ALD11, Polaribacter sp. BM10, and Plaribacter atrinae. In contrast, six species showed increased abundance under hypoxia, including Francisella tularensis, Citrobacter freundii, and Shigella sonnei.

3.5. Functional Potential Prediction of the M. chilensis Microbiome Under Normoxia and Hypoxia

Picrust2 analysis revealed alterations in the functional potential of the M. chilensis microbiome between normoxia and hypoxia conditions (Figure 6). Among the metabolic pathways analyzed, only the degradation/utilization/assimilation category showed a significant difference between conditions (Figure 6A,B). Notably, both tissues exhibited a higher proportion of sequences assigned to degradation/utilization/assimilation pathways under normoxia compared to hypoxia (Figure 6C,D).

The metabolic pathways enriched under normoxia and hypoxia differed between the gills and digestive gland (Figure 6E,F). All degradation/utilization/assimilation pathways were different for both tissues, except for the TCA cycle which was more enriched than other pathways in the digestive gland. In the gills, normoxia favored pathways related to Cofactor, Prosthetic Group, and Electron Carrier biosynthesis, as well as Fatty Acid, Lipid, and Carbohydrate biosynthesis, whereas the digestive gland exhibited a predominance of the TCA cycle, Amine and Polyamine Degradation, and Aspartate metabolism. Under hypoxia (Figure 6E), the gills showed an increase in Cofactor, Prosthetic Group, Electron Carrier biosynthesis, Carbohydrate Biosynthesis, and the TCA cycle, accompanied by a decline in Fatty Acid and Lipid Biosynthesis, Nucleoside and Nucleotide Degradation, and Secondary Metabolite Degradation. In the digestive gland (Figure 6F), hypoxia induced an upregulation of the TCA cycle and Hexuronide and Hexuronate Degradation, while pathways related to Amine and Polyamine Degradation, Aspartate metabolism, S-adenosyl-L-methionine Biosynthesis, and the Methylaspartate cycle were downregulated.

3.6. Dynamics of Bacterial Pathogens

We evaluated the presence and abundance of specific groups to explore whether changes in the microbiome due to hypoxia lead to an increase in the community of pathogenic bacteria. This objective was to assess whether hypoxia promotes bivalves as reservoirs for aquatic pathogens. Figure 7 presents a scatter plot illustrating the dynamics in the relative abundance of various fish bacterial pathogens in M. chilensis under normoxic and hypoxic conditions. It employs a dual encoding of relative abundance, utilizing both circle size and color intensity, to facilitate visualization and enable a more comprehensive interpretation of the data. The area of each circular glyph represents the relative abundance of the respective variable, quantifying its contribution to the dataset. The color intensity, in turn, reflects the frequency or prevalence of that variable within the analyzed samples. This combination of visual cues—size for quantity and color intensity for quality—allows for the discernment of patterns and relationships that would be difficult to identify with a single visual variable. Consequently, this enhances interpretability and deepens the analysis.

In the gills of M. chilensis, an increase in the relative abundance of some bacterial pathogens was observed under hypoxic conditions compared to normoxia. Pathogens that showed a significant increase included Arcobacter cryaerophilus, with smaller increases seen in Citrobacter freundii and Klebsiella pneumoniae. Conversely, some bacterial pathogens in the gills of M. chilensis disappeared under hypoxic conditions compared to normoxia. These pathogens included Flavobacterium columnare, Salmonella enterica, and Tenacibaculum ovolyticum.

Similarly, to the gills, the digestive gland of M. chilensis also experienced an increase in the relative abundance of several bacterial pathogens under hypoxia. Pathogens that showed a notable increase included Aliivibrio wodanis, Pseudomonas fluorescens, Moritella marina, and Vibrio mimicus. Additionally, in the digestive gland of M. chilensis, several bacterial pathogens disappeared under hypoxic conditions compared to normoxia. Pathogens that disappeared included Flavobacterium columnare, Tenacibaculum maritimum, and Tenacibaculum dicentrarchi. Moreover, a reduction in the relative abundance of Vibrio cholerae, Pseudomonas chlororaphis, and Citrobacter freundii was also noted.

4. Discussion

This study highlights the importance of conducting controlled hypoxia exposure experiments to understand and predict changes in the microbiome of bivalve mollusks, such as M. chilensis, under environmental stress associated with climate change [15]. Hypoxia, considered a significant environmental stressor, interacts in complex ways with other environmental factors, significantly affecting the health and performance of marine organisms [87,88]. Advancing our understanding of the microbiome’s role in the physiological response of organisms to these stressors has become essential, as it can reveal critical mechanisms of adaptation and resilience [89,90,91,92]. Expanding on previous studies, this work focused on evaluating the impact of hypoxia on the gill and digestive gland microbiota of M. chilensis [15,93,94].

In this context, our study is pioneering in demonstrating the effects of oceanic hypoxia on the gill and digestive gland microbiota of mussels at the species level. The gills, in addition to their respiratory function, play a key role in nutrition and immune defense, hosting symbiotic microbial communities that contribute to carbon and nitrogen fixation [95,96,97,98]. Our results indicate that hypoxia induces profound changes in both the composition and function of the microbiome in these tissues, primarily by eliminating bacterial groups unable to tolerate low-oxygen conditions. The observed restructuring of the microbiota under hypoxia suggests adaptive and selective effects, reinforcing the microbiome’s ability to respond to environmental stress. These changes likely result from selective pressure favoring bacterial communities with greater tolerance to hypoxic conditions, possibly mediated by adaptive mechanisms such as quorum sensing, which regulate colonization, virulence, and stress resistance in oxygen-depleted environments [96,99,100,101,102,103,104,105].

However, this adaptive capacity also presents potential risks. The selective pressure exerted by hypoxia may facilitate the proliferation of opportunistic pathogens, potentially compromising the host’s immune system and its ability to resist infections and other environmental stressors. Our results provide evidence of this selective pressure, manifested in a reduction of bacterial species richness and diversity, particularly in the gills, aligning with responses observed in other aquatic ecosystems [1,106,107].

Functional changes in the microbiota of M. chilensis, especially in lipid and fatty acid biosynthesis pathways, suggest metabolic adaptations to hypoxia. These findings corroborate previous studies in other bivalves [15]. These changes are likely strategies to optimize energy efficiency and minimize biomass production in resource-limited environments [16,108].

Alterations in the functionality of the microbiota have direct implications for host physiology [109]. Reducing the production of essential metabolites, such as vitamins and amino acids, may affect the nutrition and overall health of M. chilensis [110]. Additionally, the decrease in microbial diversity and functionality could limit the microbiome’s ability to contribute to critical metabolic processes, such as nutrient digestion and the production of immunomodulatory compounds [111].

The results reveal a decrease in the activity of degradation, utilization, and nutrient assimilation pathways in the digestive gland under hypoxic conditions. This phenomenon could affect essential cellular processes, such as proliferation and differentiation, limiting the mussel’s growth and adaptation capacity in response to additional stress factors [112,113]. The complex interaction between biosynthesis, metabolite generation, and degradation pathways under hypoxia indicates a co-evolutionary adaptation process between M. chilensis and its microbiota. This symbiotic relationship is crucial for maintaining host homeostasis; any disruption in this balance may lead to dysbiosis, increasing susceptibility to opportunistic or polymicrobial infections. In aquaculture and marine environments, dysbiosis has been associated with mass mortality events and disease outbreaks, underscoring the importance of the microbiome in host health [114].

The variation in the intestinal microbiome composition of Mytilus across individuals and populations reflects the influence of diet, host genetics, and environmental conditions [115,116,117,118,119,120]. Aquaculture practices also shape this microbial structure, as seen in the differences between the microbiota of farmed and wild mussels [121]. In the gills of M. chilensis, the increased presence of bacterial groups associated with hypoxic environments is a clear consequence of the selective pressure of low oxygen availability. These adaptations may involve the production of antioxidants and the utilization of alternative metabolic pathways, indicating a complex metabolic interaction between the bacteria and the host [122,123].

The dominant microbiota in the digestive gland of M. chilensis includes bacteria from the phyla Actinobacteria, Proteobacteria, Bacteroidetes, and Firmicutes, consistent with other mussel species and aquatic organisms [15,36,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138]. In contrast, the gill microbiota is dominated by Proteobacteria and Firmicutes, suggesting specialized and complementary roles in mussel physiology. Our findings are consistent with previous studies in other aquatic environments [126]. The prevalence of Bacteroidetes and Firmicutes in the digestive tract indicates a specialized symbiotic relationship with an enzymatic arsenal for degrading complex polysaccharides, such as cellulose and chitin, found in plankton and marine detritus [125,139,140,141]. These bacteria produce various bioactive compounds, such as antibiotics and pigments, which can benefit both the host and the microbiome [142]. Antibiotics may protect the host from pathogens, while pigments may function as antioxidants, reducing oxidative stress [143,144]. Firmicutes, with their ability to form spores, contribute to the stability of the intestinal microbiota under stress conditions [145,146,147]. Additionally, some Firmicutes produce short-chain fatty acids as fermentation byproducts, which have anti-inflammatory properties and modulate the host’s immune response [148].

Hypoxia significantly altered the microbial community structure of M. chilensis, favoring facultative anaerobes and opportunistic pathogens while reducing beneficial microbial functions. Proteobacteria in the gills of M. chilensis participate in the degradation of organic matter in marine environments [149,150,151,152]. Their metabolic versatility allows them to utilize a wide range of organic compounds, significantly contributing to nutrient cycling [153,154,155,156,157,158]. However, under hypoxic conditions, a shift toward Vibrio, Aeromonas, and Desulfovibrio was observed, indicating an increased risk of pathogenicity and altered sulfur metabolism [159,160,161]. The selective pressure exerted by low oxygen availability may enhance the survival of bacteria capable of utilizing alternative electron acceptors, potentially increasing oxidative stress in the host and disrupting metabolic homeostasis [162,163,164,165,166].

The study of the microbiota associated with M. chilensis reveals the dominance of bacteria belonging to the classes Alphaproteobacteria, Bacteroidia, Epsilonproteobacteria, and Gammaproteobacteria in the gills and digestive gland. These classes exhibit diverse and ecologically relevant metabolic functions in aquatic ecosystems, participating in biogeochemical and ecological cycles [157,167,168,169,170,171]. The presence of Epsilonproteobacteria in the gills suggests an adaptation to sulfuric niches and a potential symbiotic role in the oxidation of reduced compounds [172,173]. Alphaproteobacteria and Gammaproteobacteria are key components of the M. chilensis microbiota, particularly under normoxic conditions. Their abundance may correlate with water quality, serving as potential bioindicators of dissolved organic matter from anthropogenic sources [174,175]. Gammaproteobacteria play a crucial role in the degradation of complex organic compounds, contributing to nutrient biogeochemical cycles [176,177,178].

The order Rhodobacterales, within the class Alphaproteobacteria, dominates the gill and digestive microbiome of M. chilensis. These bacteria break down organic matter, facilitating the host’s nutrient absorption [171]. Furthermore, they may modulate the immune response, enhancing the mussel’s survival in challenging environments [179,180]. A significant aspect is the Roseobacteraceae family’s sensitivity to hypoxia. The reduction of these bacteria under low oxygen conditions suggests a disruption in microbial trophic networks, potentially affecting host health [181,182,183].

The taxonomic analysis of the digestive gland reveals the dominant presence of bacterial genera Shewanella, Aeromonas, and Vibrio as key components of the core bacterial community in the digestive gland. These bacteria are known for their pathogenic potential and ability to act as reservoirs of plasmids encoding antibiotic-resistance genes [184,185,186]. The horizontal transfer of these genes, facilitated by conjugation, transformation, and transduction mechanisms, increases the risk of antimicrobial resistance dissemination, posing a significant public health concern [187,188,189,190,191,192]. This phenomenon is not limited to clinical settings but has also been documented in natural aquatic ecosystems, suggesting that these environments can serve as global reservoirs of resistance genes [193,194,195].

The increasing antimicrobial resistance, exacerbated by climate change and the indiscriminate use of antibiotics in aquaculture, highlights the urgent need to develop more effective control strategies [196,197,198,199,200,201]. Under hypoxic conditions, a reduction in the populations of Shewanella and Vibrio was observed, while Aeromonas demonstrated a remarkable ability to adapt and maintain its presence. This adaptability may be related to its diverse genetic repertoire, which includes virulence factors and the ability to form biofilms [202,203,204,205,206,207].

Pathogenic species such as Aeromonas hydrophila and Aeromonas salmonicida cause significant losses in aquaculture and pose a global health threat due to their ability to acquire antibiotic-resistance genes [191,196,208,209,210,211,212,213,214,215,216,217]. The coevolution between Aeromonas and hosts like M. chilensis suggests a symbiotic relationship where the bacteria may facilitate digestion and offer protection against other pathogens [218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236]. However, this symbiosis can be disrupted under environmental stress, such as hypoxia, promoting Aeromonas as an opportunistic pathogen.

The decrease in Shewanella under hypoxic conditions is noteworthy, considering its ability to utilize nitrate as an electron acceptor in anaerobic environments [237,238,239]. This reduction in abundance could be due to competition for nutrients or other microbial species’ production of antimicrobial compounds. The presence of Vibrio in the core microbiota of M. chilensis suggests a potential mutualistic relationship, which is consistent with other studies [230,240,241]. However, some Vibrio species are pathogenic to both bivalves and humans [240,242,243,244,245,246,247]. Although Vibrio demonstrates ecological adaptability, its dependence on oxygen reduces its prevalence under hypoxic conditions. This decrease may benefit the health of M. chilensis by reducing Vibrio-induced diseases, but it could also increase vulnerability to other pathogenic species like Vibrio mimicus. This pathogen has been associated with disease outbreaks in various bivalve species, causing high mortality rates and significant economic losses in aquaculture [248]. V. mimicus thrives in hypoxic conditions and poses a threat due to biofilm formation, antibiotic resistance, and zoonotic potential [249,250].

An increase in the genus Acinetobacter was observed in the digestive gland under hypoxic conditions. This genus is notable for its ability to persist in diverse environments and facilitate horizontal gene transfer, including those conferring resistance to multiple clinically relevant antibiotic classes. Pathogenic strains of Acinetobacter, characterized by their multidrug resistance, represent a critical public health threat in both clinical settings and natural ecosystems [212].

Our findings confirm the absence of Psychrilyobacter in farmed M. chilensis compared to wild populations, indicating that habitat characteristics significantly influence microbial composition [121]. Anaerobic conditions in natural ecosystems support Psychrilyobacter colonization, whereas suspended aquaculture systems limit such environments [153,251,252,253,254]. Reduced microbial diversity in aquaculture could impair mussel adaptation to environmental changes and increase susceptibility to pathogens over time [255].

The presence of Aquimarina macrocephali in the gills of M. chilensis under normoxic conditions suggests an adaptation to this oxygen-rich environment. This is supported by its enzymatic profile for reducing oxidative stress and its potential role in the degradation of organic matter [256]. This bacterium may contribute to the cleaning of gill surfaces, facilitating nutrient acquisition and enhancing host health. Additionally, its ability to degrade chitin suggests a possible symbiotic interaction with the mussel, as chitin is a common structural component in plankton and microorganisms that form part of the bivalve diet [141,256]. However, A. macrocephali’s resistance to multiple antibiotics raises concerns about its role as a reservoir of resistance genes, potentially facilitating the spread of these genetic elements within the marine ecosystem and the food chain [256,257,258,259,260]. Antibiotic resistance in A. macrocephali could be linked to acquiring plasmids and exposure to subtherapeutic antibiotic doses in the aquatic environment, stemming from aquaculture and wastewater [261,262].

To effectively manipulate immune functions through microbiome-based therapies, a more personalized approach will be required, one that identifies specific microorganism-host relationships [92]. In the context of bivalves, high-throughput sequencing has been instrumental in managing diseases in aquaculture [263,264]. Prebiotics and probiotics in mussel farming offer a promising approach to enhancing larval resistance to adverse effects such as climate change [263,265]. These adaptations could support the survival of populations in an evolving marine environment. Early colonization by a complex microbiota or specific symbionts can induce lasting epigenetic modifications, promoting protective immunity and greater resilience to environmental stressors associated with climate change [266,267,268]. These early epigenetic alterations can have long-term protective effects, reducing the risk of diseases in later life stages and enhancing the organism’s adaptive capacity [266,267,268]. Additionally, molecular research on the interactions between hypoxia and other environmental factors, such as temperature, pH, and salinity, is crucial for managing the impacts of hypoxia on aquatic ecosystems [264,269]. Monitoring oxygenation conditions in marine habitats will be essential to maintain the health and robustness of key species. Identifying bacterial genera sensitive or resistant to hypoxia will provide critical insights for developing sustainable management strategies.

An integrated approach based on the hologenome concept is recommended. This approach includes developing predictive models incorporating environmental, microbial, and host factors to better understand and manage the effects of hypoxia in these ecosystems. Moreover, expanding studies on the M. chilensis hologenome, which includes viruses and microeukaryotes whose identities and functions are still in the early stages of research, is necessary [270]. Implementing sustainable management practices, such as reducing eutrophication and restoring coastal habitats, could improve water quality and protect marine biodiversity.

5. Conclusions

This study represents the first detailed report on the effects of hypoxia, a global environmental concern, on the gills and digestive glands of the mussel M. chilensis. Our findings demonstrate that hypoxia significantly impacts the structure, relative abundance, community composition, species richness, and diversity of microbial communities associated with these tissues. The hypoxia-induced alteration of the microbiota has substantial implications for mussel physiology, affecting essential processes such as digestion, nutrient absorption, and immune response. Additionally, our results indicate that hypoxia promotes the proliferation of opportunistic and pathogenic bacteria, highlighting the need for further exploration of the underlying mechanisms behind these changes and their functional implications for the organism. Consequently, maintaining optimal dissolved oxygen levels in aquaculture systems is essential to preserving the health of M. chilensis and ensuring the sustainability of aquaculture practices. Future studies are necessary to develop mitigation strategies that can minimize the negative impacts of this stressor on aquatic ecosystems.