Insights into Solea senegalensis Reproduction through Gonadal Tissue Methylation Analysis and Transcriptomic Integration

1. Introduction

Understanding reproduction is crucial for the commercial production of species. Sex determination and differentiation are vital processes in most complex organisms, with mechanisms that vary widely. Teleost fish, which boast an unparalleled diversity of species, exemplify this variation. They display all known vertebrate sex determination systems and even show differences among closely related species. These mechanisms include genetic sex determination, where sex is determined by inherited genetic factors, and environmental sex determination, where factors such as temperature, pH levels, population density, and social interactions can influence sex determination and differentiation. Additionally, some species exhibit interactions between both mechanisms [1,2,3]. Therefore, identifying and characterizing the genes and pathways involved in reproduction is a fundamental initial step.

Flatfishes are a notable example of this diversity. They undergo a significant transformation from their bilateral larval stage to their adult and characteristic flat shape adapted to a demersal lifestyle. This adaptation has led to rapid diversification and a higher rate of molecular evolution compared to their bilateral counterparts. As a result, flatfish exhibit a wide variety of sex determination mechanisms, involving different master genes and unique sex determination regions [4,5]

Solea senegalensis leads a significant expansion and investment in the European aquaculture industry in recent years. This flatfish species boasts a high growth rate with high market value, and strong demand, whose production is primarily concentrated in Spain, Iceland, and Portugal, with a continuous yearly growth [6,7]

The primary challenges in S. senegalensis aquaculture production, aside from optimizing nutrition, disease control, and addressing morphological abnormalities, include reproductive issues in captivity. Specifically, captive-bred males lack fertilization capabilities. F1 males do not exhibit the complex species-specific reproductive behaviors necessary for successful spawning, such as prespawning chasing or synchronized paired spawning. Additionally, F1 males have reduced fertilization capacity due to low sperm volume and, which hinders the implementation of genetic improvement measures through selective breeding [8,9,10,11].

Gonadal development and gametogenesis research in flatfish has primarily focused on a few species, notably Cynoglossus semilaevis, Scophthalmus maximus, and Paralichthys olivaceus. Numerous studies have aimed to identify and characterize genes homologous to sex-determining pathway genes in model species using mass sequencing technologies [12,13,14]. Genomic research on S. senegalensis includes studies analyzing transcriptomes from different tissues and culture conditions [15,16], the development of an integrated genetic map suggesting the evolution of several chromosomes in S. senegalensis as the outcome to Robertsonian fusions, pericentric inversions, and chromosomal rearrangements [17], and two recent chromosome-level genome assemblies providing evidence of the role of the fshr gene in sex determination [18,19].

Previous studies on S. senegalensis have explored the relationship between temperature and the epigenetic regulation of myogenesis. Given the interest in cultivating larger female populations, it is crucial to consider how temperature affects the sex ratio in Senegalese sole production. When larvae were exposed to daily temperatures of around 22°C and nighttime temperatures of approximately 19°C from 1 to 97 days after hatching, 70% of the population became female, with increased estradiol levels [20].

Methylation analysis of the testicular tissue in this species, focusing on genes related to reproductive traits such as sex determination and spermatogenesis, revealed higher methylation levels in F1 individuals compared to wild ones, and in immature wild fish compared to mature ones. These differences, particularly in genes linked to sexual differentiation, suggest a molecular basis for the behavioral differences observed between F1 and wild individuals, such as the absence of courtship in F1 males. The distinct cultivation conditions, with constant temperatures for F1 individuals and fluctuating temperatures for wild ones, likely contribute to these findings [21].

In this context, the quality of sperm can be significantly influenced by the methylation status of sperm DNA, and changes in DNA methylation have been linked to male infertility [12,22]. In aquaculture, fish phenotypes are shaped by various external and internal factors, such as temperature, nutrition, toxins, salinity, and population density, which often operate through distinct DNA methylation patterns.

Temporary solutions have been developed to enhance reproductive status, including artificial fertilization and nutritional strategies, as well cryopreservation protocols [23,24,25]. While hormonal induction in F1 sole has shown some success [26,27], results are not yet optimal as a solution to these reproductive issues. Despite extensive research, transcriptome data for flatfish species like S. senegalensis remain limited. Many questions about the specific pathways and genes involved in sex determination and differentiation remain unanswered, highlighting the need for further research. Additionally, the reproductive challenges faced by hatchery-born individuals complicate the implementation of precise aquaculture practices for this species.

Epigenetic changes, particularly DNA methylation, are crucial for how organisms adapt to environmental changes, influencing gene expression through the conversion of cytosines to 5′-methylcytosine by DNA methyltransferases (DNMTs). The patterns of DNA methylation are tissue-specific. Although DNA methylation is primarily found in CpG dinucleotides, these modifications can also occur at CpG, CHG, or CHH regions (where H represents A, T, or C). DNMTs facilitate methylation in these areas, while Ten-Eleven Translocation (TET) methylcytosine dioxygenases (TET1, TET2, and TET3) aid in the demethylation of the methylated cytosines [28,29,30]

These modifications influence various genetic processes, including genome stability, genomic imprinting, and transcription. While traditionally linked to gene repression, DNA methylation can also enhance gene expression depending on its location. Decreased expression has been associated with hypermethylation of CpG sites in the promoter and gene body regions. Conversely, increased expression has been linked to hypermethylation of CpG sites downstream of the gene body [31,32].

The impact of epigenetics on fish sex determination has been extensively described. Modifications in DNA methylation have been associated with temperature-induced sex differentiation and environmental sex differentiation regulation. However, while epigenetics holds great promise for improving aquaculture, practical applications remain scarce [33].

The main objective of this study was to analyze the differences in methylation patterns in the gonadal tissue of Senegalese sole (ovaries and testes) and to understand the role of methylation as an epigenetic mechanism in regulating gene expression at two different stages of maturation, under both wild and captive breeding conditions. Our results revealed significant methylation differences between sexes and rearing groups, which were also observed in gene expression comparisons. Additionally, the integration of methylome and transcriptome data showed that methylation acts as both a suppressor and enhancer of gene expression.

These findings underscore the crucial role of methylation in regulating gene expression through maturation and its potential role in sex differentiation, highlighting the complex interplay between epigenetic modifications and gene expression in this species.

2. Materials and Methods

2.1. Gonadal Tissue Sampling

