# Insights into Solea senegalensis Rreproduction through gonadal tissue methylation analysis and transcriptomic integration

Daniel Ramírez, Marco Anaya-Romero, María Esther Rodríguez, Alberto Arias-Pérez, Robert Mukiibi, Helena d'Cotta, Diego Robledo, Laureana Rebordinos

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Daniel Ramírez, Marco Anaya-Romero, María Esther Rodríguez, Alberto Arias-Pérez, Robert Mukiibi, et al.. Insights into Solea senegalensis Rreproduction through gonadal tissue methylation analysis and transcriptomic integration. Biomolecules, 2025, 15 (1), pp.54. 10.3390/biom15010054. hal-05182862

HAL Id: hal-05182862

https://hal.science/hal-05182862v1

Submitted on 23 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Article

# Insights into *Solea senegalensis* Reproduction Through Gonadal Tissue Methylation Analysis and Transcriptomic Integration

Daniel Ramírez <fnref n="1" />, Marco Anaya-Romero <fnref n="1" />,†, María Esther Rodríguez <fnref n="1" /> ID, Alberto Arias-Pérez <fnref n="1" /> ID, Robert Mukiibi <fnref n="2" />, Helena D'Cotta <fnref n="3" />, Diego Robledo <fnref n="1" />,<fnref n="2" /> and Laureana Rebordinos <fnref n="1" />,* ID

**Abstract:** Fish exhibit diverse mechanisms of sex differentiation and determination, shaped by both external and internal influences, often regulated by distinct DNA methylation patterns responding to environmental changes. In *S. senegalensis* aquaculture, reproductive issues in captivity pose significant challenges, particularly the lack of fertilization capabilities in captive-bred males, hindering genetic improvement measures. This study analyzed the methylation patterns and transcriptomic profiles in gonadal tissue DNA from groups differing in rearing conditions and sexual maturity stages. RRBS (Reduced Representation Bisulfite Sequencing) was employed to detect notable methylation variations across groups, while RNA was extracted and sequenced for differential expression analysis. Our findings suggest that DNA methylation significantly regulates gene expression, 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 observed in different CpGs neighboring specific genes linked to reproduction. Differential methylation comparisons revealed the highest number of differently methylated CpGs between maturation stages, followed by rearing conditions, and lastly between sexes. These findings underscore the crucial role of methylation in regulating gene expression and its potential role in sex differentiation, highlighting the complex interplay between epigenetic modifications and gene expression.

**Keywords:** *Solea senegalensis*; flatfish; reproduction; methylation; transcriptome

## 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,

<figure><img src="image_3.png" /></figure>

check for
updates

Academic Editor: Shicui Zhang

Received: 27 November 2024

Revised: 30 December 2024

Accepted: 2 January 2025

Published: 2 January 2025

**Citation:** Ramírez, D.; Anaya-Romero, M.; Rodríguez, M.E.; Arias-Pérez, A.; Mukiibi, R.; D'Cotta, H.; Robledo, D.; Rebordinos, L. *Insights into Solea senegalensis* Reproduction Through Gonadal Tissue Methylation Analysis and Transcriptomic Integration.
*Biomolecules* **2025**, *15*, 54. https://doi.org/10.3390/biom15010054

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">1 Área de Genética, Facultad de Ciencias del Mar y Ambientales, INMAR, Universidad de Cádiz, 11510 Cádiz, Spain; daniel.ramirez@uca.es (D.R.); alberto.arias@uca.es (A.A.-P.)</footnote>

<footnote n="2">2 The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Edinburgh EH259RG, UK</footnote>

<footnote n="3">3 UMR-Institut des Sciences de l'Evolution de Montpellier, Centre National de la Recherche Scientifique, Institut de Recherche, Pour le Développement, Ecole Pratique des Hautes Etudes, University of Montpellier, 34295 Montpellier, France</footnote>

<footnote>* Correspondence: laureana.rebordinos@uca.es</footnote>

<footnote>† Current address: Centro de Investigaciones en Ciencias de la Vida, Facultad de Ciencias Básicas y Biomédicas, Universidad Simón Bolívar, Barranquilla 080005, Colombia.</footnote>
some species exhibit interactions between both mechanisms [1–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* has led the significant expansion and investment in the European aquaculture industry in recent years. This flatfish species boasts a high growth rate, along with a high market value and strong demand. Its 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, controlling disease, 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 a reduced fertilization capacity due to their low sperm volume, which hinders the implementation of genetic improvement measures through selective breeding [8–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–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 of 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].

Preceding research on *S. senegalensis* has investigated how temperature influences the epigenetic regulation of muscle development. This is particularly relevant due to the interest in cultivating larger female populations for aquaculture purposes. Given the interest in cultivating larger female populations, a previous study considered 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 greater methylation levels in F1 individuals than wild fish, 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 different-origin 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 influenced by a variety of external and internal factors, including temperature, salinity, population density, nutrients, and
contaminants. These factors often exert their effects through specific DNA methylation patterns.

Temporary solutions have been developed to enhance the reproductive status, including artificial fertilization and nutritional strategies, as well cryopreservation protocols [23–25]. Although hormonal induction in F1 sole has achieved some success [26,27], the 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–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 primary aim of this study was to examine the variations in methylation patterns within the gonadal tissue (ovaries and testes) of Senegalese sole and to explore how methylation functions as an epigenetic mechanism in gene expression regulation 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 groups: one that was artificially bred (F1) at IFAPA and wild fish from the 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. This study included eight groups based on the previous criteria: immature wild (IW), mature wild (MW), immature F1 (IF), and mature F1 (MF). 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). Wild fish were collected using trawl nets in coastal areas. Cultured fish sample sizes were limited due to availability and high early-stage mortality rates.

Size and weight were used as indicators of sexual maturation. Mature specimens were obtained during spawning months, and their maturity was confirmed by performing gonadal tissue dissections; a portion of each sample was used for histological analysis to verify the maturation stages (Figure A1). The histological process involved standard procedures: fixation in formalin, dehydration through a graded ethanol series, paraffin embedding, microtome sectioning, and staining with hematoxylin and eosin.

F1 specimens reached sexual maturity under stable conditions. Factors like feeding patterns, population density, and oxygen levels can serve as epigenetic modulators. Even after reaching sexual maturity, F1 individuals encountered reproductive difficulties, such as lack of reproductive behavior. The fish were transferred to the laboratory in oxygenated containers and anesthetized with clove oil (90–100 mg/L). Both wild and cultured fish were maintained in oxygenated containers for 1 to 2 h before necropsy to reduce stress and ensure optimal conditions. Gonadal samples were stored in RNAlater® TissueProtect at −80 °C until DNA extraction. The experimental protocol adhered to EU guidelines (86/609/EU) and received approval from the Ethics Committee of the University of Cadiz, Spain.

## 2.2. Gonadal DNA and RNA Samples' Extraction

DNA was extracted from testicular tissue using the DNeasy Blood & Tissue Kit (Qiagen). Following extraction, the DNA samples underwent purification and size selection using the HighPrep™ PCR Clean-up System (MagBio Genomics Inc., Gaithersburg, MD, USA). A NanoDrop™ 2000C spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used for quality assessment, while concentration was measured using a Qubit 4 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).

RNA extraction was performed from 0.5 mg of gonadal tissue using TRIzol™ Reagent (Thermo Fisher Scientific, Waltham, MA, USA), following the manufacturer's protocol. Quantification and purity assessment of the extracted RNA were conducted via fluorometry using the Qubit system (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity (RIN) was determined using a 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 h using the Diagenode Premium Kit (Hologic Diagenode, Marlborough, MA, USA). Following end preparation, adaptor ligation, and size selection, sample concentrations were measured using qPCR. Samples with similar concentrations were pooled in groups of six and subjected to 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 (© Agilent Technologies, Inc., Santa Clara, CA, USA). Libraries were measured 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 remove Illumina sequencing adapters, low-quality base calls (Phred quality score < 20), and short reads (<20 bp). To ensure the quality of the fastq file data, we employed FastQC 0.12.0.

We performed RRBS alignment, through an in silico bisulfite conversion of the reference genome following the last version of the *S. senegalensis* published genome (GenBank assembly GCA_919967415.2 from NCBI). The RRBS sequence reads' alignment to the in silico converted *S. senegalensis* reference genome, as well as the subsequent methylation analysis 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, multidimensional 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. CpG sites with significant differential methylation were identified using a false discovery rate threshold of 0.01 (Benjamini-Hochberg correction).

