The role of salinity on genome‐wide DNA methylation dynamics in European sea bass gills

Epigenetic modifications, like DNA methylation, generate phenotypic diversity in fish and ultimately lead to adaptive evolutionary processes. Euryhaline marine species that migrate between salinity‐contrasted habitats have received little attention regarding the role of salinity on whole‐genome DNA methylation. Investigation of salinity‐induced DNA methylation in fish will help to better understand the potential role of this process in salinity acclimation. Using whole‐genome bisulfite sequencing, we compared DNA methylation patterns in European sea bass (Dicentrarchus labrax) juveniles in seawater and after freshwater transfer. We targeted the gill as a crucial organ involved in plastic responses to environmental changes. To investigate the function of DNA methylation in gills, we performed RNAseq and assessed DNA methylome‐transcriptome correlations. We showed a negative correlation between gene expression levels and DNA methylation levels in promoters, first introns and first exons. A significant effect of salinity on DNA methylation dynamics with an overall DNA hypomethylation in freshwater‐transferred fish compared to seawater controls was demonstrated. This suggests a role of DNA methylation changes in salinity acclimation. Genes involved in key functions as metabolism, ion transport and transepithelial permeability (junctional complexes) were differentially methylated and expressed between salinity conditions. Expression of genes involved in mitochondrial metabolism (tricarboxylic acid cycle) was increased, whereas the expression of DNA methyltransferases 3a was repressed. This study reveals novel links between DNA methylation, mainly in promoters and first exons/introns, and gene expression patterns following salinity change.


| INTRODUC TI ON
Euryhaline organisms, living in transitory habitats such as coastal lagoons or estuaries, are generally characterized by a high phenotypic plasticity. Prolonged or repeated exposure to an environmental stressor leads to changes in phenotypic responses and modification in physiological performances. Salinity strongly fluctuates in transitory habitats. Fish living in these habitats must have efficient strategies and mechanisms that operate at the gene, transcript and protein levels allowing them to respond through acclimation and to eventually adapt.
Epigenetic mechanisms like DNA methylation, histone modifications and non-coding RNA are considered as changes in chromatin without sequence changes and can either be heritable or not.
They can play an essential role in promoting phenotypic variation through the modulation of gene expression patterns within a generation (Bird, 2002), as analysed in this study. In fish, DNA methylation is the most studied epigenetic modification process in which methyl groups are transferred to the cytosines of DNA by specific DNA methyltransferases. This process potentially regulates gene expression without affecting the DNA sequence (Jones, 2012). In vertebrate genomes, DNA is methylated at a high rate, with about 60%-80% of the cytosine-phosphate guanine (CpG) dinucleotides methylated (Feng et al., 2010). However, a small subset of CpGs form clusters termed CpG islands that are often associated with genes and known to cover part of their promoter region and at least a part of one exon (Larsen et al., 1992). CpG islands are unmethylated regions which facilitate active gene transcription. DNA methylation plays a significant role in many biological functions through the regulation of gene expression (Suzuki & Bird, 2008). According to the methylated context considered, however, DNA methylation can have a different role in the regulation of gene expression with either an activation, inhibition, or, will remain without any functional effect (Jones, 2012). DNA methylation at the promoter level has often been associated with gene silencing in vertebrates (Newell-Price et al., 2000). However, recent studies suggest that DNA methylation dynamics and its regulatory role in gene expression is much more complex and depends on the cell type and genomic context (Smith et al., 2020). Promoter DNA hypermethylation has for example also been associated with high transcriptional activity by several authors (De Larco et al., 2003;Smith et al., 2020). The role of DNA methylation in gene bodies was less investigated. It could be involved in transcription elongation, alternative splicing or controlling alternative promoter usage (Jones, 2012;Maunakea et al., 2010;Suzuki & Bird, 2008). First, introns of human genes are considered as enriched in CpG islands and are thus likely involved in transcriptional regulation . In mammals, Brenet et al. (2011) have established a negative correlation between DNA methylation and gene expression in the first exon, and this correlation was stronger than between promoter DNA methylation and gene expression. Similar investigations have later been done in fish (Anastasiadi et al., 2018;Liu et al., 2022). The functional role of DNA methylation in different genomic contexts is worth considering and requires further investigations in different species, tissues and environmental conditions.
Studies in marine euryhaline fish species at the genome-wide scale with single basepair resolution methods are still lacking in order to investigate the role of DNA methylation dynamics in salinity acclimation (Metzger & Schulte, 2016).
The European sea bass Dicentrarchus labrax is a main aquaculture species in the Mediterranean area, and has recently become an important model for genetic and epigenetic studies (Vandeputte et al., 2019). In this species, environmentally-induced DNA methylation has been investigated in response to temperature (Anastasiadi et al., 2017;Navarro-Martín et al., 2011), high-density stress (Krick et al., 2021) and farming (Anastasiadi & Piferrer, 2019). D. labrax lives in coastal waters and enters estuaries and coastal lagoons that serve as feeding grounds (Dufour et al., 2009;Pickett et al., 2004).
D. labrax have also been observed in freshwater streams that are connected to the coastal lagoons or estuaries. These transitory habitats are characterized by unpredictable salinity fluctuations, with salinities ranging from 0 to over 60 ppt in Mediterranean lagoons.
The influence of salinity on D. labrax DNA methylation dynamics is still unknown. Also, the functional role of DNA methylation changes on the transcription of genes and modulation of stress response remains poorly investigated.
In this study, we provide a high-resolution analysis of DNA methylation in European sea bass using whole-genome bisulfite sequencing (WGBS) in order to address the question if a 2-week freshwater transfer affects DNA methylation patterns in the gill tissue of D. labrax juveniles. DNA methylation was analysed in different genomic regions (promoters vs gene bodies). RNAseq was performed to explore differentially expressed genes following freshwater exposure. To determine if DNA methylation has a functional role in salinity acclimation, we investigated the correlations between gene expression and DNA methylation levels. Salinity-responsive genes identified by RNAseq exhibiting differential DNA methylation patterns were highlighted in order to identify genes or gene families whose expression could be modulated by DNA methylation dynamics.

