Rearing density directly impacts fish welfare, which, in turn, affects productivity in aquaculture. Previous studies have indicated that high-density rearing during sexual development in fish can induce stress, resulting in a tendency towards male-biased sex ratios in the populations. In recent years, research has defined the relevance of the interactions between the environment and epigenetics playing a key role in the final phenotype. However, the underlying epigenetic mechanisms of individuals exposed to confinement remain elucidated. By using zebrafish (Danio rerio) the DNA methylation promotor region and the gene expression patterns of six genes: dnmt1, cyp19a1a, dmrt1, cyp11c1, hsd17b1, and hsd11b2 involved in the DNA maintenance methylation, reproduction, and stress were assessed. Zebrafish larvae were subjected to two high-density conditions (9 and 66 fish/L) during two periods of overlapping sex differentiation of this species (7 to 18 and 18 to 45 days post-fertilization, dpf). Results showed a significant masculinization in the populations of fish subjected to high densities from 18 to 45 dpf. In adulthood the dnmt1 gene was differentially hypomethylated in ovaries and its expression was significantly downregulated in the testes of fish exposed to high-density. Further, the cyp19a1a gene showed downregulation of gene expression in the ovaries of fish subjected to elevated density, as previously observed in other studies. We proposed dnmt1 as a potential testicular epimarker, and the expression of ovarian cyp19a1a as a potential biomarker for predicting stress originated from high densities during the early stages of development. These findings highlight the importance of rearing densities by long-lasting effects in adulthood conveying cautions for stocking protocols in fish hatcheries.

Keywords:

Subject: Biology and Life Sciences - Aquatic Science

1. Introduction

Rearing density holds significant importance in aquaculture as it directly impacts fish welfare, thereby influencing overall profitability. Inappropriate stocking densities can induce stress, leading to physiological disruptions that pose a threat to fish physiology (Stevens et al., 2017; Sarma et al., 2023). Persistent exposure to chronic stressors, such as high population density, can increase susceptibility to diseases, deplete energy resources, decrease growth performance, muscle and bones quality and antioxidative capacity, and, ultimately, reduce overall performance (Besson et al., 2014; Onxayvieng et al., 2021). The 'General Adaptation Syndrome' is the widely accepted concept of stress and consists of three phases (Selye, 1946; Tort et al., 1998). Firstly, there is a physiological state of 'alarm,' characterized by the production of adrenaline and cortisol. Secondly, if the stressor persists, the organism attempts to adapt and defend itself, entering a state of resistance. Finally, if the stress continues beyond the organism's coping abilities, it enters an exhaustion phase, which may lead to irreversible dysfunction or even death (Costa & Pinto, 2017).

The stress response is connected with the reproductive system by the brain–pituitary–gonadal axis in which neuro-endocrine interactions control multi-directionally the fish physiology (Besedovsky & del Rey, 1996; Weltzien et al., 2004; Acevedo-Rodriguez et al., 2018). In many fish species, the final sex of an individual is determined by genetic and environmental factors (Baroiller et al., 2009), therefore, stress needs to be into account when studying sexual development. During sex differentiation, environmental factors have the potential to influence and modify sexual phenotypes. Therefore, genes together with diverse factors in the habitat where they live are responsible to define the fate of a differentiating gonads towards an ovary or a testis. In many fish species, most well-documented abiotic factor is temperature, responsible to masculinize fish populations (Ospina-Alvarez & Piferrer, 2008; Edmands, 2021). In contrast, a limited number of studies have investigated the impact of rearing density on the final sexual phenotype. For instance, this is the case with paradise fish, Macropodus opercularus (Francis, 1984), some coral reef fish species (e.g., Centropyge potteri, Labroides dimidiatus) (Lutnesky, 1994; Kuwamura et al., 2014), the European eel, Anguilla anguilla (Roncarati et al., 1997; Krueger & Oliveira, 1999; Huertas & Cerdà, 2006), the European sea bass, Dicentrarchus labrax (Saillant et al., 2003), and zebrafish, Danio rerio (Hazlerigg et al., 2012; Ribas et al., 2017a).

Zebrafish has become a model organism in many research areas from biomedical, developmental biology, toxicology to aquaculture-related research (Dahm & Geisler, 2006; Ribas & Piferrer, 2014; Piferrer & Ribas, 2020). The advantages of zebrafish as a fish model make it an appropriate species for studying reproduction-related problems observed in fish farming. Zebrafish is a gonochoristic species of an undifferentiated type (Piferrer & Guiguen, 2008; Orban et al., 2009). The initial indication of gonadal differentiation in zebrafish occurs approximately 10 days post-fertilization (dpf), triggering the onset of a critical period during which the gonadal fate can be influenced by environmental cues. While the timing and duration of this process may vary among individuals, it is typically fully accomplished by around 50 dpf (Chen & Ge, 2012, 2013; Lawrence et al., 2012).

Wild zebrafish strains possess a chromosomal sex determination system, while some laboratory strains are characterized by a polygenic sex-determining system, both of them, environmentally influenced (Liew et al 2012; Wilson et al 2014; Valdivieso et al., 2022b). Environmental factors able to modulate sexual phenotype and consequently skew sex ratios in zebrafish are temperature (Uchida et al., 2004; Abozaid et al., 2012; Ribas et al., 2017a), hypoxia (Shang et al., 2006), nutritional resources (Lawrence et al., 2008) infection (Pradhan et al., 2012; Moraleda-Prados et al., 2021) and density (Ribas et al., 2017b, c). However, it is important to note that the effects of masculinization by density rearing were determined in experiments encompassing the entire larval development (6 to 90 dpf) (Ribas et al., 2017b, c), rather than specific windows during sex differentiation (15 to 45 dpf) of zebrafish. Thus, in order to help establish a better accurate rearing protocol to avoid masculinization by density level, more information regarding the sensible window during sex differentiation is required.

For more than a decade, aquaculture-related researchers have manifested the role of epigenetics in relation to environmental conditions to better understand fish physiology and consequently, a way to improve aquaculture production. Thus, searching for informative epimarkers is at the frontier of aquaculture-related research. The effects of temperature in altering DNA methylation levels in some canonical sex-related genes associated to the masculinization of the population have been shown in many farming fish species: in E. sea bass (Navarro-Martin et al., 2011), Nile tilapia, Oreochromis niloticus (Sun et al., 2016; Wang et al., 2017a), turbot, Scophthalmus maximus (Wang et al., 2017b), half-smooth tongue sole, Cynoglossus semilaevis (Shao et al., 2014) or Japanese flounder Paralichthys olivaceus (Wen et al., 2014).

In zebrafish, high temperature during sex differentiation revealed changes in methylation levels on genes involved in reproduction and stress and was associated with the masculinization of the populations (Valdivieso et al., 2022a). Based on these differential changes, specific epigenetic markers were used to predict the sex and the past-thermal events that occurred during the early stages of gonad development (Valdivieso et al., 2022a). Moreover, it was observed hypomethylation of the testicular epigenome in the unexposed first generation offspring derived from the heat-exposed parents (Valdivieso et al., 2020). These results confirmed that environmental disturbances may have an epigenetic memory in which molecular biomarkers can be valuable toolkits for better understanding the environmental cues that alter phenotypes during early development.

In this line, the present study consisted of identifying putative epigenetic biomarkers able to predict past-stress events by rearing density. To do this, we used a locus-specific approach to evaluate the DNA methylation of genes involved in: 1) epigenetics: DNA methyltransferase 1 (dnmt1) the enzyme responsible for re-establishing the methylation landscape in the cell (Haggerty et al., 2021); 2) reproduction: cytochrome P450, family 19, subfamily A, polypeptide 1a (gonadal aromatase, cyp19a1a), the enzyme that converts androgens to estrogens and, doublesex and mab-3 related transcription factor 1 (dmrt1), key regulator of sex determination of testis development (Herpin & Schartl, 2011); and 3) stress: cytochrome P450 family 11 subfamily b member 1 (cyp11c1), enzyme mediator in the conversion of testosterone to cortisol (Zhang et al., 2020), hydroxysteroid (17-beta) dehydrogenase 1 (hsd17b1), enzyme involved in estrogen production (Filby et al., 2010), and hydroxysteroid dehydrogenase type 2 (hsd11b2), which converts cortisol into its inactive form (Theodoridi et al., 2021). To complete the study, the expressions of these six genes were analyzed to determine their relationship with DNA methylation. In addition, in order to shed light on the density effects on the rearing protocols in the zebrafish facilities, the identification of the sensible window that skews sex ratios towards males during sex differentiation was assessed.

2. Results

2.3. Effects of elevated density on sex ratio and growth

