Skip to main content
Download PDF

A detailed transcriptome study uncovers the epigenetic characteristics associated with Aromatase inhibitor-induced masculinization in Takifugu rubripes larvae gonads

BMC Genomics volume 26, Article number: 380 (2025) Cite this article

Abstract

Background

Takifugu rubripes is an economically valuable fish species in Asia. The implementation of all-male culture for T. rubripes is highly anticipated in aquaculture. Aromatase inhibitor (AI, letrozole) treatment was found to be an efficient method to induced masculinization in T. rubripes, as reported in our previous study. Here, to further explore the underlying regulation mechanism of AI-induced masculinization, a whole-transcriptome analysis comparing was conducted between AI-induced masculinized XX (AI-XX) gonads and control (Con) gonads in T. rubripes.

Results

In Con-XX/Con-XY comparison, 1,172 differential expression (DE) mRNAs, 129 DEmiRNAs, 210 DElncRNAs, and 4 DEcircRNAs were identified. In the Con-XX/AI-XX comparison, 1,329 DEmRNAs, 174 DEmiRNAs, 6 DEcircRNAs and 280 DElncRNAs were found. Con-XX/Con-XY and Con-XX/AI-XX comparisons shared 690 DEmRNAs, 50 DEmiRNAs, 3 DEcircRNAs, and 105 DElncRNAs. The analyses of protein-protein interaction (PPI) and competitive endogenous RNA (ceRNA) network identified interactions among these shared DERNAs. Kcnh2b, trim27, cnnm2b, reln, cckbra, pkd1l2, steap4, gsg1l, hamp, and foxg1c were predicted as the top ten of hub genes. miRNAs included miRNA-27 family and miRNA-489 family showed targeting relationship with hub genes. GO and KEGG functional enrichment analysis showed that the targeted genes were mainly enriched in GO:0065008 regulation of biological quality and TGF-beta signaling pathway. qPCR validation confirmed the differential expression of selected mRNAs, and ncRNAs.

Conclusions

This research comprehensively reveals the potential regulatory effects of ncRNAs on cellular motility, fate regulation, and hormonal regulation during gonadal masculinization in T. rubripes. It may provide significant insights into the regulation mechanisms underlying sex reversal in fish.

Peer Review reports

Background

Takifugu rubripes, a carnivorous fish belonging to the order Osteichthyes, Tetraodontiformes, Tetraodontoidei, and the family Tetraodontidae, is predominantly found in warm temperate zones near the seabed, particularly along the coasts of Japan, Korean Peninsula, and China [1]. Due to its delicate flavor, fine meat texture, and significant economic value, this species is regarded as favored choice for marine cultivation. There is no significant difference in growth traits between males and females, however, the gonadosomatic index can reach 20–30% during the breeding season. Female fugu possess ovaries that contain potent toxins, rendering them inedible, while male fish have non-toxic testes, which are considered a delicacy and are highly valued in the market. Despite the market preference for males, the natural sex ratio of T. rubripes remains nearly equal. Additionally, T. rubripes shows minimal sexual dimorphism throughout its lifespan, making it difficult to distinguish between the sexes except during the breeding season [2]. Therefore, achieving all-male populations is highly desirable, as this would significantly increase testis yield and enhance the economic returns for breeding enterprises.

Sex determination and differentiation are intricate developmental processes [3]. Sex determination is the process where gonads differentiate into testes or ovaries in sexually reproducing organisms [4]. Subsequently, sex differentiation entails the maturation of gonads into testes or ovaries and the development of primordial germ cells (PGCs) into sperm or eggs, resulting in secondary sexual characteristics [5]. Estrogen is crucial for sex differentiation. In teleost fish, the specific expression of aromatase and the synthesis of endogenous estrogen during the critical period of female sex differentiation are essential for ovarian differentiation, acting as natural inducers [6, 7, 8, 9]. During this crucial phase, inhibiting aromatase activity through aromatase inhibitors (AI) or knocking out the cyp19a1a gene results in reduced estrogen levels, leading to the complete transformation of genetically female fish into phenotypic males [10]. In T. rubripes, studies have shown that AI can induce masculinization during early developmental stages [11]. The molecular mechanisms through which aromatase inhibition affects gene expression and induces sex reversal remain insufficiently understood.

Advancements in high-throughput genomic technologies have revealed that the evolution of complex developmental processes is largely driven by the regulatory potential of non-coding regions of the genome [12]. Non-coding RNAs (ncRNAs), such as circular RNAs (circRNAs), long non-coding RNAs (lncRNAs), and, microRNAs (miRNAs) are essential in gene regulation, DNA replication, and mRNA translation [13, 14, 15, 16]. Abnormal ncRNA expression significantly impacts tissue development and function. In teleosts, many sexually dimorphic expressed ncRNAs have been identified through whole transcriptome analyses. For example, 111 differentially expressed miRNAs (DEmiRNAs) were identified between the testes and ovaries of tilapia (Oreochromis niloticus) [17]. Wang et al. (2018) identified female-biased miRNAs, such as miR-21-5p, miR-21-3p, and miR-462-5p, and male-biased miRNAs, including miR-9-3p, miR-103b-3p, and miR-7b in the gonads of Pelteobagrus fulvidraco [18]. LncRNAs and circRNAs modulate gene expression by competitively interacting with miRNAs through miRNA response elements (MRE), affecting target gene expression within regulatory networks [19]. In Chinese tongue sole (Cynoglossus semilaevis), both a circRNA from the sex-determining gene dmrt1 and the lncRNA AMSDT bind to miR-196, sharing miRNA response elements with gsdf to promote testis differentiation [20]. In T. rubripes, sex-biased differences in ncRNA expression have also been observed in the gonads [2]. However, the involvement of these ncRNAs in the AI-induced masculinization process remains unclear.

This study performed whole transcriptome analysis to examine non-coding gene expression in T. rubripes gonads after AI treatment. The goal was to provide insights into the potential function of competitive endogenous RNAs (ceRNAs) in sex determination and the epigenetic mechanisms governing male gonadal differentiation in this species. This research will enhance our understanding of gonad differentiation at the molecular level and offers a theoretical basis for achieving the cultivation of all-male populations.

Results

Histological analysis and genetic sex verification

In the control group, XX individuals at 60 dat showed the presence of ovarian cavities and a large number of primary oocytes. The gonads of XY individuals displayed normally testicular seminiferous tubules and spermatogenic cells. In the AI-treated group, XX T. rubripes did not developed ovarian cavity structures but exhibited structures similar to testes (Fig. 1).

Fig. 1
figure 1

Hematoxylin-eosin stained sections of Takifugu rubripes gonad at 20, 40, and 60 dat. OG, ogonia; OCA, ovarian cavity; SL, spermatogenic cysts; SG, spermatogonia; SC, spermatocytes. Scale bar, 100 μm

Full size image

Histological observations and genetic sex identification of the individuals in the control group showed that the male ratios of T. rubripes at 20, 40, and 60 dat were 47.36%, 55.26%, and 52.63%, respectively. Notably, the male ratio in the AI-treated group was 100%, with complete masculinization of individuals (Table 1).

