The potential regulatory role of non-coding RNAs in mifepristone-induced masculinization in Takifugu rubripes gonads

Abstract

Background

The regulatory roles of non-coding RNAs (ncRNAs) during sex differentiation in teleosts have received widespread attention recently. Mifepristone (RU486, a progesterone antagonist), which acts as an endocrine disruptor, can affect reproduction and sex differentiation in teleosts.

Results

The expression of ncRNAs in the gonads of tiger puffer (Takifugu rubripes) during RU486 (500 µg/g diet) induced masculinization process was examined by RNA-sequencing. A total of 4,381 long non-coding RNAs (lncRNAs), 309 circular RNAs (circRNAs), and 1,020 microRNAs (miRNAs) were identified. The expression of 41 differentially expressed (DE) lncRNAs and 20 DE miRNAs, which showed sexual dimorphic expression patterns in genetic female gonads in the control group (C-XX) vs. genetic male gonads in the control group (C-XY), were altered in genetic female gonads in the RU486 treated group (RU-XX). The genes targeted by DE ncRNAs were mainly enriched in sex-related pathways, such as calcium signaling, ovarian steroidogenesis, and cortisol synthesis and secretion. The results of co-expression and competing endogenous RNA (ceRNA) network analysis indicated that miRNAs (e.g., miR-205-z and fru-miR-122) and lncRNAs (including XR_003890915.1 and XR_003885862.1) may have pivotal roles, and lncRNAs (including XR_003890295.1, MSTRG.11750.1, and XR_003888827.1) may act as miRNA sponges, involved in the competition between miRNAs and sex-related genes during tiger puffer masculinization process. Dual luciferase reporter assay results identified that ovarian steroidogenesis related gene hsd17b1 is a downstream target of fru-miR-122. The expression of 4 lncRNAs, 4 circRNAs, and 6 miRNAs were validated by qPCR, indicating the accuracy and dependability of RNA-Seq.

Conclusions

This study provided the evidence that ncRNAs may participate in RU486-induced masculinization in T. rubripes, and may enhance our understanding of the regulatory network of sex differentiation in fugu.

Background

Typically, sex-determining genes serve as the initial trigger to activate the differential expression of downstream genes in the gonadal differentiation, which initiate the transformation of bipotential gonads either into testes or ovaries in some vertebrates [1]. In teleosts, sex-determining genes, such as amhbY in Northern pike (Esox Lucius) [2], hsd17b1 in Seriola fishes [3], dmrt1 in Chinese tongue sole (Cynoglossus semilaevis) [4], amhy in Patagonian pejerrey (Odontesthes hatcheri) [5], and dmy in Medaka (Oryzias latipes) [6], typically expressed at early periods of gonadal development, and marked the beginning of the sexual differentiation process. Afterward, conserved sex differentiation-related genes such as gsdf, dmrt1, foxl2, and cyp19a1a were progressively expressed, and bipotential gonads began to differentiate into either ovaries or testes [1]. Differ from higher vertebrates, the sexual fate of most teleost fish is not solely determined by inherited genetic elements, but can be influenced by environmental factors (e.g., density, water temperature) and steroids [7, 8]. Among the steroid hormones, in addition to estrogen and androgen, progesterone also plays a significant role in the sex differentiation of teleosts [9, 10]. For example, treatment of fish with mifepristone (RU486), a progesterone antagonist, has been shown to affect the expression of steroid hormone receptor genes, reproduction in zebrafish (Danio rerio) [11, 12], and sex differentiation in Nile tilapia (Oreochromis niloticus) [13].

The function of non-coding RNAs (ncRNAs), including circular RNAs (circRNAs), long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and piwi-interacting RNAs (piRNAs), etc., in sex differentiation have recently attracted more attention [14,15,16,17]. For instance, the ptk2β, which may regulate cyp19a1a expression to reprogram sex differentiation, was found to be regulated by lncRNAs MSTRG.12,998 and MSTRG.38,036 in sex reversal gonads of swamp eel (Monopterus albus) [18]. MiR-34b/c, which can increase germ cell proliferation and reduce apoptosis in Amur sturgeon (Acipenser schrenckii), has an interaction with the 3′-UTR of AR mRNA [19]. miRNA sponges are RNA transcripts with various high-affinity binding sites that can bind and isolate specific miRNAs, impeding miRNAs from interacting with their target mRNAs [20]. Previous studies revealed that circRNAs and lncRNAs may regulate sex differentiation by acting as miRNA sponges. For instance, during the sex change of swamp eel, circSnd1 regulates the expression of the female-related gene foxl2 by acting as mal-miR135b/c sponge [21]. Further evidence was discovered in Chinese tongue sole, a lncRNA (AMSDT) and a circRNA (circdmrt1) both shared miRNA cse-miR-196 response elements with the male-related gene gsdf, and the testicular differentiation was facilitated by their binding to cse-miR-196 [22]. These findings suggest that the sex differentiation of vertebrates may be regulated by ncRNAs. Furthmore, several studies have shown that interfering with steroid synthesis can affect the expression of ncRNAs in vertebrates [18, 23, 24]. In mammals, dietary androgen or estrogen can affect lncRNAs expression during the development of phallus in the marsupial tammar wallaby (Macropus eugenii) [23]. In embryonic chicken, 17β-estradiol (E₂) injection caused the feminization of male gonads and reduced the expression of miR-202, which may regulate the expression of foxl2, cyp19a1a, sox9, and dmrt1 [24]. In fish, treatment with 17α-methyltestosterone (MT) inhibited the expression of lncRNAs and increased methylation in developing swamp eel gonads [18]. These findings suggest that the expression patterns of ncRNAs may be related to exogenous steroid hormone or steroid hormone antagonist-induced fish sex reversal.

The tiger puffer (Takifugu rubripes), also known as fugu, is a gonochoristic marine fish with an XX/XY sex-determining system [25], widely spread throughout the coastal areas of China, the Korean peninsula, and Japan. A missense single nucleotide polymorphism (SNP) within the anti-Müllerian hormone receptor type II gene (amhr2) has been identified as the master sex-determining polymorphism site in tiger puffer [26]. Therefore, the genetic sex can be assessed by amhr2, allowing the discrimination of genotypic and phenotypic sex to identify the sex reversal in tiger puffer. Consequently, tiger puffer represents an excellent biological species for the investigation of the mechanisms underlying sex determination and differentiation. Additionally, several genes including gsdf, dmrt1, cyp11c1, foxl2, and cyp19a1a showed sex dimorphism in the gonad of tiger puffer at early periods of sex differentiation [27]. Our previous study discovered that administering RU486 during sex differentiation can alter the expression of several sex-related genes (e.g., cyp19a1a, foxl2, dmrt1, cyp11c1, and hsd17b1) and steroid hormone receptor genes (e.g., androgen receptor AR and progesterone receptor PGR) in gonads, and induced 100% masculinization in genetic XX tiger puffer [28]. However, it remains unclear whether RU486 directly affects the expression patterns of ncRNAs during the masculinization process of tiger puffer. Due to the essential role of ncRNAs in sex differentiation, based on our previous study [28], RNA-sequencing was used here to obtain comprehensive expression patterns of lncRNAs, circRNAs, and miRNAs in the gonad of tiger puffer leave during the RU486-induced masculinization process. The study aimed to (i) illustrate their functional interactions and (ii) provide insights into the regulatory crosstalk between ncRNAs and mRNAs in the sex differentiation and sex reversal of tiger puffer.