| Sea bass husbandry and treatment
European sea bass (N = 350) hatched at Ifremer Station at Palavasles-flots (Hérault, France) at 25 ‰ and were kept in this salinity until 40 days post-hatching in order to optimize survival. After 40 days, they were reared in recirculating seawater (SW; osmolality: 1208 mOsm kg −1 , Na + : 515 mmol L −1 , Cl − : 737 mmol L −1 ) under a 12/12 h light/dark photoperiod at 20°C. At the age of 8 months (13.59 ± 0.12 cm, 32.19 ± 2.62 g), some fish were transferred to brackish water (BW; osmolality: 475 mOsm kg −1 ) for 24 h and then transferred to dechlorinated tap fresh water (FW; osmolality: 8 mOsm kg −1 , Na + : 2 mmol L −1 , Cl − : 3.5 mmol L −1 ) for 2 weeks, as in L' Honoré et al. (2019). In previous studies, juveniles at this age have been shown to be able to tolerate FW for several months (L'Honoré et al., 2019). Fish were raised in groups of 50 (in seawater) or 100 (in FW) individuals in order to form fish shoals, as European sea bass get stressed when raised individually or in too small groups. The size of the tanks was adjusted to keep the same density of individuals in each tank. Pellet food (Le Gouessant, France) was proposed daily to fish. After the 2-week salinity challenge, five fish per salinity condition were killed using a 100 ppm lethal dose of benzocaine before immediate tissue sampling.

| Tissue sampling
The four left gill arches were sampled in five fish/salinity for epigenetic analysis and stored at −80°C. In the same fish, the four right gill arches were sampled for transcriptomic analysis.

| Genomic DNA extraction
Genomic DNA (gDNA) was isolated from a pool of the four sampled gill arches of each juvenile using a Maxwell Blood®16 DNA purification kit (Promega) with the following modifications: before adding proteinase K, gill arches were homogenized with 300 μL of lysis buffer and crushed for 30 s at room temperature with a laboratory mill and stainless steel beads (Mixer Mill MM 400,Retsch,Germany) and then briefly centrifuged for 10 s. Then the manufacturer's instructions were followed for further gDNA extraction. Purity and concentration of gDNA were estimated using the Qubit™ dsDNA BR Assay Kit and the Qubit™ Fluorimeter (Thermo Scientific, Waltham, MA, USA).