Table 1 Gonadal phenotypes and sex ratio of Takifugu rubripes larvae at 20, 40 and 60 dat
Full size table

Global profiling of mRNAs from T. rubripes gonad

79,277,980 ~ 113,622,834 raw reads were obtained by whole-transcriptome RNA sequencing (Table S1). After quality control, approximately 99% of the raw reads remained as clean reads, with 94.08–95.54% successfully mapped to the reference genome. Total of 22,074 expressed genes were identified. The distribution of gene expression levels in each sample is shown in a violin plot (Fig. 2A). The result of principal component analysis (PCA) and correlation square heatmap plot is were shown in Fig. S1. 1,172 DEmRNAs were identified in the Con-XX/Con-XY comparison, including 501 upregulated and 671 downregulated in Con-XY group (Fig. 2B). These included cyp11c1 (4.70), dmrt1 (2.72), cyp19a1a (-8.11), and foxl2 (-5.25) (Table S2). 1,329 DEmRNAs were identified in the Con-XX/AI-XX comparison, including 480 upregulated and 849 downregulated in AI-XX group (Fig. 2B). These included dmrt1 (2.50), cyp19a1a (-7.61), and pgr (-1.91) (Table S3). 723 DEmRNAs were identified in the Con-XY/AI-XX comparison, including 289 upregulated and 434 downregulated in AI-XX group (Fig. 2B). These included slc12a1 (2.12), ctrb1 (1.59), stc (-2.35), and hgebf (-3.25) (Table S4).

The GO analysis revealed that DEGs in the comparisons of Con-XX/Con-XY, Con-XX/AI-XX, and Con-XY/AI-XX were primarily enriched in GO:0050896 response to stimulus and GO:0032501 multicellular organismal process under the biological processes (BP), as well as in GO:0044421 extracellular region part and GO:0005576 extracellular region under the cellular components (CC) (P < 0.05) (Fig. S2). The results of KEGG analysis suggested that the significantly enriched pathways of the DEGs in Con-XX/Con-XY, and Con-XY/AI-XX comparisons were neuroactive ligand-receptor interaction, and protein digestion and absorption. In Con-XY/AI-XX comparison, the DEGs were mainly enriched in protein digestion and absorption, and pancreatic secretion (Fig. S2).

Con-XX/Con-XY and Con-XX/AI-XX comparisons shared 690 DEGs (Fig. 2C). Most of these DEGs were all enriched in GO:0050896 response to stimulus, and GO:0032501 multicellular organismal process under BP (P < 0.05) (Fig. 2D). The significantly enriched pathways of the DEGs were TGF-beta (TGF-β) signal, protein digestion and absorption, cell adhesion molecules, and oocyte meiosis (Fig. 2E).

Fig. 2
figure 2

Bioinformatics analysis based on mRNA transcriptome sequencing results after treatment with aromatase inhibitor (AI). A: Violin plot of the expression level of genes. B: Statistic of the number of DEmRNAs. C: Venn diagram of differential expression genes in three groups. D and E: GO (D) and KEGG (E) enrichment analysis of 690 genes shared by Con-XX/Con-XY and Con-XX/AI-XX comparisons. F: PPI network of the hub genes. G: FPKM of steroid hormone synthesis enzyme

Full size image

Global profiling of miRNAs from T. rubripes gonad

In the small RNA library, 18,664,983 (Con-XX), 14,615,904 (Con-XY), and 15,849,675 (AI-XX) raw reads were generated in the Illumina sequencing. After filtering, clean tags were obtained from Con-XX (18,228,194 (97.66%)), Con-XY (14,236,079 (97.40%)), and AI-XX (15,596,209 (98.40%)) (Table S5). The result of PCA and correlation square heatmap plot is were shown in Fig. S4. Of these matched sequences, the identified miRNAs were divided into exist miRNAs (108, 11%), known miRNAs (489, 47%), novel miRNAs (434, 42%) (Fig. 3A). 129 DEmiRNAs (51 up-regulated, 78 down-regulated) were identified in Con-XX/Con-XY comparison. 174 DEmiRNAs (78 up-regulated, 96 down-regulated) were identified in Con-XX/AI-XX comparison. And, 182 DEmiRNAs (90 up-regulated, 92 down-regulated) were identified in Con-XY/AI-XX comparison (Fig. 3B).

Fig. 3
figure 3

Bioinformatics analysis based on miRNA transcriptome sequencing results after aromatase inhibitor (AI) treatment. A: Classification of the expression level of miRNA. B: Statistic of differential expression (DE) miRNAs. C: Venn diagram of DEmiRNAs in three comparison groups. D: Heat map of the DEmiRNAs shared by Con-XX/Con-XY and Con-XY/AI-XX comparisons. E and F: GO and KEGG enrichment analysis of target genes of DEmiRNAs shared by Con-XX/Con-XY and Con-XY/AI-XX comparisons

Full size image

The target genes of DEmiRNAs identified from Con-XX/AI-XX comparison were mapped into neuroactive ligand-receptor interaction, pancreatic secretion, protein digestion and absorption, and steroid hormone biosynthesis pathways (Fig. S5D). The target genes of DEmiRNAs identified from Con-XY/AI-XX comparison were mapped into protein digestion and absorption, pancreatic secretion, and cell adhesion molecules (CAMs) pathways (Fig. S5F).

Fifty DEmiRNAs were common to both Con-XX/Con-XY and Con-XX/AI-XX comparisons (Fig. 3C and D). Their target genes were most enriched in GO:0000087 mitotic M phase, and GO:0000003 reproduction under BP (P < 0.05). And they were enriched in metabolic, and aminoacyl-tRNA biosynthesis pathways (Fig. 3E, 3 F).

Trend analysis of miRNAs and their target genes

The trend analysis of mRNAs and miRNAs expression among Con-XX, Con-XY, and AI-XX identified 8 different expression patterns (Fig. S6, S7). Expression of 334 mRNA and 22 miRNAs were lower in Con-XY than that in Con-XX, and were lower in AI-treated XX than that in Con-XY (Fig. S6, S7, profile 0). Expression of 405 mRNAs and 31 miRNAs were lower in Con-XY than that in Con-XX, and no significant difference were observed between AI-treated XX and Con-XY (Fig. S4, S5, profile 1). Expression of 210 mRNA and 67 miRNAs were lower in Con-XY than that in Con-XX, and were higher during AI-treated XX than that in Con-XY (Fig. S4, S5, profile 2). Expression of 286 mRNA and 77 miRNAs were no significant difference in Con-XY than in Con-XX, and were lower during AI-treated XX than that in Con-XY (Fig. S6, S7, profile 3). Expression of 130 mRNA and 33 miRNAs were no significant difference in Con-XY than in Con-XX, and were higher AI-treated XX than that in Con-XY (Fig. S6, S7, profile 4). Expression of 191 mRNA and 25 miRNAs were higher in Con-XY than that in Con-XX, and were lower during AI-induced masculinization in XX than that in Con-XY (Fig. S6, S7, profile 5). Expression of 409 mRNAs and 30 miRNAs were higher in Con-XY than that in Con-XX, and no significant difference were observed between AI-treated XX and Con-XY (Fig. S6, S7, profile 6). Expression of 125 mRNA and 25 miRNAs were higher in Con-XY than that in Con-XX, and were higher AI-treated XX in XX than that in Con-XY (Fig. S6, S7, profile 7).