Individuals were collected from two populations: an artificially bred population (F1) at IFAPA and a wild population from the Bay of Cádiz, Andalusia. The samples were categorized based on three criteria: gender, origin, and maturity level.

Sampling occurred during the spawning season (May-June), with immature individuals collected year-round. Gonadal tissue samples were taken from 24 S. senegalensis specimens. The study included eight groups based on the previous criteria: mature wild (MW), immature wild (IW), mature F1 (MF), and immature F1 (IF). We examined 6 adult males (4 wild, 2 F1), 6 adult females (2 wild, 4 F1), 6 immature males (4 wild juveniles, 2 F1 juveniles), and 6 immature females (3 wild juveniles, 3 F1 juveniles).

Size and weight were used as indicators of sexual maturation. Mature specimens were collected during spawning months, and their maturity was confirmed through gonadal dissection; A portion of each sample was used for histological analysis to verify the maturation stages.

F1 specimens reached sexual maturity under stable conditions. Factors such as feeding patterns, population density, and oxygen levels may act as epigenetic modulators. Despite reaching sexual maturity, F1 individuals faced reproductive challenges.

Fish were transferred to the laboratory in oxygenated containers, anesthetized with clove oil (90-100 mg/L), and gonad samples were stored in RNAlater® TissueProtect at -80°C until DNA extraction. The experimental protocol followed EU guidelines (86/609/EU) and was approved by the Ethics Committee of the University of Cadiz, Spain.

2.2. DNA and RNA Extraction

Genomic DNA was extracted from the testes using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer’s instructions. Following extraction, the DNA underwent purification and size selection using the HighPrep™ PCR Clean-up System (MAGBIO). Quality assessment was performed using a NanoDrop™ 2000C spectrophotometer, while concentration was measured using a Qubit 4 fluorometer (ThermoFisher).

RNA extraction was performed from 0.5 mg of gonadal tissue using TRIzol™ Reagent (Thermo Fisher Scientific), following the manufacturer’s protocol. Quantification and purity assessment of the extracted RNA were conducted via fluorometry using the Qubit system (Invitrogen). RNA integrity (RIN) was determined using the Bioanalyzer Agilent 2100, following the manufacturer’s protocol. Samples with an RIN greater than 4 were selected for transcriptome sequencing.

2.3. RRBS and RNA Libraries Preparation and Sequencing

For RRBS library preparation genomic DNA (100 ng per sample) was digested with MspI for 12 hours using the Diagenode Premium Kit. After end preparation, adaptor ligation, and size selection, sample concentrations were quantified by qPCR. Samples with similar concentrations were pooled in groups of six and underwent bisulfite conversion, including both methylated and unmethylated control DNA. RRBS libraries were enriched by PCR, purified with AMPure XP Beads, and assessed for fragment size distribution using an Agilent 2200 Bioanalyzer. Libraries were quantified with a Qubit fluorometer, pooled, and sequenced (150 bp paired end) on the Novaseq Illumina platform by Novogene's genome sequencing service.

For mRNA library preparation mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. After fragmentation, the first strand of cDNA was synthesized with random hexamer primers, followed by the second strand synthesis. Subsequent steps included end repair, A-tailing, adapter ligation, size selection, amplification, and purification. RNA libraries were quantified using Qubit and real-time PCR, and their size distribution was assessed with a bioanalyzer. Samples with an RNA integrity number greater than 6 were selected for transcriptome sequencing. Libraries were then pooled and sequenced on Novogene's Illumina platforms based on effective library concentration and data requirements.

2.4. Quality Control and Data Processing

We used Trim Galore 0.6.61 to eliminate Illumina sequencing adapters and low-quality base calls (Phred quality score < 20), and short reads (<20 bp). To ensure the quality of the fastq files data, we employed FastQC 0.12.0.

For RRBS alignment, we performed an in-silico bisulphite conversion of the S. senegalensis reference genome (GenBank assembly GCA_919967415.2 from NCBI). The alignment of RRBS sequence reads to the converted reference genome, as well as the subsequent methylation profiling and methylation calling of each cytosine, were performed using Bismark 0.22.3 for CpG coverage analysis.

RNA paired-end reads were aligned to the genome using Hisat2 v2.0.5. Feature Counts v1.5.0-p3 was employed to count the number of reads mapping to each gene and to calculate the FPKM (Fragments Per Kilobase Million) of each gene.

2.5. Differential Analysis

For differential methylation analysis, coverage output files were processed according to Mukiibi et al. (2022) using the Bioconductor edgeR 3.42.2 R package [34]. CpG sites that were consistently methylated, unmethylated, or had low coverage (fewer than 8 reads per sample) were excluded. To explore methylation patterns, multi-dimensional scaling (MDS) analysis was performed using the R prcomp function. CpG read counts within sample libraries were normalized to the average library size for consistency. Likelihood ratio tests were used to compare methylation rates between different study groups. Significantly differentially methylated CpG sites were identified using a false discovery rate threshold of 0.01 (Benjamini–Hochberg correction).

Considering the categories based on the criteria used for sample collection (sex, maturity, and origin), 12 comparisons were conducted for both differential methylation analysis and transcriptomic profile comparisons between groups. These comparisons aimed to examine the potential effects of these factors on methylation and gene expression. The following group comparisons were performed: Mature F1 Male Vs. Mature Wild Male (MFM/MWM), Immature F1 Male Vs. Immature Wild Male (IFM/IWM), Immature Wild Male Vs. Mature Wild Male (IWM/MWM), Immature F1 Male Vs. Mature F1 Male (IFM/MFM), Mature F1 Female Vs. Mature Wild Female (MFF/MWF), Immature F1 Female Vs. Immature Wild Female (IFF/IWF), Immature Wild Female Vs. Mature Wild Female (IWF/MWF), Immature F1 Female Vs. Mature F1 Female (IFF/MFF), Mature F1 Female vs Mature F1 Male (MFF/MFM), Immature F1 Female Vs. Immature F1 Male (IFF/IFM), Mature Wild female Vs. Mature Wild Male (MWF/MWM), and Immature Wild F1 Vs. Immature Wild Male (IWF/IWM).

Methylation differences were visualized using volcano plots and heatmaps created with the R package gplots v3.1.1. Functional annotation was done using HOMER software v4.1, indicating CpG site positions relative to gene promoters, exons, introns, or transcription start/termination sites based on the Senegalese sole annotation file (GenBank assembly GCA_919967415.2 from NCBI).

