Integrated Omics Approaches Revealed the Osmotic Stress-Responsive Genes and Microbiota in Gill of Marine Medaka.

Aquatic fishes face osmotic stress continuously, and the gill is the first tissue that senses and responds to the external osmotic challenges. However, the understandings of how the gill microbiota could respond to osmotic stress and their potential host-bacterium relationships are limited. The objectives of the current study are to identify the hypotonic responsive genes in the gill cells and profile the gill microbiota communities after fresh water transfer experiment via transcriptome sequencing and 16S rRNA gene sequencing. Transcriptome sequencing identified 1,034 differentially expressed genes (DEGs), such as aquaporin and sodium potassium chloride cotransporter, after the fresh water transfer. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis further highlighted the steroid biosynthesis and glycosaminoglycan biosynthesis pathways in the gill. Moreover, the 16S rRNA gene sequencing identified Vibrio as the dominant bacterium in the seawater, which changed to Pseudomonas and Cetobacterium after the fresh water transfer. The alpha diversity analysis suggested that the gill bacterial diversity was lower in the fresh water transferred group. The KEGG and MetaCyc analysis further predicted the alteration of the glycosaminoglycan and chitin metabolisms in the gill bacteria. Collectively, the common glycosaminoglycan and chitin pathways in both the gill cells and gill microbiota suggest the host-bacterium interaction in gill facilitates the fresh water acclimation. IMPORTANCE This is the first study using the transcriptome and 16S rRNA gene sequencing to report the hypotonic responsive genes in gill cells and the compositions of gill microbiota in marine medaka. The overlapped glycosaminoglycan- and chitin-related pathways suggest host-bacterium interaction in fish gill during osmotic stress.

INTRODUCTION

Ion and water osmoregulation is critical for the maintenance of tissue and cellular functions. It is important to define cell shape, intracellular osmolality, and various cellular functions (1). Fishes have an osmoregulatory mechanism to regulate fluid and ion homeostasis to maintain a constant body osmolality during osmotic stress. The gill is the major tissue for the osmoregulatory processes. Early studies using electron microscopic analysis reported remodeling of gill cells in different salinities (2, 3). In previous decades, general characterizations of selected ion transporters and hormonal receptors, such as localization, and mRNA/protein expression levels were studied. Recent omics approaches have gained insights on the genomic or proteomic responses in different fish models (4–8). Medaka (Oryzias spp.) inhabit diverse osmotic environments worldwide (9). Due to their salinity adaptability, medaka species have been utilized in research to understand the osmoregulatory mechanisms involved in fresh water or seawater acclimation (10, 11). Marine medaka (Oryzias melastigma) is native to coastal waters. It can be found in mangroves and has the ability to adapt and spawn in freshwater (12). In the last few years, marine medaka has been used to study the environmental marine pollution (13, 14) due to its seawater habitat and well-known genome (15, 16). However, most of the salinity works in medaka species were performed in the Japanese medaka (Oryzias latipes) (10, 17–20), which lives in fresh water. The current study applied RNA sequencing to provide an overview of the hypotonic osmoregulatory mechanism in gill of marine medaka.

On the other hand, our understandings of the osmoregulatory mechanism have mainly focused on the fish itself. The effects of osmotic stress on gill bacteria are not known. Recent studies have shown the changes of taxonomic microbial compositions in gut across environmental salinity (21). The gut bacteria are suggested to play roles in nutrient absorption and immune response for host survival (22, 23). Regarding the marine medaka, a current report has demonstrated the shift of the gut bacterial communities after the seawater-to-fresh-water transfer experiment (24). Since the gill is continuously exposed to the external media, the gill bacteria must develop a mechanism to compensate for the osmotic stress. To summarize, this is the first study to integrate the transcriptomics and metagenomics approaches to understand the responses of gill and the microbiota under the hypotonic stress.

RESULTS