Methods

Fish treatment and sampling

The tiger puffer larvae at 25 days after hatching (dah) were provided by Dalian Tianzheng Industrial Co., Ltd. (Dalian, China). After a five-day acclimation period, the larvae were randomly divided into the control group and the RU486 treatment group, each group with three replicates. For the control group, larvae were fed with a commercial diet six times per day from 30 to 90 dah. For the RU486-treated group, larvae were treated with RU486 (500 µg/g diet) six times per day from 30 to 90 dah. The water temperature was kept at 21–22℃ during the experiment. At 90 dah, the gonads tissues of larvae (300 larvae per treatment) were rapidly dissected and placed in a single 1.5 ml tube with 100 µl of RNAlater reagent (Thermo Fisher Scientific, Lithuania), then the samples were deep-frozen in liquid nitrogen and then stored at -80℃ until RNA extraction. Additionally, 10 tissue samples containing the gonads (30 larvae per treatment) from each group were collected and immersion-fixed in 4% paraformaldehyde for 24 h, and then transferred to 70% ethanol. Simultaneously, a part of the fin from each larva was collected and stored at -80 ℃ for genomic DNA extraction.

Genetic sex-identification, gonads histological analysis, and RNA extraction

The genetic sex of each larva was identified using a sex-linked marker in amhr2, as described in our previous study [27]. For histological analysis, the samples in 70% ethanol were dehydrated in ascending ethanol concentrations and then embedded in paraffin. The cross-sections were cut at 4 to 6 μm thickness and subsequently stained with hematoxylin and eosin. Then, the sections were observed and captured using an optical microscope (Nikon ECLIPSE Ci-L, Japan).

For RNA extraction, the gonads of larvae with the same genetic sex were pooled together (5 gonads/ genetic sex in each tank), generating 4 samples, including C-XX, C-XY, RU-XX, and genetic female gonads in the RU486 treated group (RU-XY). The total RNA of gonads was extracted using Qiagen RNeasy Micro Kit (Qiagen, USA) according to the provided instructions. The DNAfree Kit (Qiagen, USA) was utilized to eliminate any remaining genomic DNA from the total RNA samples. Prior to placing the RNA samples in the − 80 ℃ freezer, the RNAse inhibitor (TaKaRa, Japan) was included in the RNA samples. The Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and NanoDrop ND-1000 spectrophotometer (Thermo Scientific, USA) were utilized to evaluate the quality and concentration of each RNA sample, respectively.

Libraries construction and sequencing

One microliter of RNA was used as input material for RNA-sequence. For lncRNA and circRNA library construction, rRNAs were removed from total RNA samples with a Ribo-Zero™ rRNA removal kit (Illumina, USA). The remaining ncRNAs and mRNAs were fragmented into 200 to 500 nt utilizing a fragmentation buffer. Subsequently, random primers were employed to reverse transcribe the first strand cDNA. And then, the second strand of cDNA was synthesized by dNTP (dUTP instead of dTTP), DNA polymerase I, buffer, and RNase H. Then, the QiaQuick PCR extraction kit (Qiagen, The Netherlands) was utilized for purifying the synthesized double-stranded cDNA fragments. Subsequently, the fragments were end-repaired, appended with poly (A), and connected to Illumina sequencing adapters. Next, the second-strand cDNA was treated with Uracil-N-Glycosylase. The size of the resulting fragments was then determined using agarose gel electrophoresis. The sorted fragments underwent PCR enrichment to generate 4 cDNA libraries, including C-XX, C-XY, RU-XX, and RU-XY. Sequencing of the PCR products was conducted on an Illumina HiSeqTM 4000 instrument (Gene Denovo Biotechnology Co., China). The sequencing data has been deposited in the NCBI database with the accession numbers PRJNA866845.

For miRNA library construction, RNA ranging from 18 to 30 nt was selectively enriched via polyacrylamide gel electrophoresis. The 3’ adapters were ligated to enriched 36 to 48 nt RNAs, followed by ligation of the 5’ adapters to the RNAs. The PCR amplification was performed to reverse transcribed the ligation products. The PCR products were then enriched to 140 to 160 bp to generate 4 cDNA libraries (C-XX, C-XY, RU-XX, and RU-XY), and sequenced by Illumina HiSeq Xten (Gene Denovo Biotechnology Co., Guangzhou, China). The sequencing data have been deposited in the NCBI database under the following accession numbers: PRJNA866053.

Assemble and identification of the lncRNAs, circRNAs, and miRNAs

Following the sequencing of lncRNAs and circRNAs, the raw data underwent additional filtration using fastp (version 0.18.1) to obtain high-quality clean reads [29]. Subsequently, reads that aligned with the rRNA database were eliminated. The rRNA-depleted reads from each sample were then aligned to the T. rubripes reference genome (RefSeq: GCF_000180615.1) using HISAT2 [30]. Upon alignment with the reference genome, the reads that could be mapped to the genomes were utilized for mRNA and lncRNA identification, while the reads that could not be mapped were utilized for circRNA identification.

To identify lncRNA, transcripts were reconstructed by Stringtie (version 1.3.4) [31]. Cuffcompare was then used to detect novel transcripts from all reconstructed transcripts. The protein-coding possibility of novel transcripts was assessed by CPC (version 0.9-r2) and CNCI (version 2) in default parameters [32, 33]. The novel transcripts with nonprotein-coding potential predicted by both CNCI and CPC were selected as lncRNAs. Based on their location relative to protein-coding, the lncRNAs were then classified into five classes.

For the identification of circRNAs, 20 nt were extracted from both ends of the unmapped read and then aligned with the T. rubripes reference genome to find unique anchor positions within the splice sites. CircRNAs splicing was indicated by anchor reads that aligned in reversed orientation (head-to-tail). The reads were analyzed using find_circ to identify circRNAs. To ensure the entire read aligns, the anchor alignments were extended, with breakpoints flanked by GU/AG splice sites. A circRNA candidate was considered valid if it had support from a minimum of two unique back-spliced reads in at least one sample. Finally, circRNAs were classified into six categories based on their origins.

After sequencing miRNAs, dirty reads containing adapters or low-quality bases were removed from the raw reads obtained from the sequencing machines to obtain clean tags. The clean tags were then aligned to small RNAs in the Rfam (Release 11.0) and GenBank (Release 209.0) databases to identify and remove rRNA, tRNA, snRNA, snoRNA, and scRNA. To remove tags mapped to exons, introns, and repeat sequences, the clean tags were further aligned to the T. rubripes genome. To identify existing miRNAs, the remaining clean tags were checked against the miRbase database (Release 20). The remaining clean tags also aligned to other species to identify the known miRNAs which were not included in the miRbase database. All the unannotated tags were aligned to the reference genome of T. rubripes. The novel miRNAs were identified based on their positions in the genome and their predicted hairpin structures using the mirdeep2 software.

Differential expression analysis