| Whole-genome bisulfite sequencing (WGBS)
Bisulfite-conversion, library construction and sequencing (Pairedend reads -150 bp) were performed by Génome Québec (Montréal, Canada) on an Illumina HiSeq X (Illumina, San Diego, CA, USA). Read quality analysis and alignments were performed on the local Galaxy platform (http://bioin fo.univ-perp.fr). The quality of the reads was checked using FastQC (version 0.11.8, Andrews, 2010), and MultiQC (version 1.9, Ewels et al., 2016) was used to aggregate fastQC results into a single report. Phred scores were higher than 25 for more than 90% of the reads' length for all the sequences. FastQ reads were trimmed by Trim Galore (version 0.6.3, Krueger, 2012) based on Phred score below 20 and the adaptors sequences were removed from reads. Trimmed reads were mapped with Bismark (version 0.22.1, Krueger & Andrews, 2011) and aligned to the reference genome of European sea bass (Tine et al., 2014), providing global methylation percentages per genomic context (CpG, ChG and CHH).
Methylation was called using the Bismark Methylation Extractor tool, which generated BED files.

| Genome-wide methylome analysis
After visual inspection of the quality of DNA methylation profiles using the Integrative Genomics Viewer (Thorvaldsdóttir et al., 2013), we performed the analysis of the methylome at the genome's scale. In order to evaluate the global methylation level, metagene analysis was performed using deepTools (version 2.0, Ramírez et al., 2016) command "computeMatrix" to generate read abundance from all samples over genomic regions: promoters, 5′UTR exons, coding exons, first introns, internal introns (located in non-flanking regions of genes), last introns, 3′UTR exons and Transcription End Sites (TES). This matrix was then used to create, using deepTools command "plotProfile", a metagene profile from 2 kb upstream of the Transcription Start Site (TSS) to 2 kb downstream of the TES. The same method was used to generate a profile plot of the level of methylation across all genomic regions. Mean and median DNA methylation values for each condition were compared using a non-parametric Wilcoxon test.
The methylation profiles of the samples were studied using the R package MethylKit (Akalin et al., 2012). The alignment BED files were first converted into a tabular file suitable for the MethylKit package using the methylextract2methylkit tool (version 0.1.0). In order to increase the power of the statistical tests, the samples were filtered according to read coverage. Bases that had less than 10× coverage and those that had greater than 99.9th percentile coverage in each sample were filtered out from the analysis to account for potential PCR bias. The read coverage per base was checked for each sample in order to evaluate if samples suffer from PCR bias, which seemed not to be the case. We therefore chose to not perform a PCR duplication removal step.
Using the Methylkit R package, hierarchical clustering analysis and principal component analysis were performed using the "ClusterSamples" and "PCASamples" functions, respectively. These analyses were based on similarities in the methylation patterns of the samples from each salinity condition. A distance correlation matrix was generated with the Pearson method, and the clustering was performed using the Ward method.

| Total RNA extraction
Total RNA was extracted from the pool of four gill arches sampled in each individual. Before RNA extraction, the gills were rinsed in 1× PBS to wash out previous storage traces, and a Nucleospin® RNA kit was used to extract total RNA following the manufacturer's instructions (Machery-Nagel, Germany). Purity and concentration of the total RNA were checked using the NanoDrop™ One Spectrophotometer (Thermo Scientific, Waltham, MA, USA). The integrity of RNA was assessed using the Agilent 2100 Bioanalyzer (Thermo Scientific, Waltham, MA, USA). All samples passed the quality control threshold with RNA integrity numbers >7.

| Transcriptome analysis
Library construction and sequencing were performed in Perpignan University (Bio-environment platform, Perpignan, France) using one Illumina NextSeq (Illumina, San Diego, CA, USA) high-output flow cell.
Reads were sequenced in 75 base pairs; single-end reads resulted in 25-30 million reads per sample. Read quality analysis and alignments were performed on the local Galaxy platform (http://bioin fo.univperp.fr). Raw read quality was assessed with FastQC (Andrews, 2010) and MultiQC was used to summarize FASTQC results. Phred scores were higher than 30 for more than 90% of the read's length for all the sequences. Reads were cleaned, and adaptors were removed with Cutadapt (version 1.16, Martin, 2011). Cleaned reads were mapped against the genome of European sea bass (Tine et al., 2014) using STAR (version 2.7.8a, Dobin & Gingeras, 2016). FeatureCounts (Liao et al., 2014) was used to count read-mapped per transcript.

| Correlation analysis of genome-wide DNA methylation and gene expression
In order to reveal the complex regulation of gene expression by DNA methylation, violin plots were generated with the ggplot2 R package (version 3.3.6). Correlations between DNA methylation and gene expression levels were measured in freshwater and seawater samples using Spearman's rank correlation coefficient. Gene expression levels were based on counts, normalized according to the DESeq2 method (Anders & Huber, 2010), and then divided into deciles of an equal number of genes ranging from the least (decile 1) to the most expressed (decile 10) genes. DNA methylation levels were based on the beta-value indicated in the BED files. DNA methylation levels in specific genomic regions were represented relative to the gene expression levels.

| Differential expression analysis
Estimated read counts from "featureCounts" were used as input to functions in the "DESeq2" R package (Love et al., 2014) to generate log 2 differential expression fold-difference estimates. Transcripts with less than 10 reads summed across all samples were removed from the analysis. Genes with adjusted p-value less than 5% (according to the FDR method from Benjamini & Hochberg, 1995) were declared differentially expressed. A log 2 FoldChange value greater than zero indicates an upregulation in FW compared to SW, and a log 2 FoldChange value less than 0 indicates a downregulation in FW compared to SW.
In order to explore the variability within our experiment, hierarchical clustering and principal component analysis (PCA) were performed after Variance Stabilizing Transformation (VST) of the count data.

| Gene ontology (GO) enrichment analysis
We did a gene ontology (GO) enrichment analysis in order to identify functions and processes that are overrepresented in subsets of genes with differential expression or differential methylation between salinities. A R package "" specific to sea bass annotation was generated using the "MakeOrgPackage" function from the "AnnotationForge" Package v1.40.1. A total of 23,447 genes have one or more GO annotations. This represents 16,368 GO terms for the biological process category, 20,503 GO terms for the molecular function category and 15,258 GO terms for the cellular component category. The R package "ClusterProfiler" version 4.6.0  with the "org.
Dlabr ax.eg.db" was used to analyse functional profiles of differentially methylated and expressed genes and to identify major biological functions. The functional analysis was performed using as input the DMGs, the DMPs and the differentially expressed genes (DEGs).
A hypergeometric test was performed, and enrichment p-value of gene ontology was calculated to find significantly enriched GO terms in the input list of DMGs, DMPs and DEGs. A p-value of <.05 was set as the threshold value. DMGs, DMPs and DEGs were categorized into three categories belonging to the main GO ontologies: biological process (BP), molecular function (MF) and cellular component (CC). In order to study the differences in functional annotation according to DNA methylation, we separated the DMGs and the DMPs into four groups: hypermethylated gene bodies in FW-vs SW-acclimated fish (GB hyper), hypomethylated gene bodies (GB hypo), hypermethylated promoters (PR hyper) and hypomethylated promoters (PR hypo). We proceeded in the same way for the transcriptome by indicating genes overexpressed in FW (up FW) or in SW (up SW). The 5 GO terms of each analysis with the most significant p-values have been represented. We identified the GO terms from the methylome and transcriptome enrichment analyses. Enriched GO terms between DEGs and DMRs (in genes and promoters) were extracted and used to generate three Venn diagrams, one per GO-term category: BP, MF and CC. This enabled us to identify enriched GO terms that are common to both analyses (methylome and transcriptome).

| Functional analysis using a target gene approach
In order to highlight the number of genes showing changes in methylation and expression patterns between salinity conditions, we generated a Venn diagram with the number of DEGs, DMRs (in genes and promoters) and the overlap between the two. We also identified selected functions (ion transport, integrin signalling, hormone receptor, etc.) and generated a table as well as a scatterplot displaying the expression (log 2 FC) and methylation (beta-value) changes of genes involved in these functions in order to compare their distributions. A conversion table of gene IDs between sea bass genome versions 2014 and 2021 was generated for all these genes.

| Gene expression changes in different methylation and genomic contexts
In order to determine the potential role of DNA methylation changes on gene expression patterns in different genomic contexts, the proportion (percentage) of genes according to their methylation and expression pattern was calculated and represented by horizontal bar plots.

| KEGG pathway analysis
The "enrichKEGG" function of the R package "ClusterProfiler" version 4.6.0 was applied to analyse biological pathways enriched in these four datasets (hyper-downregulated, hyper-upregulated, hypodownregulated and hypo-upregulated genes) with a p-value < .05.
ID conversion from the Hgnc symbol to Entrez was carried out by the annotation R package ("biomaRt", org.Hs.eg.db). Only the genes differentially methylated and expressed at the promoter level, 1st exon or 1st intron were kept to carry out this analysis (number of hyper_downregulated genes = 184, hyper_upregulated = 58, hypo_ downregulated = 328 and hypo_upregulated = 322). The "emapplot" function implemented in the "ClusterProfiler" package as well as the "ggplot 2" and "enrichplot" R packages were used to generate an enrichment map. The "emapplot" function organizes the enriched terms into a network with edges connecting sets of overlapping genes.
Overlapping sets of genes cluster together, making it easier to identify the functional module. The 5 KEGG pathways of each dataset with the most significant p-values have been represented. The proportion of clusters in the pie chart was determined by the number of genes.

| DNA methylation patterns in contrasted salinity treatments
Sequencing of the samples yielded 1,205,785,497 bp Illumina pairedend reads, with an average of 120,578,550 paired-end reads by library. After quality control and alignment, an average of 90,866,384 (75.4%) reads per library was uniquely mapped (Table S1), and only the uniquely aligned reads were submitted for subsequent analyses. For each sample, the bisulfite-conversion efficiency was larger than 98.7% (Table S1). For all samples, the majority of methylated cytosines were in the CpG context. The mean and median methylation levels were respectively 67.20% and 82.32% in the SW control group. In FW-exposed D. labrax, the mean and median methylation levels were 66.71% and 82.32%, respectively ( Figure S1). Even if there was no significant difference in CpG methylation levels between FW and SW samples (Wilcoxon test, p-value =.77), there seemed to be a trend of FW samples being less methylated than SW samples, which was confirmed in subsequent analyses (see subsection '3.4 The effect of freshwater transfer on DNA methylation and gene expression dynamics'). DNA methylation levels were the lowest in the promoter region at the TSS and the 5′UTR (untranslated region). Higher DNA methylation levels were observed in the exons and introns. Downstream the TES, a slightly lower DNA methylation was measured than in GB (Figure 1a

| Gene expression patterns in contrasted salinity treatments
RNAseq yielded 282,282,453 bp Illumina single-end reads, with an average of 28,228,245 single-end reads by library. After quality control and alignment, an average of 25,225,360 (89.35%) reads per library was uniquely mapped (Table S2) Figure S2 shows the global gene expression profiles among all individuals using a heatmap.  Figure S3). In the first exons, this negative correlation was even stronger (R = −.42, p-value = 2.2e-16 in FW and R = −.41, p-value = 2.2e-16 in SW) and we observed less genes with intermediate methylation levels. Most genes exhibited strong methylation levels (>60%) and had low gene expression levels (from 1st to 3rd decile). We also observed a negative correlation between gene expression and DNA methylation levels in the first introns (R = −.34, p-value = 2.2e-16 in FW and R = −.31, p-value = 2.2e-16 in SW) with a similar pattern than in promoters. Interestingly, lower methylation levels (<40%) were observed for highly expressed genes (from 7th to 10th decile). In contrast, when considering the rest of exons and introns, no significant correlation was found between DNA methylation and gene expression levels in FW (R = .01, p-value = .127). DNA methylation levels were close to 75%. In seawater-acclimated-fish, a weak positive correlation was found (R = .02, p-value = .0053) ( Figure S3).

| The effect of freshwater transfer on DNA methylation and gene expression dynamics
We identified 17,265 DMRs in D. labrax exposed to FW versus SW conditions (Table S3)

F I G U R E 2
Cluster analysis based on CpG methylation profiles in five freshwater-exposed (FW, black discs) and five seawater-exposed D. labrax (SW, white discs). (a) Hierarchical clustering of samples was performed using Ward's method based on Pearson's correlation distance for cytosine CpG methylation. (b) Principal Component Analysis (PCA) showing (PC1 × PC2) coordination plane.

F I G U R E 3
Cluster analysis based on gene expression levels in five freshwater-exposed (FW, black discs) and five seawater-exposed D. labrax (SW, white discs).  (Table S3).  Table S6.
We identified several GO terms that were associated to actin cytoskeleton organization and regulation, as well as GO-terms associated to bicellular tight junctions (black stars, Figure 6a). We also identified GO categories linked to mitochondrial metabolism as H + -transporting ATP synthase activity and oxidoreductase activity as well as proton channel activity (black stars, Figure 6a). The most significant GO terms of the DNA methylome analysis are listed in Table S4.
A differential analysis was performed to identify gene expression patterns between FW and SW exposed fish. We found 5144 DEGs, among which 2316 were upregulated and 2828 were F I G U R E 4 DNA methylation levels (in %) in different genomic context for different gene expression levels in freshwater (FW) -acclimated fish. Violin plots showing DNA methylation in promoter (n = 20,581 genes), first exon (n = 22,438 genes), first intron (n = 19,508) and the rest of exons and introns (n = 20,778) at different gene expression levels (divided into deciles based on increasing ranking of gene expression measured as log 2 -transformed normalized counts from DEseq2). Central red dots represent the median of the distribution. Correlations between DNA methylation and gene expression were measured using Spearman's rank correlation coefficient, and the significance level (p) is indicated for each plot. downregulated (Table S5). The enrichment analysis of GO terms ( Figure 6b, Tables S4 and S6) identified functional categories involving genes that were upregulated in FW linked to mitochondrial metabolism. These include genes involved in the tricarboxylic acid (TCA) cycle (as dihydrolipoamide S-succinyltransferase dlst, succinate dehydrogenases sdha, sdhd, sdhb and mitochondrial citrate synthase precursor cs), genes linked to the 'electron transport chain', 'electron transfer activity' and 'ATP biosynthetic process'. Genes linked to methionine metabolism as s-methyl-5-thioadenosine phosphorylase (mtap) were upregulated. We also identified upregulated genes involved in other metabolic processes (peptide metabolic process, cellular amino acid metabolic process) as well as ATPase and proton transporter activity (see next paragraph on ATPases).
In contrast, genes with molecular functions linked to membrane rafts and membrane microdomains were downregulated in FW ( Figure 6b). The most significant GO terms of the transcriptome analysis are listed in Table S4.

| Changes in DNA methylation and gene expression profiles are identified in specific gene pathways
We identified a number of 7120 genes and promoters with DMRs, and among them, 2035 (28.5%) also showed a differential expression between salinity conditions (Figure 7a, Table S7). This corresponds to 39.5% of all genes showing differential expression patterns between salinity conditions (5144 genes) (Figure 7a, Table S7). Enrichment analyses of gene ontology terms (GO terms) were done in both analyses (transcriptome and methylome) to F I G U R E 5 Genome-wide profile of CpG methylation. (a) Chromosomal distribution of hyper-and hypo-methylated regions in freshwatercompared to seawater-exposed D. labrax. p-value < .05 and methylation difference ≥ 25%. (b1) Distribution of the differentially methylated regions (DMRs) in different genomic contexts. (b2) Distribution of genomic features in the sea bass genome, according to Tine et al., 2014. (c) and ( Among these GO terms we identified main categories involved in active ion transport, as 'ATPase-coupled transmembrane transporter activity' (Table S8). Several genes encoding for ATPases such as different subunits of Na + /K + -ATPase (atp1a1, atp1a3b) and the V-H + -ATPase (atp6v1b2, atp6v1d) were upregulated in FW and showed differential methylation patterns between salinities (Table 1). We also identified genes encoding for other main ion transporters or channels that showed differential methylation and were upregulated (as one clcn2 paralog, clcn3, several K + channels: kcnh5, kcnk5, kcnt1) or downregulated in FW (scn4a Na + channel, atp6v1e1b subunit of V-H + -ATPase, chloride channels cftr and one of the clcn2 paralog, potassium channel kcnma1 and an ammonium channel rhcg). More generally, ion transporters seem to appear in different categories (up-or downregulated, hypo-or hypermethylated) with no specific trend ( Figure 8, Table 1). The gene encoding for Aquaporin 3, a glyceroporin, was strongly upregulated in FW (log 2 FC = 3.3) and hypomethylated at the promoter and GB (indicated in Figure 8 as 'water channel'). We also identified in the category 'cytokine receptor activity' the gene encoding for prolactin receptor (prlr), being upregulated in FW and hypomethylated F I G U R E 6 Enrichment analysis of genetic ontology (GO) terms. The five most significant GO terms for each group are represented. (a) DNA methylome. (b) RNA transcriptome. FW, fresh water; GB, gene body; PR, promoter; SW, seawater. The size of the dots represent the ratio of genes associated with the GO term and the colours indicate the p-values (threshold set at 0.05). Discussed GO terms are indicated by stars. Number of genes with GO annotation are listed in Table S6. at the promoter and GB (Table 1, Figure 8 as 'hormone receptor'). Interestingly, two genes encoding for DNA methyltransferase 3a (dnmt3a) were both downregulated in FW and showed differential methylation patterns (Table 1, Figure 8). One gene was hypomethylated, and the other was hypermethylated in GB. Methylcytosine dioxygenase 3 (tet3), a DNA demethylase, was downregulated in FW.

| Freshwater-triggered gene expression changes associated with hyper-and hypomethylation
The distribution of differentially methylated and expressed genes was examined in different genomic contexts (Figure 9). There was a balanced distribution between upregulated and downregulated genes in FW, ranging from 46% to 58% for upregulation and 42%-54% for downregulation when considering hypomethylated genes (exons and introns) and promoters (Figure 9a). In contrast, in the context of hypermethylation, 62%-81% of the genes and promoters were downregulated, and only 19%-38% were upregulated ( Figure 9b). It is worth noting that in the first introns, we found the highest proportion (81%) of downregulated genes in the context of hypermethylation.
Differential methylation was considered at promoters, first exons or first introns only. Some pathways involved in 'regulation of actin cytoskeleton', 'focal adhesion', and 'calcium signalling pathway' were statistically enriched for both DEGs and DMGs (Table S9). Genes encoding for functions related to these pathways had their expression repressed and displayed either a hypermethylation or a hypomethylation in FW ( Figure 10). Regarding the 'calcium signalling pathway' category, we identified genes encoding for the SERCA (Sarco Endoplasmic Reticulum Calcium ATPase) pump (atp2a2) and several voltage-dependent calcium channels. In several KEGG categories, we identified tropomyosin (tpm1) as being hypermethylated and downregulated, as well as several genes involved in 'integrin signalling' (itga 1, 4, 5, 9) (Figure 8).
The 'focal adhesion pathway' was enriched in downregulated genes that were either hyper-or hypomethylated. Several genes involved in integrin signalling with some of them being also hypomethylated (itga 3, 6, 7) and a transcription factor, focal adhesion kinase (ptk2), that was downregulated and hypomethylated. The pathway 'regulation of actin cytoskeleton' was significantly enriched with mainly downregulated genes that were either hyper-or hypomethylated in FW, as indicated previously. Regarding metabolism, the KEGG pathways 'glycosphingolipid biosynthesis' and 'glycerophospholipid metabolism' and 'sphingolipid signalling pathway' were also enriched (Table S9, Figure 10).
The 'tight junction' category was enriched, as in the previous analysis, and showed genes that were mainly hypomethylated and either upregulated or downregulated in FW vs SW (Figures 8, 10).
Ten genes encoding for Claudins were differently expressed and were enriched in DMRs (Table 1). Among them, eight were upregulated in FW and hypomethylated, mainly at the promoter level. In this category, we also identified one paralog of tight-junction protein 2 (tjp2), the cytoskeleton-associated cingulin-like protein (cgn), MarvelD3, a transmembrane component of tight junctions, and the junctional adhesion molecule a-like (f11r). These genes were all hypomethylated and upregulated.
In the 'mineral absorption' pathway, only two genes were highlighted, including one gene encoding for chloride channel 2 (clcn2, one of the two paralogs) and one encoding for the copper transporter ctr1 (slc31a1). Both genes appeared hypermethylated and upregulated in FW. We also noticed the upregulation of the prolactin signalling pathway with genes enriched for hypomethylation, as the prolactin receptor (prlr) ( Table S9). TA B L E 1 List of selected differentially expressed and methylated genes.  Table 1 showing differential expression and methylation in FW versus SW.

| Gene expression and DNA methylation are inversely correlated in promoters, first exons and first introns
First, we determined if there is a correlation between gene expression and DNA methylation levels in different genomic regions using gill homogenates from the same individuals for both analyses. In D. labrax, we observed a global DNA methylation all over the genome, with a depletion at the promoters and transcription start sites. This is consistent with previous observations in vertebrates (Aliaga et al., 2019;Suzuki & Bird, 2008). An inverse correlation between gene expression and DNA methylation levels was found in promoters, first introns and first exons at both salinities, with very similar correlation coefficients when considering both salinity treatments separately. This inverse correlation has previously been highlighted by Anastasiadi et al. (2018), who compared DNA methylation levels in muscles and testes of D. labrax, and later by Liu et al. (2022) in Anguilla anguilla muscles for first exons. The presence of unmethylated CpG islands in the proximity of the TSS is consistent with the role of first introns and first exons in the regulation of gene expression. Additionally to promoters, first introns and first exons could thus be involved in transcriptional regulation, as highlighted for first introns in humans . We showed that the rest of the introns and exons were highly methylated (>75%) and did not show any correlation with gene expression, as previously observed in previous studies (Anastasiadi et al., 2018).

| Freshwater transfer overall induces DNA hypomethylation in D. labrax gills
Few teleost species have been analysed regarding the effect of salinity on DNA methylation dynamics (Artemov et al., 2017;Heckwolf et al., 2020;Metzger & Schulte, 2018;Yang et al., 2023;Zhang et al., 2022). Our data showed that a freshwater transfer induced in fish gills showed an overall hypomethylation of DNA across all chromosomes. Parallelly, a downregulation of the expression of two paralogs of DNA methyltransferase 3a was observed, which is consistent with this overall decrease in DNA methylation. Predominant were detected. A salinity effect on DNA methylation level was also observed in threespine stickleback (G. aculeatus) of a marine ecotype, where overall hypomethylation was observed in whole fish that were reared for 1 month at high salinity relative to those reared at low salinity (Metzger & Schulte, 2018), which is an inverse tendency than what we observed in D. labrax gills. It therefore seems that a medium-term (2 weeks in this study) as well as a long-term salinity acclimation (1 month) can affect DNA methylation dynamics in fish tissues. It has to be noted that in the study on stickleback, salinity transfer was done from the stage of fertilization, whereas in sea bass, we transferred fish as juveniles. DNA methylation is very dynamic during development (Wang & Bhandari, 2020), and an early salinity transfer in development might affect fish differently regarding DNA methylation than a salinity transfer as juveniles. In sea bass, the distribution of DMRs across genomic features (e.g., promoters, exons, introns and intergenic regions) did not differ from the relative proportions of these features in the genome (Figure 5b2) (Tine et al., 2014) which means that FW induces methylation changes in all genomic features without targeting specific regions in the genome. Global hypomethylation has previously been shown following several stressors such as metal stress in zebrafish embryos (Bian & Gao, 2021), and salinity stress in the crustacean Daphnia magna (Jeremias et al., 2018). Global hypomethylation could be considered as a global response to a stressor, potentially playing a role in the modulation of transcription activity.

| Hypermethylation triggered by freshwater transfer overall induces downregulation of genes
Despite an overall DNA hypomethylation in the sea bass genome, we observed an interesting trend for hypermethylated genes to be silenced at a very high percentage, notably in the first introns. When considering genes that were hypermethylated in first introns and showed differential expression between salinities, 81% (n = 174) of the genes were downregulated, and only 19% (n = 42) were upregulated (Figure 9b), which is intriguing. Genes that were downregulated and hypermethylated in promoters, first exons or first introns belonged to different KEGG categories (Table S9). Among those, we identified several genes involved in cell-cell adhesion and regulation of actin cytoskeleton as tropomyosin and a change in integrin turnover. It is known that integrins assembling to the cytoskeleton is important in cell-cell adhesion (Delon & Brown, 2007). There is in fact an extensive remodelling of gills in euryhaline teleosts following salinity transfer, that involves cell proliferation and turnover pathways, leading to epithelial remodelling. In D. labrax as in other species, the density and subtypes of gill ionocytes is increased in the freshwater environment (Masroor et al., 2018). We identified several genes involved in the calcium signalling pathway. Calcium is a ubiquitous second-messenger regulating numerous cellular processes, including proliferation, cellular metabolism and cell death. In mammal studies linked to cancer research, the SERCA pump (encoded by atp2a), which sequesters Ca 2+ into the endoplasmic reticulum, as well as other calcium-regulated genes, showed altered expression patterns that coincided with increased promoter methylation (Bertocci et al., 2022). In D. labrax, the hypermethylation of the gene encoding for SERCA might contribute to its significantly lower expression in FW (log 2 FC: −0.38), which could lead to increased cytoplasmic calcium availability. We did not measure any expression difference in plasma membrane Ca 2+ ATPase (pmca) between salinity conditions, which is an important pump expressed in gill ionocytes to take up Ca 2+ from FW environments.

| Freshwater transfer affects DNA methylation in crucial genes involved in maintaining hydromineral balance
We showed that a freshwater transfer affects expression changes in genes that were differentially methylated, which suggests that a change in salinity could induce an altered pattern of DNA methylation, which in turn could have functional consequences and allow sea bass to display phenotypic variation through gene expression changes linked to hydromineral balance.
Gene expression changes in fish gills following a salinity change have been shown in numerous species (Leguen et al., 2015;Qin et al., 2022), whereas studies on DNA methylation changes in different salinity conditions are scarce. We focused in this study on genes that are involved in osmoregulation, cell volume regulation and acid-base regulation and that display gene expression as well as DNA methylation changes (Table 1). Among them, the gene encoding for water channel Aquaporin 3 (aqp3) is highly induced in FW versus SW (log 2 FC = 3.3), and hypomethylated at promoter and in the gene body. Aqp3 is expressed in ionocytes and is known to be overexpressed in gills of numerous species in FW, notably for the basolateral release of water from ionocytes to the serosal fluid to prevent cell swelling (Cutler & Cramb, 2002;Giffard-Mena et al., 2007).
In studies on mammal gastric carcinoma, aqp3 was shown to be hypermethylated at its promoter and first exon, which limited its expression (Wang et al., 2019). In fish, aqp3 gene expression is controlled by cortisol and prolactin, where prolactin induces its expression and cortisol decreases its expression in gills of Mozambique tilapia Oreochromis mossambicus (Breves et al., 2016). Interestingly, in D. labrax, the gene nr3c1 encoding for the glucocorticoid receptor was repressed and hypomethylated at its gene body (GB) ( Table 1).
Moreover, the prlra gene encoding for one of the two prolactin re- However, studies on KO mice brains with selective disruption of the dopamine D2 receptor in neurons, have shown an upregulation of prlr correlated with decreased methylation of their promoters (Brie et al., 2020).
Regarding ion channels and transporters that are involved in salt uptake, we identified chloride channels clcn2 and clcn3 that were upregulated in FW vs SW (with log 2 FC of 2.63 and 0.71 respectively). clcn2 was hypermethylated at GB level. Interestingly we found another paralog of clcn2, that was downregulated in FW (log 2 FC = −2.13) and hypomethylated in GB. Both Clcn2 and Clcn3 channels have been localized in basolateral membranes of ionocytes and are suspected to transport Cl − to the blood for its uptake (Bossus et al., 2013;Tang et al., 2010). clcn2 expression and protein changes according to salinity have been shown in several studies, although with sometimes contrasting data (Bossus et al., 2013;Root et al., 2021) that might be linked to the presence of two clcn2 paralogs that were not differentiated. Several genes encode for Na + / K + -ATPase, which is a key active ion transporting pump expressed in basolateral membranes of gill ionocytes. As expected, the major paralog (atp1a1a) encoding for subunit NKAα1 was upregulated in FW vs SW as shown previously by Blondeau-Bidet et al. (2019) in D. labrax gills. We also observed a hypomethylation at promoter and GB levels. atp1a3 (encoding for the NKAα3 subunit), which has not been investigated so far in D. labrax gills, was also upregulated and hypomethylated in promoters. In D. labrax SW-type ionocytes, the apical chloride channel CFTR and basolateral cotransporter NKCC1 are crucial proteins involved in salt secretion (Bodinier et al., 2009;Lorin-Nebel et al., 2006). Both genes encoding for these proteins were downregulated in FW and hypermethylated at GB level. (for acid-base balance) coupled to transepithelial Na + uptake (for osmoregulation). We identified three genes encoding for VHA that showed expression and methylation changes upon FW transfer (Table 1). Together, all these data point to significant methylation changes in key genes involved in hyper-and hypo-osmoregulation as well as acid-base regulation. This is consistent with the statement that salinity affects the plastic responses through DNA methylation changes, to maintain hydromineral balance, as already mentioned in other species (Heckwolf et al., 2020).

| Paracellular permeability-related genes are hypomethylated and upregulated
In clinical studies on mammal cancer where disruption of cell-cell junctions are frequently found, investigations have shown a link between claudin overexpression and DNA hypomethylation for several claudins (Li et al., 2018), including cldn4 (Kwon et al., 2011).
In fish gills, the tightening of the gill epithelium is a natural and essential process in FW in order to avoid excess water entry and limit passive ion loss. Euryhaline species have to constantly adjust gill permeability through the expression of tight or leaky junctions. A striking result of this study is the consistent promoter hypomethylation and upregulation of genes involved in paracellular permeability and the formation of tight junctions. In both analyses (KEGG pathway and GO-term enrichment analyses), tight junctions appeared among the most significant pathways. In sticklebacks, Metzger and Schulte (2018) (Tipsmark et al., 2008), killifish (Whitehead et al., 2012), and rainbow trout (Oncorhynchus mykiss) (Leguen et al., 2015). We identified two genes encoding for D. labrax Claudin 4 paralogs that were both upregulated in FW (log 2 FC: 0.87 and 1.47) and hypomethylated in promoters as well as other key genes involved in tight junction assembling (cgn, f11r, marveld3). Given the enrichment of the tight junction pathway for upregulation and hypomethylation, there is some evidence that claudin promoters, first exons or introns are a preferential target for differential methylation in changing salinity environments in euryhaline D. labrax. Studies in other teleost species on the effect of methylation on the tight junction pathway are required to confirm this trend in fish.

| Metabolism and salinity acclimation
Our analyses showed changes in the expression and methylation of genes involved in lipid, protein and carbohydrate metabolism as well as mitochondrial functions. The transition from a seawater-to a freshwater physiology requires gill epithelium turnover and active ion transport and is thus an energy-consuming process. This energy is provided by metabolites related to carbohydrates, proteins/amino acids and lipids, transported from liver stores to the blood, towards the gills (Chang et al., 2007). In addition to energy sources provided by the blood, there is a local energy supply at the gill level, mainly for carbohydrates, necessary for the modulation and stimulation of gill epithelium reorganization and ion transport mechanisms (Chang et al., 2022;. Lipid metabolism in fish gills has been widely overlooked, de-  (Lin et al., 2021). Hydrolysis of sphingomyelin due to a cell stress can lead to ceramide formation, which is implicated in numerous physiological functions. Several ion channels such as CFTR (Ramu et al., 2007) and voltage gated K + channels (Fan & Makielski, 1997) have been shown to be inhibited by ceramide. In this study, sphingomyelin synthase (sgms1) was upregulated and hypomethylated in FW.
We observed a hypomethylation and downregulation of ceramide synthase 5 (cers5) and serine palmitoyltransferase (sptlc2), which is involved in the first step of sphingolipid biosynthesis (Hanada, 2003).
Altogether, these results indicate a potential change in the ceramide species profiles (Gault et al., 2010). However, the potential functional link between methylation changes and expression changes of genes involved in sphingolipid metabolism is not clear and needs further investigation.
Studies have pointed out that the enzymes involved in methylation and demethylation have substrates that are responsive to cellular metabolism (Reid et al., 2017). Mitochondria provide key metabolites for epigenetic processes (Shaughnessy et al., 2014). The availability of these metabolites change, when fish are energetically challenged by environmental stressors.
For DNA methyltransferases, substrates and cofactors include methylthioadenosine (MTA), S-adenosylmethionine (SAM) and S-adenosylhomocysteine (SAH) for example whereas for DNA demethylases (TET-family), other metabolites are used as succinate, fumarate, … being intermediates of the TCA cycle (Reid et al., 2017). We observed changes in the expression of genes linked to methionine metabolism (ahcyl1, 2, mtap) and an overall overexpression of genes involved in the TCA cycle (dlst, sdha, sdhd, sdhb, cs). This clearly displays changes in metabolic pathways and energetic status, which could directly affect the epigenome and DNA methylation dynamics. In gills of sea bass transferred to FW, ionocytes, previously called mitochondria-rich cells due to their high abundance of mitochondria, are present at higher densities with higher energy-consuming Na + /K + -ATPase activities (Masroor et al., 2018). We also measured higher expression of paralogs encoding for the V-H + -ATPase indicating a high energy demand for transepithelial ion transport in FW. In accordance with this, we observed higher mitochondrial activities in FW. Mitochondria have a central role in energy (ATP) production but also in metabolite production in the TCA cycle and mitochondria-nuclear signalling. These processes can be linked to epigenetic regulation. As chromatin-modifying enzymes use as substrates and cofactors metabolites derived from diverse metabolic pathways including notably the TCA cycle (Lopes, 2020), salinity-driven changes in transcription of these genes might also directly affect the availability of substrates for chromatin-modifying enzymes and affect DNA methylation dynamics. Relationships between genes, environment and epigenetic marks, and the variation of those marks still require more investigations to gather a full understanding of the determinism of salinity acclimation in fish.

| CON CLUS ION
We have investigated salinity-induced DNA methylation and its role in plasticity and gene expression in gills of euryhaline European sea bass. This study highlighted that genes with low methylation levels in first exons, first introns and promoters are generally highly expressed. We also showed that FW triggers an overall hypomethylation of the genome. Our investigation showed that pathways involved in tight junctions are highly enriched in upregulated genes displaying hypomethylated promoters. We also identified other pathways as lipid metabolism, calcium signalling and regulation of actin cytoskeleton that were enriched for gene expression and DNA methylation changes in either promoters or first exons/introns. Numerous key genes involved in transepithelial ion transport of gill ionocytes also show methylation and gene expression changes. Interestingly, mitochondria metabolism is strongly activated, suggesting a modulation of metabolite availability as substrate for chromatin-modifying enzymes. We recommend further investigation of methylation dynamics in environmentally challenged fish in order to determine the role of methylation changes in phenotypic plasticity, acclimation and adaptation.

ACK N O WLE D G E M ENTS
We thank Erick Desmarais for providing the GO annotations of the

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw sequence reads (SRA) have been deposited at NCBI under BioProject PRJNA952864.