Histological evaluation of gonadal development
Figure 1 shows the results of the histological analysis of the gonads. As reported by Yan et al., sex reversed larvae were not observed in the control groups . Gonads from the C-XX group occupied the ovarian cavities, which were filled with a small number of oocytes and a large number of oogonia, closely arranged on the oviposition plate. Gonads from the C-XY group were filled with spermatogenic cells at different developmental stages (Fig. 1a and b). However, gonads from the E2-immersed group were smaller than those from the control group. In addition, exposure to E2 obviously induced the feminization of testes and an ovarian cavity was observed in all E2-treated XY torafugu (Fig. 1c and d).
Illumina Sequencing and mapping, and identification of global DEGs, in response to E2 treatment
Transcriptomic analysis generated 44,523,690 (C-XX_1), 45,264,242 (C-XX_2), 66,606,798 (C-XX_3), 43,978,518 (C-XY_1), 46,201,180 (C-XY_2), 46,409,586 (C-XY_3) 44,537,526 (E-XX_1), 40,345,586 (E-XX_2), 49,202,476 (E-XX_3), 47,352,9109 (E-XY_1), 43,506,024 (E-XY_2) and 45,728,118 (E-XY_3) raw reads from each library, respectively (Table 2). After data filtering, 43,622,952 (C-XX_1), 44,436,616 (C-XX_2), 42,850,608 (C-XX_3), 42,851,086 (C-XY_1), 45,072,134 (C-XY_2), 45,005,800 (C-XY_3), 43,442,122 (E-XX_1), 39,629,744 (E-XX_2), 48,253,210 (E-XX_3), 46,193,034 (E-XY_1), 42,700,218 (E-XY_2) and 44,804,906 (E-XY_3) clean reads were obtained from each library, respectively (Table 2).
As shown in Figs. S1, 3, only eight DEGs were identified in the C-XY versus (vs) C-XX comparison. These included three DEG that was up-regulated and one that were down‑regulated, such as aryl hydrocarbon receptor interacting protein-like 1 (aipl1), serine protease hepsin-like, retinol dehydrogenase 11-like (rdh11), and nucleoprotein TPR-like (Table 3). In the E-XX vs C-XX comparison, 514 DEGs were identified, of which 85 were up-regulated (Fig. 2a, 3-4). These included gonadotropin-releasing hormone 1 (gnrh1), cytochrome P450 aromatase (cyp19a1b), progesterone receptor (pgr), solute carrier family 6 (slc6a20) and cytochrome P450 1A1-like (cyp1a1). There were 216 down-regulated DEGs in this comparison, which included prolactin (prl), thyroid stimulating hormone (tshb), somatolactin-like (sl), glycoprotein hormones (cga) and pro-opiomelanocortin-like (pomc) (Table 3). Moreover, 362 DEGs were identified in the E-XY vs C-XY comparison, of which 52 were up-regulated (Fig. 2b, 3-4), such as vitellogenin-2-like (vtg2), pgr, gnrh1, cyp19a1b, zona pellucida sperm-binding protein 4-like (zp3) and cyp1a1. 172 down-regulated DEGs were observed in this comparison. These included potassium channel (kcnk18), WD40 repeat-containing protein (wd40), basic helix-loop-helix family (bhlhe41) and forkhead box protein O1-a (foxoa) (Table 3).
GO enrichment analysis of DEGs
In the C-XY vs C-XX, E-XX vs C-XX and E-XY vs C-XY comparisons, genes were mainly enriched in biological processes, followed by molecular function and cellular component GO terms (Fig. S2, 5). In the C-XY vs C-XX comparison, the DEGs were mainly significantly enriched in microtubule-based movement and movement of cell or subcellular component, in the biological process category. In the molecular function category, they were enriched in serine-type exopeptidase activity and exopeptidase activity (Fig. S2). In the E-XX vs C-XX comparison, DEGs were mainly enriched in response to oxygen-containing compound, response to drug and proteolysis, for biological process, in hormone activity, sequence-specific DNA binding and serine‑type peptidase activity, for molecular function, and in calcium ion binding, myosin complex and actin cytoskeleton, for cellular component. The up-regulated genes were mainly clustered in proteolysis, for the biological process category, extracellular region, for the cellular component category, and hormone activity, for the molecular function category. The down-regulated genes were mainly clustered in cell cycle arrest, for the biological process category, myosin complex, for the cellular component category, and protein kinase regulator activity, for the molecular function category (Fig. 5a). In the E-XY vs C-XY comparison, response to oxygen-containing compound, response to chemical and response to extracellular stimulus were highly represented for the biological process category. Integral component of plasma membrane, intrinsic component of plasma membrane and plasma membrane part were highly represented for the cellular component category. Sequence-specific DNA binding, heme binding and tetrapyrrole binding were highly represented for the molecular function category. The up-regulated genes were mainly clustered in response to extracellular stimulus and response to nutrient levels, in the biological process category, integral component of plasma membrane, for the cellular component category, and sequence-specific DNA binding, for the molecular function category. The down‑regulated genes were clustered in lipid transport and lipid localization in the biological process category (Fig. 5b).
KEGG enrichment analysis of DEGs
The most enriched KEGG pathways in the E-XX vs C-XX comparison were neuroactive ligand-receptor interaction, arachidonic acid metabolism, cytokine‑cytokine receptor interaction and the calcium signaling pathway (Fig. 6a). The KEGG pathways most enriched by down-regulated DEGs were neuroactive ligand‑receptor interaction, steroid hormone biosynthesis, retinol metabolism, calcium signaling pathway and GnRH signaling pathway. Eight pathways, which included neuroactive ligand-receptor interaction, notch signaling pathway, cytokine-cytokine receptor interaction, PPAR signaling pathway, steroid biosynthesis, calcium signaling pathway, metabolism of xenobiotics by cytochrome P450 and GnRH signaling pathway, were the most enriched by up-regulated DEGs. In the E-XY vs C-XY comparison, the most enriched KEGG pathways were tyrosine metabolism, phenylalanine metabolism, arachidonic acid metabolism and linoleic acid metabolism (Fig. 6b). The pathways most enriched by down-regulated DEGs were steroid hormone biosynthesis, retinol metabolism, PPAR signaling pathway, carbon metabolism, metabolism of xenobiotics by cytochrome P450, calcium signaling pathway and neuroactive ligand-receptor interaction. The tyrosine metabolism, phenylalanine metabolism, arachidonic acid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism, steroid biosynthesis, histidine metabolism, metabolism of xenobiotics by cytochrome P450, tryptophan metabolism and calcium signaling pathways were those most enriched by up-regulated DEGs.
The RNA-Seq data were completely consistent with the qPCR results (Fig. 7). In the control group, no significant difference in mRNA level of gnrh1, cyp1a1 and cyp19a1b was found between the XY and XX groups (p < 0.05). In the E2 treatment group, the expression of gnrh1, cyp1a1 and cyp19a1a, in XX and XY larvae, was significantly higher than in the control group. The level of cyp19a1b in the E2-treated XY group was significantly higher than in the E2-treated XX group. In the control group, no significant difference was observed in the expression level of prph, per1b, per3, cipc and ciart, between XY and XX larvae, whilst the level of nart1 was significantly lower in XY larvae than in XX larvae. The expression levels of bh1be, nr1d2, per1b, per3, cipc and ciart, in E2-treated larvae brains, were significantly lower than in the control group (p < 0.05), whilst the expression levels of arntl1a and cry1 were significantly higher than in the control group (p < 0.05).