The expression abundance and variations of lncRNAs were determined by StringTie software, and the fragment per kilobase of transcript per million mapped reads (FPKM) values of the transcripts of each sample were obtained. For the quantification of circRNAs, reads per million mapped reads (RPM) were calculated from back-spliced junction reads. The miRNA expression levels were calculated and normalized to transcripts per million (TPM) based on their expression in each sample.

DESeq2 [34] was used for the identification of differentially expressed lncRNAs (DE lncRNAs), while the edgeR package (version 3.12.1) [35] was used for the identification of differentially expressed circRNAs (DE circRNAs) and differentially expressed miRNAs (DE miRNAs). The DE lncRNAs and DE circRNAs were identified with a false discovery rate (FDR) parameter below 0.05 and an absolute fold change ≥ 2, while the DE miRNAs were identified with a fold change ≥ 2 and a P value < 0.05.

Targets genes prediction of lncRNAs and miRNAs, and GO and KEGG enrichment analysis of target genes

For lncRNA-mRNA association analysis, RNAplex (version 0.2) was utilized to reveal the interaction between antisense DE lncRNA and differentially expressed genes (DEGs) identified in our previous study [28]. BLAST2GO (version 4.0.7) [36] was utilized to predict the cis-target of DE lncRNAs on coding genes in less than 100 kb upstream/downstream. To analyze DE lncRNA trans-regulation, LncTar was used to identify target genes of lncRNAs by analyzing the expression correlation between lncRNAs and differentially expressed protein-coding genes [37]. The trans-target DEGs in each comparison group were then subjected to Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The target genes of DE miRNAs were predicted by TargetScan (Version 7.0) and Miranda (Version 3.3a). The target DEGs of DE miRNA were then subjected to the GO functions and the KEGG pathway enrichment analysis.