This study used marine medaka (Oryzias melastigma) to study the genome-wide changes of gene expression in gill and gill microbiota communities upon hypotonic stress. For the first part of the study, we conducted transcriptomic analysis to identify differentially expressed genes (DEGs) and the genome-wide molecular regulatory networks after fresh water transfer in marine medaka.

General sequencing information of the transcriptome.

Sequencing libraries were prepared from extracted gills of SW and FW and run in the BGISEQ-500 platform. The general sequencing information about the clean read quality is summarized in Table S1 in the supplemental material. Briefly, 24.96 million reads per sample, with an average mapping ratio with reference genome of 92.52% and 72.59% of the average mapping ratio with genes, were observed. Venn diagrams showed 25,903 identified transcripts, in which 24,704 transcripts were found commonly in both the control seawater (SW) and fresh water transferred (FW) samples (Fig. 1A). The DEGs with significant differences were used for the downstream bioinformatics analysis.

FIG 1

(A) Venn diagram showing the expressed gene between the seawater control and fresh water transferred gill samples. A total 25,903 transcripts were identified in which 24,704 transcripts were found in both groups. (B) Volcano plots of DEGs. DEGs identified in SW/FW; 568 upregulated DEGs and 466 downregulated DEGs were identified in the SW/FW group. Red points indicate the upregulated genes, while blue represents the downregulated genes. Nonsignificant transcripts were marked gray. (C) Significant DEGs in FW gill. Selected transporters are shown. The y axis presents the log2 fold change value.

The general sequencing information of the RNA-seq. Download Table S1, XLSX file, 0.01 MB.

DEG and enriched KEGG in freshwater transferred gill.

After the seawater-to-fresh-water transfer, 1,034 DEGs (up, 568; down, 466) were identified (Fig. 1B). The selected well-known seawater or fresh water transporters’ mRNA expression levels are shown in Fig. 1C. For example, seawater transporters, such as cftr, nkcc1, and nkcc2 (5, 25), were suppressed after the fresh water transfer. On the other hand, the fresh water transporters like aqp3, nhe1, and clcn2 (5, 26, 27) were upregulated. The full list of the DEGs is in Table S2. All the DEGs then underwent the GO analysis. The top enriched GO terms in biological process (BP), cellular component (CC), and molecular functions (MF) are shown in Fig. 2A. For BP, metabolic process, response to stimulus, and biological adhesion were found. Osmotic stress stimulated the cell and triggered the modifications of the adhesion proteins for osmoregulatory mechanisms. Cell migration is a suggested phenotype during the osmoregulatory progresses (28). Cells have to modulate cell–cell adhesion and gain motility to move (29, 30). Epithelial-mesenchymal transition (EMT) is a critical process to change cell adhesion and polarity (31–33). Such progress involves the increase of mesenchymal markers but decrease of epithelial markers (34, 35) that could explain the reason for identifying similar gene enrichment counts on biological adhesion in this study. For the CC and MF, membrane and transporter activity were ranked as one of the top enrichment terms that reflected the importance of the modification of gill ion transporters that were located in the cell membrane. The full list of the enriched GO terms is shown in Table S3. When we presented the BP category data as a directed acyclic graph (DAG), the chitin metabolic process was significantly enriched via amino sugar metabolism and aminoglycan metabolism (Fig. 2B). The details of the DAG are in Fig. S1.

FIG 2

(A) GO classification of DEGs between groups. GO terms identified in SW/FW group. y axis shows the number of DEGs, and the x axis indicates the GO terms that were classified into biological process, cellular component, and molecular function. Red bars indicate the upregulated genes, while blue indicates the downregulated genes. (B) Directed acyclic graph (DAG) of the selected enriched BP terms. Chitin metabolic process was significantly enriched after fresh water transfer via amino sugar metabolism and aminoglycan metabolism. (C) Pathway functional enrichment of DEGs in SW/FW group. x axis represents the enrichment factor, while y axis indicates the enriched pathways. The color indicates the q value: the blue color represents the lower value, while the white color indicates the higher value. In addition, the point size indicates the number of DEGs. The greater the rich factor, the more significant the enrichment. Glycosaminoglycan and steroid biosynthesis were found to be enriched in the SW/FW group.