Using categories based on the criteria for sample collection (sex, maturity, and origin), 12 comparisons were conducted for both differential methylation analysis and transcriptionic profile comparisons between groups. These comparisons aimed to examine the potential effects of those 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 Female vs. Immature Wild Male (IWF/IWM).

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

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/, accessed on 1 November 2024).

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 farmed juvenile and adult sample size was restricted, potentially reducing certain variability. Specifically, we could not achieve more than two F1 replicates for the male groups due to their short supply, which may have introduced some bias into the analysis. Furthermore, the mature wild female group lacked a third replicate. Despite these limitations, we believe the experiment remained scientifically valuable for studying *S. senegalensis*. 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/, accessed on 1 November 2024).

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 Patterns' Study and Differential Analysis*

A general multidimensional scaling plot of global methylation levels reveals primary differences among 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).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Multidimensional scaling (MDS) plot for methylation profiles among *S. senegalensis* samples (created using the edgeR package). Notes: The variations in methylation profiles among *S. senegalensis* samples are shown for the eight study groups—mature wild males and females, immature wild males and females, mature F1 males and females, and immature F1 males and females. Groups are each depicted in distinct colors. Individual samples are represented by colored tags, with each corresponding to its respective group. The proximity reflects the similarity among samples within groups.</figcaption></figure>

Taking into account the study groups, we conducted 12 comparisons of differential methylation profiles (Table A4). Only CpG sites with adequate 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 considered 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 overall mean 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.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Total number of DMCpGs (y-axis) across group comparisons (x-axis). Group comparison keys are shown as follows: 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 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 Female vs. Immature Wild Male (IWF/IWM).</figcaption></figure>

Focusing on the analysis of male groups, among those showing a marked reproductive dysfunction, generally a larger number of up-methylated CpGs were observed in F1 males when compared to wild males, and in immature individuals 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 exhibited 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.

When classifying samples 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 within 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 was 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. Despite its limitations, our study provides valuable insights and we hope it will inspire future research with larger sample sizes.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Volcano plots for methylation differences between *S. senegalensis* groups: (A) F1 vs. Wild; (B) Females vs. Males; (C) Immature vs. Mature. Notes: The volcano plots were generated using the R package gplots v3.1.1. The x-axis represents the fold change value, while the y-axis shows the *p*-value. Dots representing samples that exceed the established significance limits are colored green (up-methylated) and red (down-methylated).</figcaption></figure>

### 3.3. Genomic Annotation Based on Differentially Methylated Regions

The reported DMCpG sites were identified and annotated according to 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 follows: gene promoters/transcription start sites (TSS, ±1 kbp), introns, exons, transcription termination sites (TTS, ±1 kbp), or intergenic regions.

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

Upon annotation, we found that the identified DMCpG sites were located nearby 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*, *kdm6ba*).

Our findings revealed CpG methylation differences among certain groups linked to genes that play a role 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 (TSSs) as well as exons suggests potential direct impacts on gene transcription. Additionally, a GO enrichment analysis of the 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.

In the differential expression analysis, we focused on specific genes related to reproduction [36,37] identified through DMCpGs' annotation, with special attention paid to males of different origins (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, had 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 that methylation can act as a regulatory mechanism influencing gene expression
over time or under specific conditions; yet, while methylation changes are present, they
may not always translate into immediate changes in gene expression [38].

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Average methylation differences (mean LogFC) in gene-associated DMCpGs among female groups of study in *Solea senegalensis*. Notes: Group comparison keys are shown as follows: MFM/MWM (Mature F1 vs. Wild Male), IFM/IWM (Immature F1 vs. Wild Male), IWM/MWM (Immature vs. Mature Wild Male), IFM/MFM (Immature vs. Mature F1 Male), MFF/MWF (Mature F1 Female vs. Mature Wild Female), IFF/IWF (Immature F1 Female vs. Immature Wild Female), IWF/MWF (Immature Wild Female vs. Mature Wild Female), IFF/MFF (Immature F1 Female vs. Mature F1 Female).</figcaption></figure>

When we compared mature F1 and wild groups, *egr1* and *nanos1* exhibited greater 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.

When comparing immature and mature wild groups, *GATA6* showed greater methy-
lation in immature males, though it was down-regulated compared to their mature coun-
terparts (LogFC = -1.40). Conversely, *dmrt2a* showed lower methylation levels and was
up-regulated in the immature group compared to its mature counterparts (LogFC = 1.92).

Lastly, in the immature and mature F1 group comparison, *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).

3.4. *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 A2). Additionally,
the PCA was conducted independently for each sex, considering groups of different origins and maturities. When only samples of the same sex were considered, it was possible to differentiate the distributions among those of different origins, with a better grouping of those that shared the same upbringing. 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).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Counts for differentially expressed genes (y-axis) across pairwise comparisons (x-axis) in *Solea senegalensis*. Notes: Group comparison keys are shown as follows: MFM/MWM (Mature F1 Male vs. Mature Wild Male), IFM/IWM (Immature F1 Male vs. Immature Wild Male), IWM/MWM (Immature Wild Male vs. Mature Wild Male), IFM/MFM (Immature F1 Male vs. Mature F1 Male), MFF/MWF (Mature F1 Female vs. Mature Wild Female), IFF/IWF (Immature F1 Female vs. Immature Wild Female), IWF/MWF (Immature Wild Female vs. Mature Wild Female), IFF/MFF (Immature F1 Female vs. Mature F1 Female), MFF/MFM (Mature F1 Female vs. Mature F1 Male), IFF/IFM (Immature F1 Female vs. Immature F1 Male), MWF/MWM (Mature Wild Female vs. Mature Wild Male), and IWF/IWM (Immature Wild Female vs. Immature Wild Male).</figcaption></figure>

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.

To identify expression differences across general categories, samples were grouped to compare the broader categories of sex, maturity, and origin. Within 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).

### 3.5. *Analysis of Differential Gene Expression and Patterns in Males*

A comparative analysis of gene expression patterns in males showed notable differences among the groups. In the comparison between mature males with different origins, 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 1528 genes were down-regulated in mature wild males, indicating a significant shift in gene expression during maturation.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Volcano plots for differentially expressed genes between *S. senegalensis* study groups: (A) Females vs. Males; (B) Immature vs. Mature; (C) F1 vs. Wild. Notes: The volcano plots were generated using the ggplot package in R Version 3.03 on Novogene's Novomagic Plaform. The x-axis represents the fold change value, while the y-axis shows the *p*-value. Dots representing samples that exceed the established significance limits are colored red for up-regulated and blue for down-regulated, the green dots indicate values that fall below the established significance thresholds.</figcaption></figure>

Similarly, when comparing immature F1 and wild males, 169 genes were up-regulated in immature F1 males, and 56 genes were down-regulated 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 exhibited a larger shift in gene expression during maturation when compared to F1 males. Additionally, F1 males, both mature and immature, tended to show greater expression in certain genes compared to their wild counterparts. Overall, these differences suggest that both maturation and origin significantly influence gene expression patterns.

## 3.6. Analysis of Differential Gene Expression and Patterns in Females

A similar trend as that 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, 1429 genes were up-regulated in immature wild females, whereas 1830 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 a greater expression of 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, tended to show a higher expression of 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.7. 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 1008 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 A3), 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 an 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><img src="image_10.png" /><figcaption>**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 follows: MWM (Mature Wild Male), IWM (Immature Wild Male), MFF (Mature F1 Female), IWF (Immature Wild Female), IFF (Immature F1 Male), IFM (Immature F1 Male), MFM (Mature F1 Male), and MWF (Mature Wild Female).</figcaption></figure>

### 3.8. Integration of Methylation and Expression 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, 3359 sites showed a significant
correlation (p-value < 0.05) with the expression of neighboring 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|; 25% of CpGs exceeded
a |0.7| correlation, while 32% ranged between |0.4| and |0.5|.

Among the intronic CpGs, 48% showed a negative correlation, followed by intergenic (23%), exon (18%), promoter-TSS (9%), and TTS regions (2%). Overall, methylation in the putative proximal promoter-TSS regions, exons, introns, and to a lesser extent the TTS regions, 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 neighboring DMCpGs that exhibited a significant correlation between methylation and gene expression. Among the DMCpGs with a negative correlation, the 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 analyzing genes with a strong correlation (r < -0.9).

Notably, the enrichment analysis of genes neighboring DMCpGs with a positive correlation revealed fewer significantly enriched GO terms, with the most frequent being cell adhesion (GO:0007155) and regulation processes (GO:0010646, GO:0023051). Similar results were observed when analyzing 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 in the genes to which DMCpGs were associated, which showed a negative correlation (Figure A3), including some relevant to reproduction (MAPK, TGF-beta signaling pathways) [39,40].

**4. Discussion**

This study addresses a reported area of complexity by analyzing the methylome and
gonadal transcriptome of *S. senegalensis*, identifying potential interactions of 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 un-
methylated 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 an increased number of up-methylated CpGs relative 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 the genetic background and early developmental environment, which contributes to the sexual phenotypic outcome [46].

Significant differences were found in 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/TSSs 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/TSSs, exons, and TTSs. 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 a 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 their 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 neighboring 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 that its effects can differ significantly even among closely related genes.

In regard to 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 role of *GATA6* in gonadal development and reproduction in *C. semilaevis* [52] and the roles of estrogen receptors (esr) in the molecular mechanism of sex determination and differentiation in *Hippoglossus hippoglossus* [53].

The *dmrt* transcription factor's genes (double-sex and mab-3-related) are crucial for testes development and functioning [54,55]. In the flatfish *C. semilaevis*, *dmrt1*, *dmrt2*, and *dmrt3* play significant roles in the maturation of male germ cells and the differentiation of gonads. Reference [56] found that the expression of these genes is inversely related to the methylation level of their 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 cap-
tive 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 different-origin 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 an-
notation. However, integration with transcriptomic data and correlation analysis did not
show CpGs with a significant correlation (p < 0.05) neighboring 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]. This differentiation
is crucial in understanding the reproductive dysfunctions observed in *Solea senegalensis*,
which can impact the species' reproductive success and sustainability.

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 the 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 vs. F1) also play important roles, especially in wild groups.

Particularly in genes related to reproduction and sex differentiation, the gene expres-
sion 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 spermatogen-
esis. In contrast, sox19b, and hsd17b showed higher expression in females, indicating
their involvement, 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 reflected a lower expression in mature F1 males compared to wild-
type males. Previously, Sambroni et al. (2013) studied FSH control of gene expression in
fish both independently of and through steroid mediation [60]. This indicated that FSH
can regulate gene expression via multiple pathways, which could explain the observed
differences in *fshr* expression between wild and mature F1 males in *S. senegalensis*.

4.4. *Methylation and Expression Analysis Integration*

By integrating methylation coverage and expression analysis, we sought to identify
significant correlations between CpG methylation and the expression of chosen key re-
productive 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 notable differences were found in the total levels of methylation, the methylation patterns of certain sites varied among males, females, and pseudo males, with the latter showing a greater expression of *GATA6* in comparison 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 differ-
ences 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 impor-
tance in reproductive processes.

Gonadotropin-releasing hormone (GnRH) activates the MAPK cascade through the protein kinase A (PKA) signaling pathway, affecting gene expression, especially for GPα and LHβ, but not for FSHβ. This underscores the critical role of GnRH signaling in regulating reproductive hormones and processes [39,62–64]. In the DMCpG analysis of mature different-origin males, the MAPK pathway was associated with up-methylated CpGs in F1 males. Differential expression analysis showed a higher expression of *gnrh* in F1 groups, while *fshr* was more greatly expressed in mature wild males. De la Herrán et al. (2023) observed a 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. Furthermore, the GnRH signaling pathway was enriched when analyzing up-methylated genes in immature wild males versus the mature group, high-
lighting variations in methylation patterns within the MAPK pathway as fish transition
from immaturity to maturity in individuals of different origins.

Moreover, 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) found that mRNA for
olfactory receptors and other reproduction-related transcripts, including cytochrome P450,
show different expression levels 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 function
of mature females' urine in triggering the release of LH from the pituitary gland and
enhancing 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. To address reproductive issues in cultured fish, it is crucial to integrate
gene expression and methylation data with histological analysis, along with plasma hormone
levels, vitellogenin, and other relevant assessments. 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 the 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 the two maturation stages studied, followed by rearing conditions, and then be-
tween 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 differ-
ences 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 methyla-
tion patterns, maturation and origin (Wild vs. F1) also play crucial roles, particularly in
wild groups.

These findings suggest that methylation plays a crucial role in regulating gene expres-
sion 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 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.

**Supplementary Materials:** The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/biom15010054/s1, Table S1: Correlation rates between methylation coverage and gene expression in *Solea senegalensis*.

**Author Contributions:** Conceptualization, L.R.; Methodology, D.R. (Daniel Ramírez), M.A.-R. and M.E.R.; Software, D.R. (Daniel Ramírez), R.M. and A.A.-P.; Formal Analysis, D.R. (Daniel Ramírez), R.M., D.R. (Diego Robledo) and L.R.; Data Curation, D.R. (Diego Robledo) and L.R.; Writing—Original Draft Preparation, D.R. (Daniel Ramírez); Writing—Review and Editing, H.D., D.R. (Diego Robledo) and L.R.; Supervision, M.E.R., A.A.-P., D.R. (Diego Robledo) and L.R.; Project Administration, L.R.; Funding Acquisition, D.R. (Diego Robledo) and L.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Junta de Andalucía [FEDER P20-00938; FEDER CCMM-00014] and the University of Cádiz [R93REC/2019]. Diego Robledo acknowledges that he is funded by the Axencia Galega the Innovación (GAIN, Xunta de Galicia) as part of the Oportunius program, and by BBSRC Institute Strategic Grants to the Roslin Institute [BBS/E/20002172, BBS/E/D/30002275, BBS/E/D/10002070 and BBS/E/RL/230002A]. The open access fee was co-funded by the QUALIFICA Project (QUAL21-0019, Junta de Andalucía).

**Institutional Review Board Statement:** The experimental protocol followed EU guidelines (86/609/EU).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The original contributions presented in this study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

**Acknowledgments:** The authors express their gratitude to Manuel A. Merlo and Silvia Portela-Bens for their technical assistance in sample collection. They also extend their thanks to Manuel Manchado from IFAPA "El Toruño" in Cádiz, Spain, and the Servicio Central de Investigación en Cultivos Marinos (UCA) for supplying samples for this study.

**Conflicts of Interest:** The authors declare no conflicts of interest.

**Abbreviations**

CpGs (CpG Sites); DMCpGs (Differently Methylated CpG Sites); cpm (Counts Per Million); DNMTs (DNA Methyltransferases); F1 (Filial 1); GnRH (Gonadotropin-Releasing Hormone); GO (Gene Ontology); IFF (Immature F1 Female); IFM (Immature F1 Male); IWF (Immature Wild Female); IWM (Immature Wild Male); KEGG (Kyoto Encyclopedia of Genes and Genomes); LogFC (Log2 Fold Change); MAPK (Mitogen-Activated Protein Kinase); MDS (Multidimensional Scaling); MFF (Mature F1 Female); MFM (Mature F1 Male); MWF (Mature Wild Female); MWM (Mature Wild Male); PCA (Principal Component Analysis); PCR (Polymerase Chain Reaction); qPCR (Quantitative Polymerase
Chain Reaction); RRBS (Reduced Representation Bisulfite Sequencing); SD (Standard Deviation); TSS (Transcription Start Site); TTS (Transcription Termination Site).

## Appendix A. Tables

<table><caption>Table A1. Summary statistics for RRBS libraries and transcriptome sequencing data.</caption><thead><tr><th rowspan="2">Sample</th><th colspan="3">RRBS</th><th colspan="3">RNA</th></tr><tr><th>Raw Bases (Gb)</th><th>Q20 (%)</th><th>Q30 (%)</th><th>GC (%)</th><th>Raw Bases (Gb)</th><th>Q20 (%)</th><th>Q30 (%)</th><th>GC (%)</th></tr></thead><tbody><tr><td rowspan="3">Immature F1 Female</td><td>15.20</td><td>93.89</td><td>86.50</td><td>30.45</td><td>9.33</td><td>97.70</td><td>93.65</td><td>48.58</td></tr><tr><td>11.70</td><td>93.34</td><td>85.79</td><td>32.96</td><td>6.11</td><td>97.71</td><td>93.52</td><td>49.37</td></tr><tr><td>15.40</td><td>92.26</td><td>84.50</td><td>33.42</td><td>7.02</td><td>97.94</td><td>93.61</td><td>49.73</td></tr><tr><td rowspan="3">Immature F1 Male</td><td>18.70</td><td>93.90</td><td>86.29</td><td>30.61</td><td>7.08</td><td>97.40</td><td>92.86</td><td>47.62</td></tr><tr><td>11.70</td><td>93.05</td><td>85.79</td><td>32.42</td><td>7.17</td><td>97.36</td><td>92.73</td><td>47.20</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Immature Wild Female</td><td>15.30</td><td>93.49</td><td>86.17</td><td>33.69</td><td>7.00</td><td>97.80</td><td>93.30</td><td>49.22</td></tr><tr><td>7.80</td><td>94.21</td><td>86.93</td><td>29.57</td><td>8.02</td><td>97.71</td><td>93.50</td><td>49.28</td></tr><tr><td>13.10</td><td>92.26</td><td>84.47</td><td>34.07</td><td>8.85</td><td>97.78</td><td>93.37</td><td>48.94</td></tr><tr><td rowspan="3">Immature Wild Male</td><td>7.10</td><td>90.20</td><td>82.24</td><td>37.72</td><td>6.52</td><td>97.50</td><td>93.07</td><td>47.45</td></tr><tr><td>5.00</td><td>92.79</td><td>85.32</td><td>32.44</td><td>9.90</td><td>97.35</td><td>92.75</td><td>47.55</td></tr><tr><td>16.80</td><td>93.46</td><td>85.75</td><td>30.84</td><td>6.36</td><td>96.67</td><td>91.05</td><td>46.88</td></tr><tr><td rowspan="3">Mature F1 Female</td><td>16.80</td><td>93.76</td><td>86.15</td><td>30.43</td><td>7.46</td><td>97.77</td><td>93.48</td><td>49.29</td></tr><tr><td>11.90</td><td>92.37</td><td>84.77</td><td>36.16</td><td>7.92</td><td>97.72</td><td>93.11</td><td>46.11</td></tr><tr><td>13.00</td><td>92.45</td><td>84.78</td><td>31.68</td><td>8.17</td><td>97.43</td><td>92.62</td><td>48.36</td></tr><tr><td rowspan="3">Mature F1 Male</td><td>16.20</td><td>93.61</td><td>86.20</td><td>31.66</td><td>7.49</td><td>97.76</td><td>93.19</td><td>47.46</td></tr><tr><td>13.20</td><td>93.09</td><td>85.57</td><td>31.66</td><td>7.41</td><td>97.32</td><td>92.63</td><td>47.10</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Mature Wild Female</td><td>14.90</td><td>93.73</td><td>86.22</td><td>31.19</td><td>6.84</td><td>97.37</td><td>92.86</td><td>49.46</td></tr><tr><td>18.60</td><td>93.56</td><td>85.72</td><td>29.53</td><td>17.60</td><td>97.83</td><td>93.71</td><td>48.35</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Mature Wild Male</td><td>5.20</td><td>93.60</td><td>86.17</td><td>32.15</td><td>6.32</td><td>97.34</td><td>92.77</td><td>48.23</td></tr><tr><td>8.40</td><td>92.86</td><td>85.08</td><td>31.15</td><td>6.67</td><td>97.56</td><td>93.21</td><td>47.56</td></tr><tr><td>16.30</td><td>93.38</td><td>85.79</td><td>31.80</td><td>8.10</td><td>97.00</td><td>91.72</td><td>47.62</td></tr><tr><td></td><td>4.60</td><td>89.72</td><td>81.48</td><td>37.68</td><td>6.42</td><td>97.42</td><td>92.94</td><td>47.75</td></tr></tbody></table>

Notes: "F1" stands for specimens produced by aquaculture. Raw Bases (Gb) indicates the mean raw bases in gigabases per sample within groups. Q20 (%) is the proportion of sequencing reads that meet a specified quality threshold of Q20 or higher per sample within groups. Q30 (%) is the proportion of sequencing reads that meet a specified quality threshold of Q30 or higher per sample within groups. GC Content (%) is the percentage guanine-cytosine content in the sequences per sample within groups.

<table><caption>Table A2. Summary of methylation calling results.</caption><thead><tr><th rowspan="2">Sample</th><th rowspan="2">Cytosines<br>Analyzed</th><th colspan="3">Total Methylated Cs</th><th colspan="3">Total C to T Conversions</th><th colspan="3">C Methylated</th></tr><tr><th>CpG</th><th>CHG</th><th>CHH</th><th>CpG</th><th>CHG</th><th>CHH</th><th>CpG<br>Context</th><th>CHG<br>Context</th><th>CHH<br>Context</th></tr></thead><tbody><tr><td rowspan="3">Immature F1<br>Female</td><td>1,161,003,939</td><td>122,400,264</td><td>2,745,274</td><td>5,657,930</td><td>65,430,081</td><td>276,763,006</td><td>688,007,384</td><td>65.20%</td><td>1.00%</td><td>0.80%</td></tr><tr><td>965,408,146</td><td>101,336,168</td><td>2,344,612</td><td>4,820,567</td><td>55,187,119</td><td>231,465,675</td><td>570,254,005</td><td>64.70%</td><td>1.00%</td><td>0.80%</td></tr><tr><td>1,150,505,809</td><td>117,379,836</td><td>2,785,937</td><td>5,749,657</td><td>70,661,991</td><td>276,723,896</td><td>677,204,492</td><td>62.40%</td><td>1.00%</td><td>0.80%</td></tr><tr><td rowspan="3">Immature F1<br>Male</td><td>1,694,290,258</td><td>199,234,529</td><td>4,322,051</td><td>8,804,903</td><td>86,417,556</td><td>408,950,685</td><td>986,560,534</td><td>69.70%</td><td>1.00%</td><td>0.90%</td></tr><tr><td>959,664,596</td><td>115,841,877</td><td>2,463,936</td><td>4,840,611</td><td>48,808,618</td><td>232,999,065</td><td>554,710,489</td><td>70.40%</td><td>1.00%</td><td>0.90%</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Immature Wild<br>Female</td><td>1,226,494,406</td><td>128,743,569</td><td>2,969,643</td><td>6,089,415</td><td>71,282,746</td><td>294,872,229</td><td>722,536,804</td><td>64.40%</td><td>1.00%</td><td>0.80%</td></tr><tr><td>721,015,030</td><td>68,963,539</td><td>1,738,774</td><td>3,656,858</td><td>50,781,417</td><td>172,306,061</td><td>423,568,381</td><td>57.60%</td><td>1.00%</td><td>0.90%</td></tr><tr><td>1,045,075,951</td><td>110,853,128</td><td>2,576,630</td><td>5,263,851</td><td>58,901,311</td><td>251,483,299</td><td>615,997,732</td><td>65.30%</td><td>1.00%</td><td>0.80%</td></tr><tr><td rowspan="3">Immature Wild<br>Male</td><td>445,269,679</td><td>13,080,715</td><td>762,336</td><td>1,612,624</td><td>79,403,724</td><td>105,913,810</td><td>244,496,470</td><td>14.10%</td><td>0.70%</td><td>0.70%</td></tr><tr><td>393,961,163</td><td>15,535,314</td><td>825,105</td><td>1,711,055</td><td>66,333,121</td><td>94,728,423</td><td>214,828,145</td><td>19.00%</td><td>0.90%</td><td>0.80%</td></tr><tr><td>1,392,968,794</td><td>147,658,063</td><td>2,981,603</td><td>6,112,677</td><td>85,368,569</td><td>334,984,271</td><td>815,863,611</td><td>63.40%</td><td>0.90%</td><td>0.70%</td></tr><tr><td></td><td>2,494,517,598</td><td>275,885,506</td><td>6,538,261</td><td>13,028,122</td><td>159,719,419</td><td>600,354,809</td><td>1,438,991,481</td><td>63.30%</td><td>1.10%</td><td>0.90%</td></tr><tr><td rowspan="3">Mature F1<br>Female</td><td>1,744,257,380</td><td>195,337,172</td><td>4,477,104</td><td>9,050,681</td><td>94,610,073</td><td>421,568,792</td><td>1,019,213,558</td><td>67.40%</td><td>1.10%</td><td>0.90%</td></tr><tr><td>1,379,217,695</td><td>144,087,159</td><td>3,236,233</td><td>6,724,657</td><td>81,855,177</td><td>330,057,506</td><td>813,256,963</td><td>63.80%</td><td>1.00%</td><td>0.80%</td></tr><tr><td>883,508,418</td><td>96,889,871</td><td>2,147,202</td><td>4,356,367</td><td>48,382,792</td><td>213,161,844</td><td>518,570,342</td><td>66.70%</td><td>1.00%</td><td>0.80%</td></tr><tr><td></td><td>870,898,914</td><td>87,334,144</td><td>1,817,947</td><td>3,782,783</td><td>54,461,551</td><td>207,601,653</td><td>515,900,836</td><td>61.60%</td><td>0.90%</td><td>0.70%</td></tr></tbody></table>
<table><caption>**Table A2. Cont.**</caption><thead><tr><th rowspan="2">Sample</th><th rowspan="2">Cytosines<br>Analyzed</th><th colspan="3">Total Methylated Cs</th><th colspan="3">Total C to T Conversions</th><th colspan="3">C Methylated</th></tr><tr><th>CpG</th><th>CHG</th><th>CHH</th><th>CpG</th><th>CHG</th><th>CHH</th><th>CpG<br>Context</th><th>CHG<br>Context</th><th>CHH<br>Context</th></tr></thead><tbody><tr><td rowspan="2">Mature F1 Male</td><td>1,268,863,036</td><td>103,606,717</td><td>2,818,802</td><td>5,648,044</td><td>127,480,409</td><td>304,775,418</td><td>724,533,646</td><td>44.80%</td><td>0.90%</td><td>0.80%</td></tr><tr><td>937,562,250</td><td>99,915,441</td><td>2,223,903</td><td>4,271,133</td><td>64,136,781</td><td>226,344,652</td><td>540,670,340</td><td>60.90%</td><td>1.00%</td><td>0.80%</td></tr><tr><td rowspan="2">Mature Wild<br>Female</td><td>1,295,247,927</td><td>131,381,708</td><td>3,150,768</td><td>6,664,051</td><td>74,313,745</td><td>308,977,506</td><td>770,760,149</td><td>63.90%</td><td>1.00%</td><td>0.90%</td></tr><tr><td>1,570,474,527</td><td>157,303,960</td><td>3,295,521</td><td>7,050,272</td><td>93,007,036</td><td>372,067,275</td><td>937,750,463</td><td>62.80%</td><td>0.90%</td><td>0.70%</td></tr><tr><td rowspan="3">Mature Wild<br>Male</td><td>427,928,286</td><td>18,155,112</td><td>875,090</td><td>1,814,840</td><td>69,550,316</td><td>105,323,455</td><td>232,209,473</td><td>20.70%</td><td>0.80%</td><td>0.80%</td></tr><tr><td>684,891,178</td><td>24,904,912</td><td>1,321,040</td><td>2,757,633</td><td>113,021,414</td><td>166,103,962</td><td>376,782,217</td><td>18.10%</td><td>0.80%</td><td>0.70%</td></tr><tr><td>1,340,369,536</td><td>139,498,190</td><td>2,850,138</td><td>5,715,888</td><td>86,932,484</td><td>322,592,496</td><td>782,780,340</td><td>61.60%</td><td>0.90%</td><td>0.70%</td></tr><tr><td></td><td>272,155,711</td><td>8,534,487</td><td>468,304</td><td>979,457</td><td>47,311,224</td><td>65,128,833</td><td>149,733,406</td><td>15.30%</td><td>0.70%</td><td>0.60%</td></tr></tbody></table>