Notably, 127 DEmRNAs in profile 1 were predicted as the targets of 23 DEmiRNAs in profile 6 (Fig. 4A, 4B). These target genes were mainly enriched in GO:0009605 response to external stimulus, and GO:0006952 defense response (BP) (P < 0.05). And DEGs were mainly enriched in TGF-β pathway (Fig. 4C, 4D). The 145 DEmRNAs in profile 6 were predicted as the targets of 23 DEmiRNAs in profile 1 (Fig. 4E, 4F). These DEGs were mainly enriched in GO:0016021 integral component of membrane (CC), and GO:0006811 monoatomic ion transport (BP) ((P < 0.05). And they were mainly enriched in neuroactive ligand-receptor interaction pathway (Fig. 4G, 4H).

Fig. 4
figure 4

The regulatory relationship between DEmiRNAs and DEmRNAs. A: Degree distribution of DEmiRNAs in profile 6 and DEmRNAs in profile 1 interaction network. B: Interaction network diagram between DEmRNA in profile 1 and DEmiRNAs in profile 6. C and D: GO and KEGG (D) enrichment analysis of DEmRNAs in profile 1 targeted by DEmiRNAs in profile 6. E: Degree distribution of DEmiRNAs in profile 1 and DEmRNAs in profile 6 interaction network. F: Interaction network diagram between DEmRNAs in profile 6 targeted by DEmiRNAs in profile 1. G and H: GO (G) and KEGG (H) enrichment analysis of DEmiRNAs in profile 1 targeted by DEmRNAs in profile 6

Full size image

Global profiling of circRNA from T. rubripes gonad

According to their origins, the circRNAs were characterized as follows: annot_exons (152, 49%), one_exons circRNAs (30, 10%), antisense circRNAs (67, 22%), intronic circRNAs (18, 6%), intergenic circRNAs (7, 2%), and exon-intron circRNAs (35, 11%) (Fig. 5A). The result of PCA and correlation square heatmap plot is were shown in Fig. S8. Four DEcircRNAs (1 up-regulated, 3 down-regulated) were identified in Con-XX/Con-XY comparison. Six DEcircRNAs (2 up-regulated, 4 down-regulated) were identified in Con-XX/AI-XX comparison. Seven DEcircRNAs (4 up-regulated, 3 down-regulated) were identified in Con-XY/AI-XX comparison (Fig. 5B). The trend analysis of cicRNA expression among Con-XX, Con-XY, and AI-XX identified 5 different expression patterns (Fig. 5C). There were 3 circRNAs in profile 1, 1 circRNAs in profile 2, 2 circRNAs in profile 3, 3 circRNAs in profile 4, and 1 circRNAs in profile 5.

Fig. 5
figure 5

Bioinformatics analysis based on circiRNA transcriptome sequencing results. A: Classification of circRNAs. B: Statistic of DEcircRNAs. C: Statistic of DEcircRNA expression in 8 profiles. D: Venn diagram of DEcircRNAs in three comparison groups. E and F: GO enrichment analysis of target genes of DEcircRNAs identified from Con-XX/Con-XY (E) and Con-XX/AI-XXcomparison (F). G and H: KEGG enrichment analysis of target genes of DEcircRNAs identified from Con-XX/Con-XY (G) Con-XX/AI-XX comparison (H)

Full size image

Three DEcircRNAs were common to both Con-XX/Con-XY and Con-XX/AI-XX comparisons (Fig. 5D). In Con-XX/Con-XY and Con-XX/AI-XX comparison, the target genes of DEcircRNA were mainly significantly enriched in GO:0004620 phospholipase activity (MF), and GO:0007264 small GTPase mediated signal transduction (BP) (P < 0.05) (Fig. 5E, 5 H). In Con-XY/AI-XX comparison, target genes of DEcircRNA were mainly enriched in GO:0005515 protein binding (MF), and GO:0007264 small GTPase mediated signal transduction under BP (P < 0.05) (Fig. S9A). The KEGG enrichment analysis found that target genes of DEcircRNA were mainly enriched in autophagy pathway in the Con-XX/Con-XY, and Con-XX/AI-XX comparisons (Fig. 5F and G). And, the target genes of DEcircRNA were mainly enriched in cholesterol metabolism, and Retinol metabolism pathways in the Con-XY/AI-XX comparison (Fig. S9B).

Global profiling of lncRNAs from T. rubripes gonad

The study categorized 4,498 lncRNAs into several types: sense lncRNAs (109), antisense lncRNAs (1,086, 24%), bidirectional lncRNAs (367, 8%), intronic lncRNAs (125, 3%), intergenic lncRNAs (2,133, 47%), and other types (678, 15%) (Fig. 6A). The result of PCA and correlation square heatmap plot is were shown in Fig. S10. 210 DElncRNAs (133 upregulated, 77 downregulated) were identified in Con-XX/Con-XY comparison. Tow hundred and eighteen DElncRNAs (109 upregulated, 109 downregulated) identified in Con-XX/AI-XX comparison. 203 DElncRNAs (85 upregulated, 118 downregulated) were identified in Con-XY/AI-XX comparison (Fig. 6B). The trend analysis of lncRNA expression among Con-XX, Con-XY, and AI-XX identified 8 different expression patterns (Fig. 6C). There were 32 lncRNAs in profile 0, 37 lncRNAs in profile 1, 42 lncRNAs in profile 2, 62 lncRNAs in profile 3, 41 lncRNAs in profile 4, 61 lncRNAs in profile 5, 85 lncRNAs in profile 6, and 19 lncRNAs in profile 7 (Fig. 6C).

Fig. 6
figure 6

Bioinformatics analysis based on lncRNAs transcriptome sequencing results. A: Classification of lncRNAs. B: Statistic of DElncRNA. C: Statistic of DElncRNAs in 8 profiles. D: Venn diagram of DElncRNAs in three comparison groups. E and F: GO enrichment(E) and KEGG pathway (F) analysis of target genes of DElncRNAs shared by Con-XX/Con-XY and Con-XX/AI-XX comparison.

Full size image

The enriched GO categories of genes overlapping antisense lncRNAs were analyzed. The target genes of lncRNAs were mainly enriched in GO:0043269 regulation of ion transport under BP (P < 0.05) (Fig. S11A). The KEGG enrichment analysis found that target genes of lncRNAs were mainly enriched in focal adhesion, and purine metabolism pathways (Fig. S11B).

One hundred and five DElncRNAs were common to both Con-XX/Con-XY and Con-XX/AI-XX comparisons (Fig. 6D). The target genes of DElncRNAs were mainly enriched in GO:0000902 cell morphogenesis, and GO:0000278 mitotic cell cycle (CC) (P < 0.05) (Fig. 6E). The KEGG enrichment analysis found that target genes of DElncRNAs were mainly enriched in drug metabolism-other enzymes, and nucleocytoplasmic transport pathways (Fig. 6F).

CeRNA regulatory network of DEncRNAs on DEmRNAs during sex determination and differentiation

To explore the regulatory relationship of ncRNAs on mRNAs during sex determination and differentiation. The identification of ceRNAs was achieved through a filtering process from DERNAs. This network involved 180 mRNAs, 30 miRNAs, 314 lncRNAs, 6 circRNAs, forming a total of 1,265 predicted miRNA-lncRNA pairs, and 47 predicted miRNA-circRNA pairs, 418 predicted miRNA-mRNA pairs (Table S6).

Top 10 DEmRNAs (ncbi 101064060 potassium voltage-gated channel (kcnh2b), ncbi_101072370 E3 ubiquitin-protein ligase TRIM27 (trim27), ncbi 101064931 cyclin and CBS domain divalent metal cation transport mediator 2b (cnnm2b), ncbi_101074793 reelin (reln), ncbi_101076498 cholecystokinin B receptor a (cckbra), ncbi_101071016 polycystic kidney disease protein 1-like 2 (pkd1l2), ncbi_101074043 STEAP4 metalloreductase (steap4), ncbi_101076343 gsg1-like (gsg1l), ncbi_115250474 hepcidin-like (hamp), ncbi_ 101076354 forkhead box G1c (foxg1c)) for connectivity as hub genes were selected with higher degrees in networks (Fig. 7A). DElncRNAs with targeting relationships included MSTRG.11670.1, XR_008888387.1, XR_003891557.1, XR_003891558.1, XR_003886513.1, XR_003886696.1, XR_003887942.1, XR_003889001.1, XR_003889001.1, XR_003885753.1, and XR_965171.2. DEmiRNAs with targeting relationships included fru-miRNA-135a, fru-miRNA-27b, miRNA-1788-y, fru-miRNA-27e, novel-m0261-5p, miRNA-144-y, novel-m0262-5p, novel-m0303-5p, novel-m0302-5p, miRNA-27z, novel-m0279-5p, novel-m0175-5p, miRNA-31-x, novel-m0360-5p, and miRNA-489-y.

Fig. 7
figure 7

The ceRNA network construction and function annotation. A: The Sankey diagram for ceRNA network. B: GO enrichment analysis of mRNAs in ceRNA network. C: The pathway enrichment analysis of mRNAs in ceRNA network. D: Overview of TGF-β signaling pathway. E: The KEGG pathways related to TGF-β signaling pathway

Full size image

Validation of DEmRNAs, DEmiRNAs, DEcircRNAs and DElncRNAs by qPCR

As shown in Fig. 8, in the control group, individuals from Con-XX exhibit significantly lower expression levels of dmrt1, gsdf, cyp11c1, and XR_003886586 compared to Con-XY individuals. Following AI treatment, these genes were notably up-regulated in XX individuals relative to Con-XX (P < 0.05). Conversely, Con-XX individuals exhibited significantly elevated expression levels of foxl2 and cyp19a1a compared to Con-XY individuals. They were down-regulated in AI-XX individuals compared to Con-XX (P < 0.05). And there was no significance difference in the expression levels of foxl2 and cyp19a1a between Con-XY and AI-XX (P > 0.05). Furthermore, the expression levels of fru-miR-153b, fru-miR-11,980-z, and circR-00205 in Con-XX individuals were significantly higher than those in Con-XY (P < 0.05). And, the expression level was down-regulated in AI-XX individuals compared to Con-XX (P < 0.05). The levels of fru-miR-11,980-z expression in AI-XX was higher compared to Con-XY (P < 0.05). There was no significant difference in the levels of fru-miR-153b, and circR-00205 between Con-XY and AI-XX (P > 0.05).

Fig. 8
figure 8

Relative expression of randomly selected mRNAs, lncRNAs, circRNA, and miRNAs. Different lowercase indicated statistically significant differences (P < 0.05)

Full size image

Discussion

Aromatase plays a crucial role in the regulation of estrogen synthesis by catalyzing the key final step in the process and facilitating the hydroxylation of androstenedione to estrone and testosterone to estradiol [28]. Therefore, the inhibition of aromatase activity could lead to decrease the level of estrogen and induce the XX individual gonad differentiating into structures similar to testes in teleost, including T. rubripes [21]. Letrozole is a popular aromatase inhibitor with high absorption by gastrointestinal tract and specificity in inhibiting aromatase activity [22]. To better understand the molecular mechanisms of AI (letrozole)-induced masculinization in T. rubripes, whole-transcriptome data from Con-XX, Con-XY, and AI-XX individuals were compared.

There were 690 DEGs, 50 DEmiRNAs, 3 DEcircRNAs, and 105 DElncRNAs shared by Con-XX/Con-XY and Con-XX/AI-XX comparisons, indicating the involvement of ncRNAs in testis differentiation in T. rubripes. To explore potential regulatory effects for ncRNAs, ceRNA network analyses was conducted to show the interactions among these DEmRNAs, DEmiRNAs, and DElncRNAs. miRNAs can silence post-transcriptional gene by binding to partially complementary MREs on target RNAs [23, 24]. Due to shared MREs, various RNA transcripts function as ceRNAs, creating intricate regulatory networks facilitated by miRNAs [25]. Thus, the initial step in investigating the ceRNA regulatory network involves predicting the target genes of DEmiRNAs. Among the DEmiRNAs, fru-miRNA-27b, fru-miRNA-27e, miRNA-27z, novel-0252-5p, fru-miRNA-489, miRNA-489-y, novel-0317-3p, miRNA-153-x, novel-0388-3p, fru-miRNA-202, novel-0279-5p, and fru-miRNA-135a, targeted many DEmRNAs/lncRNAs in ceRNA regulatory network. The miRNA-27 family has been found to have the regulatory relationship with many DEmRNAs/lncRNAs including a2ml1, bpi, ccn2, c4, crlf1, cyp19a1, egr1, eif4ebp1, hamp, gsg1l, groa, steap4, skor2, slit2, midnb, olfm3, rab26, ptafr, pde11a, pvalb2, plppr5, pkd1l2, taar8b, trim27, rgs1, reln, rhoGAP44, kcnj10, and kcnh2. The miRNA-489 family was predicted to target hamp, kcnh2, pde11a, osbpl6, rnf23, gja3, gpr84, gsg1l, arhgap44, adm, cckar, cebpd, egr1, foxg1, furin, prkca, plppr5, sh3bp4a, slc7a8, and sybu genes. In the ceRNA regulatory network, the degree of connectivity of a specific RNA molecule (lncRNA/mRNA/circRNA) is positively correlated with the number of miRNA molecules it interacts with, indicating its potential regulatory capacity. Nodes with high connectivity, known as hub genes which are often of significant biological importance.

We focused on the top 10 DEmRNAs with the highest connectivity to identify key hub genes (kcnh2b, trim27, cnnm2b, reln, cckbra, pkd1l2, steap4, gsg1l, hamp, foxg1c) and predicted the functions in the network. Kcnh2b encodes the voltage-gated potassium channel responsible for the rectifier potassium current [26]. This channel is composed of a symmetrical tetramer structure, with each subunit consisting of voltage-sensing domains and a pore domain responsible for ion passage. In Potassium channels dysfunction can have pro-proliferative effects closely related to calcium signaling [27]. In mice, ncRNAs are produced from genomic locations where cis-regulatory elements were found associated with the expression of kcnh2b and two neighboring mRNAs, nos3 and abcb8 [28]. The protein encoded by steap4 is metalloreductase and found in the golgi apparatus. It acts as a metal reductase and is able to use NAD (+) as a receptor to reduce Fe (3+) to Fe (2+) and Cu (2+) to Cu (1+). STEAP4 is involved in iron and copper homeostasis [29]. Similarly, hamp is a key regulator in iron homeostasis [30]. Cnnm2b, encoded protein is a divalent metal cation transporter. It potentially contributes to magnesium homeostasis by facilitating epithelial Mg2+ transport [31]. GSG1L can negatively regulates single-channel conductance and calcium permeability of recombinant AMPARs [32, 33]. DElnRNAs have also been found to regulate DEmiRNAs that target these genes, but these predicted targeting relationships require further validation. These results suggested the significant involvement of ncRNAs in the regulation of ion channel protein-related genes during AI-induced masculinization in T. rubripes. Ion channels are essential in cellular signaling pathways that directly influence cell motility and fate. In vertebrates, ion channels play a pivotal role in sex differentiation, particularly when environmental factors induce this process. For instance, in Alligator mississippiensis, the TRPV4 ion channel supports a steady influx of Ca2+, which enables the Ca2+/CaM complex to initiate molecular cascades, leading to the upregulation of amh and sox9, key genes involved in testis differentiation [34]. In Trachemys scripta elegans, Ca2+ influx promotes STAT3 phosphorylation and suppresses kdm6b transcription, resulting in ovary formation [35]. Fish exhibit a remarkable plasticity in sex determination mechanisms, allowing them to adjust their sex phenotypes in response to environmental changes [36]. Environmental factors, such as temperature and social surroundings, can influence ion channel activity, which in turn regulates intracellular ion concentrations [34]. These changes can impact hormone synthesis and secretion, thereby mediating sex differentiation at a molecular level. T- and L-type Ca2+ channels, along with voltage-dependent Na+ and K+ channels, have been identified in human luteinizing granulosa cells [37]. These channels are critical for maintaining cellular excitability and facilitating steroid hormone synthesis, acting as essential physiological components in ovarian function [38, 39]. Furthermore, the interaction between ion channels and steroid hormones is evident through pregnenolone sulfate, a steroid metabolite that modulates various ion channels and transporters. This neurosteroid not only influences ion channel activity but also acts as an activating ligand for certain ion channels [40, 41, 42]. This interplay between ion channels and steroid synthesis underscores the potential of ion channels as targets for regulating sex steroid levels [38, 40]. Therefore, understanding these mechanisms offers valuable insights into environmentally induced sexual plasticity, and also presents novel approaches for inducing sex reversal in teleost fish [38, 40].

Additionally, GO/KEGG functional enrichment analysis was conducted on mRNA within the ceRNA regulatory network to identify significantly enriched gene functions and pathways. The analysis identified significant enrichment of target genes in the transforming growth factor β (TGF-β) pathway, which is essential for biological processes such as tissue and organ formation, reproductive development, and the regulation of cellular growth, proliferation, and differentiation [43, 44, 45, 46, 47]. In this study, treatment with AI induced a decrease in p15, resulting in cell-cycle arrest at the G1 phase. Moreover, the expression level of the upstream-regulated Myc proto-oncogene protein gene was up-regulated and transcription factor E2F4/5 was down-regulated. In general, the transition from the G1 phase to the S phase is essential for regulating eukaryotic cell proliferation [48]. The G1-to-S phase transition includes a ‘commitment point’, marking when a cell commits to the cell cycle and continues progression without environmental signals [49, 50]. Further investigation is needed to understand the causal link between G1 phase cell-cycle arrest and cell fate alteration. The TGF-β signaling pathway has also been found to influence Smad-independent pathways, such as Erk, SAPK/JNK, and p38 MAPK pathways [51, 52]. Activation of Rho GTPase leads to the downstream activation of target proteins, which in turn facilitate cytoskeletal molecular remodeling associated with cellular extension, growth regulation, and cytoplasmic division [53, 54]. The findings of this study indicate that ncRNAs may play a regulatory role in the expression of genes within the TGF-β signaling pathway. This result has highlighted the significance of epigenetic mechanisms in vertebrate sex determination.

Conclusion

This research provides significant insights into the molecular mechanisms of gonadal masculinization in T. rubripes, utilizing comprehensive transcriptome analysis and bioinformatics tools. The study highlights the pivotal role of ncRNAs in the complex interplay of hormonal regulation, cell signal transduction, and TGF-β signaling pathways during gonadal differentiation. Among them, the regulatory functions of ncRNAs on ion channels genes with implications for cellular motility, and fate regulation and the regulatory mechanisms of environmental factors on these ncRNAs warrant further investigation. These findings lay the groundwork for future comprehensive studies aimed at uncovering the molecular mechanisms underlying sex reversal and gonadal differentiation.

Materials and methods

Ethics statement

Tissue collection adhered to the Care and Use of Experimental Animals guidelines from the Key Laboratory of Environment Controlled Aquaculture (Dalian, China), complying with Chinese ethical standards, regulations, and laws.

Treatment of T. rubripes larvae

T. rubripes used in this study were collected from Dalian Fugu Aquatic Product Co., Ltd., Dalian, China. The fish were reared in 300 L round tanks under natural environmental conditions with the water temperature of 20 ± 2 °C from 20 day after hatching (dah) to 45 dah. Then, they were randomly assigned to the six differential tanks (400/tanks) with three replications for control and AI-treatment groups after acclimatization. From 45 to 105 dah, the larvae in AI treatment group were fed with a commercial diet incorporating AI (500 µg/g diet weight). According to previous study by Yan et al. (2021a) [11], AI (Sigma-Aldrich, St. Louis, MO, USA) was first dissolved in absolute ethanol to create a storage solution with a concentration of 5000 µg/ml. The AI solution was diluted with 75% ethanol and subsequently applied to the commercial diet at a concentration of 500 µg/g diet weight.

Tissue sampling

For the purpose of sampling collection, larvae were anesthetized by immobilizing in ice water. At 65 dah, 40 larvae were randomly chosen from each tank, and their gonad and muscle tissues were dissected and stored separately. The muscle tissues of each fish were store in a single tube, and store at -20 °C until genetic sex identification. Each gonad was preserved in 100 µl RNAlater (Thermo Fisher Scientific, Baltics, Vilnius, Lithuania), and store at -80 °C until RNA-seq. In addition, a total of 228 larvae (38 per group) were dissected to obtain trunk containing the gonads at 20, 40 and 60 day after treatment (dat). The trunks were immersed in a 4% paraformaldehyde solution for 24 h, and then transferred to a 70% ethanol solution for preservation until histological analysis.

Gonadal histologic features, sex verification and comprehensive RNA extraction

The gonads underwent dehydration through increasing ethanol concentrations before being embedded in paraffin. Subsequently, 8 μm tissue sections were prepared and stained using hematoxylin-eosin. The Imajer. A2 optical microscope from ZEISS were conducted to perform microscopic analysis and photography. The leftover paraffin-blocked trunks were employed for DNA extraction utilizing the TIANamp FFPE DNA kit (Tiangen, Beijing, China) in order to identify the genetic sex of each individual after histological evaluation.

Genetic sex identification was performed by extracting DNA from ethanol-preserved tissue using the TIANfor each individual was identified using a SNP marker on the amhr2 gene through PCRamp Marine Animals DNA kit (Tiangen, Beijing, China). Genetic sex [55]. RNA was extracted from five same-gender gonads for per sample using the Qiagen RNeasy Mini Kit (Qiagen, Valencia, CA, USA).

Whole-transcriptome sequencing analysis

The quality evaluating, and purification of total RNA were performed according to previous study [10]. Transcriptome libraries were constructed for control groups (Con-XX and Con-XY) and AI-treatment groups (AI-XX). The RNA library and a small RNA library were conducted for subsequent sequencing on the Illumina HiSeq Xten platform according to previous study [56]. To ensure data reliability, the acquisition of high-quality data was accomplished by eliminating adapters and low-quality bases. Reads aligning to rRNA were identified and removed using Bowtie2 (version 2.2.8) [57]. The remaining reads were utilized for alignment to a reference genome (https://ftp.ncbinlm.nih.govgenomes/al/GCF/901/000/725/ GCF901000725.2fTakRub1.2/) using Hisat2 v2.0.5. The transcriptome was analyzed with Stringtie software (version 1.3.4) [58]. Clean tags were compared with the miRBase database (Release 22) to identify known miRNAs for the studied species. The alignment of miRNAs with other species to identify known miRNAs. Novel miRNA candidates were identified based on their genomic locations and predicted hairpin structures using mirdeep2 software. miRNA expression levels were quantified and normalized using transcripts per million (TPM). Differential expression analysis of miRNAs was performed using edgeR, applying criteria of fold change (FC) ≥ 2 and P-value < 0.05.

CircRNAs were identified by analyzing anchor reads aligned in a reversed orientation (head-to-tail), indicative of circRNA splicing, using find_circ [59]. LncRNAs were selected by intersecting the results of both non-coding potential analyses. These lncRNAs were then categorized into five classes lncRNAs based on their relative location to protein-coding genes: intergenic, bidirectional, intronic, antisense, and sense overlapping lncRNAs. FPKM (fragments per kilobase of transcript per million mapped reads) was used to measure lncRNA expression levels and variations. On the other hand, circRNAs were scaled to RPM (Reads Per Million mapped reads). Differential expression of ncRNAs was analyzed using DESeq2, with criteria of FDR < 0.05 and absolute FC ≥ 2. The target genes of differentially expressed ncRNAs were analyzed for enrichment to assess their links to Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Protein-protein interactome network analysis

A protein-protein interaction network was constructed to investigate the mechanisms of direct or indirect interactions among DEGs, aiming to identify key regulatory genes and their roles in signaling pathways during T. rubripes gonad masculinization. DE transcript sequence IDs were retrieved from Ensemble and analyzed using the STRING interactome database (http://string-db.org) for network analysis. Cytoscape 3.9.1 [60] was used to visualize the PPI network identified in STRING v11 [61]. Hub genes were selected based on an interaction score of at least 0.40 and a minimum interaction degree of 20 in the original PPI.

ceRNA network analysis

The ceRNA network for DERNAs was constructed according to the predicted targeting relationship among DEcircRNA, DElncRNAs, DEmiRNAs, and DEmRNAs, and the correlation expressed, three conditions were followed, (a) The correlations between DEncRNAs and DEmRNAs expression levels were assessed using Pearson’s correlation coefficient (PCC), with pairs exhibiting PCC > 0.9 designated as co-expressed lncRNA/circRNA-mRNA pairs. (b) The hypergeometric cumulative distribution function test was used to identify shared miRNA sponges among target genes, retaining only gene pairs with P < 0.05 (c) The candidate ceRNA exhibits comparable levels of DEmiRNAs enrichment.

The circRNA-miRNA-mRNA network was visualized using Cytoscape software (version 3.5.1) to visualize these relationships [62]. Sankey plot, also known as Sankey diagram, was used to depict the regulatory relationships between top 10 mRNA/lncRNA and miRNA. The network of circRNA/lncRNA-miRNA-mRNA correlations was constructed using Cytoscape (v3.6.0).

qPCR

The expression validation of mRNAs, miRNA, lncRNAs, and circRNAs was performed according to previous reports [11]. Purified miRNA was reverse transcribed with Maxima Reverse Transcriptase and miRNA-specific stem-loop RT primers (Table S7). The 2X SG Fast qPCR Master Mix (High Rox) from BBI, China, was used for reverse transcription of mRNA. miRNA and gene expression were normalized using housekeeping genes U6 and ef1, respectively. Relative expression levels were determined using the 2−ΔΔCt method. Reactions were performed in triplicate, with primers designed using Primer Premier 5.0 and detailed in Table S1. Statistical analysis utilized IBM SPSS 22.0 software (IBM, Armonk, NY, USA). The expression levels of mRNAs, miRNA, lncRNA, and circRNA in the Con-XX, Con-XY, and AI-XX groups were assessed through a one-way ANOVA (P < 0.05). Subsequently, the Duncan test was performed to compare the means.

Data availability

No datasets were generated or analysed during the current study.

References

  1. Aparicio S, Chapman J, Stupka E, Putnam N, Chia JM, Dehal P, et al. Whole-genome shotgun assembly and analysis of the genome of Fugu Rubripes. Science. 2002;297:1301–10.

    Article  CAS  PubMed  Google Scholar 

  2. Yan H, Liu Q, Jiang J, Shen X, Zhang L, Yuan Z, et al. Identification of sex differentiation-related microRNA and long non-coding RNA in Takifugu rubripes gonads. Sci Rep. 2021;11:7459.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Strüssmann CA, Nakamura M. Morphology, endocrinology, and environmental modulation of gonadal sex differentiation in teleost fishes. Fish Physiol Biochem. 2002;26:13–29.

    Article  Google Scholar 

  4. Weber C, Capel B. Sex determination without sex chromosomes. Philos Trans R Soc B. 2021;376:20200109.

    Article  CAS  Google Scholar 

  5. Smolko AE, Shapiro-Kulnane L, Salz HK. The H3K9 methyltransferase SETDB1 maintains female identity in Drosophila germ cells. Nat Commun. 2018;9:4155.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Lau ESW, Zhang Z, Qin M, Zhang Q, Ge W. Knockout of zebrafish ovarian aromatase gene (cyp19a1a) by TALEN and CRISPR/Cas9 leads to all-male offspring due to failed ovarian differentiation. Sci Rep. 2016;6:36487.

    Article  Google Scholar 

  7. Zhang X, Guan G, Li M, Zhu F, Liu Q, Naruse K, et al. Autosomal gsdf acts as a male sex initiator in the fish medaka. Sci Rep. 2016;27:19738.

    Article  Google Scholar 

  8. Wu GC, Chang CF. Primary males guide the femaleness through the regulation of testicular dmrt1 and ovarian cyp19a1a in protandrous black porgy. Gen Comp Endocrinol. 2018;261:198–202.

    Article  CAS  PubMed  Google Scholar 

  9. Hou M, Feng K, Luo H, Jiang Y, Xu W, Li YM et al. Multi-locus gene editing effectively knocked out cyp19a1a and foxl2 in Monopterus albus, a hermaphroditic fish. Aquaculture. 2023;565:739130.

  10. 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 

  11. Yan H, Shen X, Jiang J, Zhang L, Yuan Z, Wu Y et al. Gene expression of Takifugu rubripes gonads during AI- or MT-induced masculinization and E2-induced feminization. Endocrinology. 2021;162.

  12. Mattick JS, Makunin IV. Non-coding RNA. Hum Mol Genet. 2006;15:17–29.

    Article  Google Scholar 

  13. Hofmann P, Boon RA. Non-coding RNA enhances cardiac development. J Mol Cell Cardiol. 2014;76:205–7.

    Article  CAS  PubMed  Google Scholar 

  14. Boon RA, Jaé N, Holdt L, Dimmeler S. Long noncoding RNAs: from clinical genetics to therapeutic targets? J Am Coll Cardiol. 2016;67:1214–26.

    Article  CAS  PubMed  Google Scholar 

  15. Rothschild G, Basu U. Lingering questions about enhancer RNA and enhancer transcription-coupled genomic instability. Trends Genet. 2017;33:143–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018;172:393–407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Tao W, Sun L, Shi H, Cheng Y, Jiang D, Fu B, et al. Integrated analysis of miRNA and mRNA expression profiles in tilapia gonads at an early stage of sex differentiation. BMC Genomics. 2016;17:1–13.

    Article  Google Scholar 

  18. Wang P, Wang L, Yang J, Luan PX, Zhang XF, Kuang YY, et al. Sex-biased miRNAs of yellow catfish (Pelteobagrus fulvidraco) and their potential role in reproductive development. Aquaculture. 2018;485:73–80.

    Article  CAS  Google Scholar 

  19. Li JQ, Yang J, Zhou P, Le YP, Gong ZH. The biological functions and regulations of competing endogenous RNA. Yi Chuan. 2015;37:756–64.

    CAS  PubMed  Google Scholar 

  20. Tang L, Huang F, You W, Poetsch A, Nóbrega RH, Power DM, et al. ceRNA crosstalk mediated by ncRNAs is a novel regulatory mechanism in fish sex determination and differentiation. Genome Res. 2022;32:1502–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Ijiri S, Kaneko H, Kobayashi T, Wang DS, Sakai F, Paul-Prasanth B, et al. Sexual dimorphic expression of genes in gonads during early differentiation of a teleost fish, the Nile tilapia Oreochromis Niloticus. Biol Reprod. 2008;78:333–41.

    Article  CAS  PubMed  Google Scholar 

  22. Lee VC, Ledger W. Aromatase inhibitors for ovulation induction and ovarian stimulation. Clin Endocrinol (Oxf). 2011;74:537–46.

    Article  CAS  PubMed  Google Scholar 

  23. Mayr C, Bartel DP. Widespread shortening of 3’UTRs by alternative cleavage and polyadenylation activates oncogenes in cancer cells. Cell. 2009;138:673–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Derda AA, Pfanne A, Bär C, Schimmel K, Kennel PJ, Xiao K, et al. Blood-based microRNA profiling in patients with cardiac amyloidosis. PLoS ONE. 2018;13:e0204235.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Sanchez-Mejias A, Tay Y. Competing endogenous RNA networks: tying the essential knots for cancer biology and therapeutics. J Hematol Oncol. 2015;8:30.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Trudeau MC, Warmke JW, Ganetzky B, Robertson GA. HERG, a human inward rectifier in the voltage-gated potassium channel family. Science. 1995;269:92–5.

    Article  CAS  PubMed  Google Scholar 

  27. Chen X, Xue B, Wang J, Liu H, Shi L, Xie J. Potassium channels: a potential therapeutic target for Parkinson’s disease. Neurosci Bull. 2018;34:341–8.

    Article  CAS  PubMed  Google Scholar 

  28. Van den Boogaard M, Van Weerd JH, Bawazeer AC, Hooijkaas IB, Van de Werken HJG, et al. Identification and characterization of a transcribed distal enhancer involved in cardiac kcnh2 regulation. Cell Rep. 2019;28:2704–14.

    Article  PubMed  Google Scholar 

  29. Scarl RT, Lawrence CM, Gordon HM, Nunemaker CS. STEAP4: its emerging role in metabolism and homeostasis of cellular iron and copper. J Endocrinol. 2017;234:123–34.

    Article  Google Scholar 

  30. Li Z, Hong WS, Qiu HT, Zhang YT, Yang MS, You XX, et al. Cloning and expression of two hepcidin genes in the mudskipper (Boleophthalmus pectinirostris) provides insights into their roles in male reproductive immunity. Fish Shellfish Immunol. 2016;56:239–47.

    Article  CAS  PubMed  Google Scholar 

  31. Franken GAC, Seker M, Bos C, Siemons LAH, Van der Eerden BCJ, Christ A et al. Cyclin M2 (Cnnm2) knockout mice show mild hypomagnesaemia and developmental defects. Sci Rep. 2021;11:8217.

  32. McGee TP, Bats C, Farrant M, Cull-Candy SG. Auxiliary subunit GSG1L acts to suppress calcium-permeable AMPA receptor function. J Neurosci. 2015;35:16171–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Gu X, Mao X, Lussier MP, Hutchison MA, Zhou L, Hamra FK, et al. GSG1L suppresses AMPA receptor-mediated synaptic transmission and uniquely modulates AMPA receptor kinetics in hippocampal neurons. Nat Commun. 2016;7:10873.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Yatsu R, Miyagawa S, Kohno S, Saito S, Lowers RH, Ogino Y, et al. TRPV4 associates environmental temperature and sex determination in the American alligator. Sci Rep. 2015;18:18581.

    Article  Google Scholar 

  35. Weber C, Zhou Y, Lee JG, Looger LL, Qian G, Ge C, et al. Temperature-dependent sex determination is mediated by pSTAT3 repression of Kdm6b. Science. 2020;368:303–6.

    Article  CAS  PubMed  Google Scholar 

  36. Robert HP, Willey LL, Jones MT, Akre TSB, King DI, Kleopfer J, et al. Is the future female for turtles? Climate change and wetland configuration predict sex ratios of a freshwater species. Glob Chang Biol. 2023;29:2643–54.

    Article  Google Scholar 

  37. Darboux I, Lingueglia E, Champigny G, Coscoy S, Barbry P, Lazdunski M. dGNaC1, a gonad-specific amiloride-sensitive na + channel. J Biol Chem. 1998;273:9424–9.

    Article  CAS  PubMed  Google Scholar 

  38. Mayerhofer A, Kunz L. Ion channels of primate ovarian endocrine cells: identification and functional significance. Expert Rev Endocrinol Metab. 2006;1:549–55.

    Article  CAS  PubMed  Google Scholar 

  39. Rotgers E, Jørgensen A, Yao HH. At the crossroads of fate-somatic cell lineage specification in the fetal gonad. Endocr Rev. 2018;39:739–59.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Harteneck C. Pregnenolone sulfate: from steroid metabolite to TRP channel ligand. Molecules. 2013;18:12012–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Morales-Lázaro SL, González-Ramírez R, Rosenbaum T. Molecular interplay between the sigma-1 receptor, steroids, and ion channels. Front Pharmacol. 2019;10:419.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Hamad M, Bajbouj K, Taneera J. The case for an estrogen-iron axis in health and disease. Exp Clin Endocrinol Diabetes. 2020;128:270–7.

    Article  CAS  PubMed  Google Scholar 

  43. Attisano L, Wrana JL. Signal transduction by the TGF-β superfamily. Science. 2002;296:1646–7.

    Article  CAS  PubMed  Google Scholar 

  44. Massagué J. TGF-β signalling in context. Nat Rev Mol Cell Biol. 2012;13:616–30.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Moustakas A, Heldin CH. The regulation of TGFbeta signal transduction. Development. 2009;136:3699–714.

    Article  CAS  PubMed  Google Scholar 

  46. Monsivais D, Matzuk MM, Pangas SA. The TGF-β family in the reproductive tract. Cold Spring Harb Perspect Biol. 2017;9:a022251.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Mullen AC, Wrana JL. TGF-β family signaling in embryonic and somatic stem-cell renewal and differentiation. Cold Spring Harb Perspect Biol. 2017;9:a022186.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Bertoli C, Skotheim JM, de Bruin RA. Control of cell cycle transcription during G1 and S phases. Nat Rev Mol Cell Biol. 2013;14:518–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Umen JG, Goodenough UW. Control of cell division by retinoblastoma protein homolog in Chlamydomonas. Genes Dev. 2001;15:1652–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Gao Q, Li L, Zhang Q. Monotropein induced apoptosis and suppressed cell cycle progression in colorectal cancer cells. Chin J Integr Med. 2024;30:25–33.

    Article  CAS  PubMed  Google Scholar 

  51. Tsoi S, Zhou C, Grant JR, Pasternak JA, Dobrinsky J, Rigault P, et al. Development of a porcine (Sus scofa) embryo-specific microarray: array annotation and validation. BMC Genomics. 2012;13:370.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Ungefroren H, Reimann J, Konukiewitz B, Braun R, Wellner UF, Lehnert H, et al. RAC1b collaborates with TAp73α-SMAD4 signaling to induce biglycan expression and inhibit basal and TGF-β-driven cell motility in human pancreatic cancer. Biomedicines. 2024;12:199.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Zhang Y, Tang HM, Liu CF, Yuan XF, Wang XY, Ma N, et al. TGF-β3 induces autophagic activity by increasing ROS generation in a NOX4-dependent pathway. Mediators Inflamm. 2019;31:3153240.

  54. Marshall-Burghardt S, Migueles-Ramírez RA, Lin Q, El Baba N, Saada R, Umar M, et al. Excitable rho dynamics control cell shape and motility by sequentially activating ERM proteins and actomyosin contractility. Sci Adv. 2024;10:eadn6858.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Kamiya T, Kai W, Tasumi S, Oka A, Matsunaga T, Mizuno N, et al. A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, Takifugu rubripes (fugu). PLoS Genet. 2012;8:e1002798.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Shen X, Yan H, Hu M, Zhou H, Wang J, Gao R, et al. The potential regulatory role of the non-coding RNAs in regulating the exogenous estrogen-induced feminization in Takifugu rubripes gonad. Aquat Toxicol. 2024;273:107022.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47:607–13.

    Article  Google Scholar 

  62. Huang BS, Yang MH, Wang PH, Li HY, Chou TY, Chen YJ. Estrogen-induced angiogenesis and implantation contribute to the development of parasitic myomas after laparoscopic morcellation. Reprod Biol Endocrin. 2016;14:64.

    Article  Google Scholar 

Download references

Acknowledgements

We are deeply grateful to Jieming Jiang, Liyu Bai and other students from our research team, who were involved in the rearing process and sampling of specimens.

Funding

This work was supported by the Youth Program of National Natural Science Foundation of China (31902347); Dalian Science and Technology Innovation Support Plan Project (2023RJ010); the General Program of Liaoning Natural Science Foundation (2022-MS-351, 2023-MSLH-001, 2023-MSLH-003); the General Project of Education Department of Liaoning Province (LJKMZ20221092); the Research and Application of Germplasm Creation and Efficient Breeding Technologies in Tiger Puffer Takifugu rubripes (2021RT07); the China Agriculture Research System of MOF and MARA (CARS-47); the Liaoning Province Science and Technology Joint Plan (Technology Research and Development Program Project) (2024011259-JH2/1026).

Author information

Authors and Affiliations

  1. College of Fisheries and Life Science, Dalian Ocean University, Dalian, Liaoning, 116023, China

    Xufang Shen, Hongwei Yan, Mingtao Hu, Huiting Zhou & Qi Zhang

  2. Key Laboratory of Environment Controlled Aquaculture (Dalian Ocean University), Ministry of Education, Dalian, 116023, China

    Xufang Shen, Hongwei Yan, Mingtao Hu, Rui Gao & Qi Liu

  3. College of Marine Science and Environment Engineering, Dalian Ocean University, Dalian, Liaoning, 116023, China

    Rui Gao & Qi Liu

  4. Dalian Tian Zheng Co., Ltd, Dalian, Liaoning, China

    Qunwen Sun

Authors
  1. Xufang Shen
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Hongwei Yan
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Mingtao Hu
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Huiting Zhou
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Qi Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Rui Gao
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Qi Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Qunwen Sun
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

S.X. dealt with the experimental materials, analyzed the data and wrote the manuscript. Y.H. designed the research and reviewed the manuscript. H.M. and Z.H. performed the experiments, and analyzed the data. Z.Q. and G.R. performed the experiments. L.Q. provided support through funding. All authors have read and approved the manuscript. S.Q. provided support through rearing technology.

Corresponding author

Correspondence to Hongwei Yan.

Ethics declarations

Ethics approval and consent to participate

Tissue collection adhered to the Care and Use of Experimental Animals guidelines from the Key Laboratory of Environment Controlled Aquaculture (Dalian, China), complying with Chinese ethical standards, regulations, and laws.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

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

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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 http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Shen, X., Yan, H., Hu, M. et al. A detailed transcriptome study uncovers the epigenetic characteristics associated with Aromatase inhibitor-induced masculinization in Takifugu rubripes larvae gonads. BMC Genomics 26, 380 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11375-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11375-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us