Skip to main content

Comparative transcriptome analysis of heat-induced domesticated zebrafish during gonadal differentiation



The influence of environmental factors, especially temperature, on sex ratio is of great significance to elucidate the mechanism of sex determination. However, the molecular mechanisms by which temperature affects sex determination remains unclear, although a few candidate genes have been found to play a role in the process. In this study, we conducted transcriptome analysis of the effects induced by high temperature on zebrafish during gonad differentiation period.


Totals of 1171, 1022 and 2921 differentially expressed genes (DEGs) between high temperature and normal temperature were identified at 35, 45 and 60 days post-fertilization (dpf) respectively, revealing that heat shock proteins (HSPs) and DNA methyltransferases (DNMTs) were involved in the heat-exposed sex reversal. The Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway that were enriched in individuals after heat treatment included Fanconi anemia (FA) pathway, cell cycle, oocyte meiosis and homologous recombination.


Our study provides the results of comparative transcriptome analyses between high temperature and normal temperature, and reveals that the molecular mechanism of heat-induced masculinization in zebrafish is strongly related to the expression of HSPs and DNMTs and FA pathway during gonad differentiation.

Peer Review reports


Differentiation at high temperature is one of the most critical issues in developmental biology, and plays important roles in evolution, behavior and gonadal development in animals [1,2,3]. Fish exist a variety of sex determination strategies, including genetic sex determination (GSD), temperature-dependent sex determination (TSD) and GSD + TE (GSD plus temperature effects) [4, 5]. In species with GSD, sex is determined after fertilization and depends on the genetic constitution of an individual, while in TSD species, gonads remain undifferentiated until exposed to environmental temperature at the sensitive period. Interestingly, in GSD + TE species, primary sex determined by the genetic factors, can be changed by extreme environmental temperature [6,7,8,9].

In order to reveal the molecular mechanism underlying GSD + TE, the effects of temperature on sex reversal have been investigated in different species, and the following hypotheses have been proposed. Sex hormone has been suspected to play a critical role in temperature influencing sex determination. The cyp19a1 gene, also known as “aromatase”, which encodes an aromatase that converts androgens into estrogens, is inhibited after exposed to high temperature, resulting in masculinization of individuals [7, 8]. Another typical biological response to temperature is epigenetic modification. Past studies have reported that temperature-specific patterns of DNA methylation and histone modifications of sexual development genes, such as the histone demethylase KDM6B in red-eared slider turtles that regulates transcriptomic of the male sex-determining gene dmrt1 and leads to sex reverse [4, 9,10,11]. Findings in Nile tilapia and American alligator suggested that the heat shock proteins also play an important role in the influence of temperature on sex determination [12, 13].

Zebrafish is a small freshwater teleost that has become an important model to study developmental biology, environmental toxicology and medicine in vertebrates. Zebrafish, all individuals develop a “juvenile ovary” prior to the final differentiation on mature ovaries or testes. Ovaries start their differentiation before at 14 dpf and are fully differentiated at 60 dpf [14, 15], on the other hand, testes begin to differentiate at 30 dpf and are fully differentiated about 60 dpf. What’s more, coexistence of degenerating diplotene oocytes and spermatogonia was evident at 26 dpf, and this structure finally transformed into testes at 34 dpf [14]. All males go through the “juvenile ovary-to-testes” transformation process, during which large numbers of apoptotic oocytes have been observed and were thought to play a crucial role in gonadal development [16, 17]. Several signaling pathways have been reported to be involved in the transformation, such as Tp53-apoptosis, NF-κB and canonical Wnt [18,19,20].

Zebrafish lost sex determining genes during laboratory domestication, and the current studies did not detect heteromorphic sex chromosomes [21, 22]. In the latest years, several molecular studies support that zebrafish has a polygenic sex determination [23,24,25]. The effects of elevated temperature have been examined in zebrafish [26,27,28]. In agreement with the general temperature responsive pattern in fish, elevated temperatures result in a higher proportion of males. However, most of the previous experiments have focuses on cyp19a1, dmrt1, and other genes involved in sex determination or differentiation in zebrafish are known, and little attention has been paid to other genes, such as dnmt1, which is involved in DNA methylation modification [26,27,28]. In addition, the time-course transcriptome analysis of GSD + TE process is minimal by far [28].

Thus, zebrafish provide an excellent opportunity to investigate the influences of temperature on sex determination. In this study, we combined RNA-seq and qPCR approaches to characterize the sex ratio response to temperature in this model, and reveal the genes network responsible for heat-induced masculinization.


Effect of rearing temperature on phenotypic sex ratios in zebrafish