List of DEGs after fresh water transfer. Download Table S2, XLSX file, 0.2 MB.

List of enriched GO terms after fresh water transfer. Download Table S3, XLSX file, 0.1 MB.

Directed acyclic graph (DAG) of the enriched BP category after fresh water transfer. Download FIG S1, PDF file, 0.03 MB.

KEGG pathway analysis was performed to show the functional enrichment from the DEGs. The full lists of the enriched terms (level 2) are shown in Table S4. Pathways such as fatty acid metabolism, glycosphingolipid biosynthesis, glycosaminoglycan biosynthesis, and steroid biosynthesis were enriched (Fig. 2C). The KEGG-DEG relationship network was drafted to have a more detailed understanding of the whole regulatory system in both groups. The network clearly showed that the metabolic pathway was linked with the mentioned glycosphingolipid and steroid biosynthesis. Readers can refer to the PDF file (Fig. S2).

List of enriched KEGG pathways after fresh water transfer. Download Table S4, XLSX file, 0.05 MB.

The KEGG-DEG network after fresh water transfer. Download FIG S2, PDF file, 1.7 MB.

Hypotonic stress alters the gill microbial diversity in marine medaka.

The 16S rRNA gene sequencing in gill of marine medaka was performed. A total of 399 sequences (246 in SW, 196 in FW) were identified by the amplicon sequence variant (ASV) method. Among them, 43 sequences showed overlap between SW and FW (Fig. 3A). Alpha diversity was used to analyze the complexity of species diversity in the medaka fish gill (36). Various indexes are shown in Fig. 3B. For example, the Shannon index reflects species diversity of the community, which is affected by both species richness and evenness. Our results showed that FW could lead to a significant reduction of gill microbial diversity compared to SW (Fig. 3B). Moreover, the calypso analysis was performed to obtain a higher-level analysis of microbial community composition (37). The network relation is shown in Fig. 3C. The control SW is shown in red, while the FW data are shown in blue. The result demonstrated that the microbiota communities were changed after the transfer, in which Vibrio was spotted in SW and Pseudomonas in FW. Cetobacterium was found in both SW and FW. Phylogenetic tree diagram showed that certain microbiotas in FW overlapped with SW (Fig. 3C, right bottom).

FIG 3

(A) Venn diagram of the ASV detected in the gill of SW control group and the FW group. Different colors represent different groups; the intersection represents the set of ASV commonly present in the counterpart groups. Of the 399 ASV recognized (246 in SW, 196 in FW), 43 were found in common in SW and FW. (B) Various alpha diversity measurements of the two groups. Results indicated that the transfer of fish from SW to FW could lead to significant reduction of gill microbial diversity. (C) The network relationship among the two groups was obtained from Calypso analysis. Vibrio was mainly found in SW (red spot), while Pseudomonas was in FW (blue spot). Cetobacterium was found in both SW and FW. Tree diagram at the right bottom shows that the FW (blue spot) shares some microbiota with SW (red spot).

Osmotic stress triggers changes in the fish gill microbial taxonomic composition.