The R package DESeq2 (1.20.0) was employed for differential expression analysis between two conditions/groups. We analyzed the correlation of gene expression levels between samples using the FPKM values of all genes in each sample and calculated the correlation coefficients within and between groups. These coefficients were visualized in heatmaps to display the differences between groups and the consistency within groups. For differential expression analysis with two biological replicates per condition, Principal Component Analysis (PCA), volcano plots, and heatmaps were also generated using the ggplot2 package in R (Version 3.0.3) within Novogene’s NovoMagic Tool (https://eu-magic.novogene.com/).

To integrate methylation and transcription data, a correlation analysis was conducted to identify significant associations (P-value < 0.05) between methylation levels and gene expression. Differentially methylated regions and differentially expressed genes (DEGs) were identified and compared to find overlapping genes. Pearson’s correlations between the percentage of DNA methylation at CpG sites and the counts per million (cpm) of their overlapping genes were computed.

Due to the difficulty in obtaining specimens we were unable to acquire the desired number of samples. Consequently, our sample size for farmed juveniles and adults was limited, potentially reducing certain variability. Specifically, we could not achieve more than two replicates for the F1 male groups due to their limited availability, which may have introduced some bias into the analysis. Furthermore, the mature wild female group lacks a third replicate. Despite these limitations, we believe the experiment remains scientifically valuable for studying the Senegalese sole. However, it is important to consider these constraints in the analysis and interpretation of the results. Further studies with larger sample sizes are necessary to validate and expand upon these findings.

2.6. GO and KEGG Enrichment Analysis

The clusterProfiler R package was used to perform Gene Ontology (GO) and KEGG enrichment analyses on differentially methylated and expressed genes using Novogene’s Novomagic Platform (https://eu-magic.novogene.com/).

GO terms with a corrected P-value below 0.05 were deemed significantly enriched by these genes. Additionally, the clusterProfiler R package was used to test the statistical enrichment of differentially expressed genes in KEGG pathways. Associations with GO and KEGG terms were visualized in bar and scatter plots.

3. Results

3.1. Data Quality Control of Sequencing Information

RRBS libraries from gonadal tissue sequencing yielded a total of 2,180,046,558 raw paired-end reads (n = 12). The Q20 and Q30 values ranged from 93.04 ± 1.17% to 85.46 ± 1.47% (Mean ± S.D), respectively, across samples (Table A1). After filtering, 2,133,878,766 reads remained. The mapping efficiency was 62.33 ± 4.54%. Following quality control, a total of 2,623,861,381 cytosines in the CpG context were reported, showing some level of methylation, which accounted for 53.63 ± 0.20% of the total number of analyzed cytosines (Table A2).

RNA sequencing of gonadal samples yielded a total of 1,239,301,554 clean reads. The Q20 and Q30 values ranged from 97.52 ± 0.28% to 92.99 ± 0.60%, respectively, across samples (Table A1). On average, 88.57 ± 1.00% of the sequences aligned with the S. senegalensis reference genome; the mapping efficiency was 94.60 ± 0.01%. Among the mapped reads, 82.48 ± 4.88% corresponded to exons, 5.94 ± 1.94% to introns, and 11.57 ± 3.20% to intergenic regions. Overall, these results indicate that the sequencing experiment was successful and yielded high-quality data.

3.2. Methylation Profiling and Differential Analysis

A general multi-dimensional scaling plot of overall methylation levels reveals primary differences between males and females. Within each sex, F1 and wild groups can be distinguished, particularly among male samples, with F1 males exhibiting the most significant methylation levels. Independent comparisons of males and females also revealed differences between the F1 and wild groups. Notably, significant differences were observed between groups of different maturity (immature and wild) in both sexes (Figure 1).

Considering the groups of study, we conducted 12 differential methylation profile comparisons (Table A3). Only CpG sites with sufficient coverage (>8 counts) and variation in methylation levels across samples with a false discovery rate below 0.01, adjusted using the Benjamini–Hochberg method, were deemed significantly differentially methylated. The analysis revealed a total of 19,117 CpG sites with significant differential methylation across the performed comparisons (Figure 2), which were considered for the genetic annotation step.

Generally, in terms of differently methylated CpGs (DMCpGs), the comparisons revealed that the largest number of DMCpGs were observed between different maturation states, followed by comparisons between different rearing groups, and lastly between sexes.

The general average methylation difference, measured as log2 fold change (logFC), was 1.50 ± 0.20. While the average fold change values remained close among comparisons (ranging from 1.04 ± 0.22 to 1.75 ± 0.31) the highest differences among the studied groups were observed when comparing groups of different maturity, followed by comparisons between different sex groups, and finally groups differentiated by their rearing.

Focusing on the analysis of male groups, those showing a marked reproductive dysfunction, generally a larger number of up-methylated CpGs were observed towards F1 males when compared to wild-type males, and in immature individuals when compared to their mature counterparts. This trend was also observed in females apart from the comparison between immature and mature F1 female groups, where the mature group showed a larger number of DMCpGs. When comparing groups of males and females, the latter consistently exhibited higher levels of methylation. However, an exception was observed in the immature F1 groups, where males showed a larger number of up-methylated CpGs.

While having samples classified into specific categories, they were grouped to compare the general categories of sex, maturity, and origin. This approach allowed for the identification of methylation differences between these broad categories, considering the total number of individuals. Consequently, comparisons were made between males and females, mature and immature individuals, and F1 and wild individuals. Significant differences were observed in all comparisons, with higher methylation levels in females, immature individuals, and F1 individuals when compared to their respective counterparts (Figure 3).

Despite our efforts to include more than three samples in each studied group, the limited number of F1 samples allowed us to include only two replicates for the immature and mature F1 male groups, which might have caused some bias in the analysis. Additionally, a third replicate in the mature wild female group is missing. However, RRBS offers extensive coverage of CpG-rich regions in the genome, which can often offset the absence of biological replicates [35].

The high-resolution data from RRBS can yield robust and reliable results, even with a limited number of samples. We believe that our study, despite its limitations, makes a valuable contribution to the field and hope that it will pave the way for future research with larger sample sizes.

3.2. Genomic Annotation Based on Differentially Methylated Regions

The reported DMCpG sites were annotated based on their locus and proximity to genes, following the annotation of the Senegalese sole genome (GenBank assembly GCA_919967415.2 from NCBI), and further categorized based on their location regarding the nearest gene, as within gene promoters/transcription start sites (TSS, ±1kbp), exons, introns, intergenic regions, or transcription termination sites (TTS, ±1kbp).

These DMCpGs were annotated to 4,581 different genes. Most DMCpGs were found in intronic and intergenic regions. The distribution was as follows: 8.00% in promoters/transcription start sites (TSS, ±1kbp), 15.15% in exons, 40.45% in introns, 31.53% in intergenic regions, and 4.87% in transcription termination sites (TTS, ±1kbp). Additionally, 13.08% were in first intron regions.

Upon annotation, the identified DMCpG sites were found to be situated near or within various sox family transcription factors, genes related to the steroidogenic pathway (dmrt2ab, hsd3b7), and genes involved in methylation and demethylation processes (tet3, kdm1b, kdm5b, and kdm6ba).

Our findings revealed CpG methylation differences among certain groups associated with genes involved in sex determination and reproduction, gonadal differentiation and sperm production, and germ cell regulation and development (fshr, dmrt, sox9a, bcar1, hsd17b1, hox).

Methylation changes in intronic and intergenic regions in these areas could significantly influence gene expression, while the presence of DMCpGs in promoters and transcription start sites (TSS) as well as exons suggests potential direct impacts on gene transcription. Additionally, a GO enrichment analyses for annotated genes, highlighted several enriched regulatory processes, such as the regulation of signaling (GO:0023051) and response to stimuli (GO:0050896), particularly in the mature F1/Wild males’ comparison. These findings are especially relevant given the challenges associated with farming this species.

Upon differential expression analysis, we focused on specific genes related to reproduction [36,37] identified through DMCpGs annotation in both F1 and wild males (Figure 4) to study whether they exhibited differential expression. Some of the genes related to reproduction, such as those in the sox transcription factor family and the estrogen receptor genes, showed differentially methylated CpG sites (DMCpGs) in their vicinity according to gene annotation. However, these genes did not exhibit expression differences that could be directly associated with the observed methylation differences. This suggests methylation can act as a regulatory mechanism that may influence gene expression over time or under specific conditions; while methylation changes are present, they may not always translate into immediate changes in gene expression [38].

When comparing mature F1 and wild groups, egr1 and nanos1 showed higher methylation in the mature F1 group compared to its wild counterpart and were also down-regulated (LogFC > -2). However, when comparing immature F1 and wild-type groups, despite several genes showing methylation differences between the groups, no significant differences in the expression of these genes were observed.

In the comparison between immature and mature wild groups, GATA6 showed greater methylation in immature males while being down-regulated compared to their mature counterparts (LogFC = -1.40). Conversely, dmrt2a showed lower methylation levels and was up regulated in the immature group compared to their mature counterparts (LogFC = 1.92).

Lastly, when comparing immature and mature F1 groups, dmrt2b, esr2a, hsd3b7, and sox1a showed greater methylation levels and lower expression in the immature group compared to their mature counterparts, although the expression differences were less significant (LogFC < 1.5).

Figure 4. Average methylation differences (mean LogFC) in gene-associated DMCpGs among male groups of study in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFM/MWM (Mature F1 Versus Wild male), IFM/IWM (Immature F1 Versus Wild Male), IWM/MWM (Immature Versus Mature Wild Male), IFM/MFM (Immature Versus Mature F1 Male).

Figure 4. Average methylation differences (mean LogFC) in gene-associated DMCpGs among male groups of study in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFM/MWM (Mature F1 Versus Wild male), IFM/IWM (Immature F1 Versus Wild Male), IWM/MWM (Immature Versus Mature Wild Male), IFM/MFM (Immature Versus Mature F1 Male).

Figure 4. Average methylation differences (mean LogFC) in gene-associated DMCpGs among female groups of study in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFF/MWF (Mature F1 Female Versus Mature Wild Female), IFF/IWF (Immature F1 Female Versus Immature Wild Female), IWF/MWF (Immature Wild Female Versus Mature Wild Female), IFF/MFF (Immature F1 Female Versus Mature F1 Female).

Figure 4. Average methylation differences (mean LogFC) in gene-associated DMCpGs among female groups of study in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFF/MWF (Mature F1 Female Versus Mature Wild Female), IFF/IWF (Immature F1 Female Versus Immature Wild Female), IWF/MWF (Immature Wild Female Versus Mature Wild Female), IFF/MFF (Immature F1 Female Versus Mature F1 Female).

3.3. Outline of Transcriptome Profiles and Differential Expression Analysis

A principal component analysis (PCA) of the overall samples revealed significant differences between the sexes. The analysis showed that the samples formed two distinct clusters, one for each sex, indicating a clear sexual dimorphism (Figure A1). Additionally, the PCA was conducted independently for each sex, considering groups of different origins and maturity. When only samples of the same sex were considered, it was possible to differentiate a different distribution among those of different origin, and a better grouping of those that shared up bringing. Within each sex, F1 and wild groups were distinguishable, particularly among female samples. These differences were also reflected in the differential gene expression analysis (Figure 5).

The comparisons between samples of different sexes revealed the highest number of differentially expressed genes (DEGs). This was followed by differences in maturation status, and lastly, by differences between wild-type and F1 individuals.

Figure 5. Counts for differentially expressed genes (y-axis) across pairwise comparisons (x-axis) in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFM/MWM (Mature F1 Male Versus Mature Wild Male), IFM/IWM (Immature F1 Male Versus Immature Wild Male), IWM/MWM (Immature Wild Male Versus Mature Wild Male), IFM/MFM (Immature F1 Male Versus Mature F1 Male), MFF/MWF (Mature F1 Female Versus Mature Wild Female), IFF/IWF (Immature F1 Female Versus Immature Wild Female), IWF/MWF (Immature Wild Female Versus Mature Wild Female), IFF/MFF (Immature F1 Female Versus Mature F1 Female), MFF/MFM (Mature F1 Female vs Mature F1 Male), IFF/IFM (Immature F1 Female Versus Immature F1 Male), MWF/MWM (Mature Wild female Versus Mature Wild Male), and IWF/IWM (Immature Wild F1 Versus Immature Wild Male).

Figure 5. Counts for differentially expressed genes (y-axis) across pairwise comparisons (x-axis) in Solea senegalensis. Notes: Group comparisons keys are shown as following: MFM/MWM (Mature F1 Male Versus Mature Wild Male), IFM/IWM (Immature F1 Male Versus Immature Wild Male), IWM/MWM (Immature Wild Male Versus Mature Wild Male), IFM/MFM (Immature F1 Male Versus Mature F1 Male), MFF/MWF (Mature F1 Female Versus Mature Wild Female), IFF/IWF (Immature F1 Female Versus Immature Wild Female), IWF/MWF (Immature Wild Female Versus Mature Wild Female), IFF/MFF (Immature F1 Female Versus Mature F1 Female), MFF/MFM (Mature F1 Female vs Mature F1 Male), IFF/IFM (Immature F1 Female Versus Immature F1 Male), MWF/MWM (Mature Wild female Versus Mature Wild Male), and IWF/IWM (Immature Wild F1 Versus Immature Wild Male).

To identify expression differences across general categories, samples were grouped to compare the broader categories of sex, maturity, and origin. Considering the total number of individuals, comparisons were made between males and females, mature and immature individuals, and F1 and wild individuals. Significant expression differences were observed in all comparisons, with a larger number of down-regulated genes in males, mature individuals, and wild-type specimens compared to their respective counterparts (Figure 6).

Figure 6. Volcano plots for differentially expressed genes between S. senegalensis study groups: A) Females Versus Males B); C) Immature Versus Mature; C) F1 Versus Wild. Notes: The volcano plots were created using the ggplot package R Version 3.03 implemented in Novogene’s Novomagic Plaform. Fold change value and p-value are represented in the x-axis and y-axis, respectively. When the established significance limits are exceeded, the samples represented by dots are displayed coloured red (up-regulated) and blue (down-regulated).