Results showed that the effect of the density was significantly dependent on the exposure period in which larvae was subjected. Comparing the effect of density levels (9 and 66 fish/L) in the 7–18 dpf group, we did not find differences in proportion of males. However, the impact of the 66 fish/L density level was significant when the treatments were performed during 18– 45 dpf (Figure 1A, and Table S3). The overall percent of males from the six families increased from (mean ± SD) 55.22 ± 9.10% at 9 fish/L to 67.54 ± 14.62 at 66 fish/L treatment in the 18–45 dpf group (Figure 1B and Table S3). However, when analyzed the effects of masculinization by density at individual family, only Family#4 and #6 showed a significant (χ² = 6.86 and 3.96, both P < 0.05, respectively) increase of males (Figure 1B), confirming previous results of family-specifics response.

Because we did not have equal representation of all the families tested in control and 18–45 dpf group, and the latter only in two of six families the effects of masculinization by density were substantially evident, we evaluated whether the masculinization effects of elevated density through a GLMM. Results showed that density had a significant effect on the masculinization rate (P < 0.001) (Table 1). Based on the prediction fit of the GLMM analysis, we plotted the expected masculinization in the populations of zebrafish according to density level (Figure 1C).

Growth was inversely related to rearing density, with sex-related differences (Figure 2 and Table S4). Weight (mean ± SEM) in females remained unaffected (0.40 ± 0.02 and 0.39 ± 0.01 for 9 and 66 fish/L, respectively) while those males exposed to 66 fish/L showed a significant (H = 44.38, d.f. = 1, N = 412, P < 0.001) reduction (0.31 ± 0.01) when compared to males from control treatment (0.23 ± 0.01) (Figure 2A and B). Elevated density during 18 to 45 dpf significantly decreased length in both sexes (H = 4.60, d.f. = 1, N = 224 and H = 44.37, d.f. =1, N = 412 with P < 0.05 and < 0.01 for females and males, respectively) but not during 7–18 dpf (Figure 2C and D).

2.4. Methylation patterns in mature gonads

Based on the results on sex ratio and the effects of elevated density on growth, for methylation analysis we selected ten males and females of the family #4 from the 18–45 dpf group (Table S5). The total number of raw paired reads obtained from sequencing was 4,245.896 with a mean of 106.15 ± 48.33 per sample (Table S5). After trimming, the number of reads was 501.24 and 171.52 for females and males from 9 fish/L and 262.97 and 109.02, for females and males from 66 fish/L density levels, respectively (Table S5). Mapping efficiency was 68.07 ± 8.67% and efficiency conversion 99.55 ± 0.08% from all samples (Table S5).

We examined the DNA methylation in the promoter region of genes coding for four different steroidogenic enzymes (cyp19a1a, cyp11c1, hsd17b1, and hsd11b2), one for a transcription factor in sex-related development (dmrt1) and one epigenetic gene (dnmt1) and how the effects of elevated density during the period 18-45 dpf affect their methylation profile. Results showed the effects of the high-density treatment (66 fish/L) caused an effect on the methylation of a single gene, dnmt1, with low methylation levels in the ovary (~40%) compared to the ovary of control fish (9 fish/L) (~60%). The mean DNA methylation levels were significantly higher in testes for cyp19a1a, cyp11c1, hsd17b1, and hsd11b2 and lower DNA methylation levels for dmrt1 (P < 0.001) when compared to ovaries at 90 dpf. (Figure 3). These differences were due to sexual dimorphism but the effects of elevated density did not affect.

2.5. Gene expression in mature gonads

Dnmt1 was downregulated under the high-density treatment for both ovary and testes, showing significant (P < 0.05) differences only in the latter (fold change of -2.5) compared to the control (Figure 4). The gene expressions of dmrt1 and hsd11b2 were not significantly altered for any of the sexes under the treatment, although a decrease of gene expression of both genes appeared in the testes subjected to the 66 fish/L density. By contrast, cyp19a1a gene was significant downregulated (P < 0.001) in ovaries in the stressed fish (six-fold changes) than females of control group. In addition, no changes in the expression of cyp11c1 were observed for any of the gonadal types at different densities.

2.6. Correlation DNA methylation vs. gene expression in mature gonads

We carried out correlation analyses between DNA methylation and gene expression levels for all of the genes selected (Figure 5). Negative trends in the correlation of methylation and gene expression were observed under both densities (9 fish/L and 66 fish/L) for dmrt1, cyp19a1a and dnmt1 in ovaries, despite non-significant correlation (Figure 5). Hsd11b2, however, showed a positive correlation trend in ovaries under the high-density treatment (ρ = 0.75) although without statistical significance. In testes, negative trends in the correlation of methylation and gene expression were observed for dmrt1, hsd11b2 and dnmt1, although these correlations were not statistically significant. Despite the lack of statistical signification, Cyp11c1 showed a positive correlation trend in males under the high-density treatment both control (ρ = 0.11) and treatment (ρ = 0.72), and cyp19a1a in males of the high-density group (ρ = 0.24).

3. Discussion

Stress is defined as a coordinated series of behavioral and physiological responses to any appreciable challenge to homeostasis or allostasis (Li et al., 2021). As other vertebrates, fish respond to environmental challenges with a series of adaptive neuroendocrine adjustments that are collectively termed ‘the stress response’, which in the short-term could be beneficial but a prolonged activation of the stress response may lead to immunosuppression, reduced growth, and reproductive dysfunction (Braithwaite & Ebbesson, 2014). These deregulations in the reproductive axis are key in the development of the gonads in species in which the environment influences the sexual determination (Castañeda Cortes & Fernandino, 2021).

Subjecting zebrafish larvae during 18 to 45 dpf to high rearing densities (66 fish/L) showed a sex ratio bias towards males. The beginning of the sexual differentiation in zebrafish into an ovary or a testis occurs from 20 to 25 dpf (Chen et al., 2017; Kossack & Draper, 2019; Ye et al., 2022). Nevertheless, Pradhan and Olson (2016) described the first sign of gonadal differentiation at 10 dpf, when the primordial germ cells show signs of oogenesis, although they confirmed that in the period between 20–25 dpf, either the ovary maturation or transformation with progressive degeneration of oocytes to initiate testis differentiation occurs. In our study, the masculinization by effects of elevated density was not observed when applied in earlier periods of development (7–18 dpf). In a similar manner, treatments with DNA methyltransferase inhibitor (5-aza-dC) altered sex ratios when larvae were exposed at 20–30 dpf, with no significant effects at 10–20 dpf (Ribas et al., 2017d). However, high temperature during 7–21 dpf was able to masculinize populations of zebrafish (Ribas et al., 2017a). Here we showed that the window of sensitivity to density factor occurs approximately halfway through gonad formation and should be considered when rearing zebrafish in the husbandry protocols.