Distinct diversity patterns between the control marine medaka (SW, red spot), and the fresh water transferred medaka (FW, blue spot) can be seen after the canonical correspondence analysis (CCA) (Fig. 4A). At the phylum level, the control marine medaka contained Proteobacteria and Fusobacteriota. After the fresh water transfer, these two bacteria could still be identified in the samples (Fig. 4B). At the genus level, Vibrio (pink spot) was the dominant bacteria in the control seawater gill, and the Pseudomonas (blue spot) was found mainly in the FW gill. Cetobacterium (light green spot) could be found in both SW and FW (Fig. 4C). Regarding the abundances, Vibrio (∼55%) was the dominant bacteria at the genus level in the control seawater gill. When the fish was transferred to fresh water, Vibrio was reduced to about 5%, and Pseudomonas was increased to 17% in FW. Cetobacterium maintained its abundance at around 24% in SW and 37% in FW (Fig. 4D). The raw data of the samples can be found in Table S5. To understand the origins of rising FW bacteria during the transfer, we collected the rearing water samples that underwent the sequencing. The principal coordinate analysis (PCoA) result indicated that the SW gill microbiota (orange) and FW gill (green) were located apart from the seawater (yellow) and the fresh water (purple) (Fig. 4E and Fig. S3). Such result suggested that the changes of external aquatic microbiota composition were not the major factor contributing to the shift of gill bacteria. Such a result matched the notion that the new external aquatic microbiotas were not able to complete with the well-established microbes (38). Lastly, volcano plots were further used to show the microbiota composition differences between SW and FW. At genus level, 14 bacteria had significant differences (Fig. 5A). The hypotonic stress resulted in the decrease of Vibrio (red underlined) and increase of Pseudomonas (blue underlined) abundances (Fig. 5B). The full list and statistical data can be referred to (Table S6). Lastly, the KEGG and MetaCyc analyses were performed to identify the possible functional differences between microbiota in SW and FW. In the KEGG analysis, only the glycosaminoglycan degradation was identified in the SW/FW group (Fig. 5C). It should be noted that all the enriched pathways were identified in the SW group, which indicated the gill microbiota after hypotonic stress have relatively lower metabolic activities. Pathways such as chitin derivative degradation and methionine biosynthesis were found (Fig. 5D). The full list with static data can be referred to (Table S7).

FIG 4

(A) Canonical correspondence analysis (CCA) was used to show the microbiota communities under the two conditions. Distinct compositions were found in the control marine medaka (SW, red spot) and progressive-transfer medaka (FW, blue spot). (B) The microbiota community distribution at phylum level. Proteobacteria was mainly found in the samples. (C) The microbiota community distribution at genus level. Vibrio (pink spot) was mainly found in SW, while Pseudomonas (blue spot) was in FW. Cetobacterium (light green spot) was found in both SW and FW. (D) The relative abundance of selected gill bacteria at genus level. Vibrio was highly present in the SW gill. When the fish was transferred to fresh water (FW), Vibrio was reduced and Pseudomonas was increased. (E) PCoA of rearing water and gill samples. The SW (orange) and FW (green) gill microbiota were located apart from the seawater (yellow) and the fresh water (purple).

FIG 5

(A) Comparison of microbiota at genus level between the gill of SW control group and the FW group. y axis shows the genus, and the x axis represents the log2 fold value. Vibrio was mainly found in SW, while Pseudomonas was dominant in FW gill. (B) Volcano plot of the microbiota. Red indicates the significant changes in abundances, while green refers to changes without statistical significance. SW was used as the reference. Vibrio was reduced (red underline) after the fresh water transfer, while Pseudomonas (green underline) increased its abundance. (C) KEGG bioinformatics analysis on the gill microbiota. Only the glycosaminoglycan degradation was enriched significantly in SW group. (D) MetaCyc analysis showed the 20 enriched pathways with significant changes. The bar plot shows mean proportions of differential MetaCyc pathways. The difference in proportions between the groups is shown with 95% confidence intervals. Methionine biosynthesis and chitin derivative degradation were identified that may play roles with the gill cell in osmoregulation.

Raw data of the abundances of microbiota at each level. Download Table S5, XLSX file, 0.2 MB.

Bar chart showing the relative abundance of microbiota of water samples and gill samples at genus level. Download FIG S3, PDF file, 0.3 MB.

Raw data of the significant changes between SW/FW. Download Table S6, XLSX file, 0.02 MB.

Raw data of the MetaCyc analysis. Download Table S7, XLSX file, 0.02 MB.

DISCUSSION

Fish gill is the first and major tissue that contributes to osmoregulation. Marine medaka has the ability to acclimate in both seawater and fresh water, which indicates its gill undergoes rapid osmoregulatory mechanisms to control ion loss and water gain during hypotonic stress.