Notes: Summary of methylation calling results in *S. senegalensis* for eight study groups. “F1” stands for specimens produced by aquaculture. The reported cytosine counts are presented for the CpG, CHG, and CHH contexts (where H represents C, T, or A).

<table><caption>**Table A3.** Performed comparisons between groups in *Solea senegalensis* for differential methylation analysis.</caption><thead><tr><th>GROUP</th><th>IFM</th><th>IWF</th><th>IWM</th><th>MFF</th><th>MFM</th><th>MWF</th><th>MWM</th></tr></thead><tbody><tr><td>IFF</td><td>IFF/IFM</td><td>IFF/IWF</td><td></td><td>IFF/MFF</td><td></td><td></td><td></td></tr><tr><td>IFM</td><td></td><td></td><td>IFM/IWM</td><td></td><td>IFM/MFM</td><td></td><td></td></tr><tr><td>IWF</td><td></td><td></td><td>IWF/IWM</td><td></td><td></td><td>IWM/MWF</td><td></td></tr><tr><td>IWM</td><td></td><td></td><td></td><td></td><td></td><td></td><td>IWM/MWM</td></tr><tr><td>MFF</td><td></td><td></td><td></td><td></td><td>MFF/MFM</td><td>MFF/MWF</td><td></td></tr><tr><td>MFM</td><td></td><td></td><td></td><td></td><td></td><td></td><td>MFM/MWM</td></tr><tr><td>MWF</td><td></td><td></td><td></td><td></td><td></td><td></td><td>MWF/MWM</td></tr></tbody></table>

Notes: Group comparisons keys are shown as follows: MFM/MWM (Mature F1 Male vs. Mature Wild Male), IFM/IWM (Immature F1 Male vs. Immature Wild Male), IWM/MWM (Immature Wild Male vs. Mature Wild Male), IFM/MFM (Immature F1 Male vs. Mature F1 Male), MFF/MWF (Mature F1 Female vs. Mature Wild Female), IFF/IWF (Immature F1 Female vs. Immature Wild Female), IWF/MWF (Immature Wild Female vs. Mature Wild Female), IFF/MFF (Immature F1 Female vs. Mature F1 Female), MFF/MFM (Mature F1 Female vs. Mature F1 Male), IFF/IFM (Immature F1 Female vs. Immature F1 Male), MWF/MWM (Mature Wild Female vs. Mature Wild Male), and IWF/IWM (Immature Wild Female vs. Immature Wild Male). To maintain consistency and order in the comparisons, the following groups were established: F1 vs. Wild, Immature vs. Mature, and Females vs. Males. While the order of the groups does not affect the results, it is important to consider this order when interpreting the findings.

<table><caption>**Table A4.** Gene expression levels measured as gene counts for reproduction-related genes across the studied groups in *Solea senegalensis*.</caption><thead><tr><th>Gene</th><th>MWM</th><th>IWM</th><th>MFM</th><th>IFM</th><th>MWF</th><th>IWF</th><th>MFF</th><th>IFF</th><th>SD</th></tr></thead><tbody><tr><td>gsdf</td><td>20,585.25</td><td>22,900.00</td><td>28,574.00</td><td>32,425.00</td><td>57.50</td><td>56.33</td><td>49.50</td><td>50.00</td><td>14,371.75</td></tr><tr><td>sox19b</td><td>28.50</td><td>43.75</td><td>65.50</td><td>105.00</td><td>9985.50</td><td>7506.67</td><td>6473.50</td><td>7789.00</td><td>4320.61</td></tr><tr><td>fshr</td><td>9677.00</td><td>6814.75</td><td>6331.50</td><td>6229.00</td><td>7.50</td><td>9.67</td><td>5.25</td><td>11.00</td><td>4021.85</td></tr><tr><td>dmrt1</td><td>1537.25</td><td>1178.25</td><td>1220.50</td><td>1422.50</td><td>14.00</td><td>5.67</td><td>6.50</td><td>8.67</td><td>720.01</td></tr><tr><td>sox9</td><td>1513.75</td><td>1072.75</td><td>1294.00</td><td>1076.50</td><td>13.00</td><td>38.67</td><td>51.50</td><td>40.33</td><td>657.87</td></tr><tr><td>hsd17b</td><td>147.00</td><td>151.25</td><td>111.50</td><td>116.00</td><td>1338.50</td><td>1198.00</td><td>1197.25</td><td>1438.00</td><td>625.73</td></tr><tr><td>GATA6</td><td>1458.50</td><td>737.50</td><td>992.50</td><td>927.00</td><td>19.50</td><td>8.33</td><td>112.00</td><td>17.67</td><td>566.51</td></tr><tr><td>dmrt3a</td><td>648.25</td><td>461.00</td><td>853.00</td><td>472.50</td><td>208.50</td><td>131.00</td><td>146.50</td><td>133.67</td><td>271.87</td></tr><tr><td>gnrh</td><td>113.75</td><td>221.50</td><td>363.50</td><td>307.50</td><td>2.50</td><td>1.00</td><td>4.50</td><td>2.00</td><td>151.01</td></tr><tr><td>anth</td><td>233.75</td><td>252.25</td><td>310.00</td><td>401.00</td><td>148.50</td><td>41.67</td><td>48.00</td><td>56.00</td><td>134.40</td></tr><tr><td>cyp1b1</td><td>287.50</td><td>172.25</td><td>168.50</td><td>222.50</td><td>3.00</td><td>2.33</td><td>65.00</td><td>2.67</td><td>112.00</td></tr><tr><td>nanos1</td><td>351.25</td><td>167.00</td><td>220.50</td><td>211.50</td><td>73.50</td><td>37.00</td><td>151.50</td><td>31.00</td><td>108.29</td></tr><tr><td>dmrt2</td><td>113.75</td><td>137.75</td><td>348.00</td><td>229.00</td><td>178.50</td><td>41.00</td><td>62.50</td><td>37.67</td><td>106.60</td></tr><tr><td>sox3</td><td>168.75</td><td>260.25</td><td>198.00</td><td>175.00</td><td>42.50</td><td>32.67</td><td>22.75</td><td>25.67</td><td>94.87</td></tr><tr><td>sox8</td><td>250.00</td><td>32.50</td><td>35.00</td><td>46.50</td><td>32.50</td><td>10.33</td><td>40.50</td><td>12.33</td><td>78.82</td></tr><tr><td>sox6</td><td>30.50</td><td>23.25</td><td>10.00</td><td>18.50</td><td>4.50</td><td>3.67</td><td>6.00</td><td>1.67</td><td>10.58</td></tr></tbody></table>

Notes: Group keys are shown as follows: 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); Standard Deviation (SD).
Appendix B. Figures

<figure><img src="image_11.png" /><figcaption>**Figure A1.** Histological analysis of gonadal maturation stages in *Solea senegalensis*. (A) IWM: Immature Wild Male, (B) MWM: Mature Wild Male, (C) IWF: Immature Wild Female, (D) MWF: Mature Wild Female, (E) IFM: Immature F1 Male, (F) MFM: Mature F1 Male, (G) IFF: Immature F1 Female, and (H) MFF: Mature F1 Female. Notes: Meaning of abbreviations, spt: spermatozoa, st: seminiferous tubule, s: spermatogonia, str: stroma, o: oogonia (germ cell), pvo: pre-vitellogenic oocyte, vo: vitellogenic oocyte (cell in maturation phase), vd: vitello drops, and m: membrane (formed by the theca, granulosa cells, and pellucid membrane). Scale bar in (A–C,E–G) (10×): 50 µm; and scale bar in (D) and (H) (4×): 100 µm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A2.** Principal Component Analysis (PCA) plot based on gene expression profiles among *S. senegalensis* samples (generated using Novogene's Novomagic Platform). Notes: PC1 (Principal Component 1) and PC2 (Principal Component 2) show the amount of variability in the data in different directions, representing distinct gene expression profiles. The variations in expression profiles among *S. senegalensis* samples are shown for the eight study groups—mature wild males and females, immature wild males and females, mature F1 males and females, and immature F1 males and females. Groups are each depicted in distinct colors. Individual samples are represented by colored tags, with each corresponding to its respective group. The closeness indicates the degree of similarity among the *S. senegalensis* samples within each group.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure A3.** Differential expression Venn diagrams illustrating the differential genes of the compared groups in *Solea senegalensis*. Notes: The sum of all the numbers in the circle represents the total number in the compared groups, and the overlapping area indicates the number of differential genes shared between the groups in different comparison sets. Group comparison keys are shown as follows: MFM/MWM (Mature F1 Male vs. Mature Wild Male), IFM/IWM (Immature F1 Male vs. Immature Wild Male), IWM/MWM (Immature Wild Male vs. Mature Wild Male), IFM/MFM (Immature F1 Male vs. Mature F1 Male), MFF/MWF (Mature F1 Female vs. Mature Wild Female), IFF/IWF (Immature F1 Female vs. Immature Wild Female), IWF/MWF (Immature Wild Female vs. Mature Wild Female), IFF/MFF (Immature F1 Female vs. Mature F1 Female), MFF/MFM (Mature F1 Female vs. Mature F1 Male), IFF/IFM (Immature F1 Female vs. Immature F1 Male), MWF/MWM (Mature Wild Female vs. Mature Wild Male), and IWF/IWM (Immature Wild Female vs. Immature Wild Male).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A4.** KEGG enrichment analysis scatter plot for genes neighboring DMCpGs that exhibited a significant negative correlation between methylation and gene expression in male gonadal tissue in *S. senegalensis*. Notes: The enrichment analysis scatter plot was generated using Novogene's Novomagic platform. The Padj value, displayed on a color scale from 0 to 1, indicates the significance of enrichment. The gene ratio reflects the number of genes counted in relation to the total number of genes analyzed.</figcaption></figure>
## References

1. Luckenbach, J.A.; Borski, R.J.; Daniels, H.V.; Godwin, J. *Sex Determination in Flatfishes: Mechanisms and Environmental Influences*; Elsevier Ltd.: Amsterdam, The Netherlands, 2009. [CrossRef]

2. Holleley, C.E.; O'Meally, D.; Sarre, S.D.; Marshall Graves, J.A.; Ezaz, T.; Matsubara, K.; Azad, B.; Zhang, X.; Georges, A. Sex reversal triggers the rapid transition from genetic to temperature-dependent sex. *Nature* **2015**, *523*, 79–82. [CrossRef] [PubMed]

3. Kitano, J.; Ansai, S.; Takehana, Y.; Yamamoto, Y. Diversity and Convergence of Sex-Determination Mechanisms in Teleost Fish. *Guest* **2024**, *12*, 233–259. [CrossRef] [PubMed]

4. Nagahama, Y.; Chakraborty, T.; Paul-Prasanth, B.; Ohta, K.; Nakamura, M. Sex determination, gonadal sex differentiation, and plasticity in vertebrate species. *Am. Physiol. Soc.* **2021**, *101*, 1237–1308. [CrossRef] [PubMed]

5. Wang, J.; Tao, W.; Kocher, T.D.; Wang, D. Sex chromosome turnover and biodiversity in fishes. *J. Genet. Genom.* **2024**, *51*, 1351–1360. [CrossRef] [PubMed]

6. Morais, S.; Aragão, C.; Cabrita, E.; Conceição, L.E.; Constenla, M.; Costas, B.; Dias, J.; Duncan, N.; Engrola, S.; Estevez, A.; et al. *New Developments and Biological Insights into the Farming of Solea Senegalensis Reinforcing Its Aquaculture Potential*; Wiley-Blackwell: Hoboken, NJ, USA, 2016. [CrossRef]

7. APROMAR (Asociación Empresarial de Acuicultura de España). Aquaculture in Spain 2023. *APROMAR*. 2023. Available online: https://apromar.es/wp-content/uploads/2023/10/Aquaculture_in_Spain_2023_APROMAR.pdf (accessed on 1 September 2023).

8. Riesco, M.F.; Valcarce, D.G.; Martínez-Vázquez, J.M.; Martín, I.; Calderón-García, A.Á.; Gonzalez-Nunez, V.; Robles, V. Male reproductive dysfunction in *Solea senegalensis*: New insights into an unsolved question. *Reprod. Fertil. Dev.* **2019**, *31*, 1104–1115. [CrossRef] [PubMed]

9. Martín, I.; Carazo, I.; Rasines, I.; Rodríguez, C.; Fernández, R.; Martínez, P.; Norambuena, F.; Cherenguini, O.; Duncan, N. Reproductive performance of captive *Senegalese sole*, *Solea senegalensis*, according to the origin (Wild or cultured) and gender. *Span. J. Agric. Res.* **2019**, *17*, e0608. [CrossRef]

10. Fatsini, E.; González, W.; Ibarra-Zatarain, Z.; Napuchi, J.; Duncan, N.J. The presence of wild Senegalese sole breeders improves courtship and reproductive success in cultured conspecifics. *Aquaculture* **2020**, *519*, 734922. [CrossRef]

11. González-López, W.Á.; Ramos-Júdez, S.; Duncan, N.J. Reproductive behaviour and fertilized spawns in cultured Solea senegalensis broodstock co-housed with wild breeders during their juvenile stages. *Gen. Comp. Endocrinol.* **2024**, *354*, 14546. [CrossRef] [PubMed]

12. Liu, J.; Liu, X.; Jin, C.; Du, X.; He, Y.; Zhang, Q. Transcriptome profiling insights the feature of sex reversal induced by high temperature in tongue sole cynoglossus semilaevis. *Front. Genet.* **2019**, *10*, 522. [CrossRef] [PubMed]

13. Martínez, P.; Robledo, D.; Taboada, X.; Blanco, A.; Moser, M.; Maroso, F.; Hermida, M.; Gómez-Tato, A.; Álvarez-Blázquez, B.; Cabaleiro, S.; et al. A genome-wide association study, supported by a new chromosome-level genome assembly, suggests sox2 as a main driver of the undifferentiated ZZ/ZW sex determination of turbot (*Scophthalmus maximus*). *Genomics* **2021**, *113*, 1705–1718. [CrossRef]

14. Zou, C.; Wang, W.; Shu, C.; Liang, S.; Zou, Y.; Wang, L.; Wu, Z.; Liu, Y.; You, F. Expression characteristics of Hsd3b7 in the gonads of Paralichthys olivaceus. *Comp. Biochem. Physiol. B Biochem. Mol. Biol.* **2023**, *266*, 110848. [CrossRef] [PubMed]

15. Benzekri, H.; Armesto, P.; Cousin, X.; Rovira, M.; Crespo, D.; Merlo, M.A.; Mazurais, D.; Bautista, R.; Guerrero-Fernández, D.; Fernandez-Pozo, N.; et al. De novo assembly, characterization and functional annotation of Senegalese sole (Solea senegalensis) and common sole (Solea solea) transcriptomes: Integration in a database and design of a microarray. *BMC Genom.* **2014**, *15*, 1–8. [CrossRef]

16. Chauvigné, F.; Fatsini, E.; Duncan, N.; Ollé, J.; Zanuy, S.; Gómez, A.; Cerdà, J. Plasma levels of follicle-stimulating and luteinizing hormones during the reproductive cycle of wild and cultured Senegalese sole (Solea senegalensis). *Comp. Biochem. Physiol. A Mol. Integr. Physiol.* **2016**, *191*, 35–43. [CrossRef] [PubMed]

17. Merlo, M.A.; Portela-Bens, S.; Rodríguez, M.E.; García-Angulo, A.; Cross, I.; Arias-Pérez, A.; García, E.; Rebordinos, L. A comprehensive integrated genetic map of the complete karyotype of solea senegalensis (Kaup 1858). *Genes* **2021**, *12*, 49. [CrossRef] [PubMed]

18. Guerrero-Cózar, I.; Perez-Garcia, C.; Benzekri, H.; Sánchez, J.J.; Seoane, P.; Cruz, F.; Gut, M.; Zamorano, M.J.; Claros, M.G.; Manchado, M. Development of whole-genome multiplex assays and construction of an integrated genetic map using SSR markers in Senegalese sole. *Sci. Rep.* **2020**, *10*, 21905. [CrossRef]

19. de la Herrán, R.; Hermida, M.; Rubiolo, J.A.; Gómez-Garrido, J.; Cruz, F.; Robles, F.; Navajas-Pérez, R.; Blanco, A. A chromosome-level genome assembly enables the identification of the follicle stimulating hormone receptor as the master sex-determining gene in the flatfish *Solea senegalensis*. *Mol. Ecol. Resour.* **2023**, *23*, 886–904. [CrossRef]

20. Blanco-Vives, B.; Vera, L.M.; Ramos, J.; Bayarri, M.J.; Mananos, E.; Sánchez-Vázquez, F.J. Exposure of larvae to daily thermocycles affects gonad development, sex ratio and sexual steroids in Solea senegalensis. *J. Exp. Zool. Part A Ecol. Genet. Physiol.* **2011**, *315*, 162–169. [CrossRef] [PubMed]
21. Ramírez, D.; Rodríguez, M.E.; Mukiibi, R.; Peñaloza, C.; D'Cotta, H.; Robledo, D.; Rebordinos, L. Methylation profile of the testes of the flatfish *Solea senegalensis*. *Aquac. Rep.* **2024**, *39*, 102405. [CrossRef]

22. Woods, L.C., III; Li, Y.; Ding, Y.; Liu, J.; Reading, B.J.; Fuller, S.A.; Song, J. DNA methylation profiles correlated to striped bass sperm fertility. *BMC Genom.* **2018**, *19*, 224. [CrossRef]

23. Agulleiro, M.J.; Anguis, V.; Cañavate, J.P.; Martínez-Rodríguez, G.; Mylonas, C.C.; Cerdà, J. Induction of spawning of captive-reared Senegal sole (*Solea senegalensis*) using different administration methods for gonadotropin-releasing hormone agonist. *Aquaculture* **2006**, *257*, 511–524. [CrossRef]

24. Guzmán, J.M.; Ramos, J.; Mylonas, C.C.; Mañanós, E.L. Spawning performance and plasma levels of GnRHa and sex steroids in cultured female Senegalese sole (*Solea senegalensis*) treated with different GnRHa-delivery systems. *Aquaculture* **2009**, *291*, 200–209. [CrossRef]

25. Riesco, M.F.; Oliveira, C.; Soares, F.; Gavaia, P.J.; Dinis, M.T.; Cabrita, E. *Solea senegalensis* sperm cryopreservation: New insights on sperm quality. *PLoS ONE* **2017**, *12*, e0186542. [CrossRef] [PubMed]

26. Chauvigne, F.; Olle, J.; Gonzalez, W.; Duncan, N.; Gimenez, I.; Cerda, J. Toward developing recombinant onadotropin-based hormone therapies for increasing fertility in the flatfish Senegalese sole. *PLoS ONE* **2017**, *12*, e0174387. [CrossRef] [PubMed]

27. Chauvigné, F.; González, W.; Ramos, S.; Ducat, C.; Duncan, N.; Giménez, I.; Cerdà, J. Seasonal-and dose-dependent effects of recombinant gonadotropins on sperm production and quality in the flatfish *Solea senegalensis*. *Comp. Biochem. Physiol. A Mol. Integr. Physiol.* **2018**, *225*, 59–64. [CrossRef]

28. Firmino, J.; Carballo, C.; Armesto, P.; Campinho, M.A.; Power, D.M.; Manchado, M. Phylogeny, expression patterns and regulation of DNA Methyltransferases in early development of the flatfish, *Solea senegalensis*. *BMC Dev. Biol.* **2017**, *17*, 1–14. [CrossRef]

29. Law, P.P.; Holland, M.L. DNA Methylation at the Crossroads of Gene and Environment Interactions; Portland Press Ltd.: London, UK, 2019. [CrossRef]

30. Zhang, S.; Cheng, Y.; Věchtová, P.; Boryshpolets, S.; Shazada, N.E.; Alavi, S.M.; Cosson, J.; Linhart, O. Potential Implications of Sperm DNA Methylation Functional Properties in Aquaculture Management; John Wiley and Sons Inc.: Hoboken, NJ, USA, 2023. [CrossRef]

31. Qin, W.; Scicluna, B.P.; van der Poll, T. The Role of Host Cell DNA Methylation in the Immune Response to Bacterial Infection. *Front. Immunol.* **2021**, *12*, 696280. [CrossRef] [PubMed]

32. Salem, M.; Al-Tobasei, R.; Ali, A.; Kenney, B. Integrated Analyses of DNA Methylation and Gene Expression of Rainbow Trout Muscle under Variable Ploidy and Muscle Atrophy Conditions. *Genes* **2022**, *13*, 1151. [CrossRef] [PubMed]

33. Piferrer, F.; Anastasiadi, D.; Valdivieso, A.; Sánchez-Baizán, N.; Moraleda-Prados, J.; Ribas, L. The model of the conserved epigenetic regulation of sex. *Front. Genet.* **2019**, *10*, 857. [CrossRef] [PubMed]

34. Mukiibi, R.; Peñaloza, C.; Gutierrez, A.; Yáñez, J.M.; Houston, R.D.; Robledo, D. The impact of *Piscirickettsia salmonis* infection on genome-wide DNA methylation profile in Atlantic Salmon. *Genomics* **2022**, *114*, 110503. [CrossRef] [PubMed]

35. Smyth, G.K.; Chen, Y.; Pal, B.; Visvader, J.E. Differential methylation analysis of reduced representation bisulfite sequencing experiments using edgeR. *F1000Research* **2017**, *6*, 2055. [CrossRef]

36. Laing, L.V.; Viana, J.; Dempster, E.L.; Webster, T.U.; Van Aerle, R.; Mill, J.; Santos, E.M. Sex-specific transcription and DNA methylation profiles of reproductive and epigenetic associated genes in the gonads and livers of breeding zebrafish. *Comp. Biochem. Physiol. A Mol. Integr. Physiol.* **2018**, *222*, 16–25. [CrossRef] [PubMed]

37. Anaya, M.; Ramírez, D.; Arias-Pérez, A.; Rodríguez, M.E.; Robledo, D.; Rebordinos, L. Comparative transcriptomics of gonadal sex differentiation mechanisms in *Solea senegalensis*. (In preparation). Manuscript Draft-Powered by Editorial Manager® and ProduXion Manager® from Aries Systems Corporation.

38. Beemelmanns, A.; Ribas, L.; Anastasiadi, D.; Moraleda-Prados, J.; Zanuzzo, F.S.; Rise, M.L.; Gamperl, A.K. DNA Methylation Dynamics in Atlantic Salmon (*Salmo salar*) Challenged With High Temperature and Moderate Hypoxia. *Front. Mar. Sci.* **2021**, *7*, 604878. [CrossRef]

39. Kohli, G.; Hu, S.; Clelland, E.; Di Muccio, T.; Rothenstein, J.; Peng, C. Cloning of transforming growth factor-β1 (TGF-β1) and its type II receptor from zebrafish ovary and role of TGF-β1 in oocyte maturation. *Endocrinology* **2003**, *144*, 1931–1941. [CrossRef] [PubMed]

40. Zheng, W.; Xu, X.W.; Zechen, E.; Liu, Y.; Chen, S. Genome-wide identification of the MAPK gene family in turbot and its involvement in abiotic and biotic stress responses. *Front. Mar. Sci.* **2022**, *9*, 1005401. [CrossRef]

41. Blondeau-Bidet, E.; Banousse, G.; L'Honoré, T.; Farcy, E.; Cosseau, C.; Lorin-Nebel, C. The role of salinity on genome-wide DNA methylation dynamics in European sea bass gills. *Mol. Ecol.* **2023**, *32*, 5089–5109. [CrossRef] [PubMed]

42. Feng, S.; Cokus, S.J.; Zhang, X.; Chen, P.Y.; Bostick, M.; Goll, M.G.; Hetzel, J.; Jain, J.; Strauss, S.H.; Halpern, M.E.; et al. Conservation and divergence of methylation patterning in plants and animals. *Proc. Natl. Acad. Sci. USA* **2010**, *107*, 8689–8694. [CrossRef] [PubMed]

43. Jiang, N.; Wang, L.; Chen, J.; Wang, L.; Leach, L.; Luo, Z. Conserved and divergent patterns of DNA methylation in higher vertebrates. *Genome Biol. Evol.* **2014**, *6*, 2998–3014. [CrossRef]
44. Marina, E.K.; Aurélien, B.; Amin, S.; Delphine, L.; Catherine, L.; Audrey, L. Analysis of DNA methylation in rainbow trout spermatozoa: The strengths and limitations of RRBS. *bioRxiv* **2023**. [CrossRef]

45. Rodriguez Barreto, D.; Garcia de Leaniz, C.; Verspoor, E.; Sobolewska, H.; Coulson, M.; Consuegra, S. DNA Methylation Changes in the Sperm of Captive-Reared Fish: A Route to Epigenetic Introgression in Wild Populations. *Mol. Biol. Evol.* **2019**, *36*, 2205–2211.
[CrossRef] [PubMed]

46. Anastasiadi, D.; Esteve-Codina, A.; Piferrer, F. Consistent inverse correlation between DNA methylation of the first intron and gene expression across tissues and species. *Epigenetics Chromatin* **2018**, *11*, 1–17. [CrossRef]

47. Wei, L.; Yang, C.; Tao, W.; Wang, D. Genome-wide identification and transcriptome-based expression profiling of the Sox gene family in the nile tilapia (*Oreochromis niloticus*). *Int. J. Mol. Sci.* **2016**, *17*, 270. [CrossRef] [PubMed]

48. Aoki, Y.; Nakamura, S.; Ishikawa, Y.; Tanaka, M. Expression and syntenic analyses of four nanos genes in medaka. *Zool. Sci.* **2009**, *26*, 112–118. [CrossRef] [PubMed]

49. Zhu, W.; Wang, T.; Zhao, C.; Wang, D.; Zhang, X.; Zhang, H.; Chi, M.; Yin, S.; Jia, Y. Evolutionary conservation and divergence of Vasa, Dazl and *Nanos1* during embryogenesis and gametogenesis in dark sleeper (*Odontobutis potamophila*). *Gene* **2018**, *672*, 21–33.
[CrossRef] [PubMed]

50. Lai, F.; Zhou, Y.; Luo, X.; Fox, J.; King, M.L. *Nanos1* functions as a translational repressor in the Xenopus germline. *Mech. Dev.* **2011**, *128*, 153–163. [CrossRef] [PubMed]

51. Lai, F.; Singh, A.; King, M.L. Xenopus *nanos1* is required to prevent endoderm gene expression and apoptosis in primordial germ cells. *Development* **2012**, *139*, 1476–1486. [CrossRef]

52. Liu, J.; Zhang, W.; Sun, Y.; Wang, Z.; Zhang, Q.; Wang, X. Molecular characterization and expression profiles of GATA6 in tongue sole (*Cynoglossus semilaevis*). *Comp. Biochem. Physiol. B Biochem. Mol. Biol.* **2016**, *198*, 19–26. [CrossRef]

53. Van Nes, S.; Andersen, Ø. Temperature effects on sex determination and ontogenetic gene expression of the aromatases cyp19a and cyp19b, and the estrogen receptors esr1 and esr2 in Atlantic Halibut (*Hippoglossus Hippoglossus*). *Mol. Reprod. Dev.* **2006**, *73*, 1481–1490. [CrossRef]

54. Portela-Bens, S.; Merlo, M.A.; Rodriguez, M.E.; Cross, I.; Manchado, M.; Kosyakova, N.; Liehr, T.; Rebordinos, L. Integrated gene mapping and synteny studies give insights into the evolution of a sex proto-chromosome in *Solea senegalensis*. *Chromosoma* **2017**, *126*, 261–277. [CrossRef]

55. Cross, I.; García, E.; Rodríguez, M.E.; Arias-Pérez, A.; Portela-Bens, S.; Merlo, M.A.; Rebordinos, L. The genomic structure of the highly conserved dmrt1 gene in Solea senegalensis (Kaup, 1868) shows an unexpected intragenic duplication. *PLoS ONE* **2020**, *15*, e0241518. [CrossRef]

56. Zhu, Y.; Cui, Z.; Yang, Y.; Xu, W.; Shao, C.; Fu, X.; Li, Y.; Chen, S. Expression analysis and characterization of dmrt2 in Chinese tongue sole (*Cynoglossus semilaevis*). *Theriogenology* **2019**, *138*, 1–8. [CrossRef]

57. Liu, Z.; Zhou, T.; Gao, D. Genetic and epigenetic regulation of growth, reproduction, disease resistance and stress responses in aquaculture. *Front. Genet.* **2022**, *13*, 994471. [CrossRef] [PubMed]

58. Navarro-Martín, L.; Galay-Burgos, M.; Piferrer, F.; Sweeney, G. Characterisation and expression during sex differentiation of Sox19 from the sea bass *Dicentrarchus labrax*. *Comp. Biochem. Physiol. B Biochem. Mol. Biol.* **2012**, *163*, 316–323. [CrossRef]

59. Hsu, C.W.; Chung, B.C. Evolution, expression, and function of gonadal somatic cell-derived factor. *Front. Cell Dev. Biol.* **2021**, *9*, 684352. [CrossRef] [PubMed]

60. Sambroni, E.; Lareyre, J.J.; Le Gac, F. Fsh Controls Gene Expression in Fish both Independently of and through Steroid Mediation. *PLoS ONE* **2013**, *8*, e76684. [CrossRef] [PubMed]

61. Li, L.; Lin, W.; Wang, Z.; Huang, R.; Xia, H.; Li, Z.; Deng, J.; Ye, T.; Huang, Y.; Yang, Y. Hormone Regulation in Testicular Development and Function. *Int. J. Mol. Sci.* **2024**, *25*, 5805. [CrossRef] [PubMed]

62. Gur, G.; Bonfil, D.; Safarian, H.; Naor, Z.; Yaron, Z. GnRH receptor signaling in tilapia pituitary cells: Role of mitogen-activated protein kinase (MAPK). *Comp. Biochem. Physiol. B Biochem. Mol. Biol.* **2001**, *129*, 517–524. [CrossRef]

63. Durán-Pastén, M.L.; Fiordelisio, T. GnRH-induced Ca2+ signaling patterns and gonadotropin secretion in pituitary gonadotrophs. Functional adaptations to both ordinary and extraordinary physiological demands. *Front. Endocrinol.* **2013**, *4*, 127. [CrossRef]

64. Wei, X.; Zhang, Y.; Li, C.; Ai, K.; Li, K.; Li, H.; Yang, J. The evolutionarily conserved MAPK/Erk signaling promotes ancestral T-cell immunity in fish via c-Myc-mediated glycolysis. *J. Biol. Chem.* **2020**, *295*, 3000–3016. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.