The target genes were mapped to the Gene Ontology database (http://www.geneontology.org/) for GO enrichment analysis. The number of genes included in each of the terms was calculated. The T. rubripes genome was used as a reference set to define significantly enriched GO terms via hypergeometric tests. Under the T. rubripes genome background, the target genes were mapped to the KEGG database [38] to identify significantly enriched signal transduction pathways or metabolic pathways. Calculated P-values underwent FDR correction, with FDR ≤ 0.05 as a threshold for significant GO terms and KEGG pathways.

Construction of the CeRNA interaction network

All DE lncRNAs, DE miRNAs, and eight sex-related DEGs were used for ceRNA interaction network construction. The Spearman rank correlation coefficient (SCC) was utilized to assess the correlation of expression between pairs of lncRNA-miRNA, circRNA-miRNA, or mRNA-miRNA. Pairs with SCC < -0.7 were considered co-expressed pairs, and mRNA, lncRNA, and circRNA were identified as miRNA target genes. The Pearson correlation coefficient (PCC) was utilized to assess the correlation of expression between lncRNA-mRNA or circRNA-mRNA pairs. Pairs with a PCC greater than 0.9 were considered co-expressed. Both mRNA and lncRNA in this pair were targeted and exhibited negative co-expression with a common miRNA, or both the mRNA and circRNA in this pair were targeted and exhibited negative co-expression with a common miRNA. Consequently, only the gene pairs with a P-value less than 0.05 were selected. The ceRNA network was visualized using Cytoscape [39].

Validation of lncRNAs, circRNAs, and miRNAs

To validate the RNA-sequence data, 4 lncRNAs, 4 circRNAs, and 6 miRNAs were further randomly selected for qPCR. The elongation factor 1-alpha (ef1α) was utilized as a reference gene for lncRNAs and circRNAs, while U6 was utilized as a reference gene for miRNAs. The primers for lncRNAs, circRNAs, miRNAs, and reference genes were shown in Table S1. The PCR mixture contained 5 µl 2XSYBR Green qPCR Master Mix (TaKaRa Bio, Dalian, China), 0.1 µl upstream primers, 0.1 µl downstream primers, and 10 µl cDNAs. All reactions were carried out with three technical replicates. PCR was carried out at 95℃ for 10 min, followed by 40 cycles of 95℃ for 15 s and 60℃ for 30 s, followed by signal acquisition. The melting curve was generated in the following way: 95℃ for 15 s, followed by 60 ℃ for 1 min, and then heating temperature from 60 ℃ to 95 ℃ by 1 ℃ every 5 s, and followed by signal acquisition. The relative expression of lncRNAs, miRNAs, and circRNAs was calculated by 2^−ΔΔCT method. Statistical significance of relative expression between control and treatment groups was assessed using One-way ANOVA followed by Duncan’s test (SPSS, version 22.0; IBM, Chicago, IL, USA), with a P value less than 0.05 was defined to be significant.

Dualluciferase reporter assay

The targeting relationship between fru-miR-122 and hsd17b1 was performed using dual luciferase reporter assay. The fru-miR-122 mimics (Sense 5’-3’ UGGAGUGUGACAAUGGUGUUUG, Antisense 5’-3’ AACACCAUUGUCACACUCCAUU) and the negative control (mimics-NC) (Sense 5’-3’ UUCUCCGAACGUGUCACGUTT, Antisense 5’-3’ ACGUGACACGUUCGGAGAATT) were designed and synthesized from GenePharma (Shanghai, China). The hsd17b1 3′-UTRs wild-type (WT) or mutant-type (MUT) sequence containing the predicted fru-miR-122 binding site was synthesized, and then ligated into the dual-luciferase pmirGLO reporter vector which was then cotransfected into HEK-293T cells together with fru-miR-122 mimics or mimics-NC using Lipofectamine 2000 (Invitrogen, USA). Dual-luciferase activities was measured after 48 h of transfection using the Dual-Luciferase Assay Kit (Promega, Madison, USA). The mock-vehicle (MV) group and the positive control (PC) group were used to test the reliability of the experiment. Two-way ANOVA followed by Sidak’s test (SPSS, version 22.0; IBM, Chicago, IL, USA) was used to assess the statistical significance of relative luciferase activity, with a P value less than 0.05 was defined to be significant.

Results

Paraffin section observation of gonadal tissues

In control group, the gonad of larvae contained ovary, testis, and undifferentiation gonad tissues accounted for 26.7%, 50%, and 23.3% (n = 30), respectively (Fig. 1A). The ovaries exhibited plenty of oogonia and a few oocytes (Fig. 1B), and the testis in the control group exhibited proliferating spermatogonia (Fig. 1C). In RU486 treated group, the gonad of larvae contained ovary, testis, and undifferentiation gonad tissues accounted for 0, 80%, and 20%, respectively (Fig. 1A). The testis-like structure, similarly to the control testis, was observed in the gonad of larvae in RU486-treated group (Fig. 1D).

Fig. 1

Histological analysis of gonad in Takifugu rubripes at 90 dah and classification of non-coding RNA. (A) Percentage of gonadal phenotype (n = 30). (B) The gonad of control tiger puffer with testis. (C) The gonad of control tiger puffer with ovary. (D) The gonad of RU486 treated tiger puffer with ovary. (E) Classification of lncRNAs, circRNAs, and miRNAs. OG: oogonia; OC: oocyte; OCA: ovarian cavity; SG: spermatogonium; SL: spermatogenic cysts. bar = 50 μm

Identification and characterization of lncRNAs, circRNAs, and miRNAs

Based on the RNA-sequence library, a total of 82,979,962, 75,444,460, 79,686,400, and 77,746,290 raw data from the lncRNA and circRNA libraries were obtained from C-XX, C-XY, RU-XX, and genetic male gonads in the RU486 treated group (RU-XY) groups, respectively. After filtering out low-quality reads, adapters, and reads with over 10% unidentified nucleotides, 82,824,060 (C-XX), 75,298,904 (C-XY), 79,515,874 (RU-XX), and 77,576,450 (RU-XY) clean data were obtained, respectively. The clean data with removed rRNA were aligned to the reference genome with mapping rates ranging from 96.32% (RU-XX) to 97.11% (C-XX) (Table S2). Overall, from the small RNA library, a total of 9,205,295, 11,049,535, 9,384,791, and 13,760,955 clean reads were acquired from C-XX, C-XY, RU-XX, and RU-XY groups, respectively. After removing reads with 5’ connectors, reads without 3’ connectors, low-quality reads, reads without insertion fragments, reads with insertion fragments less than 18 nt, and reads containing polyA (> 70% of the bases in a read are A), 9,017,776 (97.9629%), 10,836,424 (98.0713%), 9,171,818 (97.7307%), and 13,418,129 (97.5087%) clean tags were obtained from the C-XX, C-XY, RU-XX, and RU-XY groups, respectively (Table S3).

In total, 4,142 known lncRNAs and 239 novel lncRNAs were found in four libraries from tiger puffer gonads. The 4,381 lncRNAs were classified into different categories according to their relative location to nearby protein-coding genes. These categories include sense lncRNAs (65, 1.48%), antisense lncRNAs (1,080, 24.65%), intronic lncRNAs (122, 2.78%), bidirectional lncRNAs (359, 8.19%), intergenic lncRNAs (2,097, 47.87%), and other undefined lncRNAs (658, 15.02%). According to origins of circRNAs, 309 circRNAs identified in four libraries were classified into six categories, including annot-exons circRNAs (182, 58.90%), antisense circRNAs (48, 15.53%), exon-intron circRNAs (30, 9.71%), intergenic circRNAs (12, 3.88%), intronic circRNAs (8, 2.60%), and one-exon circRNAs (29, 9.39%). A total of 1,020 miRNAs were identified based on small RNA-sequence data, including 439 known miRNAs, 109 existing miRNAs, and 472 novel miRNAs (Fig. 1E).

Identification of DE lncRNAs, DE circRNAs, and DE miRNAs

In total, 344 DE lncRNAs, 14 DE circRNAs, and 322 DE miRNAs in the gonads of tiger puffer were revealed in comparisons between six pairs of groups. Compared to the C-XX group, the expression of 55 DE lncRNAs, 2 DE circRNAs, and 55 DE miRNAs were up-regulated, and the expression of 51 DE lncRNAs, 2 DE circRNAs, and 26 DE miRNAs were down-regulated in C-XY group. Compared to the C-XX group, the expression of 94 DE lncRNAs, 0 DE circRNAs, and 101 DE miRNAs were up-regulated, and the expression of 50 DE lncRNAs, 2 DE circRNAs, and 33 DE miRNAs were down-regulated in RU-XX group. Additionally, the expression of 62 DE lncRNAs, 2 DE circRNAs, and 62 DE miRNAs were up-regulated, and the expression of 54 DE lncRNAs, 4 DE circRNAs, and 43 DE miRNAs were down-regulated in RU-XY group compared to C-XX group (Fig. 2). Furthermore, compared to the C-XY group, the RU-XY group exhibited 112 DE lncRNAs, 3 DE circRNAs, and 98 miRNAs, while the RU-XX group had 137 DE lncRNAs, 2 DE circRNAs, and 151 DE miRNAs. A total of 120 DE lncRNAs, 5 DE circRNAs, and 122 DE miRNAs were identified in the RU-XY group compared to the RU-XX group. The expression level of 41 lncRNAs and 20 miRNAs, which showed sexual dimorphic expression patterns in C-XX vs. C-XY groups, were altered in RU486-induced masculinization XX gonads of T. rubripes (Fig. 2B and F).

Fig. 2

DE lncRNAs, DE circRNAs, and DE miRNAs in C-XX, C-XY, RU-XX, and RU-XY group. The pairwise comparison of the numbers of DE lncRNAs (A), DE circRNAs (C), and DE miRNAs (E) between different groups. The Venn diagram displaying the number of DE lncRNAs (B), DE circRNAs (D), and DE miRNAs (F) identified between the comparison groups of C-XX vs. C-XY, C-XX vs. RU-XX, and C-XY vs. RU-XX

Series test of the cluster of DE lncRNAs, DE circRNAs, and DE miRNAs

In the series test of the clusters of DE lncRNAs (Fig. 3A), profile 18 and profile 11 were two significant profiles (P < 0.05). In the C-XX vs. C-XY and C-XX vs. RU-XX comparison groups, the different expression levels of DE lncRNAs in profile 18, profile 17, profile 3, profile 0, profile 19, profile 16, profile 1, and profile2. There were 16, 43, 23, 7, 7, 16, 6, and 9 DE lncRNAs in profile 18, profile 17, profile 3, profile 0, profile 19, profile 16, profile 1, and profile2, respectively (Table S4). No significant profile was found in the series test of the clusters of DE circRNAs (Fig. 3B). In the series test of the clusters of DE miRNAs (Fig. 3C), profile 11, profile 18, profile 12 were identified as three significant profiles (P < 0.05). A total of 108 miRNAs in profile 0- profile 3 and profile 16- profile 19 may be crucial for XX tiger puffer masculinization (Table S5).

Fig. 3

DE lncRNAs (A), DE circRNAs (B), and DE miRNAs (C) expression profiles

Target DEGs of DE lncRNAs, and DE miRNAs

The associations between the above DE lncRNAs and DE mRNA were predicted and analyzed, encompassing antisense, cis, and trans interactions, to elucidate the putative functional roles of lncRNA through mRNA targeting. Eleven pairs of antisense lncRNA-mRNA were generated with 11 DE lncRNAs and 10 mRNAs (Table S6). A total of 97 cis-lncRNA-mRNA pairs were generated to investigate the role of DE lncRNA in local gene regulation. Notably, the dataset comprised 77 DE lncRNAs and 72 mRNAs (Table S7). In total 74,438 trans-lncRNA-mRNA pairs, comprising 344 DE lncRNAs and 3,595 mRNAs, were generated according to the expression level analysis. The target genes of several DE lncRNAs were found to be associated with sex-related genes. Specifically, a total of 10, 18, 21, 22, 26, and 38 DE lncRNAs targeted cyp11c1, dmrt1, dmrt3, gsdf, amhr2, and nanos2, respectively. Additionally, 20, 20, and 19 DE lncRNAs were found to target cyp19a1a, foxl2, and hsd17b1 respectively (Table S8).

Based on trans interactions, the co-expression network between the DE lncRNA in profile 0 - profile 3 and profile 16 - profile 19 (mentioned in 3.4 part) and 27 sex-related genes (cyp11c1, dmrt1, gsdf, cyp19a1a, foxl2, etc.) (Table S9) was constructed, and 279 mRNA-lncRNA pairs which including 92 lncRNAs and 23 mRNAs were generated (Figure S1, Table S10). Among the DE lncRNAs, DE lncRNAs in profile2 (XR_003885863.1, XR_003888756.1, XR_003888884.1, XR_003889260.1, XR_003889849.1, XR_003890064.1, XR_003890751.1, XR_003890915.1, and XR_003891384.1) all targeted female-related genes (cyp19a1, foxl2, hsd17b1, esrrg, zp3, and zp4), only expressed in C-XX gonads (Table S4). Six lncRNAs (MSTRG.13468.1, MSTRG.17819.1, XR_003886055.1, XR_003886586.1, XR_003887080.1, and XR_003887947.1) in profile 17 that were highly expressed in C-XY and RU-XX gonads than in C-XX gonads, were predicted to target dmrt1. Four lncRNAs (MSTRG.3717.2, XR_003885862.1, XR_003888314.1, and XR_003890080.1) that were highly expressed in C-XY and RU-XX gonads than in C-XX gonads, were all predicted to target gsdf and amhr2. Three DE lncRNA, including XR_966292.2 (LOC105419695), XR_172564.3 (LOC101071485), and XR_003887473.1 (LOC115248715), all targeted cyp11c1 and pgr. Four DE lncRNA, including XR_003889315.1 (LOC105418423), MSTRG.1846.7, XR_003890318.1 (LOC105417209), and MSTRG.8300.1, all targeted pgr, foxl2, cyp19a1a, and hsd17b1.

A total of 83 potential miRNA-mRNA interactions between 17 sex-related mRNAs and 63 miRNAs were detected (Table S11). The results show that cyp19a1a were targeted by novel-m0182-5p, novel-m0183-5p, miR-499-x, and miR-27-z, while foxl2 were targeted by miR-14-y, miR-93-z, miR-205-z, and miR-451-z. Six miRNAs, including miR-205-x, novel-m0327-3p, novel-m0326-3p, novel-m0340-5p, novel-m0341-5p, and novel-m0342-5p targeted gsdf, one miRNA, miR-730-x targeted dmrt1. Four miRNAs, including miR-8159-x. fru-miR-217, novel-m0083-3p, and miR-8159-z, was predicted to target esr1, while esr2 was targeted by miR-19-y. In addition, fru-miR-122 and miR-122-x were predicted to target hsd17b1 (Fig. 4A).

Fig. 4

GO and KEGG enrichment analysis of the differentially expressed genes (DEGs) targeted by DE lncRNAs for C-XX versus C-XY (A, C), and C-XX versus RU-XX (B, D). GO and KEGG enrichment analysis of DEGs targeted by DE miRNAs for C-XX versus C-XY (E, G), and C-XX versus RU-XX (F, H)

To validate the target relationships, the fru-miR-122 and hsd17b1 were selected for dual luciferase reporter assay. The result showed that fru-miR-122 mimics significantly reduced luciferase activity in HEK-293T cells co-transfected with hsd17b1 WT (P < 0.05), while no significant difference was found in hsd17b1 MUT cells co-transfected with fru-miR-122 mimics or mimics-NC group (P < 0.05), suggesting that hsd17b1 is a downstream target of fru-miR-122 (Fig. 4B).

GO and KEGG enrichment analyses of target DEGs of DE lncRNAs and DE miRNAs

In C-XX vs. C-XY comparison group, results of GO enrichment analysis showed that the trans-target DEGs of DE lncRNAs were primarily enriched in hydrolase activity (GO: 0016787), extracellular region (GO: 0044421), extracellular region (GO:0005576), and hormone metabolic process (GO:0042445) (Fig. 5A). In C-XX vs. RU-XX comparison group, the DEGs were primarily enriched in multicellular organismal process (GO: 003250), single-multicellular organism process (GO: 0044707), developmental process (GO: 0032502), and receptor activity (GO: 0004872) (Fig. 5B). In C-XY vs. RU-XX comparison group, the DEGs were primarily enriched in multicellular organismal process (GO: 003250), response to external stimulus (GO: 0009605), extracellular region (GO:0005576), and reproductive process (GO: 0022414) (Figure S2A). In the C-XX vs. RU-XY comparison group, the DEGs were primarily enriched in response to chemical (GO: 0042221), multicellular organismal process (GO: 003250), receptor activity (GO:0004872), and molecular transducer activity (GO: 0060089) (Figure S2B). In C-XY vs. RU-XY, the DEGs were primarily enriched in response to endogenous stimulus (GO:0009719), response to external stimulus (GO:0009605), and extracellular region (GO: 0005576) (Figure S2C). In RU-XX vs. RU-XY, the DEGs were primarily enriched in muscle contraction (GO: 0006936), muscle system process (GO: 0003012), contractile fiber (GO:0043292), and modification of morphology or physiology of other organisms (GO: 0035821) (Figure S2D).

Fig. 5

Schematic diagram of interactions between DE miRNAs and DE mRNAs (A), and the relative luciferase activity of HEK-293T cells co-transfected with mimic fru-miR-122/NC (B). Wild-type (WT); Mutants (MUT); Mock-vehicle (MV); positive control (PC); * significance level at 0.05

In C-XX vs. C-XY comparison group, results of KEGG enrichment analysis showed that the trans-target DEGs of DE lncRNAs were primarily enriched in steroid biosynthesis, steroid hormone biosynthesis, neuroactive ligand-receptor interaction, and ovarian steroidogenesis pathways (Fig. 5C). In C-XX vs. RU-XX comparison group, the target DEGs were primarily enriched in calcium signaling pathways, cortisol synthesis and secretion, neuroactive ligand-receptor interaction, steroid hormone biosynthesis, and ovarian steroidogenesis pathways (Fig. 5D). Importantly, 12 DEGs in this comparison group, such as cyp17a1, igf1, hsd3b, hsd17b1, and cyp19a1a were screened in ovarian steroidogenesis pathway. In the C-XY vs. RU-XX comparison group, DEGs were primarily enriched in calcium signaling, neuroactive steroid biosynthesis, steroid biosynthesis, cortisol synthesis and secretion, and ligand-receptor interaction pathways (Figure S2E). In the C-XX vs. RU-XY comparison group, the DEGs were primarily enriched in circadian rhythm, steroid hormone biosynthesis, neuroactive ligand-receptor interaction, ovarian steroidogenesis, and cortisol synthesis and secretion pathways (Figure S2F). In the C-XY vs. RU-XY comparison group, the DEGs were primarily enriched in circadian rhythm, neuroactive ligand-receptor interaction, ovarian steroidogenesis, and steroid biosynthesis pathways (Figure S2G). In the RU-XX vs. RU-XY comparison group, the DEGs were primarily enriched in systemic lupus erythematosus, cell adhesion molecules, and estrogen signaling pathways (Figure S2H).

GO and KEGG enrichment of target DEGs of miRNAs was analyzed to recognize the molecular functions of DE miRNAs. In the comparison group of C-XX vs. C-XY, the DEGs mainly enriched in molecular transducer activity (GO:0060089), transmembrane receptor activity (GO:0099600), signal transducer activity (GO:0004871), and intrinsic component of membrane (GO:0031224) (Fig. 5E). In C-XX vs. RU-XX and C-XY vs. RU-XX comparison groups, the DEGs primarily enriched in the plasma membrane (GO:0005886), developmental process (GO:0032502), system development (GO:0048731), and multicellular organismal process (GO:0032501) (Fig. 5F, Figure S3A). In C-XX vs. RU-XY and C-XY vs. RU-XY comparison groups, the DEGs were primarily enriched in receptor activity (GO:0004872), molecular transducer activity (GO:0060089), circadian rhythm (GO:0007623), and signaling receptor activity (GO:0038023) (Figure S3B-S3C). In the RU-XX vs. RU-XY comparison group, the DEGs were primarily enriched in actomyosin structure organization (GO:0031032), defense response to bacterium (GO:0042742), signal transducer activity (GO:0004871), and receptor activity (GO:0004872) (Figure S3D).

The target DEGs of DE miRNAs in the C-XX vs. C-XY comparison group were primarily enriched in steroid hormone biosynthesis, neuroactive ligand-receptor interaction, ovarian steroidogenesis, and steroid biosynthesis pathways (Fig. 5G). In C-XX vs. RU-XX comparison group, the DEGs were primarily enriched in ovarian steroidogenesis, calcium signaling, cortisol synthesis and secretion, neuroactive ligand-receptor interaction, steroid hormone biosynthesis, and TGF-beta signaling pathways (Fig. 5H). Importantly, 11 DEGs in this comparison group, such as igf1, hsd3b, hsd17b1, cyp2j6, cyp19a1, and cyp17a1, were involved in ovarian steroidogenesis pathway. 12 DEGs, such as hsd3b, plcb1, and cyp11c1, were involved in cortisol synthesis and secretion pathway (Table S12). And, 16 DEGs, such as hamp, nog, and fsta, were involved in the TGF-beta signaling pathway. In the C-XY vs. RU-XX comparison group, the DEGs were primarily enriched in calcium signaling, circadian entrainment, neuroactive ligand-receptor interaction, TGF-beta signaling, and circadian rhythm pathways (Figure S3E). Both in the comparison groups of C-XY vs. RU-XX and C-XX vs. RU-XY, the DEGs were primarily enriched in circadian rhythm, neuroactive ligand-receptor interaction, ovarian steroidogenesis, cortisol synthesis and secretion, and circadian entrainment pathways (Figure S3F-S3G). In the RU-XX vs. RU-XY comparison group, the DEGs were primarily enriched in TGF-beta signaling, estrogen signaling, and cytokine and cytokine receptor pathways (Figure S3H).

qPCR

Among the four lncRNAs, the expression of XR_003890876 and XR_003889089 was significantly higher in the C-XX gonads than in the C-XY gonads, whereas the expression of XR_003886586 and XR_172564 was significantly higher in the C-XY gonads than that in the C-XX gonads (P<0.05). And, the expression of XR_003886586 was significantly increased in RU-XX gonads compared to that in C-XX gonads (P<0.05). Among the four circRNAs, circ_000134, circ_000244, and circ_000097_3 were significantly higher expressed in C-XX gonads, while circ_000245 was significantly higher expressed in C-XY gonads (P<0.05). The expression of circ_000245 was significantly increased in RU-XX gonads than that in C-XX gonads (P<0.05) (Fig. 6). Among the six miRNAs, novel-m0200-5p was significantly higher expressed in C-XX gonads, while the other miRNAs, including miR_749_z, fru_miR_122, m0001_5p_2, m0134_3p, miR_2187_y were significantly higher expressed in male gonads (P<0.05). The expression of novel-m0200-5p was down-regulated, while m0001_5p_2 was up-regulated in the gonads of RU-XX than that in C-XX gonads (P<0.05) (Fig. 6).

Fig. 6

The expression profiles of four lncRNAs, four circRNAs and six miRNAs in the gonads of C-XX, C-XY, and RU-XX groups. C-XX: control genetic female; C-XY: control genetic male; RU-XX: RU486-treated female. Each value represents the mean ± SEM of three measurements, and P < 0.05 shows significant difference

ceRNA regulatory networks

A sex differentiation ceRNA regulatory network including eight mRNAs, 23 miRNAs, and 48 lncRNAs was constructed (Fig. 7A). The results showed that one lncRNA could target several miRNAs, and one mRNA could be targeted by several miRNAs. It was found that gsdf was targeted by six miRNAs (novel-m0340-5p, novel-m0341-5p, novel-m0342-5p, miR-205-x, novel-m0326-3p, and novel-m0327-3p), while novel-m0340-5p, novel-m0341-5p, novel-m0342-5p, and miR-205-x were targeted by three lncRNAs (XR_003886319.1, MSTRG.6561.1, and MSTRG.9559.1), miR-205-x, novel-m0326-3p, and novel-m0327-3p were targeted by six lncRNAs (XR_003885989.1, XR_003886814.1, XR_003887458.1, XR_003888383.1, XR_003889060.1, and XR_003890136.1). One miRNA, miR-730-x, was predicted to regulate dmrt1 expression, was targeted by XR_003890295.1, XR_003887080.1, XR_003887012.1, and MSTRG.5738.1. Three lncRNAs, including MSTRG.13054.4, XR_003888827.1, and XR_003888884.1, were predicted to target novel-m0182-5p, and novel-m0183-5p, which regulate the gene cyp19a1a expression. Five lncRNAs, including MSTRG.11750.1, MSTRG.12471.3, XR_003891719.1, XR_003887833.1, and XR_003887090.1, were predicted to target miR-19-y, which regulate the expression of foxl2. The expression of hsd17b1 may regulated by two miRNAs (miR-122-x and fru-miR-122), while these two miRNAs were targeted by three lncRNAs (XR 00388331.1, XR 003891384.1, and XR 003889260.1). In the ceRNA regulatory networks, the expression patterns of most ncRNAs were in the opposite direction in RU-XX groups compared with C-XX groups (Fig. 7B).

Fig. 7

The ceRNAs and four genes with the characteristic expression pattern related to sex determination and differentiation. (A) The ceRNA regulatory networks of four sex-related genes and relative ncRNAs in the gonads. The red, purple, and green nodes represent mRNAs, miRNAs, and lncRNAs, respectively. (B) Heat map showing the mRNAs, miRNAs, and lncRNAs expression patterns in the ceRNA networks in C-XX, C-XY, RU-XX, and RU-XY

Discussion

ncRNAs have been demonstrated to be crucially important in sex differentiation, gonadal development, and gametogenesis in animals [40]. In teleost, sexual dimorphism expression patterns of ncRNAs were also been found in some species, such as swamp eel [21], golden pompano (Trachinotus blochii) [41, 42] Nile tilapia [43] and tiger puffer [44]. However, limited reports exist on the expression of ncRNAs in the gonads during the sex reversal process in teleosts. In this study, the potential functional roles of ncRNAs were further investigated during the RU486-induced masculinization process in T. rubripes gonads. Here, the qPCR results were in accordance with the RNA-sequence data, indicating the accuracy and dependability of our findings. This study gives insight into the epigenetic regulation of fish masculinization induced by RU486.

Although the sex-determining genes are diverse in teleosts, the downstream genes involved in gonadal sex differentiation appear to be conserved, such as dmrt1and gsdf were confirmed to be involved in testicular differentiation, while foxl2 and cyp19a1 were involved in ovarian differentiation [1]. The homozygous mutation of dmrt1 XY gonads developed into ovaries in medaka [45], zebrafish [46], and Nile tilapia [47]. The mutation of gsdf and amhr2 resulted in male-to-female sex reversal in medaka [48, 49] and Nile tilapia [50, 51]. The foxl2 or cyp19a1a homozygous mutation can lead to masculinization in some teleosts [52,53,54]. In tiger puffer, the expression of dmrt1 and gsdf was much higher in male gonads, while the expression of foxl2 and cyp19a1a were much higher in female gonads at early periods of tiger puffer sex differentiation [27]. Furthermore, the expression of dmrt1 and gsdf were increased, and foxl2 and cyp19a1a were decreased during aromatase inhibitor (AI) induced tiger puffer masculinization [55]. These findings indicate the importance of dmrt1, gsdf, foxl2, and cyp19a1a in sex differentiation and sex reversal of tiger puffer. Recently, increasing studies on the epigenetic regulation mechanism have indicated that lots of ncRNAs have potential roles in regulating sex-related genes in teleosts [40]. lncRNAs can regulate sex determination, sex differentiation, and sex reversal through positive and negative regulation of sex-related gene expression by cis or trans-acting regulatory mechanisms [16, 56]. In this study, many DE lncRNAs were predicted to trans target sex-related DEGs, such as XR_003885863.1, XR_003888756.1, and XR_003888884.1.etc were found may target female-related genes (cyp19a1, foxl2, and hsd17b1. etc.), MSTRG.13468.1, MSTRG.17819.1, and XR_003886055.1. etc. were found may target dmrt1, suggesting that lncRNAs may be involved in regulating the sex differentiation of tiger puffers. Similar results were also found in other vertebrates. In humans, it has been demonstrated that TCONS_00025195 and TCONS_00025196 might target sox9 in testes [15]. In Chinese soft-shell turtle (Pelodiscus sinensis), dmrt1, gata4, and cyp19a1a were key sex differentiation-related genes, and were cis or trans-targeted by several lncRNAs [57]. LncRNA MSTRG.24,346 positively regulated the expression of sex-determining candidate gene dmrt1 in large yellow croaker (Larimichthys crocea) [58]. Sex-related genes, including cyp19a1a, and hsd17β1 were regulated by the lncRNA TCONS_00083175 in golden pompano, indicating the importance of lncRNAs in golden pompano sex differentiation [41]. In swamp eel, lncRNA MSTRG.12,998 and MSTRG.38,036 were potentially involved in sex reversal, were highly expressed in the ovotestis than that in the ovary or testis, and were significantly upregulated during 17α‑methyltestosterone induced masculinization [18].

miRNAs play a crucial role in regulating the expression of mRNA by binding to complementary sequences in the 3’-untranslated regions (3’-UTRs) of target mRNAs. This binding can either inhibit or destabilize the expression of target genes [59,60,61]. More and more evidence has proved that miRNAs play an important role in gonadal differentiation and development in teleosts, such as in Atlantic halibut [62], Nile tilapia [63], and yellow catfish (Pelteobagrus fulvidraco) [64]. In this study, some DE miRNAs, such as the miR-202 family (include fru-miR-202, miR-202-x, miR-202-y, and miR-202-z), miR-205-z, miR-122-x, and fru-miR-122 were highly expressed in C-XY and RU-XX gonads than that in C-XX gonads, while some DE miRNAs, such as novel-m0326-3p, and novel-m0327-3p were highly expressed in C-XX gonads than that in C-XY and RU-XX gonads. Furthermore, some DE miRNAs were predicted to target sex-related genes, (including miR-205-z/foxl2 and ctnna2, miR-122, and fru-miR-122/hsd17b1). In vertebrates, the sequence of miR-202 is conserved [65], was identified expressed in adult testes of Atlantic halibut (Hippoglossus hippoglossus) [62], African clawed frogs (Xenopus) [66], mouse [67], and human [68], and fetal testes of chicken [69]. In the mouse, the primary transcript of miR-202 was strongly expressed in the testis during the gonadal differentiation [65]. In chicken gonads, miR-202 was expressed in a sexually dimorphic manner during gonadal sex differentiation, with a high expression in the testis [69]. The expression level of miR-202 decreased during E₂-induced feminization, whereas it increased during aromatase inhibitor-induced masculinization in embryonic chicken gonads [24]. This suggests that increased miR-202 expression is associated with testicular differentiation in the gonads of embryonic chicken [24]. In zebrafish, miR-202-5p was found to be a constituent of germplasm and a potential primordial germ cells (PGCs) marker, and it exhibited specific expression in the gonad of embryo zebrafish, with a greater level of expression observed in the testis compared to the ovary [70], and the same result was found in the Japanese flounder (Paralichthys olivaceus) [71]. In medaka, miR-202-5p was germ cell specific and identified as a key candidate for male differentiation and development [72]. All of these suggest the importance of the miR-202 family on animal sex differentiation. Previous studies identified that miR-205 was related to gonads development in some teleosts [63, 73, 74]. In common carp (Cyprinus carpio), miR-205 was predicted to regulate the female gonad development-related gene inhibin beta A chain and TGF-β signaling-related gene pdk1. And, the miR-205 was down-regulated in the primordial and juvenile gonads during the process of atrazine exposure-induced femininization in this species [73]. In tilapia gonads, miR-205b-3p was predicted to target cyp19a1a [63]. In Chinese longsnout catfish (Leiocassis longirostris), miR-205-3, which was predicted to have the potential to regulate genes in the ovary, was observed to be highly expressed in the testes [74]. In this study, miR-205-z was up-regulated in the C-XY and RU-XX gonads than that in the C-XX, and was predicted to target foxl2, and a female gonad development-related gene, ctnna2, indicating the potential roles of miR-205-z in tiger puffer masculinization. The miR-122 was identified as related to sex steroid hormone biosynthesis in humans and Nile tilapia [75, 76]. Additionally, it was discovered that miR-122 has the potential to regulate the mRNA binding protein of the luteinizing hormone receptor in rat ovaries [77]. Furthermore, a recent study revealed that miR-122 can regulate the expression of a female sex-related gene, femla, in Chinese longsnout catfish [74]. In common carp, miR-122, which was predicted to target esr1 and esr2, was down-regulated during the process of atrazine-induced femininization [73]. In this study, miR-122 and fru-miR-122 were predicted to regulate the E₂ synthesizing enzyme gene, hsd17b1. And these two miRNAs were highly expressed in C-XY and RU-XX gonads than that in C-XX gonads, which is consistent with the found in chicken embryonic and common carp [73, 78]. Additionally, our previous study found that the expression of hsd17b1 was down-regulated in RU-XX gonads than in C-XX gonads [28]. Therefore, we speculated that high levels of miR-122 and fru-miR-122 may suppress hsd17b1 expression during RU486-induced masculinization. Furthermore, our present study also observed lots of potential interactions of novel miRNAs and sex-related genes, such as novel-m0326-3p, novel-m0327-3p, novel-m0340-5p, novel-m0341-5p, and novel-m0342-5p were predicted to target gsdf, novel-m0182-5p and novel-m0183-5p were predicted to target cyp19a1a. These results may provide a good starting point for investigating the regulatory mechanism of miRNA on teleost sex differentiation.

Several researchers revealed that the complex mechanism of fish sex differentiation requires a coordinated network of ncRNAs and mRNAs, such as in Chinese tongue sole [22], Amur sturgeon [19], swamp eel [18, 21] A recent study identified that lncRNAs and circRNAs can act as sponges for miRNAs, thereby protecting target mRNAs from repression in fish sex determination and differentiation [22]. In the present study, we found that the expression of most lncRNAs that targeted gsdf, dmrt1, cyp17a1, esr1, or esr2 was up-regulated in C-XY and RU-XX groups compared to that in C-XX groups, which was consistent with the expression pattern of their target genes but contrary to the miRNAs that targeted these genes. In contrast, the expression of most lncRNAs targeting foxl2, cyp19a1a, or hsd17b1 was down-regulated in C-XY and RU-XX groups compared to that in C-XX groups, which was consistent with the expression pattern of their target genes but contrary to the miRNAs targeting foxl2, cyp19a1a, and hsd17b1. All these suggest that these lncRNAs may be involved in the competition between miRNAs and sex-related genes in tiger puffer sex differentiation and sex reversal.

In this study, the KEGG analysis showed that the DEGs targeted by DE lncRNAs or DE miRNAs were significantly enriched in sex differentiation-related pathways in RU-XX vs. C-XX, such as calcium signaling, ovarian steroidogenesis, and cortisol synthesis and secretion pathways. The calcium signaling pathway is associated with a variety of biological processes [79]. The environment influenced gender regulation by changing calcium ions in the red-eared slider turtle, and high temperatures caused an increase in Ca²⁺ influx in the somatic cells of gonads, thereby triggering the stat3 signaling factor and inhibiting the expression of kdm6b, leading to the activation of the female pathway [80]. Moreover, during early sexual development, exposure to letrozole resulted in female-to-male sex reversal of XX Nile tilapia, and at least ten transcripts associated with the calcium signaling pathway in the gonads were up-regulated [81]. Furthermore, our previous study revealed that several DEGs were involved in the calcium signaling pathway during E₂-induced feminization and MT- and AI-induced masculinization in tiger puffer [55]. All these indicate a key role of the calcium signaling pathway in sex differentiation and sex reversal in vertebrates. Intracellular Ca²⁺ is widely recognized as a second messenger in the rapid effects of most known steroids [82]. In fathead minnow (Pimephales promelas), cyclotrimethylenetrinitramine (RDX) exposure altered the expression level of at least three transcripts related to calcium transport, binding, and signaling in brain tissue [83]. In this study, the expression of at least 20 DEGs (including atp2a1, chrm3, chrm5, and chrna7) which targeted by DE lncRNAs and/or DE miRNAs, associated with calcium signaling pathway were upregulated in C-XX vs. RU-XX groups. It is speculated that miRNAs and lncRNAs may modulate the expression of DEGs in the calcium signaling pathway under RU486 treatment, mediating intracellular calcium homeostasis, and modulating multiple pathways. This may result in the activation of the male pathway in the genetic female tiger puffer.

Sex steroids are relatively conserved in sex determination and differentiation in teleosts [84]. Changes in sex steroid hormone concentration during appropriate developmental stages can cause sex reversal [85]. Sex steroidogenesis commences with the transportation of cholesterol into the mitochondria, facilitated by steroidogenic acute regulatory protein (StAR) [86]. Within the mitochondria, cholesterol is transformed into pregnenolone, the initial precursor in the steroidogenic cascade [87]. The downstream enzymes of the synthesis pathway, encoded by cyp19a1a, cyp17a1, cyp11a1, cyp11c1, cyp17a2, and so on, conserved the pregnenolone to E₂ or 11-KT in the gonads [84]. The ovary produces steroid hormones that are essential for maintaining the female phenotype and normal ovarian processes, such as sex differentiation, follicle growth, oocyte maturation, and ovulation [84, 88]. In the process of ovarian steroidogenesis, the enzymes encoded by cyp19a1a, cyp17a1, and star2, are involved in the production of estrogen and exclusively expressed in the gonads. These genes mutation leads to the activation of male pathways and female-to-male sex reversal in fish [84]. Notably, a previous study found that ovarian steroidogenesis was regulated by miRNAs in mice. Specifically, miR-132 attenuates steroidogenesis by repressing star expression and inducing 20α-hydroxysteroid dehydrogenase (20α-HSD) via inhibition of methyl-CpG binding protein 2 (MeCP2) to generate a biologically inactive 20α-hydroxyprogesterone (20α-OHP) [89]. In this study, the expression of at least ten DEGs (including cyp19a1, cyp17a1, hsd17b1, hsd3b, igf1, cyp2j6, and so on), which are targeted by DE lncRNAs and/or DE miRNAs, involved in ovarian steroidogenesis were altered in the XX gonads treated by RU486. Among the DEGs, cyp19a1 and hsd17b1 were significantly decreased in the RU-XX group. Therefore, we hypothesized that the alteration of cyp19a1a and hsd17b1 mRNA expressions may be associated with the expression of miRNAs and lncRNAs that targeted these two genes during the process of RU486-induced tiger puffer masculinization (Fig. 8).

Fig. 8

The potential regulatory role of non-coding RNAs in mifepristone-induced masculinization in T. rubripes gonads

Cortisol, the dominant glucocorticoid in fish, is directly associated with environmental stress. Recent studies have confirmed that environmental stressors such as high density, high temperature, and bright background color, as well as cortisol administration, can induce the masculinization of genetically female fish, which may be associated with increased cortisol levels [90]. In the C-XX vs. RU-XX comparison group of this study, the DEGs targeted by DE lncRNAs and/or DE miRNAs were significantly enriched in cortisol synthesis and secretion pathway. Among the DEGs, at least ten DEGs (cyp11c1, cyp17a1, hsd3b, adcy8, plcb1, and so on) were up-regulated, suggesting that RU486 may act as an environmental endocrine disruptor to promote the synthesis and secretion of cortisol. Furthermore, previous studies have identified that RU486 is a diverse endocrine, and acts as both a progesterone and cortisol antagonist [91]. RU486 treatment increased the whole-body cortisol concentration in summer flounder, and increased plasma cortisol levels in Rainbow trout [92]. Therefore, it is hypothesized that the altered expression patterns of genes involved in cortisol synthesis and secretion may be related to the expression of miRNAs and lncRNAs, and the high expression of these genes may activate a series of downstream female-to-male pathways in RU-XX gonads in the present study.

Conclusions

Our results provide the evidence that ncRNAs may participate in RU486-induced masculinization in T. rubripes, and may enhance our understanding of the regulatory network of sex differentiation in fugu. Also, further research should be conducted to elucidate the functional significance of those identified ncRNAs.