Modification of selected ion transporters/water channel mRNA expression level upon hypotonic stress.

Numerous genes were changed upon the hypotonic stress in gills. Our data clearly showed that the two well-known fresh water channels, chloride channel (clcn2) and water channel (aqp3), were significantly induced in FW gill. It is known that there are two major gill cells in fishes, which are the pavement cells and the mitochondrion-rich cells (ionocytes) (39). Ionocyctes are believed to play major ion regulatory roles during osmoregulation (40). Fresh water acclimation requires the activation of sodium uptake and acid excretion, chloride uptake, and base secretion (1, 17, 40, 41). Clcn2 belongs to the CLC family of chloride channel/transport proteins. It involves the chloride efflux pathway and takes part in chloride uptake (42), and clcn2 in eel ionocytes has a higher mRNA expression in fresh water (17, 40). Aqp3 is the water channel and was found to be induced in eel gill after fresh water transfer (25). Moreover, sodium hydrogen exchanger 1 (nhe1) was also upregulated in FW, which is involved in sodium ion and proton regulation, known to be highly expressed in fresh water environments (25, 27). In addition, it plays important roles in gaseous exchanges and the acid-base balance in fishes (43–46). On the other hand, seawater transporters like sodium potassium chloride cotransporter, nkcc1 (slc12a2), and nkcc2 (slc12a1) decreased their expression about 2- and 4-fold in FW, respectively. In the seawater environment, nkcc1 and cftr were highly expressed in the seawater-type ionocytes (1, 17, 40, 41, 47). In eel gill, nkcc1 was found upregulated after seawater transfer (25, 27), while nkcc2 was suggested to play a transition role in hyperosmotic transfer in red drum (48). Furthermore, cftr decreased its mRNA expression in FW and matched other publications on fishes (5, 14, 25). Collectively, our findings in whole gill matched those of other publications on fishes and confirmed the reliability of our sequencing data.

Enriched hormonal and metabolic pathways upon hypotonic stress.

Steroid biosynthesis and glycosphingolipid biosynthesis were highly enriched after fresh water transfer, indicating their possible chronic roles in fresh water environments. A study in Pacific white shrimp also identified these two pathways upon chronic low-salinity stress (49). Steroid hormones are involved in various cellular functions, such as control metabolism, inflammation, immune functions, and salt and water balance (50). For example, season-dependent changes in immune status and activities of immune cells in fishes were suggested to correlate with changes in the levels of circulating sex hormones (51, 52). Enrichment of steroid hormone biosynthesis pathways during osmoregulation was found in other aquatic organisms, such as cetaceans and white shrimp (49, 53). Moreover, it has been suggested that the salinity tolerance capacities were related to differential expression in immune response genes (54). Various transcriptomics studies further demonstrated the immune-related genes are altered in fresh water environments (55, 56).

On the other hand, glycerophospholipids are glycerol-based phospholipids that are the main component of biological membranes (57). They play roles in protecting the cell surface by maintaining the stability of the plasma membrane via modifying the plasma membrane lipid bilayer (58, 59). To maintain and modify the membrane structure, it is suggested that extra energy is needed, and the identified enriched fatty acid biosynthesis could provide such energy for osmoregulation compensatory processes (60, 61).

Relationship between fish and bacteria.

Fish have unique and direct interaction with the surrounding water. They contact microorganisms throughout their lifetime, a relationship that may be beneficial or pathogenic. The gut microbes have been extensively studied, and researchers identified the host-microbial relationship that contributed to nutritional provisioning, metabolic homeostasis, and immune response (62, 63). Diet composition and osmotic stress were the factors that alter the gut microbiota composition (24, 64, 65). On the other hand, studies of gill bacteria were extremely limited and were mainly focused on the pathological infection issues (66, 67). Currently, a study of gills of reef fish has suggested that the gill microbiome composition differed significantly from that of the gut for both adults and juveniles across 15 teleost fish families (68). However, to our knowledge, there is no study on the effects of osmotic stress in gill microbiota. The adaptive microbial community shifts concomitant with the host habitat change may contribute the necessary physiological changes required for host survival (22).