First, we conducted experiments to determine whether high temperature affected individual masculinization. The offspring of 21 adult fish from the same family was obtained to carried out heat treatment experiment during temperature sensitive period (25–35 dpf). The data shown that the male ratio of the normal temperature group was 23.17 ± 3.05% (n = 95) (Fig. 1). And the average male ratio of the temperature-treated group (35 °C, 25–35 dpf) was 82.16 ± 1.79%(n = 104), a significant increase in the masculinization rate was observed compared with the normal temperature group (p<0.001) (Fig. 1).

Fig. 1
figure 1

Early treatment (25–35 dpf) with high temperature (35 °C) increased the proportion of males. The mean ± SEM of three biological replicates for control (28 °C; n = 95) and treated group (35 °C; n = 104) is shown. Statistically significant differences were found between the two sex ratios (***P < 0.001; T-test)

Phenotypic identification and analysis of adult zebrafish gonads

To explore the heat treatment effect on gonadal development, the gonadosomatic index (GSI) and histological analysis were performed for adult fish at 90 dpf. The average GSI of adult female in high temperature group was 24.42% and that in control temperature group was 17.19%, indicating a significant difference (P < 0.05) (Fig. 2a). The GSI analysis of males compared between high and normal temperature-treated also carried out, but there was not statistical difference (Fig. 2a). With regard to the terms of developmental stages of oocytes in zebrafish, there were appeared a series of oocyte maturation stages, including stage IA, stage IB, stage II and stage III in the ovaries of the control group. Among the 20 females from thermal treatment group, 16 of them were observed with several stage IA, stage IB and few stage II oocytes in their ovaries (Fig. 2b).

Fig. 2
figure 2

Phenotypic identification and analysis of temperature-treated fish. (a) the comparative analysis of GSI between high temperature (35 °C) and normal temperature (28 °C) (T-test, *P < 0.05). (b) Histological analysis of gonad tissue. Stage IA and stage IB, primary growth stage; stage II, cortial alveolus stage; stage III, vitellogenesis stage; SC, spermatocytes; SG, spermatogonia; SZ, spermatozoa

Transcriptomic analysis between fish exposed to high and normal temperature

To investigate the effects of temperature during gonad differentiation, RNA-seq on the BGISEQ-500 platform was used to sequence the cDNA library from the whole body of zebrafish at 25, 35, 45 and 60 dpf respectively. A mean of 45.51 M filtered clean reads were obtained from each sample. 82.18% of which mapped to the Danio rerio genome (see Additional file 1). Traditionally, qRT-PCR is used to validate the expression levels quantified by high throughput technology RNA-seq. Therefore, 20 candidate genes associated with cell cycle, oocyte meiosis and homologous recombination were selected from hundreds of differently expressed transcripts between the high and normal temperature treated groups (see Additional file 2 and Additional file 3). Then, the Person’s correlation coefficient indicated that the expression levels were found to be highly reliable for genes that are determined to be significantly differentially expressed by RNA-seq (r = 0.73, p-value< 0.0001) (see Additional file 5).

To explore differences in gene expression between high temperature and normal temperature treatments during gonadal differentiation, we conducted comparative analysis of different expression genes on the three transcriptomic groups, T35D vs. C35D, T45D vs. C45D and T60D vs. C60D respectively (Fig. 3a). Venn diagram was shown (Fig. 3a). To determine the biological function of DEGs during gonadal differentiation after heat-induced, GO classification and KEGG pathway enrichment analysis were performed. (Additional file 6).

Fig. 3
figure 3

Analysis of DEGs between high temperature and normal temperature treatments at 35dpf, 45dpf and 60dpf. Number of DEGs identified from each comparison groups (a), shown by Venn diagram (b), and list of pathways [29, 30] involved in the gonad development in juvenile fish (c)

Comparing T35D and C35D, there were 1171 DEGs, including 681 up-regulated and 490 down-regulated genes (Fig. 3a). GO analysis of DEGs between T35D and C35D was conducted, and the results revealed that the subclasses double−strand break repair, DNA recombination (in biological process category), DNA polymerase activity, nuclease activity (in molecular function category), chromosomal region, replication fork (in cellular component category) were the aspects significantly regulated by heat-induced (Additional file 6a). These results suggest that DNA damage repair and DNA replication were mainly involved in heat-induced of zebrafish. The results of KEGG enrichment analysis revealed that top 20 pathways were significantly enriched(q < 0.05), including FA pathway, cell cycle, oocyte meiosis, homologous recombination and progesterone-mediated oocyte maturation, which were up-regulated (Fig. 3c and Additional file 6d). Especially in FA pathway, 13 genes including fancf, fancg, fanci and fanco were up-regulated in heat-treated fish compared to normal control fish (Fig. 3c and Additional file 3). In addition, the expression level of genes coding DNA methyltransferase family members (dnmt1, dnmt3bb.2) and heat shock proteins (hsp90a, hsp70l, and hspbp1) were also significant upregulated in high-temperature group except hsp4l gene (Additional file 3 and Additional file 5).

