Experimental fish were collected from aquaculture facilities located in Xiangshan Bay, Zhejiang Province, China. Before initiating the experiments, 120 individuals were randomly selected and allowed to acclimate in a tank (3 m diameter) with 25‰ salinity seawater for 30 days. Next, 60 similarly-sized healthy fish (body weight 69.49 ± 15.20 g, body length 168.82 ± 12.99 mm) were randomly selected and transferred to six tanks (1 m diameter; 10 individuals per tank with 500 L 25‰ seawater) to adapt to the experimental conditions for 7 days. During this period, the fish were fed daily with commercial mixed feed and the water was exchanged once per day to eliminate residual feed and excrement. After this secondary adaptation period, all of the fish from each tank were randomly assigned to two salinity treatments [5‰ (S-group) and 25‰ (C-group)] in triplicate. During the experiment, the temperature and DO were 21–23°C and 8.3–8.9 mg L⁻¹, respectively. The fish were fed once per day with commercial feed and 1/3 of the water was exchanged daily. After 7 days of exposure, 5 individuals from each tank were randomly selected and anesthetized with 50 mg/L MS-222, after which the gill tissues were dissected and immediately flash-frozen in liquid nitrogen. The samples were divided into five sections (technical triplicates were obtained for each section) to conduct gene expression, methylation sensitive amplification polymorphism (MSAP) analysis, RNA-seq, WGBS, and Mass Array analysis.

Genomic DNA from gill tissues was isolated using the TIANamp Marine Animals DNA Kit (TIANGEN, Beijing, China) according to the manufacturer’s instructions. DNA concentration and integrity were examined using a NanoDrop 8000 spectrophotometer and 1.5% agarose gel electrophoresis, respectively.

MSAP experiments were conducted as described by Xiong et al. (1999). Gill DNA samples were double digested by (EcoRI + HpaII) or (EcoRI + MspI). Each reaction was conducted using the following components: genomic DNA, 400 ng; 10× T4 Buffer, 2 μl; 20 μM EcoRI, 0.4 μl; HpaII (or MspI), 0.4 μl; 0.1% BSA, 0.2 μl; sterile water was added to a final volume of 20 μl. The mixture was incubated at 16°C for 12 h. Afterward, the ligated DNA was pre-amplified in a 20 μl reaction system as follows: digested product, 2 μl; 10 mM dNTPs, 0.4 μL; 10× Buffer, 2 μL; 5 U/μL Taq polymerase, 0.2 μL; 10 μM pre-amplification E-A primer, 0.5 μL; 10 μM pre-amplification HM-T primer, 0.5 μL; sterile water, 14.6 μL. PCR was conducted under the following thermal profile: 94°C denaturation for 5 min; 26 cycles of 94°C for 30 s, 56°C for 1 min, and 72°C for 1 min; 72°C extension for 10 min. The pre-amplification products were 20-fold diluted with sterile water and used as the template for the selective amplification reactions. The selective amplification reactions system had the following components: 20-fold diluted pre-amplification product, 2 μL; 10 mM dNTPs, 0.4 μL; 10× Buffer, 2 μL; 5 U/μL Taq Polymerase, 0.2 μL; 10 μM E-primer, 0.5 μL; 10 μM M-primer, 0.5 μL; sterile water was added to reach a 20 μL volume. The reactions were conducted using the program described by Li et al. (2017). The sequences of the adaptors and pre-amplification primers are summarized in Table 1 . The selective amplification products were separated and examined using capillary electrophoresis. The GeneMarker software (version 2.2) was used to analyze and visualize gill MSAP data. Data analyses were performed as described by Li et al. (2017).

The expression of genes related to DNA methylation (dnmt1, dnmt3bb, tet1 and tet2) was analyzed through qRT-PCR. The primers were designed using the Primer3 program ( Table 1 ). qRT-PCR was performed using the TB Green Premix Ex Taq II Kit (TaKaRa, Dalian, China). The amplification program was set as indicated by the manufacturers. The 18S rRNA gene was used as a reference.

The total RNA of gill tissues was extracted from six large yellow croaker individuals from the S-group and C-group in triplicate (i.e., three fish per group). The concentration and integrity of the total RNA samples were measured using Qubit^® RNA Assay Kit in a Qubit^® 2.0 Fluorometer and 1.5% agarose gel electrophoresis, respectively. Afterward, a strand-specific RNA library was constructed using the NEB Next Ultra™ RNA Library Prep Kit (NEB, MA, USA), and the obtained library reads were purified using the AMPure XP system (Beckman Coulter, CA, USA). After library preparation, 150 bp paired-end sequencing was performed on a NovaSeq 6000 sequencer. Differential expression analysis was performed using DESeq2 (Love et al., 2014), with a false discovery rate (FDR) ≤ 0.05, and |log₂(fold change)| > 1. GO and KEGG pathway enrichment analyses of DEGs were performed using ‘GOseq’ R package (Ashburner et al., 2000) and KOBAS software (Xie et al., 2011), respectively.