In the first part of the study, we suggested that the enriched fatty acid biosynthesis provides extra energy for osmoregulation in gill cell. A study showed that the enhanced sphingolipid synthesis could improve osmotic tolerance in bacteria (69). Thus, it is reasonable to link the host-bacterium interaction in the fish gill. In the second part of the study, we aimed to identify the effects of hypotonic stress on gill microbiota communities and further determine the common enriched pathways that were shared by the gill cells and its bacteria.

The core microbial habitat in fish gill and its osmosensing mechanism.

Studies in fish guts have identified that Proteobacteria is a major phylum in numerous species and were found to compose 90% of the gut microbiota (22, 70–72). Vibrio was found to be the most dominant reported microbiota in marine fish, while Pseudomonas and Cetobacterium were found in the fresh water fishes (24, 63, 73, 74). Here, our result showed that the dominant microbiota in the external tissue (gill) was similar to that of the internal organ (gut). Vibrio, as the dominant bacterium in the gill of marine medaka, was suggested to have an evolutionary association with marine fish (75). It has been shown to produce hydrolytic enzymes for breaking down dietary carbohydrates and lipids (76, 77). On the other hand, Cetobacterium isolates under the phylum Fusobacteria were increased after fresh water transfer (FW). It has been detected in different freshwater fish guts (71, 73) and was shown to produce vitamin B12 (78, 79).

Osmoregulatory changes for acclimation can be divided into two major stages, which are the adaptive period with changes in osmotic parameters and the chronic regulatory period for achieving homeostasis (80). During the adaptive phase, changes in external osmolality could induce water fluxes across the cytoplasmic membrane that lead to the modification of intracellular environments such as cellular hydration, molecular crowding, magnitude of turgor, and cellular integrity (81). For the prolonged osmoregulatory mechanism, metabolic pathway activation plays a major role. The enriched methionine metabolism in this study is linked to folate cycle, and the one-carbon metabolism of the histone-like nucleoid structuring protein is a key regulator in sensing the change in osmolality (82).

Shared enriched chitin-related pathways in gill DEGs and gill bacteria.

The study tried to identify the shared pathways that were enriched in both gill cells and gill bacteria. The shared glycosaminoglycan and chitin metabolism (Fig. 2B and C and 5C and D) may provide hints for the bacterium-host relationship in gill during osmotic stress. Chitin is a cellulose-like biopolymer that is mainly found in marine invertebrates, insects, fungi, and yeasts (83). Chitinase-like proteins have been found to be involved in immune response (84). A previous study in crab demonstrated that the gills could participate in chitin degradation and may prevent pathogen infection (85). However, its function in fish gill is not known. On the other hand, glycosaminoglycan and aminoglycan metabolic processes were shown to be enriched in various cancers (86, 87), similar to chitin, their roles in fish gill osmoregulation are poorly understood. It is not possible to confirm the osmoregulatory functions of these pathways by the current bioinformatics analysis. Nevertheless, the study here suggested chitin and glycoaminoglycan metabolism are important in fish gill osmoregulatory events. Further investigation on this group of genes in the future might discover their new roles in fish osmoregulation.

Conclusions.

This is the first study using the transcriptome and 16S rRNA gene sequencing to report the hypotonic responsive genes in gill cells and the compositions of gill microbiota in marine medaka. The overlapping glycosaminoglycan- and chitin-related pathways suggest the host-bacterium interaction in fish gill during osmotic stress.

MATERIALS AND METHODS

Fish maintenance and experimental setup.