Figure 6. Volcano plots for differentially expressed genes between S. senegalensis study groups: A) Females Versus Males B); C) Immature Versus Mature; C) F1 Versus Wild. Notes: The volcano plots were created using the ggplot package R Version 3.03 implemented in Novogene’s Novomagic Plaform. Fold change value and p-value are represented in the x-axis and y-axis, respectively. When the established significance limits are exceeded, the samples represented by dots are displayed coloured red (up-regulated) and blue (down-regulated).

3.4. Analysis of Differential Gene Expression and Patterns in Males

Comparative analysis of gene expression patterns in males revealed significant differences across groups. In the comparison between mature F1 and wild males, 154 genes were up regulated in mature F1 males, while 374 genes were down regulated in mature wild males, highlighting differences between these groups.

When comparing immature to mature wild males, 129 genes were up regulated in immature wild males, while 1,528 genes were down regulated in mature wild males, indicating a significant shift in gene expression during maturation.

Similarly, when comparing immature F1 and wild males, 169 genes were up regulated in immature F1 males, and 56 genes were downregulated in immature wild males, suggesting higher gene expression in F1 males.

Lastly, the comparison between immature and mature F1 males showed 302 genes up regulated in immature F1 males, and 179 genes down regulated in mature F1 males, indicating significant changes in gene expression during maturation.

