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

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.