In the present study, inter-family variation in sex ratios was observed. In fact, although the overall percent of males increased under high-density in all tested families, only two out of the six families (#4 and #6) showed a significant increase of males. The inter-family variation in zebrafish was described by Liew et al. (2012), reporting the importance of the interaction between genetic and the environment (GxE) factors in this fish species. Later temperature and density experiments have corroborated the existence of these inter-family variation in zebrafish (Delomas & Dabrowski, 2017; Ribas et al 2017a, b, c). Moreover, the importance of the ‘natural’ strain with genetic sex determination using different families has been highlighted, being responsible of the sensitivity to the environment (Valdivieso et al., 2022b). Therefore, although the sex ratio of all families was not altered by the environment, it is important to take into consideration the environmental conditions to which the fish are exposed during gonadal development.

Growth was inversely related to stocking density during 18 to 45 dpf. This result was observed by Hazlerigg et al. (2012) in zebrafish, and in other cultured species (Björnsson, 1994; El-Sayed, 2002; Onxayvieng et al., 2021). In particular, the males showed lower growth than the females. The presence of the sexual dimorphism in the growth response of those fishes subjected to high-density might be related to the faster growth of the females during development in zebrafish (Lawrence et al., 2008). This effect might counteract with the environmental cues during the sensitive window of the sex differentiation.

In the last decade in aquaculture, understanding the epigenetic events under environmental cues has been expanded. The most studied epigenetic event is DNA methylation. Nevertheless, efforts toward identifying the epigenetic marks that affect the expression of genes and how they persist throughout life still need years of research. After different temperature regimes during early development in the E. sea bass, the methylation patterns of 70 genes were identified as potential biomarkers as diagnosis and prognosis being on a gene that fulfilled all the criteria, the keratin-associated protein 10–4 (krtap10–4) (Valdivieso et al., 2023). Recently in zebrafish, we have identified biomarkers able to predict past-thermal effects (Valdivieso et al., 2022a). By using machine learning strategies, differential methylation levels on CpGs in the promoter region of key genes (i.e., cyp19a1a and foxl2 for the females and amh for males) were enough to predict whether an adult zebrafish had suffered high temperature regimes during early gonadal development. With the aim of deciphering the underlying epigenetic mechanisms triggered by high-density insults, in the present study, we selected six genes involved in epigenetics (dnmt1), reproduction (cyp19a1a, dmrt1), and stress (cyp11c1, hsd17b1, hsd11b2). Our results showed that the main alteration in methyl groups of the studied genes was in dnmt1 in ovaries, an epimarker that were significantly hypomethylated in the ovaries almost two months after the density treatment disappeared. DNA methylation changes occur due to the activity of DNA methyltransferases (dnmts) (Jeltsch & Jurkowska, 2014) and thus, they play a central role in the DNA methylation mechanisms. In fish, many studies have shown the alteration of the expression of dnmts, although less data exists on the alteration of the DNA methylation patterns. Most of the studies revealed dnmt alterations of fish subjected to toxic insults in the water, for example, by lead or atrazine treatments (Campos et al., 2012; Wirbisky-Hershberger et al., 2017) and in much less extent their relation with the gonads. For example, dnmt1 was significantly downregulated in ovaries and testes associated with bisphenol exposure concentrations in zebrafish and cyprinid rare minnow, Gobiocypris rarus, which was associated with a reduced of global DNA methylation in the gonads (Liu et al., 2014; Laing et al., 2016). Transgenerational disturbances after bisphenol were observed in some steroidogenic genes in the ovary of rare minnow (Zhu et al., 2021).

Together with the DNA methylation changes in the ovaries, our data showed that the exposure to the rearing stressor was contributing to decrease dnmt1 transcript in female oocytes, but more sharply in the testicular tissue. A relationship between the masculinization effects of high temperature treatments at early stages and dnmts alterations have been observed, for example, in Nile tilapia and zebrafish (Li et al., 2014; Ribas et al., 2017a). The masculinizing steroid 17α-methyltestosterone downregulated the expression of genes involved in chromatin histone modification and DNA methylation pathways in epigenetics, including dnmt1, in 18 to 19 dpf zebrafish gonads (Lee et al., 2017). Zebrafish testes development occurs via programmed cell death from the initial undifferentiated ovary-like gonads (Liew & Orban, 2014). The apoptosis events required a cascade of genes involved in several pathways like p53, wnt signaling pathways and B-cell lymphoma/leukemia-2 (Bcl-2) family (Rodríguez-Marí & Postlethwait, 2011). In mammals, the upregulation of dnmt1 expression under oxidative stress-induced apoptosis via the hypermethylation of Bcl-2 family (Zeng et al.., 2020), thus might alter the sex differentiation in the gonads. In this sense, our observations conclude that changes in the dnmt1 methylation together with its expression in the gonad during sex differentiation were triggered by changes in the environment, i.e., high density.

Based on our data, no more significant changes in the DNA methylation after rearing treatments of the other studied genes (cyp19a1a, dmrt1, cyp11c1, hsd17b1, and hsd11b2) were observed. This might be explained by the fact that the stressor was not severe enough to cause permanent DNA methylation changes in the studied genes as that occurred by temperature in a similar experiment in which six out of ten genes showed significant changes in the DNA methylation patterns (Valdivieso et al., 2022a). Another explanation of our data is the limitation of the performed method. Although it has been useful in several studies (Anastasiadi et al., 2018b; Valdivieso et al., 2022a), the candidate gene approach method is limited to just a small part of the genome, nearby the coding region of the genes and thus wide-genome strategies would have been required to full assess modifications in the epigenome after rearing conditions.

In contrast, clear sexual dimorphism in the DNA methylation patterns were observed in our data. With the emergence of epigenetic-related data in the last years, evidences highlighted the importance of the sexual dimorphic methylation patterns in fish. This is the case, for example in Nile tilapia muscles (Wan et al., 2016) and zebrafish brains (Chatterjee et al., 2016). Sexual dimorphism was present in the methylation pattern of two immune in the zebrafish gonads, indicating that the immune system responses differed between fish sexes (Caballero-Huertas et al., 2020). In fact, the methylation levels of certain genes were linked to a particular sex defining the Essential Epigenetic Marks (EEM) as those informative epigenetic marks that were essential for a particular sexual phenotype (Piferrer et al., 2019). Two genes, cyp19a1a and dmrt1, fulfilled the EEM criteria, while in others, it was not as clear, indicating that more information regarding sexual difference in the fish epigenome is needed. In the present data, we added more information regarding the importance of the sexual differences of the DNA methylation of key genes related to reproduction and stress that might be potential EEM.

High-density during early development was able to downregulate the expression level of cyp19a1a in the ovaries. Our data is in accordance with that observed by Valdivieso et al. (2019), after subjecting fish to even higher densities (i.e., 74 fish/L). Other environmental factors were able to decrease the levels of cyp19a1a expression in fish, as occurred with high temperature (Navarro-Martin et al., 2011; Li et al., 2014; Ribas et al., 2017a) and hypoxia (Shang et al., 2006; Hala et al., 2012). Therefore, based on all the available data, the lower expression of cyp19a1a in the ovaries of fish subjected to environmental stressors might be considered a good biomarker. The other studied genes did not show a significant regulation of their expression after the stressor, data in accordance with no more changes in the DNA methylation patterns.

Correlation analysis of methylation and gene expression in the six studied genes indicated no significant differences neither in ovaries nor testes. Most of the studied genes generally indicated negative relationships following the classical dogma pattern in which low DNA methylation of CpG-rich promoters is associated with the activation of the gene transcription machinery (Jones & Takai, 2001; Deaton & Bird, 2011). In contrast, in our data we found positive correlations between DNA methylation and gene expression in hsd11b2 in the ovary, and cyp19a1a and cyp11c1 in the testis. In a similar manner, positive correlations were observed in hsd11b2 and cyp11c1 in both testis and ovaries of zebrafish exposed to high temperature (Valdivieso et al., 2022a), corroborating our findings. In fact, it is currently known that the DNA methylation patterns are more complex than originally thought as other genomic elements rather than the promoters, such as gene body or introns, can contribute to transcriptional regulation (Brenet et al., 2011; Blattler et al., 2014; Anastasiadi et al., 2018a).

4. Materials and methods

4.1. Animal rearing conditions and facility

Zebrafish (AB strain) were housed in facility at the ‘Institut de Ciències del Mar’ (ICM-CSIC) in Barcelona. Fish were reared in a commercial water rack system (Aquaneering, USA). Husbandry and water conditions were carried out as previously described (Ribas et al., 2017b). Fish were fed ad libitum three times per day with zebrafish commercial food (Aqua Schwarz Gmbh, Germany) and supplemented with brine shrimp (Artemia nauplii). To prevent any potential negative effects of masculinization prior to commencing the density experiments, the larvae were initially reared at a density of 9 fish per liter (fish/L) (Ribas et al., 2017b). The experimental density of 66 fish/L was selected as a potential treatment to induce a skewed sex ratio towards males while ensuring minimal mortality that could hinder the representativeness of the population (Ribas et al., 2017b).

4.2. Experimental design

Larvae offspring from eight unrelated families (#1–8) originated from independent pairs were used to assess the effects of density during sex differentiation. A total of three and six families were used for the first (7–18 dpf) or the second (18– 45 dpf) studied periods, respectively. For each family, in order to optimize the number of survived larvae and to accomplish the desired level density treatment in each tank, we calculated the proper volume to confine the larvae using a micro-perforated barrier in Aquaneering tanks (2.8 L volume, ZT280). The barrier was positioned along the x-axis position while the y- and z-axis remained fixed (Figure S1A). Three groups were created: control, 7–18 dpf, and 18–45 dpf (Figure S1B). For the control group, tanks started with 9 fish/L at 7 dpf and the density level remained unchanged until 90 dpf. During the experiment, fish survival was recorded and the barrier was adjusted as needed to maintain density at 66 fish/L. Water quality parameters were monitored routinely in order to avoid any other environmental factor that could interfere with fish health and development. In the 7–18 dpf group, tanks started with larvae confined with 66 fish/L at 7 dpf, and larvae remained confined until 18 dpf. At this point, the physical barrier was removed to create a non-masculinization density effect that would not lead to masculinization (specifically below 16 fish/L) (Ribas et al 2017b). In the 18–45 dpf group, the larvae were initially maintained at a density of 9 fish/L until 18 dpf, after which they were confined at 66 fish/L until 45 dpf. From 45 dpf onwards, the barrier was removed and fish were kept at a density below 16 fish/L until 90 dpf. If it was not possible to achieve the target density threshold of 16 fish/L after removing the barrier, larvae were evenly distributed into new tanks. At 90 dpf, fish were euthanized using iced water. Body weight (BW ± 0.05 g) and standard length (SL ± 0.01 cm) were recorded. Fish gonads were recorded for sex ratio analysis and gonads were kept at -80ºC for further molecular analysis.

4.3. DNA and RNA extractions

Genomic DNA and total RNA were isolated from the same gonad sample to allow comparisons of DNA methylation and gene expression in the same individual. For DNA extraction, gonads were treated overnight at 60°C with a digestion buffer containing 1 μg of proteinase K (P2308, Sigma-Aldrich, St. Louis, Missouri). Then, a standard phenol-chloroform-isoamyl alcohol protocol (PCI 25:24:1, v/v/v) with 0.5 μg ribonuclease A (12091021, PureLink RNase A, Life Technologies, Carlsbad, California) was performed to eliminate RNA traces. Five hundred nanograms of DNA per sample were bisulfite-converted using the EZ DNA Methylation-Direct™ Kit (ZymoResearch; D5023). For RNA extraction, samples were carried out by Trizol Reagent (T9424, Sigma-Aldrich, St. Louis, Missouri) according to the manufacturer’s instructions and RNA. Quality and quantity of DNA and RNA were measured by NanoDrop (ND-1000) spectrophotometer and Qubit (Thermo Fisher Scientific, Waltham, Massachusetts).

4.4. Methylation Bisulfite Sequencing analysis

Gonadal DNA methylation levels were studied by Methylation Bisulfite Sequencing (MBS) technique following the procedures described elsewhere (Anastasiadi et al., 2018b; Moraleda-Prados et al., 2021; Valdivieso et al., 2022a). The studied genes were epigenetic-related (dnmt1), reproduction-related (dmrt1, cyp19a1a), and stress-related (cyp11c1, hsd17b1, and hsd11b2). The selected region for each gene, included the promotor, the first exon, and the first intron extended as much as possible. The targeted portion of 500 bp approximately was amplified with PCR by using designed primers from Valdivieso et al., 2023) (Table S1). Adaptor sequences for 16 S metagenomic library preparation (Illumina) were added to the 5ʹ ends of the primers designed: Forward-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG and Reverse-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG. PCR products were indexed by Nextera XT index Kit Set A (Illumina; FC-131–2001) according to Illumina’s protocol for 16 S metagenomic library preparation and were pooled in an equimolar manner to obtain a single multiplexed library which was sequenced in a MiSeq (Illumina, San Diego, California) using the paired-end (PE) reads 250 bp protocol at the National Center of Genomic Analysis (CNAG, Barcelona). Raw sequencing data were submitted in Gene Expression Omnibus with accession number GSE134646.

4.5. Bioinformatic analysis

The adapters were removed using Trim galore! software (v. 0.4.5) (Krueger, 2015) and those reads with low-quality filtered (Phred score < 20) were discarded. To ensure that the adapters were cut off correctly, pre- and post-trim quality control was carried out using MultiQC and FastQC (Andrews, 2010; Ewels et al., 2016). Bismark (v.20.0) was used to generate an in silico bisulfite converted zebrafish genome (GCF_000002035) using bismark_genome_preparation function. Deduplication and methylation extraction at CpG_context was carried out using the deduplicate_bismark and bismark_methylation_extractor functions, respectively. Bisulfite conversion efficiency was calculated with a minimum threshold of 99.0% to accept the sample for methylation analysis. The coordinates of the CpG sites of the targeted genes were assessed through ‘bedr’ package (v. 1.0.4) (Haider et al., 2016). From the targeted CpG sites, only those that showed coverage > 5 times were retained for the methylation analysis.

4.6. Gene expression by quantitative PCR

For each sample, 200 ng of RNA were treated with DNAse I, Amplification Grade (Thermo Fisher Scientific Inc., Wilmington, DE, USA) and retrotranscribed to cDNA with SuperScript III RNase Transcriptase (Invitrogen, Spain) with Random hexamer (Invitrogen, Spain). Quantitative PCR (qPCR) was carried out in technical triplicates for each sample with the SYBR Green chemistry (Power SYBR Green PCR Master Mix; Applied Biosystems). The conditions in the thermocycler were: 50ºC for 2 min, 95ºC for 10 min, followed by 40 cycles of 95ºC for 15 s and 60ºC for 1 min in a 384-well plate (CFX-386, Touch BioRad). The qPCRs were run in optically clear 384-well plates. Finally, a temperature- determining dissociation step was performed at 95ºC for 15 s (sec), 60ºC for 15 s and 95ºC for 15 s at the end of the amplification phase. Dissociation step, primers efficiency curves and PCR product sequencing confirmed the specificity for each primer pair. The six qPCR primers used together with the two reference genes are shown in Table S2.

4.7. Statistical analysis

All statistical analyses were performed using R software (R Development Core Team, 2018). Data were expressed as mean ± SEM and the differences were considered significant when P < 0.05. Graphs were generated using the ‘ggplot2 package (v. 3.1.0) (Wickham et al., 2016). Sex ratio analysis was calculated using the Chi-square test with the application of the Yates correction (Yates, 1934) for each family. The overall analysis of density was assessed by Generalized Linear Mixed Model (GLMM). Density level (continuous variable) was considered fixed effect factors and family was random effect factors where the response variable was the proportion of males and females. For biometry data, normality was checked with the Kolmogorov–Smirnov test and logarithmic transformations were applied when necessary. The homoscedasticity of variances was checked with Levene's test. Means were compared by one-way analysis of variance with a Tukey's post hoc multiple-range test. Significant differences were accepted when P ≤ 0.05.

To obtain the DNA methylation across all the CpG sites in the studied gene region, the DNA methylation level for each site was averaged from all samples in each group. The effects of density on DNA methylation were evaluated by two-way ANOVA followed by a Tukey's HSD test. Normality of the residuals was checked by the Shapiro–Wilk test, and a Levene's test for homogeneity of variance was used.

Data obtained from qPCR were collected by SDS 2.3 and RQ Manager 1.2 software. For each sample, the relative quantity (RQ) values of the genes of interest were used to normalize against the geometric mean value of two references genes validated for zebrafish (Tang et al., 2007) and the fold change was calculated using the 2ΔΔCt method (Schmittgen & Livak, 2008). Two-way ANOVA was used to detect differences in gene expression between treatments (sex and density), with the previous check of homoscedasticity by Levene's test for every single group, as well as normality by the Shapiro-Wilk test for each group. When normality was not followed, a Kruskal-Wallis’ test was performed. Tukey's test was used to perform post hoc multiple comparisons. DNA methylation and gene expression correlation analyses were carried out by the application of Spearman's rank correlation (ρ) test using the ‘corrplot; package (v. 0.84).

4.8. Ethics statement

The experimental protocol was approved by the Spanish National Research Council (CSIC) Ethics Committee within the project AGL2015-73864-JIN and licensed by the Bioethical Committee of the Generalitat de Catalunya under reference code 9977. European regulations of animal welfare (ETS N8 123, 01/01/91) were respected regarding fish maintenance. Likewise, ICM facilities were validated for animal experimentation by the Ministry of Agriculture and Fisheries (certificate number 08039-46-A) in accordance with Spanish law (R.D. 223 of March 1988).

5. Conclusions

This study confirms the impact of continued exposure to stressors during gonadal differentiation on DNA methylation, gene expression, and final sexual phenotypes. Present data revealed that the sensible window for high-density occurred between 18 to 45 dpf, and not earlier, skewing sex ratios towards males. Among six studied genes related to epigenetics, reproduction, and stress, dnmt1 showed significant differences in the DNA hypomethylation acting as an epimarker in the male fish subjected to stressors in the two previous months. In a similar manner, the inhibition of cyp19a1a gene expression could be considered as a promising biomarker in the ovaries of fish subjected to rearing densities during sex differentiation. This analysis shed light on the characteristic of fish reproductive system and the environment's role in modulating epigenetic and gene expression molecular events. Thus, caution in the rearing protocols must be considered in the fish facilities, particularly those that use zebrafish as an animal model for research.

Supplementary Materials

Supplementary tables, Supplementary table S1. Information of the sex ratio obtained for the experiment. Supplementary table S2. Biometry of fish after density treatments at 90 days post fertilization. Supplementary table S3. Information of the primer sequences, the coordinate location in the Danio rerio genome (GRCz11) and the number of CpG sites of three immune related genes used for DNA methylation analysis. Abbreviation: base pairs, bp; transcription start site, TSS. Supplementary table S4. Information of primer sequences used for quantitative (qPCR) gene expression. Abbreviation: base pairs, bp. Supplementary table S5. Information of quality raw data, the number of reads before and after trimming process, alignment efficiency and the bisulfite conversion efficiency. Supplementary figures legends, Figure S1. Experimental design to induce stress by density. A) Tanks used for this experiment (2.8 L volume, ZT280, Aquaneering) used and how the animals were confined to control or high-density conditions. B) Experimental design which consisted of control, 7–18 days post fertilization (dpf), and 18–45 dpf, in which for control group, tanks started with 9 fish/L at 7 dpf and the density level remained unmodified until 90 dpf; in 7–18 group, tanks started with 66 fish/L at 7 dpf and fish remained confined until 18 dpf; in 18–45 group, larvae were set up at 9 fish/L until 18 dpf, and then larvae were confined at 66 fish/L until 45 dpf and set up at 9 fish/L until 90 dpf.