Six-month-old marine medaka (O. melastigma) was maintained in seawater at 26°C. Thirty fish were transferred to 50% seawater for 7 days, followed by fresh water for another 7 days before gill isolation for RNA sequencing and 16S rRNA gene sequencing. The transfer experiment was based on our previously published protocol (24). Seawater-to-seawater transfer was performed in another 30 fish as the control group (SW). Fish were kept in one tank for each condition after the transfer. Gill samples were extracted after the transfer. All experimental protocols were approved by the ethics committee of Kyushu University, Japan (A19-165-1).

Library construction and Illumina RNA-seq.

Total RNA from seawater and fresh water transferred gills was extracted for transcriptome sequencing (RNA-seq) by TRIzol reagent (Life Technologies, CA, USA). Five gills (one side) were pooled as one biological sample, with each group containing four replicates for library construction. The sequencing method was based on our previously published protocol (5). Briefly, RNA concentrations were measured using a Qubit RNA assay kit on a Qubit 2.0 fluorometer (Life Technologies, CA, USA); 300 ng total RNA with an RNA integrity number of >8 was used for library construction (Agilent Technologies, CA, USA). The Agilent 2100 Bioanalyzer system was used for qualification. cDNA libraries were prepared using the TruSeq stranded mRNA LT sample prep kit (Illumina, San Diego, USA) per the protocol. Index codes were ligated to identify the individual samples. mRNA was purified from total RNA using poly(T) oligonucleotide-attached magnetic beads (Illumina, San Diego, USA) before fragmentation. cDNAs were then synthesized by using random oligonucleotides and SuperScript II with DNA polymerase I and RNase H treatment. Overhangs were blunted by treatment with exonuclease/polymerase followed by 3′-end adenylation and ligation to Illumina PE adaptor oligonucleotides. DNA fragments with adaptor molecules on both ends were enriched by using the Illumina PCR primer cocktail for 15 PCR cycles. Products were purified using the AMPure XP system and quantified by the Agilent Bioanalyzer 2100 system and then sequenced by the BGISEQ-500 platform. Sequence reads were filtered by SOAPunke software (v1.5.2) to remove reads with adaptors, >0.1% unknown bases (N), and low-quality reads (i.e., percentage of bases with quality of less than 20 is greater than 50% in a read).

Sequencing reads were mapped to the reference genome using Bowtie2 (v2.2.5) (88) with the parameter settings “-q –phred64 –sensitive –dpad 0 –gbar 99999999 –mp 1,1 –np 1 –score-min L,0,-0.1 -p 16 -k 200,” and then we calculated gene expression level with the RSME software package (89) with default settings to estimate gene and isoform expression levels from RNA-seq data. The DEGs were detected by PossionDis with a fold change of ≥2.0 and false discovery rate (FDR) of ≤0.001 (90). GO and KEGG analysis were performed by using phyper, a function of R.

16S rRNA metagenomic sequencing.

Five gill samples (the other side of the same fish) were pooled as one sample, with each group containing four replicates for metagenomic 16S rRNA gene sequencing. Rearing seawater and fresh water were collected at the end of the experiment. Bacterial DNA was extracted by using the DNeasy blood and tissue kit (Qiagen, Hilden, Germany). Bacterial genomic DNA was then collected and quantified using the Qubit dsDNA HS assay kit (Life Technologies, Carlsbad, CA, USA) as previously described (91).

A 30-ng sample of genomic DNA was subjected to amplicon PCR reaction to amplify the DNA fragment flanking the V3 and V4 regions of the 16S rRNA gene as previously described (91). Primer sequences were the following: forward primer, 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACG GGN GGC WGC AG-3′; reverse primer, 5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTA CHV GGG TAT CTA ATC C-3′. PCR was performed as initial denaturation at 95°C for 3 min, followed by 25 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, elongation at 72°C for 30 s, and a final elongation step at 72°C for 5 min. The 16S V3 and V4 amplicon was purified from free primers and primer-dimer species using Ampure XP beads (Agencourt Bioscience, Beverly, MA, USA). The library was quantified using real-time quantitative PCR and quality checked using an Agilent 2100 bioanalyzer instrument (EvaGreen; Santa Clara, CA, USA). The normalized library was then sequenced by the BGI sequencer platform.