Genes expression changes in individuals maintained at high temperature conditions were identified by comparing T45D and C45D. Total of 1022 DEGs, including 436 up-regulated and 586 down-regulated genes (Fig. 3a). The GO analysis results showed that the enriched GO terms were visual perception, synaptic vesicle localization and transport (in biological process category), structural constituent of eye lens, aminoacyl-tRNA ligase activity (in molecular function category), photoreceptor cell cilium, transmembrane transporter complex (in cellular component category) (Additional file 6b). These results suggest that increased temperature affects light perception and synaptic vesicle. In addition, there were top 20 enriched KEGG pathways(q < 0.05), including synaptic vesicle cycle, glutamatergic synapse, fatty acid elongation and dopaminergic synapse (Additional file 6e).

Genes that were differentially expressed in high-temperature group and control group after sex differentiation were evaluated by identifying DEGs in T60D vs C60D. There were 2921 DEGs, including 2362 up-regulated and 559 down-regulated genes (Fig. 3a). The GO analysis results of significantly terms were double−strand break repair, ribosome biogenesis (in biological process category), methyltransferase activity, S-adenosylmethionine-dependent methyltransferase activity (in molecular function category), replication fork, chromosomal region (in cellular component category) (Additional file 6c). From the above results, it is suggested that ribosome biogenesis, transferase activity and DNA replication process were mainly involved in heat-induced of zebrafish. KEGG enrichment analysis identified top 20 pathways(q < 0.05), including metabolic pathway, phototransduction and steroid biosynthesis (Additional file 6f).


Temperature represents one of abiotic factors essential for the development and growth of fish species [31, 32]. The sex reversal effects of high temperature were carried out in zebrafish or other fish species. However, previous researches mainly focused on early embryonic period or different expression genes related to sex differentiation in adult [7, 26,27,28, 33]. While, we performed transcriptome gene expression analysis during the sex differentiation period in zebrafish to investigate the time-course effects of temperature on sex determination.

Numerous studies have shown that DNA damage repair system, heat shock proteins and other factors may effectively respond to heat stress, also known as heat shock response (HSR), which is a defensive adaptive response characterized by changes in gene expression under environmental stress [34,35,36,37]. In our study, GO enrichment analysis results at 35dpf and 60dpf showed that processes such as double-strand break repair, DNA polymerase activity and nuclease activity were significantly regulated by heat induction, suggesting that the DNA damage repair system contributes to HSR. In addition, most HSPs are rapidly produced via heat shock factor (HSF)-regulated gene expression in stressed cells during HSR [37, 38]. However, it is worth mentioning that HSP70 and HSP47 were induced by thermal stress, whereas HSP90A, HSP90B and HSF1 were expressed constitutively and not induced by heat stress in zebrafish [31, 39,40,41]. Interestingly, we observed that exposure to thermal increased the transcript levels of hsp90a, hsp70l and hspbp1, and decreased hsp4l, suggesting that HSPs are also involved in other biological processes.

Previous evidence suggests that pathways such as HSR, FA pathway and cytochrome P450 are critical in gonadal differentiation [7, 8, 39, 42,43,44]. First of all, HSR-related factors such as HSP90A, HSPBP1, HSP4l and HSF5 are involved in spermatogenesis and fertilization in mammals [45,46,47,48]. Particularly, HSP90 is a major molecular chaperones that even mediated estrogen receptor α (ERα) methylation and repressed the ERα activity [49]. Our data support the works in Nile tilapia and American alligator that the HSPs may play multiple roles in TSD, but further research is needed [12, 13].

Furthermore, several gonadal development-related pathways were upregulated in heat-exposed zebrafish compared to normal group, including FA pathway, oocyte meiosis, progesterone-mediated oocyte maturation, cell cycle and homologous recombination. We focused on the FA pathway and discovered that the expression of uhrf1, fancf and other 11 genes were raised in heat-temperature compared to normal temperature. Surprisingly, we detected a significant increase in the expression level of FA signaling pathway after heat stress, which is reported for the first time. What’s more, previous studies have shown that there is crosstalk between FA pathway and cell cycle and homologous recombination pathway, even playing an important role in gonads sterility and sex-reversion [42,43,44]. Therefore, we hypothesize that the FA signaling pathway may play an important role in heat-exposed masculinization in collaboration with other signaling pathways.