Funding

This study was supported by the Spanish Ministry of Science and Innovation grant AGL2015-73864-JIN “Ambisex” and 2PID2020- 113781RB-I00 “MicroMet” and by the Consejo Superior de Investigaciones Científicas (CSIC) grant 02030E004 “Interomics” to LR. We thank the lab technician Sílvia Joly for her essential assistance in our team and Gemma Fusté for her assistance in fish facilities. This study was supported by funding from the Spanish government through the ‘Severo Ochoa Centre of Excellence accreditation (CEX2019-000928-S).

Competing interests

The authors declare no competing or financial interests.

Credit author statement

JMP and MCH performed the experiments and sampling data. JMP and AV carried out bioinformatic and statistical analyses. AV, MCH, JMP, FP, and LR wrote the article. LR conceived the study. All authors read and approved the final manuscript.

References

Abozaid, H.; Wessels, S.; Hörstgen-Schwark, G. Elevated Temperature Applied during Gonadal Transformation Leads to Male Bias in Zebrafish (Danio rerio). Sex. Dev. 2012, 6, 201–209. [Google Scholar] [CrossRef] [PubMed]
Acevedo-Rodriguez, A.; Kauffman, A.S.; Cherrington, B.D.; Borges, C.S.; Roepke, T.A.; Laconi, M. Emerging insights into hypothalamic-pituitary-gonadal axis regulation and interaction with stress signalling. J. Neuroendocr. 2018, 30, e12590. [Google Scholar] [CrossRef] [PubMed]
Ambrosi, C.; Manzo, M.; Baubec, T. Dynamics and Context-Dependent Roles of DNA Methylation. J. Mol. Biol. 2017, 429, 1459–1475. [Google Scholar] [CrossRef] [PubMed]
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. [Google Scholar] [CrossRef] [PubMed]
Anastasiadi, D.; Vandeputte, M.; Sánchez-Baizán, N.; Allal, F.; Piferrer, F. Dynamic epimarks in sex-related genes predict gonad phenotype in the European sea bass, a fish with mixed genetic and environmental sex determination. Epigenetics 2018, 13, 988–1011. [Google Scholar] [CrossRef]
Andrews, S. FastQC: a quality control tool for high throughput sequence data. 2010.
Baroiller, J.; D’cotta, H.; Saillant, E. Environmental Effects on Fish Sex Determination and Differentiation. Sex. Dev. 2009, 3, 118–135. [Google Scholar] [CrossRef] [PubMed]
Besedovsky, H.O.; del Rey, A. Immune-neuro-endocrine interactions: facts and hypotheses. Endocrine reviews 1996, 17, 64–102. [Google Scholar] [CrossRef]
Besson, M.; Komen, H.; Aubin, J.; de Boer, I.J.M.; Poelman, M.; Quillet, E.; Vancoillie, C.; Vandeputte, M.; van Arendonk, J.A.M. Economic values of growth and feed efficiency for fish farming in recirculating aquaculture system with density and nitrogen output limitations: a case study with African catfish (Clarias gariepinus). Journal of Animal Science 2014, 92, 5394–5405. [Google Scholar] [CrossRef]
Best, C.; Ikert, H.; Kostyniuk, D.J.; Craig, P.M.; Navarro-Martin, L.; Marandel, L.; Mennigen, J.A. Epigenetics in teleost fish: From molecular mechanisms to physiological phenotypes. Comp. Biochem. Physiol. Part B: Biochem. Mol. Biol. 2018, 224, 210–244. [Google Scholar] [CrossRef]
Bestor, T.H. The DNA methyltransferases of mammals. Hum. Mol. Genet. 2000, 9, 2395–2402. [Google Scholar] [CrossRef]
Bird, A. DNA methylation patterns and epigenetic memory. Minerva Anestesiol. 2002, 16, 6–21. [Google Scholar] [CrossRef]
Björnsson, B. Effects of stocking density on growth rate of halibut (Hippoglossus hippoglossus L.) reared in large circular tanks for three years. Aquaculture 1994, 123, 259–270. [Google Scholar] [CrossRef]
Blattler, A.; Yao, L.; Witt, H.; Guo, Y.; Nicolet, C.M.; Berman, B.P.; Farnham, P.J. Global loss of DNA methylation uncovers intronic enhancers in genes showing expression changes. Genome Biol. 2014, 15, 1–16. [Google Scholar] [CrossRef]
Braithwaite, V.; Ebbesson, L. Pain and stress responses in farmed fish. Rev Sci Tech 2014, 33, 245–253. [Google Scholar] [CrossRef] [PubMed]
Brenet, F.; Moh, M.; Funk, P.; Feierstein, E.; Viale, A.J.; Socci, N.D.; Scandura, J.M. DNA Methylation of the First Exon Is Tightly Linked to Transcriptional Silencing. PLoS ONE 2011, 6, e14524. [Google Scholar] [CrossRef]
Caballero-Huertas, M.; Moraleda-Prados, J.; Joly, S.; Ribas, L. Immune genes, IL1β and Casp9, show sexual dimorphic methylation patterns in zebrafish gonads. Fish Shellfish. Immunol. 2020, 97, 648–655. [Google Scholar] [CrossRef]
Campos, C.; Valente, L.M.; Fernandes, J.M. Molecular evolution of zebrafish dnmt3 genes and thermal plasticity of their expression during embryonic development. Gene 2012, 500, 93–100. [Google Scholar] [CrossRef] [PubMed]
Castañeda-Cortés, D.C.; Fernandino, J.I. Stress and sex determination in fish: from brain to gonads. Int. J. Dev. Biol. 2021, 65, 207–214. [Google Scholar] [CrossRef] [PubMed]
Cavalieri, V.; Spinelli, G. Environmental epigenetics in zebrafish. Epigenetics Chromatin 2017, 10, 1–11. [Google Scholar] [CrossRef]
Chatterjee, A.; Lagisz, M.; Rodger, E.J.; Zhen, L.; Stockwell, P.A.; Duncan, E.J.; Horsfield, J.A.; Jeyakani, J.; Mathavan, S.; Ozaki, Y.; et al. Sex differences in DNA methylation and expression in zebrafish brain: a test of an extended ‘male sex drive’ hypothesis. Gene 2016, 590, 307–316. [Google Scholar] [CrossRef]
Chen, W.; Ge, W. Ontogenic Expression Profiles of Gonadotropins (fshb and lhb) and Growth Hormone (gh) During Sexual Differentiation and Puberty Onset in Female Zebrafish1. Biol. Reprod. 2012, 86, 73. [Google Scholar] [CrossRef]
Chen, W.; Ge, W. Gonad differentiation and puberty onset in the zebrafish: Evidence for the dependence of puberty onset on body growth but not age in females. Mol. Reprod. Dev. 2013, 80, 384–392. [Google Scholar] [CrossRef]
Chen, W.; Liu, L.; Ge, W. Expression analysis of growth differentiation factor 9 (Gdf9/gdf9), anti-müllerian hormone (Amh/amh) and aromatase (Cyp19a1a/cyp19a1a) during gonadal differentiation of the zebrafish, Danio rerio. Biology of Reproduction 2017, 96, 401–413. [Google Scholar] [CrossRef] [PubMed]
Costa, B.; Pinto, I.C. Stress, burnout and coping in health professionals: A literature review. Journal of Psychology and Brain Studies 2017, 1, 1–8. [Google Scholar]
Dahm, R.; Geisler, R. Learning from Small Fry: The Zebrafish as a Genetic Model Organism for Aquaculture Fish Species. Mar. Biotechnol. 2006, 8, 329–345. [Google Scholar] [CrossRef]
Deaton, A.M.; Bird, A. CpG islands and the regulation of transcription. Minerva Anestesiol. 2011, 25, 1010–1022. [Google Scholar] [CrossRef]
Delomas, T.A.; Dabrowski, K. The importance of controlling genetic variation-remarks on ‘Appropriate rearing density in domesticated zebrafish to avoid masculinization: links with the stress response’. J. Exp. Biol. 2017, 220, 4078–4078. [Google Scholar] [CrossRef]
Dimitriadi, A.; Beis, D.; Arvanitidis, C.; Adriaens, D.; Koumoundouros, G. Developmental temperature has persistent, sexually dimorphic effects on zebrafish cardiac anatomy. Sci. Rep. 2018, 8, 8125. [Google Scholar] [CrossRef] [PubMed]
Edmands, S. Sex Ratios in a Warming World: Thermal Effects on Sex-Biased Survival, Sex Determination, and Sex Reversal. J. Hered. 2021, 112, 155–164. [Google Scholar] [CrossRef]
Edwards, J.R.; Yarychkivska, O.; Boulard, M.; Bestor, T.H. DNA methylation and DNA methyltransferases. Epigenetics Chromatin 2017, 10, 23. [Google Scholar] [CrossRef]
El-Sayed, A.F.M. Effects of stocking density and feeding levels on growth and feed efficiency of Nile tilapia (Oreochromis niloticus L.) fry. Aquaculture research 2002, 33, 621–626. [Google Scholar] [CrossRef]
Ewels, P.; Magnusson, M.; Lundin, S.; Käller, M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 2016, 32, 3047–3048. [Google Scholar] [CrossRef]
Filby, A.L.; Paull, G.C.; Bartlett, E.J.; Van Look, K.J.; Tyler, C.R. Physiological and health consequences of social status in zebrafish (Danio rerio). Physiol. Behav. 2010, 101, 576–587. [Google Scholar] [CrossRef]
Haggerty, C.; Kretzmer, H.; Riemenschneider, C.; Kumar, A.S.; Mattei, A.L.; Bailly, N.; Gottfreund, J.; Giesselmann, P.; Weigert, R.; Brändl, B.; et al. Dnmt1 has de novo activity targeted to transposable elements. Nature structural & molecular biology 2021, 28, 594–603. [Google Scholar]
Haider, S.; Waggott, D.; Lalonde, E.; Fung, C.; Liu, F.-F.; Boutros, P.C. A bedr way of genomic interval processing. Source Code Biol. Med. 2016, 11, 1–7. [Google Scholar] [CrossRef]
Hala, D.; Petersen, L.H.; Martinovic, D.; Huggett, D.B. Constraints-based stoichiometric analysis of hypoxic stress on steroidogenesis in fathead minnows, Pimephales promelas. J. Exp. Biol. 2012, 215, 1753–1765. [Google Scholar] [CrossRef]
Hazlerigg, C.R.E.; Lorenzen, K.; Thorbek, P.; Wheeler, J.R.; Tyler, C.R. Density-Dependent Processes in the Life History of Fishes: Evidence from Laboratory Populations of Zebrafish Danio rerio. PLoS ONE 2012, 7, e37550. [Google Scholar] [CrossRef]
Herpin, A.; Schartl, M. Dmrt1 genes at the crossroads: a widespread and central class of sexual development factors in fish. FEBS J. 2011, 278, 1010–1019. [Google Scholar] [CrossRef]
Hsu, C.W.; Pan, Y.J.; Wang, Y.W.; Tong, S.K.; Chung, B.C. Changes in the morphology and gene expression of developing zebrafish gonads. General and comparative endocrinology 2018, 265, 154–159. [Google Scholar] [CrossRef] [PubMed]
Jeltsch, A.; Jurkowska, R.Z. New concepts in DNA methylation. Trends Biochem. Sci. 2014, 39, 310–318. [Google Scholar] [CrossRef] [PubMed]
Jones, P.A. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 2012, 13, 484–492. [Google Scholar] [CrossRef]
Jones, P.A.; Takai, D. The Role of DNA Methylation in Mammalian Epigenetics. Science 2001, 293, 1068–1070. [Google Scholar] [CrossRef] [PubMed]
Kossack, M.E.; Draper, B.W. Genetic regulation of sex determination and maintenance in zebrafish (Danio rerio). Current topics in developmental biology 2019, 134, 119–149. [Google Scholar]
Krueger, F. Trim galore. A wrapper tool around Cutadapt and FastQC to consistently apply quality and adapter trimming to FastQ files. 2015.
Kuwamura, T.; Kadota, T.; Suzuki, S. Testing the Low-density Hypothesis for Reversed Sex Change in Polygynous Fish: Experiments in Labroides dimidiatus. Sci. Rep. 2014, 4, 4369. [Google Scholar] [CrossRef] [PubMed]
Laing, L.V.; Viana, J.; Dempster, E.L.; Trznadel, M.; Trunkfield, L.A.; Uren Webster, T.M.; van Aerle, R.; Paull, G.C.; Wilson, R.J.; Mill, J.; Santos, E. M. Bisphenol A causes reproductive toxicity, decreases dnmt1 transcription, and reduces global DNA methylation in breeding zebrafish (Danio rerio). Epigenetics 2016, 11, 526–538. [Google Scholar] [CrossRef] [PubMed]
Lawrence, C.; Best, J.; James, A.; Maloney, K. The effects of feeding frequency on growth and reproduction in zebrafish (Danio rerio). Aquaculture 2012, 368–369, 103–108. [Google Scholar] [CrossRef]
Lawrence, C.; Ebersole, J.P.; Kesseli, R.V. Rapid growth and out-crossing promote female development in zebrafish (Danio rerio). Environ. Biol. Fishes 2007, 81, 239–246. [Google Scholar] [CrossRef]
Lee, S.L.J.; Horsfield, J.A.; Black, M.A.; Rutherford, K.; Fisher, A.; Gemmell, N.J. Histological and transcriptomic effects of 17α-methyltestosterone on zebrafish gonad development. BMC Genom. 2017, 18, 557. [Google Scholar] [CrossRef]
Li, L.; Shen, Y.; Yang, W.; Xu, X.; Li, J. Effect of different stocking densities on fish growth performance: A meta-analysis. Aquaculture 2021, 544. [Google Scholar] [CrossRef]
Li, C.G.; Wang, H.; Chen, H.J.; Zhao, Y.; Fu, P.S.; Ji, X.S. Differential expression analysis of genes involved in high-temperature induced sex differentiation in Nile tilapia. Comp. Biochem. Physiol. Part B: Biochem. Mol. Biol. 2014, 177-178, 36–45. [Google Scholar] [CrossRef] [PubMed]
Liew, W.C.; Bartfai, R.; Lim, Z.; Sreenivasan, R.; Siegfried, K.R.; Orban, L. Polygenic Sex Determination System in Zebrafish. PLoS ONE 2012, 7, e34397. [Google Scholar] [CrossRef]
Liew, W.C.; Orbán, L. Zebrafish sex: a complicated affair. Briefings Funct. Genom. 2013, 13, 172–187. [Google Scholar] [CrossRef] [PubMed]
Liu, Y.; Yuan, C.; Chen, S.; Zheng, Y.; Zhang, Y.; Gao, J.; Wang, Z. Global and cyp19a1a gene specific DNA methylation in gonads of adult rare minnow Gobiocypris rarus under bisphenol A exposure. Aquat. Toxicol. 2014, 156, 10–16. [Google Scholar] [CrossRef]
Lutnesky, M.M.F. Density-dependent protogynous sex change in territorial-haremic fishes: models and evidence. Behav. Ecol. 1994, 5, 375–383. [Google Scholar] [CrossRef]
Luzio, A.; Monteiro, S.M.; Rocha, E.; Fontaínhas-Fernandes, A.A.; Coimbra, A.M. Development and recovery of histopathological alterations in the gonads of zebrafish (Danio rerio) after single and combined exposure to endocrine disruptors (17α-ethinylestradiol and fadrozole). Aquat. Toxicol. 2016, 175, 90–105. [Google Scholar] [CrossRef]
Mhanni, A.A.; McGowan, R.A. Variations in DNA (cytosine-5)-methyltransferase-1 expression during oogenesis and early development of the zebrafish. Development genes and evolution 2002, 212, 530–533. [Google Scholar] [CrossRef]
Moraleda-Prados, J.; Caballero-Huertas, M.; Valdivieso, A.; Joly, S.; Ji, J.; Roher, N.; Ribas, L. Epigenetic differences in the innate response after immune stimulation during zebrafish sex differentiation. Dev. Comp. Immunol. 2020, 114, 103848. [Google Scholar] [CrossRef]
Navarro-Martín, L.; Viñas, J.; Ribas, L.; Díaz, N.; Gutiérrez, A.; Di Croce, L.; Piferrer, F. DNA Methylation of the Gonadal Aromatase (cyp19a) Promoter Is Involved in Temperature-Dependent Sex Ratio Shifts in the European Sea Bass. PLoS Genet. 2011, 7, e1002447. [Google Scholar] [CrossRef]
Onxayvieng, K.; Piria, M.; Fuka, M.M.; Gavrilović, A.; Liang, X.; Liu, L.; Tang, R.; Li, L.; Li, D. High stocking density alters growth performance, blood biochemical profiles, and hepatic antioxidative capacity in gibel carp (Carassius gibelio). Fish Physiol. Biochem. 2021, 47, 203–212. [Google Scholar] [CrossRef]
Orban, L.; Sreenivasan, R.; Olsson, P.-E. Long and winding roads: Testis differentiation in zebrafish. Mol. Cell. Endocrinol. 2009, 312, 35–41. [Google Scholar] [CrossRef]
Ospina-Álvarez, N.; Piferrer, F. Temperature-Dependent Sex Determination in Fish Revisited: Prevalence, a Single Sex Ratio Response Pattern, and Possible Effects of Climate Change. PLoS ONE 2008, 3, e2837. [Google Scholar] [CrossRef]
Pagès, H. BSgenome: Software infrastructure for efficient representation of full genomes and their SNPs. R package version 1.66.3. 2021.
Park, K.; Han, E.J.; Ahn, G.; Kwak, I.-S. Effects of thermal stress-induced lead (Pb) toxicity on apoptotic cell death, inflammatory response, oxidative defense, and DNA methylation in zebrafish (Danio rerio) embryos. Aquat. Toxicol. 2020, 224, 105479. [Google Scholar] [CrossRef] [PubMed]
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. [Google Scholar] [CrossRef] [PubMed]
Piferrer, F.; Guiguen, Y. Fish Gonadogenesis. Part II: Molecular Biology and Genomics of Sex Differentiation. Rev. Fish. Sci. 2008, 16, 35–55. [Google Scholar] [CrossRef]
Piferrer, F.; Ribas, L. The use of the zebrafish as a model in fish aquaculture research. Fish Physiology 2020, 38, 273–313. [Google Scholar]
Pradhan, A.; Khalaf, H.; Ochsner, S.A.; Sreenivasan, R.; Koskinen, J.; Karlsson, M.; Karlsson, J.; McKenna, N.J.; Orbán, L.; Olsson, P.-E. Activation of NF-κB protein prevents the transition from juvenile ovary to testis and promotes ovarian development in zebrafish. Journal of Biological Chemistry 2020, 287, 37926–37938. [Google Scholar] [CrossRef]
Pradhan, A.; Olsson, P.-E. Regulation of zebrafish gonadal sex differentiation. AIMS Mol. Sci. 2016, 3, 567–584. [Google Scholar] [CrossRef]
Quinlan, A.R.; Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26, 841–842. [Google Scholar] [CrossRef]
R Development Core Team. R: A language and environment for statistical computing; R Found. Stat. Comput.: Vienna, Austria, 2018. [Google Scholar]
Rai, K.; Nadauld, L.D.; Chidester, S.; Manos, E.J.; James, S.R.; Karpf, A.R.; Cairns, B.R.; Jones, D.A. Zebra Fish Dnmt1 and Suv39h1 Regulate Organ-Specific Terminal Differentiation during Development. Mol. Cell. Biol. 2006, 26, 7077–7085. [Google Scholar] [CrossRef]
Ribas, L.; Liew, W.C.; Díaz, N.; Sreenivasan, R.; Orbán, L.; Piferrer, F. Heat-induced masculinization in domesticated zebrafish is family-specific and yields a set of different gonadal transcriptomes. Proceedings of the National Academy of Sciences 2017, 114, E941–E950. [Google Scholar] [CrossRef]
Ribas, L.; Piferrer, F. The zebrafish (Danio rerio) as a model organism, with emphasis on applications for finfish aquaculture research. Rev. Aquac. 2013, 6, 209–240. [Google Scholar] [CrossRef]
Ribas, L.; Valdivieso, A.; Díaz, N.; Piferrer, F. Response to “The importance of controlling genetic variation – remarks on ‘Appropriate rearing density in domesticated zebrafish to avoid masculinization: links with the stress response’”. J. Exp. Biol. 2017, 220, 4079–4080. [Google Scholar] [CrossRef] [PubMed]
Ribas, L.; Valdivieso, A.; Díaz, N.; Piferrer, F. Response to “The importance of controlling genetic variation-remarks on ‘Appropriate rearing density in domesticated zebrafish to avoid masculinization: links with the stress response’”. Journal of Experimental Biology 2017, 220, 4079–4080. [Google Scholar] [CrossRef]
Ribas, L.; Vanezis, K.; Imués, M.A.; Piferrer, F. Treatment with a DNA methyltransferase inhibitor feminizes zebrafish and induces long-term expression changes in the gonads. Epigenetics Chromatin 2017, 10, 1–16. [Google Scholar] [CrossRef] [PubMed]
Rodríguez-Marí, A.; Postlethwait, J.H. The role of Fanconi anemia/BRCA genes in zebrafish sex determination. In Methods in cell biology; Academic Press, 2011; Volume 105, pp. 461–490.
Santos, D.; Luzio, A.; Coimbra, A.M. Zebrafish sex differentiation and gonad development: A review on the impact of environmental factors. Aquat. Toxicol. 2017, 191, 141–163. [Google Scholar] [CrossRef] [PubMed]
Sarma, O.S.; Frymus, N.; Axling, F.; Thörnqvist, P.-O.; Roman, E.; Winberg, S. Optimizing zebrafish rearing−Effects of fish density and environmental enrichment. Front. Behav. Neurosci. 2023, 17, 1204021. [Google Scholar] [CrossRef] [PubMed]
Segner, H. Zebrafish (Danio rerio) as a model organism for investigating endocrine disruption. Comp. Biochem. Physiol. Part C Toxicol. Pharmacol. 2009, 149, 187–195. [Google Scholar] [CrossRef]
Selye, H. The general adaptation syndrome and the diseases of adaptation1. J. Clin. Endocrinol. Metab. 1946, 6, 117–230. [Google Scholar] [CrossRef]
Shang, E.H.; Yu, R.M.; Wu, R.S. Hypoxia affects sex differentiation and development, leading to a male-dominated population in zebrafish (Danio rerio). Environmental science & technology 2006, 40, 3118–3122. [Google Scholar]
Shao, C.; Li, Q.; Chen, S.; Zhang, P.; Lian, J.; Hu, Q.; Sun, B.; Jin, L.; Liu, S.; Wang, Z.; et al. Epigenetic modification and inheritance in sexual reversal of fish. Genome Res. 2014, 24, 604–615. [Google Scholar] [CrossRef]
Siegfried, K.R. In search of determinants: gene expression during gonadal sex differentiation. J. Fish Biol. 2010, 76, 1879–1902. [Google Scholar] [CrossRef]
Silva, P.; Rocha, M.J.; Cruzeiro, C.; Malhão, F.; Reis, B.; Urbatzka, R.; Monteiro, R.A.; Rocha, E. Testing the effects of ethinylestradiol and of an environmentally relevant mixture of xenoestrogens as found in the Douro River (Portugal) on the maturation of fish gonads—A stereological study using the zebrafish (Danio rerio) as model. Aquat. Toxicol. 2012, 124-125, 1–10. [Google Scholar] [CrossRef] [PubMed]
Smith, J.; Sen, S.; Weeks, R.J.; Eccles, M.R.; Chatterjee, A. Promoter DNA Hypermethylation and Paradoxical Gene Activation. Trends Cancer 2020, 6, 392–406. [Google Scholar] [CrossRef] [PubMed]
Stevens, C.H.; Croft, D.P.; Paull, G.C.; Tyler, C.R. Stress and welfare in ornamental fishes: what can be learned from aquaculture? J. Fish Biol. 2017, 91, 409–428. [Google Scholar] [CrossRef]
Straussman, R.; Nejman, D.; Roberts, D.; Steinfeld, I.; Blum, B.; Benvenisty, N.; Simon, I.; Yakhini, Z.; Cedar, H. Developmental programming of CpG island methylation profiles in the human genome. Nat. Struct. Mol. Biol. 2009, 16, 564–571. [Google Scholar] [CrossRef] [PubMed]
Sun, L.-X.; Wang, Y.-Y.; Zhao, Y.; Wang, H.; Li, N.; Ji, X.S. Global DNA Methylation Changes in Nile Tilapia Gonads during High Temperature-Induced Masculinization. PLoS ONE 2016, 11, e0158483. [Google Scholar] [CrossRef]
Takahashi, H. Juvenile hermaphroditism in the zebrafish, Brachydanio rerio. Bulletin of Fisheries Sciences, Hokkaido University 1977, 28, 57–65. [Google Scholar]
Theodoridi, A.; Dinarello, A.; Badenetti, L.; Pavlidis, M.; Valle, L.D.; Tsalafouta, A. Knockout of the hsd11b2 Gene Extends the Cortisol Stress Response in Both Zebrafish Larvae and Adults. Int. J. Mol. Sci. 2021, 22, 12525. [Google Scholar] [CrossRef]
Tort, L.; Rotllant, J.; Rovira, L. Immunological suppression in gilthead sea bream Sparus aurata of the North-West Mediterranean at low temperatures. Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology 1998, 120, 175–179. [Google Scholar]
Uchida, D.; Yamashita, M.; Kitano, T.; Iguchi, T. Oocyte apoptosis during the transition from ovary-like tissue to testes during sex differentiation of juvenile zebrafish. Journal of Experimental Biology 2002, 205, 711–718. [Google Scholar] [CrossRef]
Uchida, D.; Yamashita, M.; Kitano, T.; Iguchi, T. An aromatase inhibitor or high water temperature induce oocyte apoptosis and depletion of P450 aromatase activity in the gonads of genetic female zebrafish during sex-reversal. Comp. Biochem. Physiol. Part A: Mol. Integr. Physiol. 2003, 137, 11–20. [Google Scholar] [CrossRef]
Valdivieso, A.; Anastasiadi, D.; Ribas, L.; Piferrer, F. Development of epigenetic biomarkers for the identification of sex and thermal stress in fish using DNA methylation analysis and machine learning procedures. Mol. Ecol. Resour. 2022, 23, 453–470. [Google Scholar] [CrossRef] [PubMed]
Valdivieso, A.; Ribas, L.; Monleón-Getino, A.; Orbán, L.; Piferrer, F. Exposure of zebrafish to elevated temperature induces sex ratio shifts and alterations in the testicular epigenome of unexposed offspring. Environ. Res. 2020, 186, 109601. [Google Scholar] [CrossRef]
Valdivieso, A.; Sánchez-Baizán, N.; Mitrizakis, N.; Papandroulakis, N.; Piferrer, F. Development of epigenetic biomarkers with diagnostic and prognostic value to assess the lasting effects of early temperature changes in farmed fish. Aquaculture 2023, 563. [Google Scholar] [CrossRef]
Valdivieso, A.; Wilson, C.A.; Amores, A.; da Silva Rodrigues, M.; Nóbrega, R.H.; Ribas, L.; Postlethwait, J.H.; Piferrer, F. (). Environmentally-induced sex reversal in fish with chromosomal vs. polygenic sex determination. Environmental Research 2022, 213, 113549. [Google Scholar] [CrossRef] [PubMed]
Wan, Z.Y.; Xia, J.H.; Lin, G.; Wang, L.; Lin, V.C.L.; Yue, G.H. Genome-wide methylation analysis identified sexually dimorphic methylated regions in hybrid tilapia. Sci. Rep. 2016, 6, 35903. [Google Scholar] [CrossRef] [PubMed]
Wang, Y.Y.; Sun, L.X.; Zhu, J.J.; Zhao, Y.; Wang, H.; Liu, H.J.; Ji, X.S. Epigenetic control of cyp19a1a expression is critical for high temperature induced Nile tilapia masculinization. J. Therm. Biol. 2017, 69, 76–84. [Google Scholar] [CrossRef]
Wang, H.; Wang, B.; Liu, X.; Liu, Y.; Du, X.; Zhang, Q.; Wang, X. Identification and expression of piwil2 in turbot Scophthalmus maximus, with implications of the involvement in embryonic and gonadal development. Comparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology 2017, 208, 84–93. [Google Scholar] [CrossRef]
Wei, T.; Simko, V.; Levy, M.; Xie, Y.; Jin, Y.; Zemla, J. Package ‘corrplot’. Statistician 2017, 56, e24. [Google Scholar]
Weltzien, F.A.; Andersson, E.; Andersen, Ø.; Shalchian-Tabrizi, K.; Norberg, B. The brain–pituitary–gonad axis in male teleosts, with special emphasis on flatfish (Pleuronectiformes). Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology 2004, 137, 447–477. [Google Scholar]
Wen, A.Y.; You, F.; Sun, P.; Li, J.; Xu, D.D.; Wu, Z.H.; Ma, D.Y.; Zhang, P.J. CpG methylation of dmrt1 and cyp19a promoters in relation to their sexual dimorphic expression in the Japanese flounder Paralichthys olivaceus. Journal of Fish Biology 2014, 84, 193–205. [Google Scholar] [CrossRef]
Wickham, H.; Chang, W.; Wickham, M.H. Package ‘ggplot2’. Create elegant data visualisations using the grammar of graphics. Version 2016, 2, 1–189. [Google Scholar]
A Wilson, C.; High, S.K.; McCluskey, B.M.; Amores, A.; Yan, Y.-L.; A Titus, T.; Anderson, J.L.; Batzel, P.; Carvan, M.J.; Schartl, M.; et al. Wild Sex in Zebrafish: Loss of the Natural Sex Determinant in Domesticated Strains. Genetics 2014, 198, 1291–1308. [Google Scholar] [CrossRef] [PubMed]
Wirbisky-Hershberger, S.E.; Sanchez, O.F.; Horzmann, K.A.; Thanki, D.; Yuan, C.; Freeman, J.L. Atrazine exposure decreases the activity of DNMTs, global DNA methylation levels, and dnmt expression. Food Chem. Toxicol. 2017, 109, 727–734. [Google Scholar] [CrossRef]
Wu, T.P.; Wang, T.; Seetin, M.G.; Lai, Y.; Zhu, S.; Lin, K.; Liu, Y.; Byrum, S.D.; Mackintosh, S.G.; Zhong, M.; et al. DNA methylation on N6-adenine in mammalian embryonic stem cells. Nature 2016, 532, 329–333. [Google Scholar] [CrossRef] [PubMed]
Yates, F. The analysis of multiple classifications with unequal numbers in the different classes. Journal of the American Statistical Association 1934, 29, 51–66. [Google Scholar] [CrossRef]
Ye, D.; Tu, Y.-X.; Wang, H.; He, M.; Wang, Y.; Chen, Z.; Chen, Z.-X.; Sun, Y. A landscape of differentiated biological processes involved in the initiation of sex differentiation in zebrafish. Water Biol. Secur. 2022, 1. [Google Scholar] [CrossRef]
Zeng, H.; Li, T.; He, X.; Cai, S.; Luo, H.; Chen, P.; Chen, Y. Oxidative stress mediates the apoptosis and epigenetic modification of the Bcl-2 promoter via DNMT1 in a cigarette smoke-induced emphysema model. Respir. Res. 2020, 21, 1–14. [Google Scholar] [CrossRef]
Zhang, Q.; Ye, D.; Wang, H.; Wang, Y.; Hu, W.; Sun, Y. Zebrafish cyp11c1 Knockout Reveals the Roles of 11-ketotestosterone and Cortisol in Sexual Development and Reproduction. Endocrinology 2020, 161. [Google Scholar] [CrossRef]
Zhu, L.; Liu, Y.; Xue, X.; Yuan, C.; Wang, Z. BPA’s transgenerational disturbance to transcription of ovarian steroidogenic genes in rare minnow Gobiocypris rarus via DNA and histone methylation. Science of The Total Environment 2021, 762, 143055. [Google Scholar] [CrossRef] [PubMed]

