Epigenetic regulation of gene expression is involved in the sexual development of gonochoristic fish with different types of sex-determining mechanisms, as well as in driving the process of sex change in different types of hermaphrodites.

Here, we searched the published literature in fish and collected information on the DNA methylation of genes related to sexual development. A WGBS was used in the half-smooth tongue sole, (Shao et al., 2014), while a MBS was used in the European sea bass (Anastasiadi et al., 2018b). These are the exceptions because in the rest of studies carried so far, which concern around 15 different species, just one or two genes have been analyzed in each case ( Table 1 ). DNA methylation at a single CpG is of a binary nature, since a given CpG can be either methylated or unmethylated. However, mean percent DNA methylation can, theoretically, fall in any value between 0 and 100%. This applies regardless of whether one considers the promoter or the first intron (Anastasiadi et al., 2018a) or other genomic features in a predefined window of a given length. Information drawn from the primary literature shows that DNA methylation levels are more or less evenly distributed across five arbitrarily defined methylation classes (0–20%, 21–40%, 41–60%, 61–80%, and 81–100%), perhaps with a higher preponderance in the 0–20% class, regardless of other considerations such as method of analysis, targeted genomic feature, sex, species, etc. Thus, these preliminary data indicate that there are no preferred or typical DNA methylation values for the sex-related genes as a whole ( Figure 3A ). Again, it should be remembered that DNA methylation values represent the combined values resulting from the different cell types making up the gonads. Thus, the correlation with gene expression, if present, should take this into account.

Gonadal aromatase (cyp19a1a) was the first gene shown to be under epigenetic regulation during sexual development in a vertebrate, the European sea bass (Navarro-Martín et al., 2011). This is not surprising because it is the only steroidogenic enzyme responsible for the balance between androgens and estrogens and because estrogens are needed for ovarian differentiation in all nonmammalian vertebrates (Guiguen et al., 2010). Since then, the DNA methylation of only few genes has been studied in more than two species: cyp19a1a just cited earlier (studied in 10 species), doublesex- and mab-3-related transcription factor 1 (dmrt1) (6 species), anti-Müllerian hormone or Müllerian-inhibiting hormone (amh) (4 species), and the member of the winged helix/forkhead group (foxl2a) (3 species) ( Figure 3B ). From the analysis of the published data and our own unpublished data, we found that mean DNA methylation levels of cyp19a1a were typically <50% in ovaries (mean: 46.2%, sd: 15.98) and >75% in testes (mean: 77.0%, sd: 24.89) (t-test: -4.0439; df = 28, p = 0.00037) in a fairly consistent manner across species (see list of species in Table 1 ). This finding was in accordance with the constitutive higher expression of cyp19a1a in ovaries when compared with testes (Piferrer and Blázquez, 2005; Guiguen et al., 2010). Likewise, mean DNA methylation levels of dmrt1 were ∼30% in ovaries (mean: 32.54%, sd: 15.98) and <10% in testes (mean: 5.54%, sd: 5.21) (t-test: 4.54; df = 14, p = 0.00046), also in accordance with the higher constitutive expression of dmrt1 in testes when compared with ovaries (Herpin and Schartl, 2011) ( Figure 3B ). Therefore, these two important genes for sex differentiation, which have been used as sex markers in some fish species, e.g., turbot, Scophthalmus maximus (Ribas et al., 2016), and medaka, Oryzias latipes (Herpin and Schartl, 2011), do indeed conform to the CERS predictions, since there is an inverse relationship between DNA methylation and gene expression with clear sex-specific differences.

This inverse relationship does not seem apparent when two other well-known genes with sex-biased expression in fish are considered: amh and foxl2a ( Figure 3B ). Amh is a member of the TGF-β superfamily of growth and differentiation factors involved in sex differentiation from mammals to fish (Piferrer and Guiguen, 2008). Relatively low and equal levels of amh expression are detected in gonads prior to the appearance of sex-specific differences. However, once sex differentiation is underway, higher amh levels are typically associated with testis differentiation in several species analyzed (reviewed in Pfennig et al., 2015). Here, we found that mean DNA methylation levels of amh were 54.05% (sd: 25.26) in ovaries and 80.24% (sd: 12.74) in testes, a difference that did not reach statistical significance with the data available so far (t-test: -2.07; df = 8, p = 0.07211). In the same way, foxl2a is expressed at higher levels in the ovary when compared with the testis (reviewed in Bertho et al., 2016), like cyp19a1a. On the other hand, foxl2a is actually one of the earliest transcriptional activators of cyp19a1a that co-localizes in the granulosa cells (Wang et al., 2004). However, DNA methylation levels were clearly not different (mean = 3.08% and sd = 3.88 in ovaries and mean = 2.59% and sd = 3.3 in testes) (t-test: 0.1656; df = 4, p = 0.8765). Therefore, unlike cyp19a1a and dmrt1, and with the information available so far, data suggest that amh and foxl2a do not seem to conform to CERS predictions or that, in these genes, the relationship between DNA methylation and gene expression is positive ( Figure 3B ), although, clearly, further research is needed.

There are other genes related to sex differentiation at different degrees for which it may be premature to attempt any sort of generalizations. These genes include amhr2, cyp11a, hsd3b2, nr3c1, sox9, vasa, and gsdf ( Figure 3C ). Here, it is worth noting that allelic diversification of amhr2 in Takigugu rubripes results in a dominant master sex-determining gene, while allelic diversification of gsdf has given rise to the sex-determining gene in some fish species, including Oryzias luzonensis and Anoplopoma fimbria (reviewed in Piferrer, 2018; Guiguen et al., 2019). DNA methylation levels of amhr2 in the European sea bass were ∼50% without sex-related differences (Anastasiadi et al., 2018b), while in the half-smooth tongue sole DNA methylation levels are higher in females (Shao et al., 2014). Similarly, in the latter species, the only species where gsdf DNA methylation values have been determined, these values are clearly lower in males, in accordance with the higher expression of gsdf in males (Shao et al., 2014) ( Figure 3C ).

Except in the half-smooth tongue sole (Shao et al., 2014), where WGBS was used, in the rest of the studies reported in Table 1 and used to draw Figure 3 , targeted approaches were utilized to query the DNA methylation status of the target genes. For these studies, an average of ∼9 late juvenile or adult fish per sex was used. Typically, amplicons spawn ∼450 bp and usually include ∼15 CpG located around the transcription start site, although the latter figure may vary considerably among species. It is interesting to note that while sex-specific differences involve change in DNA methylation of several CpGs in some genes, in contrast, in other genes, sex-differences involve only a low number of CpGs ( Figure 4 ). A more comprehensive picture will emerge when genome-wide DNA methylation techniques such as WGBS or RRBS will be employed in lieu of the targeted approaches used so far in most studies.

For the rest of the genes, cyp11a, hsd3b2, and nr3c1, there are only preliminary data gathered in our lab with the European sea bass (Anastasiadi et al., unpublished) and zebrafish (Danio rerio) (Valdivieso et al., unpublished; Moraleda-Prados et al., unpublished). In the European sea bass, methylation values of hsd3b2 are higher in females. This is in agreement with the expression of this gene that is male-skewed in the developing gonads of Nile tilapia, Oreochromis niloticus (Ijiri et al., 2008). On the other hand, DNA methylation values of nr3c1 and sox9 were quite different between the two species.

In many species, sex determination has an environmental component. Hence, it is worth mentioning that an environmental factor such as temperature or population density may be connected to sex through epigenetic mechanisms. DNA methylation changes in sex-related genes is the type of epigenetic modification most commonly studied so far. Temperature can affect DNA methylation of many genes, as shown by MeDIP-seq in the Nile tilapia (Sun et al., 2016), although the exact mechanism is not known yet. In the European sea bass, elevated temperature induces hypermethylation in the promoter of cyp19a1a, and this prevents the binding of cyp19a1a transcriptional activators such as sf1 and foxl2a (Navarro-Martín et al., 2011). Other epigenetic modifications can also be involved in the connection between environmental factors and sex. Thus, temperature increases the transcription of lysine-specific demethylase 6B (kdm6b), a chromatin modifier gene in the red-eared slider turtle, Trachemys scripta. Kdm6b eliminates the trimethylation of H3K27 in the promoter of dmrt1, leading to upregulation of its expression and male development (Ge et al., 2018; Georges and Holleley, 2018).

We would like to mention three considerations for further testing the CERS model. First, what species are worth testing? Obviously, fish, due to their great diversity of sexual systems and sex determining systems, which can vary even in closely related species. Reptiles can also provide very relevant information. Many reptiles possess temperature-dependent sex determination and thus offer the opportunity to test whether DNA methylation in key genes do correlate with gene expression and phenotypic sex under different incubation temperatures during the thermosensitive period. Thus, in the red-eared slider turtle cyp19a1 DNA methylation levels conformed to CERS predictions (Matsumoto et al., 2013). The same is true in the alligator, Alligator mississippiensis, for cyp19a1 and sox9 (Parrott et al., 2014) and in the sea turtle, Lepidochelys olivacea, for sox9 (Venegas et al., 2016). In birds and mammals, sexual development is strongly canalized (Capel, 2017), and therefore there is little or no room for sexual plasticity. Nevertheless, in such canalized systems, it would also be interesting to determine to what extent DNA methylation of key genes correlates with expression and whether this is established before the completion of gonadal differentiation. In any case, and regardless of the species of choice, testing the role of epigenetic regulation on the expression of key sex-related genes during the process of sex differentiation should involve, in our opinion, the analysis of at least three different time points. The first one, ideally, should be prior to any morphological sign of sex differentiation, the second around the middle of the process, and the third towards the end or after the completion of sex differentiation.

Second, what other genes can be targeted? In our view, the genes to be tested should include at least the ones that consistently follow or not the predictions of the CERS model, namely, cyp19a1a, dmrt1, amh, and foxl2a, as shown in this paper. However, other genes with known functions in sexual development in vertebrates, including mammals, birds, and reptiles, should also be studied. We propose here a list of some of the most relevant genes found in the literature ( Table 2 ). Information on DNA methylation of additional genes during gonadal differentiation and any possible sex-related differences will help to better understand the epigenetic regulation of sexual development.

Third, what other approaches can be used? Gene-editing techniques such as CRISPR/Cas9 or the more recently developed technique to edit the methylome in the mammalian genome by Liu et al. (2016b) can be very useful. To date, knockout mutants of sex-related genes in fish have been mostly developed for some model species, e.g., in zebrafish: cyp19a1a (Lau et al., 2016), amh, and dmrt1 (Lin et al., 2017), and in medaka: estrogen receptor 1 (esr1) (Tohyama et al., 2017), gnrh family genes (Marvel et al., 2018), and cyp19a1a knockout (Nakamoto et al., 2018). Lau et al. (2016) found that all knockout mutants of cyp19a1a were males, supporting the view that aromatase plays an essential role in ovarian differentiation and development. Yet, Lin et al. (2017) found that dmrt1 and amh knockout zebrafish mutants displayed female-biased sex ratios, but the development of abnormal testes was still possible. Dmrt1 was suggested to be necessary for the maintenance, self-renewal, and differentiation of male germ cells, and amh was proposed to control the balance between proliferation and differentiation of these cells. Therefore, it would be interesting to analyze the DNA methylation of dmrt1 and other male-biased genes in amh knockout mutants and vice versa, the DNA methylation of amh and other male-biased genes of the network in dmrt1 knockout mutants.

There are some aspects worth discussing regarding future studies of the involvement of DNA methylation on the regulation of sexual development. First, one aspect concerns the genomic feature on which one should focus when the goal is to associate DNA methylation with gene expression levels. Determination of the expression should be accurate and consistent for each gene assuming that the method of measurement, e.g., qPCR, is properly employed according to the appropriate standards (primers, reference genes, etc.). In contrast, DNA methylation levels can vary across different genomic features of the same gene. In most studies, the promoter region has been typically targeted. However, methylation of other genomic regions has been found to be equally or even better associated with gene silencing. Indeed, it was shown that the first exon is tightly linked to transcriptional silencing (Brenet et al., 2011). Furthermore, in a systematic study aimed at addressing this question, it was found that the first intron, more than the promoter and the first exon, is tightly related to gene silencing. This seems to be conserved across vertebrate species since it was observed in fish (Japanese puffer, Takifugu rubipres, and the European sea bass), frog (Xenopus), and humans (Anastasiadi et al., 2018a). Thus, for the epigenetic regulation of sex, as well as for sex-related development of biomarkers, it is better to focus around the transcription start site and to prioritize the CpGs localized in the first intron, first exon, and promoter regions, in the order mentioned. Furthermore, gene expression can also be positively correlated to tissue-specific DNA methylation, and this should be kept in mind (Lokk et al., 2014; Wan et al., 2015; Anastasiadi et al., 2018a).

Another aspect concerns the possible effect of genetic variation on DNA methylation levels and how to account for it in the data analysis (Lea et al., 2017; Anastasiadi et al., 2018b). This is related to the number of samples to be analyzed per treatment in studies of DNA methylation, which has been discussed elsewhere (Bock et al., 2016). Also, it would be desirable to overcome the noise induced by the cell heterogeneity of the gonadal tissue. In this regard, recent technological advances allow to determine the epigenome of single cells (Farlik et al., 2015). Efforts toward such type of measurements would definitively help in obtaining more robust measurements of DNA methylation.

Furthermore, DNA methylation and gene expression levels discussed throughout this paper refer to the gonads. DNA methylation is known to be tissue-specific. However, it cannot be ruled out that the methylation patterns of the gonads could be replicated in other tissues. This could be the case of tissues involved in the control of reproduction (e.g., the hypothalamus) or that present sex dimorphism (e.g., secondary sexual characters) because they are under the control of hormonal steroids. To the best of our knowledge, this information does not still exist despite increasing evidence of sex-related differences in DNA methylation for many genes in nonreproductive tissues, such as the muscle or the liver (Davegårdh et al., 2019; Grimm et al., 2019).

Another major challenge will be to determine the sexual phenotype just by the DNA methylation levels of selected EEMs before it can be determined by other means (e.g., by analyzing transcriptomic or histological changes). This would be achievable if demonstrated that the epigenetic modifications precede changes in gene expression. In this case, the EEMs and the CERS model can be foreseen as having potentially useful applications. For example, a defined set of EEMs could be used to predict the sexual phenotype in species with marked sexual growth dimorphism (Parker, 1992; Wang et al., 2019). EEMs could allow to predict the sex ratio in a subsample of a clutch before gonadal differentiation. This would aid in the stock management and in the selection of future broodstock. The same principle could be applied in ornamental fish culture, where the secondary sexual characteristics of males make them usually more desirable than females (Piferrer and Lim, 1997). Another case would be to aid in selection of broodstock fish with a certain epigenetic profile that is suitable to withstand, for example, a masculinization environment due to elevated density or temperature. In reptiles, the use of EEMs combined with temperature manipulations could aid in the research toward our understanding of the underlying molecular mechanism of temperature-dependent sex determination.

Finally, epigenetic modifications can recapitulate past environmental influences (Turner, 2009; Vogt, 2017). Taking advantage of this, EEMs could help to determine whether animals in the wild were exposed to altered environmental conditions such as, for example, exposure to pollutants or elevated temperatures. These EEMs could therefore be useful in conservation programs aimed at determining the environmental hazards to which natural populations may have been previously exposed. As an example along these lines, Guillette et al. (2016) identified epigenetic biomarkers to assess the environmental exposures and health impacts on populations of alligators from lakes contaminated with endocrine-disrupting compounds. The effects of endocrine-disrupting compounds on DNA methylation in the field of aquatic toxicology and biodiversity conservation have recently been reviewed by Tubbs and McDonough, (2018). This approach would allow determining whether a wild population was subjected to a sex-altering condition in the past. To our knowledge, this type of applications has not been fully developed to date, but several studies have started to identify biomarkers with this aim in mind.

All datasets analyzed for this study are included in the manuscript and the supplementary files.

FP conceived the study, supervised the collection of data, coined the concepts of EEM and CERS, and wrote the paper. DA developed the technique of MBS and wrote the paper. AV, NS-B, JM-P, and LR collected data on different species and wrote the paper. All authors approved the final version of the manuscript.

This study was supported by the Spanish Ministry of Science grants AGL2016–787107-R “Epimark” to FP and AGL2015-73864-JIN “Ambisex” to LR. DA was supported by an Epimark contract, AV and NS-B were supported by Spanish government scholarships (BES-2014-069051 and BES-2017-079744, respectively); LR and JM-P were supported by Ambisex contracts.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.