Genomic DNA of six gill tissues from the S-group and C-group were extracted using the TIANamp Marine Animals DNA Kit (TIANGEN, Beijing, China). The concentration and integrity of genomic DNA was examined using a Nanodrop ND2000 spectrophotometer (Thermo Scientific, MA, USA) and 1% agarose gel electrophoresis, respectively. Genomic DNA samples were fragmented to 200–300 bp using a Covaris S220 sonicator (Covaris, MA, USA). Afterward, the DNA samples were bisulfite-treated according to the EZ DNA Methylation-Gold kit (Zymo Research, CA, USA) instructions. After bisulfite treatment, un-methylated cytosine became uracil, and methylated cytosine remained unchanged. After end repair and adenylation reaction, the bisulfite DNA fragments were subjected to PCR amplification to construct the sequencing library. After quantifying the library and examining the inserted DNA fragments using a Qubit 2.0 fluorometer and an Agilent 2100 Bioanalyzer, qRT-PCR was conducted to determine the final concentration of the library (library concentration >2n M). The qualified DNA libraries were then sequenced using an Illumina HiSeq 2500 sequencer.

The raw reads were quality checked and filtered to obtain clean reads. Then, the sequencing reads were compared to the reference genome of large yellow croaker (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/972/) using the Bismark software (version 0.12.5) (Krueger and Andrews, 2011). After bisulfite treatment, the unmethylated cytosine (C) converted into thymine (T), and the methylated C remain unchanged; thus, the reads cannot be located accurately, and the methylated sites cannot be identified in the reference genome. However, using in silico bisulfite conversion algorithm of the Bismark software, we can comparatively analyze the mapping result of the bisulfite treated reads and no bisulfite treated reads to definite the methylated C sites by changing the entire C base in the reference genome into T base. Next, the methylated sites and differentially methylated regions (DMRs) identified using the MethylKit software (version 0.99.2) (Akalin et al., 2012) were compared with the results obtained with the Bismark software. In order to analyze the methylation status of different gene elements, we carried out a profiling analysis. The analytical procedure is dividing the gene elements into 20 bins first, and then comparing the methylation level of each bin to reflect the methylation difference of each gene elements of each gill samples of S-group and C-group. Moreover, the number of methylated and demethylated DMRs in different gene elements and chromosomes were statistically analyzed to reflect the methylation dynamics in gill under low-salinity stress. Upon characterizing the distribution of DMRs, differentially methylated genes (DMGs) were identified by detecting which genes overlapped with the DMRs. GO and KEGG analyses of the identified DMGs were performed using the ‘cluster Profiler’ R package (Yu et al., 2012). Significant differences were identified via Student’s t-test (p < 0.05).

To evaluate the correlation between DNA methylation and gene expression in the gills of large yellow croaker submitted to low-salinity stress, Pearson correlation analyses were conducted to screen the transcriptome data and WGBS data, and identify the mutual DMGs among the DEGs in RNA-seq and the DMR target genes in WGBS. Afterward, GO and KEGG enrichment analyses of the obtained DMGs were conducted using the ‘GOseq’ R package and the KOBAS software, respectively.

10 DEGs (arhgef37, gnb4, grtp1, LOC104929349, LOC104928056, LOC104918888, LOC104935555, pdhb, slc2a4, and slc9a3) and 5 DMGs (gnb4, LOC104929349, LOC104928056, slc2a4, and slc9a3) were selected to verify the reliability of the RNA-seq and WGBS data, respectively. The validation of the DEGs and DMRs (which were significantly enriched in ion channels, energy metabolism, cytoskeleton, and other key gene pathways) was conducted using qRT-PCR and Mass Array methylation analysis. qRT-PCR was carried out as described above. For the Mass Array procedure, genomic DNA was isolated using the QIAamp DNA Mini Kit (QIAGEN, Dusseldorf, Germany), and bisulfite modification of the genomic DNA was carried out using an EpiTect Bisulfite Kit (QIAGEN, Dusseldorf, Germany). Quantitative methylation analyses of the 5’ promoter of the DMGs were performed using the Sequenom Mass Array platform (Oebiotech, Shanghai, China) according to the manufacturer’s instructions. The primers used for qRT-PCR and Mass Array methylation analysis are summarized in Table 1 .

Gene expression differences between the S-group and C-group were identified via two-way ANOVA followed by Bonferroni’s post hoc test. The differences in the methylation levels between the S-group and C-group were identified via one-way ANOVA with Duncan’s multiple range tests. A p-value <0.05 was considered statistically significant.

To clarify the DNA methylation status in the gills of large yellow croaker under low-salinity stress, we first detected the expression of DNA transmethylase genes (dnmt1 and dnmt3bb) and demethyltransferase genes (tet1 and tet2). The expression levels of the four genes were significantly increased in the S-group compared to the C-group ( Figure 1A ). Moreover, we examined the methylation status in gill tissues using the MSAP technique. The hemi-methylation, full-methylation, and total methylation levels of gill genomic DNA in the C-group were 26.92%, 23.68%, and 50.60%, whereas the levels of the S-group were 29.19%, 27.41%, and 56.60%, respectively ( Figure 1B ). All three methylation indices were significantly different between the S-group and C-group. The analyses of the genomic DNA CpG methylated sites indicated that 33.36% CpG sites remained unchanged, whereas 32.83% CpG sites were methylated and 33.81% were demethylated in the S-group compared to the C-group ( Figure 1C ).

After filtering low quality reads, a total of 263 million clean reads were generated, and at least 88.75% clean reads were uniquely mapped to the reference genome ( Table 2 ). To clarify the effects of low-salinity stress on gene expression in the gills, principal component analysis (PCA) was conducted between the S-group and C-group individuals. As expected, the results showed a high similarity in gene expression among individuals under the same salinity treatment, and significant differences between the S-group and C-group ( Figure 2A ). Moreover, comparative analysis of FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values indicated that the level of gene expression in the S-group individuals was significantly lower than that of the C-group individuals (p < 0.05) ( Figure 2B ). Furthermore, the gene expression levels in each of the gill samples were analyzed based on the range of the FPKM value ( Figure 2C ), and the results showed that gene expression was higher in the S-group compared to the C-group. Overall, these results demonstrated that low salinity stress significantly affected the gene expression patterns of gills.

To determine the transcriptome dynamics in gill under low-salinity stress, a total of 581 differential expressed genes (DEGs) were identified based on a ∣log₂(Fold Change)∣> 1.5 and q-value < 0.05 threshold. Among the identified DEGs, 330 were down-regulated and 281 were up-regulated in the S-group compared to the C-group ( Figure 3A ). Functional enrichment analysis of the DEGs showed that these DEGs were significantly enriched in GO terms involved in “metabolic process,” “stress response related biological regulation,” among other pathways ( Figure 3B , Table 3 ). KEGG pathway analysis indicated that the up- and down-regulated DEGs were significantly enriched in the “cell adhesion molecules (CAMs),” “cytokine-cytokine receptor interaction,” “cellular process,” and “tight junction” gene pathways, among others ( Figure 3C , Table 4 ).

A total of 10 DEGs in the gill transcriptome were selected for expression validation using qRT-PCR. Our results confirmed that the expression of arhgef37, gnb4, LOC104929349, LOC104928056, and LOC104918888 were significantly decreased, whereas the expression of grtp1, LOC104935555, pdhb, slc2a4, and slc9a3 were significantly increased. These findings were consistent with our RNA-seq results, thus confirming the reliability of the transcriptome sequencing data ( Figure 4 ).

To reveal the effects of low-salinity stress on the gill tissues of large yellow croaker, bisulfite-treated genomic DNA libraries were constructed using samples from the gill tissues of three individuals each from the S-group and the C-group. After quality control, 490.87 million and 472.16 million clean reads were obtained for the C-group and the S-group, respectively, and at least 67.04% of the clean reads were uniquely mapped to the reference genome. Therefore, the sequencing depth and density were deemed sufficient to conduct genome-wide methylation analysis. Furthermore, the bisulfite conversion rate of each sample exceeded 99%, thus ensuring reliable results for the WGBS analysis ( Table 5 ). The genome-wide cytosine methylation levels of three sequence contexts (CG, CHG, and CHH; where “H” represents either A, T, or C) were also calculated. The average methylation levels of CG, CHG, and CHH in the C-group were 75.72 ± 0.76%, 10.20 ± 0.56%, and 10.44 ± 0.61%, respectively, whereas the levels for the S-group were 75.48 ± 1.88%, 10.03 ± 0.07%, and 10.26 ± 0.14%.

To identify the genomic DNA methylation dynamics between the S-group and the C-group, we next analyzed the genome-wide gill methylation. Our findings indicated that exon and intron regions showed higher methylation levels than promoter and downstream regions. Furthermore, the S-group exhibited higher methylation levels than the C-group among different genomic elements ( Figures 5A, B ), and similar CG-type methylation patterns were observed at the chromosome level ( Figure 5C ). Moreover, differential methylation analysis showed that there were fewer hyper-methylated regions than hypo-methylated regions among different gene regions ( Figure 5D ).

According to previous studies, the CG context accounts for the largest proportion of genome-wide DNA methylation (Dodge et al., 2002; Lister et al., 2009; Lister et al., 2013), which was consistent with our findings. Therefore, downstream analyses of the relationship between DNA methylation and gene expression were conducted based on the CG context methylation. WGBS analyses elucidated a total of 16,027 DMRs and 3,858 DMPs (differentially methylated promoters). The GO and KEGG analysis results of the DMR target genes are summarized in Table 6 and Table 7 , respectively. Functional analysis indicated that low salinity induced changes in DNA methylation were primarily involved in ion exchange, energy material transport, energy metabolism, cytoskeleton system, nervous system, signal transduction, and other biological processes in the gills of large yellow croaker.

The integrative analysis of RNA-seq with WGBS indicated that there were 422 and 109 DMGs among the identified DMRs and DMPs, respectively ( Figure 6A ). The expression/methylation correlations of the DMR and DMP target DMGs are illustrated in Figures 6B, C , respectively. Furthermore, the differences in the expression level of the DEGs and the methylation status of the entire DMRs in gill between the S-group and C-group are illustrated in Figure 6D .

Functional analysis indicated that the 422 target DMGs of the DMRs were significantly enriched in 228 sub-categories of three major GO categories, including 128 biological process (BP) sub-categories, 23 cellular component (CC) sub-categories, and 77 molecular function (MF) sub-categories with a p-value < 0.01. The top 30 enriched GO terms were mainly involved in osmotic regulation, energy metabolism, gene transcription, and signal transduction and are summarized in Table 8 . Additionally, the top 20 enriched KEGG pathways were consistent with those elucidated through GO analysis and are summarized in Table 9 .

Previous studies reported that promoter methylation has a negative correlation with gene transcription and expression (Law and Jacobsen, 2010). Therefore, we also analyzed the target DMGs of DMPs. GO analysis indicated that 23 hypo-methylated and up-regulated DMGs were mainly enriched in plasma membrane (GO: 0005886), GTPase activator activity (GO: 0005096), signal transduction (GO: 0007165), and 22 other GO categories with a p-value < 0.01 ( Figure 7A ). Moreover, these DMGs were also enriched in the cAMP signaling pathway and other 6 KEGG pathways (p-value < 0.01) ( Figure 7B ). Additionally, 28 hyper-methylated and down-regulated DMGs were enriched in cell junction (GO: 0030054), plasma membrane (GO: 0005886), actin cytoskeleton organization (GO: 0030036), and 29 other sub-categories (p-value < 0.01) ( Figure 7C ), and were also enriched in axon guidance and 21 other KEGG pathways (p-value < 0.01) ( Figure 7D ). The functional analysis of these 51 DMGs (23 hypo-methylated + 28 hyper-methylated) indicated that DNA methylation might play a critical role in the regulation of ion exchange, energy metabolism, adaptive changes in cytoskeleton and nervous system, and other biological processes in gills during low-salinity adaption.

This work selected 10 DMGs from the integrative analysis to verify the correlation between promoter methylation status and gene expression in gills. The result indicated that the hyper-methylated status in the promoters of arhgef37, gnb4, LOC104929349, LOC104928056, and LOC104918888 significantly decreased the expression of these genes in gill under low salinity stress. In contrast, the hypo-methylation of the promoters in grtp1, LOC104935555, pdhb, slc2a4, and slc9a3 notably increased the expression of these genes ( Figure 8 ).

We selected 5 DMGs (gnb4, LOC104929349, LOC104928056, slc2a4, and slc9a3) for validation of the methylation level. One DMP region of each DMG was selected to comparatively validate the methylation levels in gill tissues among S-group and C-group individuals. As shown in Figure 9 , gnb4 and LOC104928056 promoter regions exhibited higher methylation levels in the S-group than in the C-group; conversely, the LOC104929349, slc2a4, and slc9a3 promoter regions exhibited lower methylation levels in the S-group. These results confirmed that the methylation status in the Mass Array test was consistent with the WGBS data.