Figure 1. A) Overall percent of males at different exposure windows: 7–18 and 18–45 days post-fertilization (dpf) at two stocking densities (9 and 66 fish/L); B) Proportion of males in each individual family (six families, #F1-6; symbols) and the mean (boxplot) found at 90 days post fertilization of fish exposed to 9 fish/L (control) and 66 fish/L (high density) during 18–45 dpf; C) Plot of the expected masculinization in the populations of zebrafish according to density level based on the prediction fit of a generalized linear mixed model (GLMM). Sex ratio analysis between density treatments was analyzed by χ² test. Abbreviations: N.S. = no significant; *** = P < 0.001. The total number of fishes used for sex ratio was 1,003; 369 (control), 221 (7–18 dpf group), and 413 (18–45 dpf group).

Figure 2. A and B) Body mass (weight in g) and B and C) length (cm) of females and males at different exposure windows: 7–18 and 18–45 days post-fertilization (dpf) at two stocking densities (9 fish/L and 66 fish/L). The total number of fishes used was 1,003; 369 (control), 221 (7–18 dpf), and 413 (18–45 dpf).

Figure 3. Mean DNA methylation levels in the promoter region of A) dnmt1, B) dmrt1, C) cyp19a1a, D) cyp11c1, E) hsd17b1, and F) hsd11b2 in mature gonads of females and males (90 days post-fertilization, dpf) exposed to two rearing densities (9 and 66 fish/L) from the 18–45 dpf group. The number of fishes analysed in each groupis n = 10. Two-way ANOVA followed by a post hoc Tukey test was applied. The P-values for the factor effects of sex (S), density (D) and the interaction of both factors (S x D) are reported for each gene. A robust no parametric two-way ANOVA with trimmed means was applied when data did not follow normality. Data is shown as mean ± S.E.M.