Bioinformatics analysis, data processing, and ASV prediction.

Bioinformatics analysis was carried out by following a workflow that includes (i) quality assessment and preprocessing of raw sequence reads with adapter, primer, and base quality trimming, (ii) custom reference database generation, (iii) amplicon sequence variant (ASV) prediction and taxonomy assignment, (iv) statistical analysis, and (v) visualization. Unless stated otherwise, data handling, statistics, and visualization were performed using custom-made Python, Perl, and R scripts.

The quality of raw sequence data was assessed based on the presence of adapters, primers, and low-quality bases. The presence of these in the data was addressed by trimming the bases off from both ends of the reads using either standalone or a combination of FastQC v0.11.8, Cutadapt v2.10, and Trim Galore v0.6.6 while maintaining a minimum Phred quality score of 20 and length of 150 bp for paired-end (PE) reads. A taxonomy classifier based on 341F-805R universal primer (92), which corresponds to the V3-V4 hypervariable region in the 16S small subunit ribosomal DNA, was built from Silva 138.1 SSU NR99 reference database (Silva 138) (93–95) using RESCRIPt plugin in Qiime v2 (Qiime2), 2020.6 release (96). In brief, Silva 138 was processed starting with trimming of low-quality sequences, filtering sequences by length and taxonomy (900, 1,200, and 1,400 bases for archaea, bacteria, and eukaryotes, respectively), dereplicating sequences and taxonomy by lowest common ancestor (LCA), and extracting amplicon-specific region classifier by Qiime2 feature-classifier extract-reads function. Meanwhile, amplicon sequences were processed and analyzed by Qiime2 pipeline involving steps that import, merge, denoise, and classify features from PE reads and associated metadata into a collection of Qiime2 artifact and visualization files. While —p-min-fold-parent-over-abundance was set to 2, no further trimming and truncation of sequences were carried out during denoising. A quality control step was carried out on denoised sequences to address un- or poorly characterized features using Qiime2 “quality-control exclude seqs” module at 97% identity threshold and 95% query alignment, searching against full-length Silva 138 with vsearch. Features classified as chloroplast or mitochondria were excluded from the data set. Phylogenetic analysis was performed on ASV using Qiime2 alignment and phylogeny functions, generating a midpoint-rooted tree. Features represented as a table, sequence, distance matrix, phylogenetic tree, and biom files were exported individually and combined into a phyloseq object using qiime2R for downstream analysis.

Taxonomic and functional analysis.

Taxonomic profiling was carried using phyloseq package (97) in R 4.0.2. The resulting phylogenetic tree was visualized using the plot_tree function in the phyloseq package and was annotated with ggtree package (98, 99) in R. Significant difference on alpha diversity indices between SW and FW groups was tested using one-way analysis of variance followed by post hoc Tukey honestly significant difference at a 95% confidence level. Likewise, significant differentially abundant taxa at various levels between treatments were identified using DESeq2 (100). Inferences derived from DESeq2 were based on Wald test for two groups with a parametric fit for dispersion and Benjamini-Hochberg correction for multiple testing at a 95% confidence level. On the other hand, functional profiling was performed using PICRUSt2 v2.0.0-b with built-in EC/MetaCyc and KEGG/KO databases (101–103). MetaCyc and KEGG pathway abundances were further analyzed using statistical analysis of taxonomic and functional profiles (STAMP) v2.1.3 (104). Statistical significance in pathway enrichments was calculated between treatments based on White’s nonparametric t test (two-sided) with bootstrapping at 95% confidence intervals and Storey’s FDR for multiple testing at the 5% significance level. Finally, cooccurrence network analysis at the family and genus levels was done using Calypso v8.84 with default parameters (37).

Data accessibility.

The sequencing data from this study have been submitted to the NCBI BioProject (https://www.ncbi.nlm.nih.gov/bioproject) under the accession numbers PRJNA702883 and PRJNA588335.