In many organisms, genome-wide methylation level in testis (or ovaries) are significantly higher in the hyperthermia group than in the normal group, and epigenetic remodeling is widely believed to play a crucial role in heat-induced masculinization [4, 5, 9, 10, 50]. However, these works focused on differentially expressed genes in sexually mature individuals, and no significant difference was detected in dnmts gene expression. Nevertheless, in this study, we tracked DEGs during gonadal differentiation and found the expression level of the DNA methyltransferase genes dnmt1 and dnmt3bb.2 were up-regulated under heat temperature at 35, 45 and 60 dpf, consistent with the results of Dorts et al [51] Strangely, the results of Ribas et al. showed that dnmt3bb.2 was downregulated after high temperature treatment [28], opposite to our results, which may be related to acute heat treatment or chronic heat treatment. At present, DNMT3A and DNMT1 have been found the effect on fertility by programming the methylation patterns of germ cells via cell cycle process in mammal [52,53,54]. So, what’s the role of dnmts in heat-induced masculinization in zebrafish? Is there an (direct or indirect) interaction between DNA methyltransferase and histone methylation in high-temperature sex reversal? Whole Genome Bisulfite Sequencing approaches combined with classical reverse genetic methods are needed to reveal these questions.


In our study, we present the RNA-seq analysis to learn the mechanism of heat-induced masculinization in zebrafish. Our data demonstrated elevated temperature during the temperature sensitive period leads to an increase in male ratio in domesticated zebrafish, supporting a GSD + TE sex determination. We also show that normal ovarian development is somewhat inhibited in individuals exposed to high temperature. Strikingly, we detected DNA methyltransferase genes (dnmt1 and dnmt3bb.2) and FA pathway were significantly upregulated after heat treated for the first time. These results provide insights into heat-induced sex reversal and help us further understanding the mechanisms of vertebrata sex determination.


Fish and ethics statement

The AB strain were maintained in a glass aquarium with a circulating water system at 28 ± 0.5 °C, with photoperiods of 14 L-10D. All procedures were approved by the Animal Care and Use Committee of Hunan University of Science and Technology. Fish were maintained and conducted in compliance with the ARRIVE guideline [27].

Embryo collection and thermal treatment

In order to avoid the influence of different families on the sex ratio, 21 fish (14 females and 7 males) from the same family mating took place spontaneously in spawning tanks. Fertilized eggs were collected and separated in Petri dishes and maintained at 28 ± 0.5 °C. Eggs hatch at 3 dpf and total 1200 larvae were divided into 2 groups as follows. The control group was kept at 28 °C, whereas the thermal treatment group was exposed to 35 °C for 10 days during 25–35 dpf. And each group had three biological replicates. In order to avoid the effect of high density on gonad differentiation, the number of fish in each 5-liter tank was 50–55. After the treatment, fish were returned at 28 °C until the fish reached to sexual maturity (90 dpf).

The GSI and histological analysis in fish

Dissected gonads from 20 adult individuals of each group at 90 dpf were analyzed. First, their gonads were removed and weighed immediately. The GSI was calculated as the ratio of gonad (Wg) to the total fish weight (W), expressed as a percentage, GSI = [(Wg /W) × 100]. Then, gonads were fixed in 4% paraformaldehyde for 24 h at room temperature, washed in phosphate buffered saline, dehydrated in an ascending gradient of ethanol concentrations, cleared in xylene and embedded. 14 μm thick sections were cut and flatted on glass slides, stained with haematoxylin and eosin and sealed with neutral resin. Sections were examined with a light microscope and images were captured.

Total RNA extraction

5–15 fish were collected at 25, 35, 45 and 60 dpf. Three independent biological replicates were processed at each time node and temperature. Thus, 21 samples labeled as C25D-1 ~ 3, C35D-1 ~ 3, C45D-1 ~ 3, C60D-1 ~ 3, T35D-1 ~ 3, T45D-1 ~ 3, and T60D-1 ~ 3 were acquired. Total RNA for each sample was extracted from the whole body of fish using Trizol (Invitrogen, USA) according to manual instruction. About 60 mg of tissues were ground into powder by liquid nitrogen in a 2 mL tube, followed by being homogenized for 2 minutes and rested horizontally for 5 minutes. The mix was centrifuged for 5 minutes at 12,000×g at 4 °C, then the supernatant was transferred into a new EP tube with 0.3 mL chloroform/isoamyl alcohol (24:1). The mix was shacked vigorously for 15 s, and then centrifuged at 12,000×g for 10 minutes at 4 °C. After centrifugation, the upper aqueous phase where RNA remained was transferred into a new tube with equal volume of supernatant of isopropyl alcohol, then centrifuged at 13600 rpm for 20 minutes at 4 °C. After deserting the supernatant, the RNA pellet was washed twice with 1 mL 75% ethanol, then the mix was centrifuged at 13600 rpm for 3 minutes at 4 °C to collect residual ethanol, followed by the pellet air dry for 5–10 minutes in the biosafety cabinet. Finally, 25–100 μL of DEPC-treated water was added to dissolve the RNA. Subsequently, total RNA was qualified using a Nano Drop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, USA).