Figure 4. Gene expression of dnmt1, dmrt1, cyp19a1a, cyp11c1, and hsd11b2 in A) ovary and B) testis of 90 days post-fertilization (dpf) after the exposure to two rearing densities (9 and 66 fish/L) from 18 to 45 dpf. Data shown as mean ± SEM. Fold change values of control group (9 fish/L) was set up at 1 as reference. Sample size n = 10 ovaries and n = 10 testes each group. Significant differences between sexes were analyzed by Student's t-test. Abbreviations: * = P < 0.05; *** = P < 0.001; N.S. = no significant.

Figure 5. Correlations of DNA methylation of the promoter regions and gene expression levels for dmrt1, cyp19a1a, cyp11c1, hsd11b2, and dnmt1 in the gonads of females and males (90 days post-fertilization, dpf) exposed to 9 and 66 fish/L during 18–45 dpf. Spearman’s rank correlation coefficients (ρ) are shown. The direction of the long axis of the ellipses and the color indicate the type of correlation: negative is shown in red and positive is in blue shades. The short axis of the ellipse and the intensity of the color are proportional to the correlation coefficients. Significant correlations are considered when P < 0.05.

Table 1. Generalized linear mixed model (GLMM) with logit link function to test effects of elevated density in masculinization of zebrafish during 18–45 days post fertilization (dpf). The density factor (9 and 66 fish/L) was set as fixed whereas family as random factor effects.

Fixed effects	Coefficient	S. E. ^a	Z-value ^b	Pr(>\|Z\|) ^c
Intercept	0.102	0.237	0.430	0.667
Density	0.010	0.003	3.244	0.001

Density 9 fish/L set as references (intercept). Number of observations = 12 and groups = 6. Intercepts of the random effects of the variance and standard deviation for #Family = 0.176 and 0.419, respectively. ^a Standard Error of parameter estimate. ^b Z-value estimate to standard error ratio. ^c Pr (>|Z|) statistic for Z-value.

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.

© 2023 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 (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Exploring the Effects of Rearing Densities on Epigenetic Modifications in the Zebrafish Gonads

Abstract