Wild males exhibit a larger shift in gene expression during maturation when compared to F1 males. Additionally, F1 males, both mature and immature, tend to show greater expression in certain genes compared to their wild counterparts. Overall, these differences suggest both maturation and origin significantly influence gene expression patterns.

3.5. Analysis of Differential Gene Expression and Patterns in Females

A similar trend as the observed in male specimens was also observed in females. Specifically, in the comparison between mature wild and F1 females, 161 genes were up-regulated in mature F1 females, while 62 genes were down-regulated in mature wild females.

When comparing immature and mature wild females, 1,429 genes were up-regulated in immature wild females, whereas 1,830 genes were down-regulated in mature wild females.

Similarly, when comparing immature F1 and wild females, 706 genes were up-regulated in immature F1 females, and 528 genes were down-regulated in immature wild females.

Lastly, the comparison between immature and mature F1 females showed 27 genes up-regulated in immature F1 females and 889 genes down-regulated in mature F1 females.

As in males, wild females exhibit a larger shift in gene expression during maturation compared to F1 females. However, mature and immature F1 females showed greater expression in certain genes when compared to the wild groups.

Furthermore, the comparative analysis of gene expression patterns revealed substantial differences between males and females, highlighting the significant influence of sex on gene expression. Females, both immature and mature, tend to show greater higher expression in certain genes compared to their male counterparts, with more pronounced differences observed in wild groups compared to F1 groups. Additionally, the differences in gene expression between immature and mature individuals within the same sex (for both wild and F1 groups) suggest that maturation significantly influences gene expression, with sex-specific changes.

3.6. Annotation and Enrichment Analysis of Differentially Expressed Genes

Following the gene expression level analysis, a total of 26,174 genes were annotated. The GO enrichment analysis revealed 1,008 enriched GO terms, including those related to reproduction and reproductive processes (GO:0000003, GO:0022414). As previously described, 12 comparisons were conducted among groups (Figure A2), all of which showed significant differences.

The results indicate that gene expression varies significantly between sexes and developmental stages, particularly in genes related to reproduction and sex differentiation (Figure 7 and Table A3). Genes such as gsdf, fshr, and amh exhibited higher expression levels in both immature and mature male groups compared to female groups, underscoring their crucial roles in male reproductive functions. This pattern was further seen in several transcription factors showing elevated expression in males, particularly in mature wild and F1 males. Conversely, sox19b and hsd17b genes were more prominently expressed in females, indicating their involvement in female reproductive processes.

Figure 7. Bars graph representation of gene expression levels measured as gene counts for reproduction related genes across studied groups in Solea senegalensis. Notes: Group keys are shown as following: MWM (Mature Wild Male), IWM (Immature Wild Male), MFF (Mature F1 Female), IWF (Immature Wild Female), IFF (Immature F1 Female), IFM (Immature F1 Male), MFM (Mature F1 Male) and MWF (Mature Wild Female).

Figure 7. Bars graph representation of gene expression levels measured as gene counts for reproduction related genes across studied groups in Solea senegalensis. Notes: Group keys are shown as following: MWM (Mature Wild Male), IWM (Immature Wild Male), MFF (Mature F1 Female), IWF (Immature Wild Female), IFF (Immature F1 Female), IFM (Immature F1 Male), MFM (Mature F1 Male) and MWF (Mature Wild Female).

3.6. Integration of Methylation Data

We assessed the correlation between methylation levels and gene expression to study methylation's potential impact on transcriptional regulation. We focused on mature male sole groups to understand the relationship between DNA methylation and gene expression. By analyzing the methylation percentage in the reported cytosines from the coverage files and the expression levels as counts per million (cpm) of the annotated genes in the RNA-seq data, we aimed to investigate the potential repression of gene expression due to increased methylation.

Out of 74,523 CpG sites with variable methylation, 3,359 sites showed significant correlation (p-value < 0.05) with the expression of neighbouring genes. The mean correlation values remained consistent across different contexts (intro, exon, intergenic, promoter-TSS regions, TTS region), with a mean correlation value of |0.5 ± 0.28|; The 25% of CpGs exceeded a |0.7| correlation, while 32% ranged between |0.4| and |0.5|.

Among the intronic CpGs, 48% showed negative correlation, followed by intergenic (23%), exon (18%), promoter-TSS regions (9%), and TTS region (2%). Overall, methylation in the putative proximal promoter-TSS regions, exon, intron, and to a lesser extent the TTS region, exhibited the strongest negative correlations with gene expression. No significant correlation values differences were observed between first-intron/exon regions. Approximately 45% of the methylated sites negatively correlated with gene expression, while 55% showed a positive correlation (Table S1).

A Gene Ontology (GO) functional enrichment analysis was conducted on genes neighbouring DMCpGs that exhibited a significant correlation between methylation and gene expression. Regarding those DMCpGs with a negative correlation, most significant terms (p-value <0.05) included developmental processes, transcriptional regulation, and stimulus response pathways (GO:0032502, GO:0140110, GO:0042221, GO:0009719). Additionally, calcium ion binding (GO:0005509) was significantly enriched when exclusively analysing genes with a strong correlation (r <-0.9).

Notably, the enrichment analysis of genes neighbouring DMCpGs with a positive correlation revealed fewer significant enriched GO terms, with the most frequent being cell adhesion (GO:0007155) and regulation processes (GO:0010646, GO:0023051). Similar results were observed when analysing genes with a strong (r >0.9) positive correlation, including transcription regulator activity (GO:0140110) and developmental processes (GO:0032502).

Additionally, The KEGG enrichment analysis revealed significantly (p<0.05) different pathways regarding those genes to which DMCpGs were associated that showed a negative correlation (Figure A3), including some relevant in reproduction (MAPK, TGF-beta signaling pathways) [39,40].

4. Discussion

This study addresses this complexity by analyzing the methylome and gonadal transcriptome of S. senegalensis, identifying potential interactions on specific genes related to reproduction, and examining the effects of sex, maturity, and rearing conditions on these mechanisms.

4.1. RRBS Methylation Analysis of Gonadal Tissue

The sequencing of RRBS libraries from gonadal tissue revealed that over 2.6 billion cytosines in the CpG context showed methylation, accounting for about 53.63% of the total analyzed cytosines. While CpG islands, often associated with genes and their promoter regions, may exhibit different methylation patterns [41], this percentage is slightly lower than the typical 60%–80% methylation observed in vertebrate genomes [42]. This discrepancy could be attributed to the species-specific genomic and epigenetic landscape, including differences in the distribution and density of CpG sites, as well as the presence of unmethylated CpG islands. Additionally, the specific tissue type (gonadal tissue) and the developmental stage of the samples analyzed could also influence the overall methylation levels. Overall, while the percentage is lower than the typical range, it remains within a biologically plausible range given these factors [43].

The mean mapping efficiency was 62%, which falls towards the lower end of the typically accepted range (>60%), likely due to biases from the bisulfite conversion process, the repetitive characteristics of CpG regions targeted by RRBS, and possible discrepancies with the reference genome [44].

The analysis of DMCpG levels revealed significant differences when comparing sex, rearing conditions, and maturity groups. Primary differences in overall methylation were observed between males and females, indicating a substantial role of sex in methylation patterns. This could be due to the specificity of ovarian and testicular gonadal tissues, among other factors [36]. Grouping samples by sex, maturity, and origin allowed for the identification of methylation differences across these broad categories, with significant differences observed in all comparisons. As further discussed, higher methylation levels were found in females, immature individuals, and F1 individuals compared to their counterparts.

Distinct methylation patterns were observed within each sex, distinguishing F1 and wild groups. F1 males exhibited the highest methylation levels, particularly those with reproductive dysfunction, who showed a higher number of up-methylated CpGs compared to wild-type males. This suggests that captivity may significantly impact methylation patterns and potentially affect reproductive abilities [21].

Previous studies suggest that changes in gene expression caused by raising salmonids in captivity can affect the fitness of the first generation of hatchery fish, highlighting the potential role of methylation after just one generation [45]. Raising organisms in captivity often results in stable conditions, with minimal fluctuations in factors like temperature, salinity, pH, and daylight. In European seabass, the role of temperature in the methylation of genes related to reproduction has been highlighted, revealing a complex epigenetic layer influenced by both genetic background and early developmental environment, which contributes to the sexual phenotypic outcome [46].

Significant differences were found between immature individuals compared to their mature counterparts in both sexes, with females generally exhibiting higher methylation levels than males. These findings suggest a potential sex-specific role of methylation patterns influencing maturity. Similarly, in the model organism Danio rerio, sex-specific epigenetic regulation and transcription of genes involved in reproduction have been described when analyzing gonadal tissue, showing both positive and negative correlations [36].

While the limited number of F1 samples and lack of replicates must be considered, the RRBS data analysis revealed significant differences among methylation profiles. Overall, differential methylation comparisons showed the largest number of DMCpG sites between different maturation statuses, followed by rearing conditions, and finally between sexes. These findings suggest a potential role of environmental factors in shaping methylation profiles in S. senegalensis, which play a crucial role in maturation [47].

4.2. Annotation and Distribution of Differentially Methylated CpG Sites

Methylation changes in intronic and intergenic regions could significantly influence gene expression by affecting regulatory elements, while the presence of DMCpGs in promoters/TSS and exons suggests potential direct impacts on gene transcription. In our study, the majority of DMCpGs were in intronic and intergenic regions, with fewer in promoters/TSS, exons, and TTS. This distribution suggests that methylation changes are more prevalent in non-coding regions, which could influence gene regulation indirectly [38,46].

The analysis of DMCpGs in testicular samples from mature male Senegalese sole revealed intricate relationships between methylation and gene expression. For instance, although the sox transcription factor gene family was identified through DMCpG annotation by proximity, this did not result in significant expression differences in the differential expression analysis. Previous findings in Oreochromis niloticus showed that sox genes have stage-specific and sexually dimorphic expression, highlighting their complex roles in sex-linked functions. The overlapping distribution of Sox genes across tissues, their multifunctionality, and functional overlap, along with species-specific genetic mechanisms and diverse sex-determining systems, contribute to the complexity and uncertainty in their functions [47].

On the other hand, the observed differences in methylation between mature F1 and wild groups, particularly in genes like nanos1 and egr1, suggest that epigenetic mechanisms play a significant role in regulating gene expression. These mechanisms can both repress and enhance gene activity, highlighting their importance in the modulation of genetic functions. The higher methylation and down-regulation of these genes in the mature F1 group indicate a potential link between methylation and gene suppression. Nanos1 showed up to six intergenic DMCpGs: one displaying a strong positive correlation (r = 0.9) and the others showing weaker negative correlations ranging from -0.2 to -0.6. Similarly, several DMCpG sites in both intronic and intergenic regions were associated with err2a and esrrb, showing positive (r = 0.5) and negative (r = -0.3) correlations, respectively.

Regarding the observed differences in methylation and expression between immature and mature wild groups, GATA6 showed greater methylation and down-regulation in immature males, while dmrt2a exhibited lower methylation levels and was up-regulated in the same group. The varying correlation values between methylation and expression for GATA6, dmrt2b, and dmrt3b underscore the complexity and specificity of methylation's impact. For example, GATA6 had four neighbouring DMCpGs with varying negative correlation values (r = -0.3 ~ -0.7), while intergenic DMCpGs associated with dmrt2b and dmrt3b exhibited strong negative correlations (r = -0.95 and r = -0.8). These findings suggest that methylation can influence gene expression in a gene-specific manner, and its effects can differ significantly even among closely related genes.

In regard of the implication of the studied genes, nanos1 has been identified in various fish species and plays a crucial role in germ cell differentiation [48]; A deficiency in this gene is associated with a reduction in primordial germ cells [49], suggesting its essential role in the migration of PGCs and preserving germline identity [50,51]. Additionally, previous studies have described the roll of GATA6 in gonadal development and reproduction in C. semilaevis [52] and the of roles of estrogen receptors (esr) in the molecular mechanism of sex determination and differentiation in Hippoglossus hippoglossus [53].

The dmrt genes (double-sex and mab-3-related transcription factor) are crucial for testes development and functioning [54,55]. In the flatfish C. semilaevis, dmrt1, dmrt2, and dmrt3 play significant roles in male germ cell maturation and gonadal differentiation. [56] found that the expression of these genes is negatively correlated with the methylation level of their associated promoters, indicating that methylation may regulate their expression and, consequently, their role in reproductive processes.

The methylation differences between two differently reared groups suggest that captive conditions may play a role in shaping methylation patterns. The significant correlation between methylation and expression indicates that these differences may affect certain biological processes and how they unfold in wild versus captive-bred individuals. Our results also highlight complex dynamics regarding dmrt genes and their methylation/expression, the influence of other regulatory factors as well as differences in how maturity affects the expression of these genes in wild-type and F1 individuals.

Similar differential methylation analysis performed in females resulted in different genes related to reproduction to which different DMCpGs were associated by gene annotation. However, integration with transcriptomic data and correlation analysis did not show CpGs with significant correlation (p < 0.05) neighbouring these genes. These results highlight the distinct epigenetic mechanisms between males and females, potentially due to tissue specificity and the involvement of other regulatory factors [57].

4.3. Differential Gene Expression in Sex Differentiation and Reproductive Development

Significant differences in both methylation and gene expression were observed when comparing individuals of different sexes. This was evident in both the differentially expressed PCA and the MDS based on methylation differences measured in fold change. Among samples of the same sex, distinct separation was observed between those with different rearing conditions. This separation was most notable in males when analyzing methylation differences, and in females when analyzing expression differences. The highest number of differentially expressed genes (DEGs) was observed between sexes, followed by differences in maturation status, and lastly, by differences between wild-type and F1 individuals. This pattern suggests that sex is the most influential factor in gene expression differences, while maturation and origin also play significant roles.

The differential gene expression analysis further highlighted substantial differences between males and females, with females generally showing higher expression compared to males. These differences were more pronounced in wild groups compared to F1 groups. Additionally, maturation significantly influenced gene expression within the same sex, with sex-specific changes. Wild individuals exhibited a larger shift in gene expression during maturation compared to F1 individuals. While sex seems like a primary factor when comparing gene expression and methylation patterns, with significant differences observed between males and females, maturation and origin (wild Versus F1) also play important roles, especially in wild groups.

Particularly in genes related to reproduction and sex differentiation, the gene expression analysis revealed significant sex-specific differences; Genes such as gsdf, fshr, and amh exhibited higher expression levels in male groups both immature and mature compared to female groups. This was further observed in dmrt and sox genes, indicating their crucial roles in male reproductive functions, such as testicular development and spermatogenesis. In contrast, sox19b, and hsd17b showed higher expression in females, indicating their involvement in such as in steroidogenesis, ovarian development and general female reproductive processes [19,58].

The gsdf gene, which is vital for sex differentiation in fish, showed higher expression levels in immature males (both wild and F1) compared to mature groups, indicating its importance in the early stages of male development, and was nearly absent in females. This highlights its significant role in the sex differentiation processes of S. senegalensis. Additionally, the differential expression of the fshr and amh genes suggests their predominant role in the differentiation of male gonads. In teleosts, gsdf functions mainly in gonad differentiation, with its expression, regulation, and function varying significantly across species. It is usually expressed in Sertoli cells surrounding spermatogonia, suggesting its involvement in self-renewal, proliferation, and differentiation of spermatogonia [59]. This aligns with our findings, where gsdf showed higher expression in immature males, underscoring its role in early male development.

While the pathways and genes involved in sex differentiation are still not well understood, and the reproductive challenges faced by hatchery-born individuals remain unresolved, the fshr gene was identified as playing a key role in sex determination in this species [19]. Although no significant differences in methylation or expression related to fshr were detected with the statistical criteria chosen for the group comparison analysis, genetic quantification reflects a lower expression in F1 mature males compared to wild-type males. Previously, Sambroni et al. (2013) studied FSH control on gene expression in fish both independently of and through steroid mediation [60]. This indicates that FSH can regulate gene expression via multiple pathways, which could explain the observed differences in fshr expression between wild and F1 mature males in S. senegalensis.

4.4. Methylation Integration

By integrating methylation coverage and expression analysis, we sought to identify significant correlations between CpG methylation and the expression of chosen key reproductive genes. Upon integrating methylation and expression data, significant positive correlations in CpGs neighboring dmrt1 (r = 0.99), gnrh (r = 0.95) and cyp1b1 (r = 0.74) were observed in mature females. Only transcription factor GATA6 showed a significant negative correlation (r = -0.68). These correlations were observed in intergenic and intronic regions, respectively.

In C. semilaevis, Liu et al. (2016) described the role of GATA6 in gonadal development and reproduction [52]. While no significant differences were found in overall methylation levels, the patterns of methylation sites varied among males, females, and pseudo males with the latter showing higher expression of GATA6 compared to females. This study underscores the significant role of methylation in gene expression regulation. Similarly, our results highlight the crucial impact of methylation on gene expression, suggesting that epigenetic mechanisms are essential in regulating reproductive traits in Senegalese sole.

The absence of these genes when previously analyzing significant methylation differences between mature females suggests that while these genes are regulated by methylation, the differences in methylation levels between individuals may not be large enough to be detected as significant in a differential methylation analysis.

In contrast, the lack of significant correlations (r > 0.5) in male testes for selected genes indicates the influence of other regulatory factors in certain reproductive mechanisms, such as hormonal interactions [61]. These findings underscore the complex regulatory mechanisms governing reproduction in this species.

Our study further focused on genes with significant negative methylation/expression correlations to understand the sexual dysfunction observed in male sole. The enrichment analysis revealed significant enrichment of the MAPK, TGF-beta, and calcium signaling pathways. These pathways are crucial for cell differentiation, gonadal development, and signal transduction in maturation and cell activation processes, highlighting their importance in reproductive processes.

Gonadotropin-releasing hormone (GnRH) activates the MAPK cascade via the protein kinase A (PKA) signaling pathway, influencing gene expression, particularly for GPα and LHβ, but not for FSHβ. This underscores the critical role of GnRH signaling in regulating reproductive hormones and processes [39,62,63,64]. In the DMCpGs analysis of mature wild versus F1 males, the MAPK pathway was associated with up-methylated CpGs in F1 males. Differential expression analysis showed higher expression of gnrh in F1 groups, while fshr was more expressed in mature wild males. De la Herrán et al. (2023) observed consistent overexpression of the fshr gene in males at all developmental stages, from undifferentiated primordium to mature adults, indicating its crucial role in male reproductive development [19]. These results highlight the differences in reproductive gene regulation between wild and F1 males, with potential implications for their reproductive capabilities and success. This suggests differences in hormonal regulation between these groups, which could be due to environmental factors and the implication of potential epigenetic modifications. Additionally, the GnRH signaling pathway was enriched when analysing up-methylated genes in immature wild males versus the mature group, highlighting differences in methylation patterns within the MAPK pathway during the transition from immaturity to maturity in wild versus F1 individuals.

Furthermore, when studying the methylation/expression correlation in mature male sole, receptors belonging to the OR1G1 family showed a significant positive correlation (r > 0.7). When analyzing this correlation between F1 and wild-type males independently, no significant differences were found. Chauvigné et al. (2016) reported that mRNA for olfactory receptors and other reproduction-related transcripts, such as cytochrome P450, are differentially expressed in the upper olfactory epithelia of wild-caught and F1 sole [16]. This suggests that urine and intestinal fluid may play crucial roles in chemical communication during reproduction in sole, with the regulation at the olfactory epithelium level being sensitive to both the sex and maturity of the receiver. The physiological role of mature female urine in stimulating the release of LH from the pituitary gland and increasing spermatogenesis and milt production in males as previously described suggests that methylation may play a role in optimizing reproductive success.

Functional enrichment analyses further support the involvement of methylation in key biological processes and pathways linked to reproduction. These insights contribute to our understanding of the complex interplay between epigenetic modifications and gene expression in this species. Future studies could focus on tracking methylation and gene expression changes over time in the same individuals, providing insights into the dynamics of these processes. Additionally, comparative studies across different environmental conditions, such as temperature changes, which have been shown to have an epigenetic influence, could elucidate how these factors affect methylation.

In the absence of further validation, we strictly considered those genes with significant differences and variations in methylation alongside a significant correlation. However, previous studies have shown that minor alterations in DNA methylation levels can act as drivers for the epigenetic regulation of gene expression. For instance, in Atlantic salmon, small changes in DNA methylation at CpG sites were found to correlate with transcriptional changes, suggesting that even minimal methylation variations can serve as regulatory epimarkers. Overall, our findings highlight the intricate nature of epigenetic regulation in reproductive mechanisms and emphasize the need for further research to fully understand the interplay between methylation and other regulatory factors in gene expression [38].

5. Conclusions

This study explores the methylome and gonadal transcriptome of S. senegalensis, focusing on the interactions of specific genes related to reproduction and the effects of sex, maturity, and rearing conditions. Our findings suggest that DNA methylation significantly regulates gene expression in S. senegalensis, acting as a mechanism that can both repress and enhance gene expression depending on the genomic context. The complexity of this epigenetic mechanism is evident from the varying levels of methylation and correlation rates, both positive and negative, observed in different DMCpGs associated with specific genes linked to reproduction. This variability could explain the differences in methylation and expression when comparing groups, as many genes show different associated DMCpGs with varying correlation rates.

Differential methylation comparisons revealed the highest number of DMCpG sites between both maturation stages studied, followed by rearing conditions, and then between sexes. Additionally, the analysis of differential gene expression showed substantial differences between males and females, with females generally exhibiting higher expression levels. These differences were more pronounced in wild groups compared to F1 groups. The comparison between maturation stages revealed significant gene expression differences within the same sex, with wild individuals showing a larger shift during maturation compared to F1 individuals. While sex is a major factor in gene expression and methylation patterns, maturation and origin (wild Versus F1) also play crucial roles, particularly in wild groups.

These findings suggest that methylation plays a crucial role in regulating gene expression in processes such as sexual differentiation and maturation. However, the influence of other factors on this regulation must also be considered. Our findings indicate that while rearing conditions do affect methylation patterns and gene expression, these effects must be evaluated alongside other factors, such as sex-specific mechanisms during maturation. This underscores the complex interaction between epigenetic modifications and gene expression in Senegalese sole.