mRNA library construction

Oligo (dT)-attached magnetic beads were used to purified mRNA. Purified mRNA was fragmented into small pieces with fragment buffer at appropriate temperature. Then first-strand cDNA was generated using random hexamer-primed reverse transcription, followed by a second-strand cDNA synthesis. Afterwards, A-Tailing Mix and RNA Index Adapters were added by incubating to end repair. The cDNA fragments obtained from previous step were amplified by PCR, and products were purified by Ampure XP Beads, then dissolved in EB solution. The product was validated on the Agilent Technologies 2100 bioanalyzer for quality control. The double stranded PCR products from previous step were heated denatured and circularized by the splint oligo sequence to get the final library. The single strand circle DNA (ssCir DNA) was formatted as the final library. The final library was amplified with phi29 to make DNA nanoball (DNB) which had more than 300 copies of one molecular, DNBs were loaded into the patterned nanoarray and single end 50 bases reads were generated on BGISEQ-500 platform (BGI-Shenzhen, China).

RNA-sequencing data analysis

The sequencing data was filtered with SOAPnuke (v1.5.2) [55] by removing reads containing sequencing adapter; removing reads whose low-quality base ratio (base quality less than or equal to 5) is more than 20%; removing reads whose unknown base (‘N’ base) ratio is more than 5%, afterwards clean reads were obtained and stored in FASTQ format. The clean reads were mapped to the reference genome using HISAT2(v2.0.4) pipeline [56]. Bowtie2 (v2.2.5) [57] was applied to align the clean reads to the reference coding gene set (GCF_000002035.6_GRCz11), then expression level of gene was calculated by RSEM (v1.2.12) [58]. Essentially, differential expression analysis was performed using the DESeq2(v1.4.5) [59] with Q value ≤0.05. To take insight to the change of phenotype, GO ( and KEGG ( enrichment analysis of annotated DEGs were performed by Phyper ( Hypergeometric_distribution) based on Hypergeometric test. The significant level of terms and pathways were corrected by Q value with a rigorous threshold (Q value≤0.05) by Bonferroni [60].

Validation of transcriptome data by real-time PCR

To verify accuracy of the transcriptome data, qRT-PCR analyses were performed using samples from 35, 45 and 60 dpf. Twenty candidate genes with key functions in cell cycle, oocyte meiosis and homologous recombination were randomly selected for qRT-PCR analysis. Total RNAs were reverse transcribed into cDNAs with the Primescript™ RT reagent kit (Takara, Japan). The PCR reaction system comprised the following ingredients: 10 μl 2 × SYBR Green PCR master Mix (Biorad, China), 1 μl of each primer, 2 μl of template cDNA, and 6 μl of RNase-free ddH2O. The cycling program was carried out under the following protocol: 95 °C for 1 min, 40 cycles with each cycle for 15 s at 95 °C, 15 s at 60 °C, and 30 s at 72 °C, before finally ramping from 65 to 95 °C at 0.5 °C per 5 s to generate a melting curve demonstrated no dimer primers, nor unwanted additional PCR products. The zebrafish eef1a1l1 and β2m genes were selected as the control genes. Each qRT-PCR experiments was repeated three times. The relative expression levels of 20 genes in different samples were calculated using the 2-ΔΔCT method among three biological replicates. The primers employed were listed in the Additional file 4.

Statistical analyses

GraphPad Prism 8 was used to analyze all statistical data and make the graphs. Data are presented as the mean ± SEM. Differences between experimental group were made by independent sample T-test. Then the Pearson’s correlation analysis for the data of RNA-seq and qRT-PCR was performed. The statistical significance was declared at P < 0.05, P < 0.01 or P < 0.001.

Availability of data and materials

The complete clean data have been uploaded to NCBI’s Gene Expression Omnibus (GEO) with accession number GSE201748 (



Differentially expressed genes


Days post-fertilization


Heat shock proteins


DNA methyltransferases


Gene Ontology


Kyoto Encyclopedia of Genes and Genomes


Fanconi anemia


Genetic sex determination


Temperature-dependent sex determination


GSD plus temperature effects


Gonadosomatic index


Heat shock response


Heat shock factor


Estrogen receptor α


  1. Shen ZG, Wang HP. Molecular players involved in temperature-dependent sex determination and sex differentiation in teleost fish. Genet Sel Evol. 2014;46(1):26.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Devlin RH, Nagahama Y. Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture. 2002;208:191–364.

    Article  CAS  Google Scholar 

  3. Sun LX, Teng J, Zhao Y, et al. Gonad transcriptome analysis of High-temperature-treated females and High-temperature-induced sex-reversed Neomales in Nile Tilapia. Int J Mol Sci. 2018;19(3).

  4. Shao C, Li Q, Chen S, et al. Epigenetic modification and inheritance in sexual reversal of fish. Genome Res. 2014;24(4):604–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Fan Z, Zou Y, Jiao S, et al. Significant association of cyp19a promoter methylation with environmental factors and gonadal differentiation in olive flounder Paralichthys olivaceus. Comp Biochem Physiol A Mol Integr Physiol. 2017;208:70–9.

    Article  CAS  PubMed  Google Scholar 

  6. Diaz N, Piferrer F. Lasting effects of early exposure to temperature on the gonadal transcriptome at the time of sex differentiation in the European sea bass, a fish with mixed genetic and environmental sex determination. BMC Genomics. 2015;16:679.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Uchida D, Yamashita M, Kitano T, et al. 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 A Mol Integr Physiol. 2004;137(1):11–20.

    Article  PubMed  CAS  Google Scholar 

  8. Navarro-Martin L, Vinas J, Ribas L, et al. 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(12):e1002447.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Ge C, Ye J, Weber C, et al. The histone demethylase KDM6B regulates temperature-dependent sex determination in a turtle species. Science. 2018;360(6389):645–8.

    Article  CAS  PubMed  Google Scholar 

  10. Sun LX, Wang YY, Zhao Y, et al. Global DNA methylation changes in Nile Tilapia gonads during High temperature-induced masculinization. PLoS One. 2016;11(8):e0158483.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. Lai YS, Shen D, Zhang W, et al. Temperature and photoperiod changes affect cucumber sex expression by different epigenetic regulations. BMC Plant Biol. 2018;18(1):268.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Kohno S, Katsu Y, Urushitani H, et al. Potential contributions of heat shock proteins to temperature-dependent sex determination in the American alligator. Sex Dev. 2010;4(1–2):73–87.

    Article  CAS  PubMed  Google Scholar 

  13. Li CG, Wang H, Chen HJ, Zhao Y, Fu PS, Ji XS. Differential expression analysis of genes involved in high-temperature induced sex differentiation in Nile tilapia. Comp Biochem Physiol B: Biochem Mol Biol. 2014;177–178:36–45.

    Article  CAS  Google Scholar 

  14. Hsiao CD, Tsai HJ. Transgenic zebrafish with fluorescent germ cell: a useful tool to visualize germ cell proliferation and juvenile hermaphroditism in vivo. Dev Biol. 2003;262(2):313–23.

    Article  CAS  PubMed  Google Scholar 

  15. Maack G, Segner H. Morphological development of the gonads in zebrafish. J Fish Biol. 2010;62(4):895–906.

    Article  Google Scholar 

  16. Uchida D, Yamashita M, Kitano T, et al. Oocyte apoptosis during the transition from ovary-like tissue to testes during sex differentiation of juvenile zebrafish. J Exp Biol. 2002;205(Pt 6):711–8.

    Article  PubMed  Google Scholar 

  17. Tzung KW, Goto R, Saju JM, et al. Early depletion of primordial germ cells in zebrafish promotes testis formation. Stem Cell Reports. 2015;4(1):61–73.

    Article  CAS  PubMed  Google Scholar 

  18. Pradhan A, Khalaf H, Ochsner SA, et al. Activation of NF-κB protein prevents the transition from juvenile ovary to testis and promotes ovarian development in zebrafish. J Biol Chem. 2012;287(45):37926–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Rodríguez-Marí A, Cañestro C, Bremiller RA, et al. Sex reversal in zebrafish fancl mutants is caused by Tp53-mediated germ cell apoptosis. PLoS Genet. 2010;6(7):e1001034.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Sreenivasan R, Jiang J, Wang X, et al. Gonad differentiation in zebrafish is regulated by the canonical Wnt signaling pathway. Biol Reprod. 2014;90(2):45.

    Article  PubMed  CAS  Google Scholar 

  21. Wilson CA, High SK, McCluskey BM, Amores A, Yan YL, Titus TA, et al. Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains. Genetics. 2014;198(3):1291–308.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Amores A, Postlethwait JH. Banded chromosomes and the zebrafish karyotype. Methods Cell Biol. 1999;60:323–38.

    Article  CAS  PubMed  Google Scholar 

  23. Bradley KM, Breyer JP, Melville DB, et al. An SNP-based linkage map for zebrafish reveals sex determination loci. G3 (Bethesda). 2011;1(1):3–9.

    Article  CAS  Google Scholar 

  24. Liew WC, Bartfai R, Lim Z, et al. Polygenic sex determination system in zebrafish. PLoS One. 2012;7(4):e34397.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Liew WC, Orban L. Zebrafish sex: a complicated affair. Brief Funct Genomics. 2014;13(2):172–87.

    Article  PubMed  Google Scholar 

  26. Abozaid H, Wessels S, Horstgen-Schwark G. Elevated temperature applied during gonadal transformation leads to male bias in zebrafish (Danio rerio). Sex Dev. 2012;6(4):201–9.

    Article  CAS  PubMed  Google Scholar 

  27. Abozaid H, Wessels S, Hörstgen-Schwark G. Effect of rearing temperatures during embryonic development on the phenotypic sex in zebrafish (Danio rerio). Sex Dev. 2011;5(5):259–65.

    Article  CAS  PubMed  Google Scholar 

  28. Ribas L, Liew WC, Diaz N, et al. Heat-induced masculinization in domesticated zebrafish is family-specific and yields a set of different gonadal transcriptomes. Proc Natl Acad Sci U S A. 2017;114(6):E941–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Mottola G, Nikinmaa M, Anttila K. Hsp70s transcription-translation relationship depends on the heat shock temperature in zebrafish. Comp Biochem Physiol A Mol Integr Physiol. 2020;240:110629.

    Article  CAS  PubMed  Google Scholar 

  32. Long Y, Li L, Li Q, He X, Cui Z. Transcriptomic characterization of temperature stress responses in larval zebrafish. PLoS One. 2012;7(5):e37209.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Hosseini S, Brenig B, Tetens J, et al. Phenotypic plasticity induced using high ambient temperature during embryogenesis in domesticated zebrafish Danio rerio. Reprod Domest Anim. 2019;54(3):435–44.

    Article  PubMed  Google Scholar 

  34. Kantidze OL, Velichko AK, Luzhin AV, et al. Heat stress-induced DNA damage. Acta Nat. 2016;8(2):75–8.

    Article  CAS  Google Scholar 

  35. Houston BJ, Nixon B, Martin JH, et al. Heat exposure induces oxidative stress and DNA damage in the male germ line. Biol Reprod. 2018;98(4):593–606.

    Article  PubMed  Google Scholar 

  36. Lopes AR, Figueiredo C, Sampaio E, et al. Impaired antioxidant defenses and DNA damage in the European glass eel (Anguilla anguilla) exposed to ocean warming and acidification. Sci Total Environ. 2021;774:145499.

    Article  CAS  PubMed  Google Scholar 

  37. Chien LC, Wu YH, Ho TN, et al. Heat stress modulates nucleotide excision repair capacity in zebrafish (Danio rerio) early and mid-early embryos via distinct mechanisms. Chemosphere. 2020;238:124653.

    Article  CAS  PubMed  Google Scholar 

  38. Dubrez L, Causse S, Borges Bonan N, Dumétier B, Garrido C. Heat-shock proteins: chaperoning DNA repair. Oncogene. 2020;39(3):516–29.

    Article  CAS  PubMed  Google Scholar 

  39. Keller JM, Escara-Wilke JF, Keller ET. Heat stress-induced heat shock protein 70 expression is dependent on ERK activation in zebrafish (Danio rerio) cells. Comp Biochem Physiol A Mol Integr Physiol. 2008;150(3):307–14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Feugere L, Scott VF, Rodriguez-Barucg Q, Beltran-Alvarez P, Wollenberg Valero KC. Thermal stress induces a positive phenotypic and molecular feedback loop in zebrafish embryos. J Therm Biol. 2021;102:103114.

    Article  CAS  PubMed  Google Scholar 

  41. Murtha JM, Keller ET. Characterization of the heat shock response in mature zebrafish (Danio rerio). Exp Gerontol. 2003;38(6):683–91.

    Article  CAS  PubMed  Google Scholar 

  42. Ramanagoudr-Bhojappa R, Carrington B, Ramaswami M, et al. Multiplexed CRISPR/Cas9-mediated knockout of 19 Fanconi anemia pathway genes in zebrafish revealed their roles in growth, sexual development and fertility. PLoS Genet. 2018;14(12):e1007821.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Tsui V, Crismani W. The Fanconi Anemia pathway and fertility. Trends Genet. 2019;35(3):199–214.

    Article  CAS  PubMed  Google Scholar 

  44. Ceccaldi R, Sarangi P, D'andrea AD. The Fanconi anaemia pathway: new players and new functions. Nat Rev Mol Cell Biol. 2016;17(6):337–49.

    Article  CAS  PubMed  Google Scholar 

  45. Grad I, Cederroth CR, Walicki J, et al. The molecular chaperone Hsp90alpha is required for meiotic progression of spermatocytes beyond pachytene in the mouse. PLoS One. 2010;5(12):e15770.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Held T, Paprotta I, Khulan J, et al. Hspa4l-deficient mice display increased incidence of male infertility and hydronephrosis development. Mol Cell Biol. 2006;26(21):8099–108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Rogon C, Ulbricht A, Hesse M, et al. HSP70-binding protein HSPBP1 regulates chaperone expression at a posttranslational level and is essential for spermatogenesis. Mol Biol Cell. 2014;25(15):2260–71.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Saju JM, Hossain MS, Liew WC, et al. Heat shock factor 5 is essential for spermatogenesis in zebrafish. Cell Rep. 2018;25(12):3252–3261.e4.

    Article  CAS  PubMed  Google Scholar 

  49. Obermann WMJ. A motif in HSP90 and P23 that links molecular chaperones to efficient estrogen receptor alpha methylation by the lysine methyltransferase SMYD2.J. Biol Chem. 2018;293(42):16479–87.

    Article  CAS  Google Scholar 

  50. Skinner MK, Guerrero-Bosagna C, Haque M, et al. Environmentally induced transgenerational epigenetic reprogramming of primordial germ cells and the subsequent germ line. PLoS One. 2013;8(7):e66318.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Dorts J, Falisse E, Schoofs E, Flamion E, Kestemont P, Silvestre F. DNA methyltransferases and stress-related genes expression in zebrafish larvae after exposure to heat and copper during reprogramming of DNA methylation. Sci Rep. 2016;6:34254. Published 2016 Oct 12.

  52. Ren W, Gao L, Song J. Structural basis of DNMT1 and DNMT3A-mediated DNA methylation. Genes (Basel). 2018;9(12).

  53. Kibe K, Shirane K, Ohishi H, Uemura S, Toh H, Sasaki H. The DNMT3A PWWP domain is essential for the normal DNA methylation landscape in mouse somatic cells and oocytes. PLoS Genet. 2021;17(5):e1009570.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Li Y, Zhang Z, Chen J, et al. Stella safeguards the oocyte methylome by preventing de novo methylation mediated by DNMT1. Nature. 2018;564(7734):136–40.

    Article  CAS  PubMed  Google Scholar 

  55. Li R, Li Y, Kristiansen K, et al. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008;24(5):713–4.

    Article  CAS  PubMed  Google Scholar 

  56. Kim D, Langmead B, Salzberg S. L.HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Li B, Dewey C N.RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics,2011, 12: 323.

  59. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  60. Abdi H. The Bonferonni and Šidák corrections for multiple comparisons. Encycloped Measurement Stat. 2007;1:1–9.

    Google Scholar 

Download references


We thank Wei Hu for his valuable comments on the manuscript from Institute of Hydrobiology, Chinese Academy of Sciences. We are also grateful to China Zebrafish Resource Center (CZRC) for the zebrafish material.


This work was supported by the Hunan Province science and technology Department (2018JJ3149), and Education Department of Hunan Province(20B235).

Author information

Authors and Affiliations



Xiaojuan Cui, Yuandong Sun and Yifei Zhang: initial conceptual and experimental design of the study. Chenchen Wang, Xuhuai Chen and Yu Dai: performed the experiment, interpretation of data, key discussions on principle findings. Chenchen Wang and Xiaojuan Cui: wrote and edited the manuscript. All authors read and approved the final version of the manuscript.

Corresponding author

Correspondence to Xiaojuan Cui.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Animal Care and Use Committee of Hunan University of Science and Technology. Fish were maintained and conducted in compliance with the ARRIVE guidelines. This manuscript does not involve in the use of any human material or human data, and all procedures were performed in accordance with relevant guidelines.

Consent for publication

Not applicable.

Competing interests

I declare that the authors have no competing interests as defined by BMC, or other interests that might be perceived to influence the results and discussion reported in this paper.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

Summary statistics of the transcriptome sequencing of D. rerio.

Additional file 2.

Type of heat effect on gene in juveniles at 35, 45 and 60 dpf.

Additional file 3.

RNA-Seq and qRT-PCR validation results.

Additional file 4.

Gene symbols, Refseq IDs and primer sequences for all genes used in qPCR.

Additional file 5 Fig. S1.

The person’s correlation coefficient between RNA-Seq and qRT-PCR.

Additional file 6 Fig. S2.

Enriched GO terms and KEGG pathways.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, C., Chen, X., Dai, Y. et al. Comparative transcriptome analysis of heat-induced domesticated zebrafish during gonadal differentiation. BMC Genom Data 23, 39 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: