- Open Access
Transcriptome analysis revealed gene expression feminization of testis after exogenous tetrodotoxin administration in pufferfish Takifugu flavidus
BMC Genomics volume 23, Article number: 553 (2022)
Tetrodotoxin (TTX) is a deadly neurotoxin and usually accumulates in large amounts in the ovaries but is non-toxic or low toxic in the testis of pufferfish. The molecular mechanism underlying sexual dimorphism accumulation of TTX in ovary and testis, and the relationship between TTX accumulation with sex related genes expression remain largely unknown. The present study investigated the effects of exogenous TTX treatment on Takifugu flavidus. The results demonstrated that exogenous TTX administration significantly incresed level of TTX concentration in kidney, cholecyst, skin, liver, heart, muscle, ovary and testis of the treatment group (TG) than that of the control group (CG). Transcriptome sequencing and analysis were performed to study differential expression profiles of mRNA and piRNA after TTX administration of the ovary and testis. The results showed that compared with female control group (FCG) and male control group (MCG), TTX administration resulted in 80 and 23 piRNAs, 126 and 223 genes up and down regulated expression in female TTX-treated group (FTG), meanwhile, 286 and 223 piRNAs, 2 and 443 genes up and down regulated expression in male TTX-treated group (MTG). The female dominant genes cyp19a1, gdf9 and foxl2 were found to be up-regulated in MTG. The cyp19a1, whose corresponding target piRNA uniq_554482 was identified as down-regulated in the MTG, indicating the gene expression feminization in testis after exogenous TTX administration. The KEGG enrichment analysis revealed that differentially expressed genes (DEGs) and piRNAs (DEpiRNAs) in MTG vs MCG group were more enriched in metabolism pathways, indicating that the testis produced more metabolic pathways in response to exogenous TTX, which might be a reason for the sexual dimorphism of TTX distribution in gonads. In addition, TdT-mediated dUTP-biotin nick end labeling staining showed that significant apoptosis was detected in the MTG testis, and the role of the cell apoptotic pathways was further confirmed. Overall, our research revealed that the response of the ovary and testis to TTX administration was largely different, the ovary is more tolerant whereas the testis is more sensitive to TTX. These data will deepen our understanding on the accumulation of TTX sexual dimorphism in Takifugu.
Tetrodotoxin (TTX), one of the most potent neurotoxins which is specific to voltage-gated sodium channels (VGSCs) on excitable membranes of muscle and nerve tissues, was long believed to occur exclusively in TTX-bearing pufferfish, mainly the genus Takifugu [1, 2]. TTX is produced primarily by marine bacteria and accumulated extremely high in pufferfish through the food chain with tissue-specific and seasonal changes [3,4,5]. It is suggested that TTX may function as a chemical defense against predators and as a pheromone during spawning [6, 7].
The accumulation of TTX in pufferfish has been studied extensively [1, 2, 8,9,10,11]. The distribution of TTX in pufferfish bodies appears to be species-specific [2, 3, 5, 7, 12,13,14]. In marine species, the liver and ovary generally have the highest toxicity, followed by the intestines and skin. Muscles and/or testis are non-toxic or weakly toxic . From a seasonal perspective, whole body TTX content was significantly higher during the maturation/spawning period than in the ordinary period . Accumulation of TTX also demonstrated sex dimorphism, ovary generally accumulated a higher level of TTX than the testis [2, 3, 5, 7, 12, 13]. Interestingly, through all seasons, TTX was detected in the skin or ovary of females and the skin or liver of males in Takifugu. poecilonotus and Takifugu niphobles [4, 11, 15], the difference in TTX accumulation between females and males was particularly evident during the spawning period. These tissue-specific seasonal changes in TTX content suggest that TTX accumulation might be affected by sexual maturation .
The larvae of Takifugu pufferfish were protected by maternal TTX which had accumulated in the eggs during their development in the ovary [16, 17]. In addition, studies have shown that TTX as a toxin or pheromone accumulated in the ovary protects pufferfish larvae from predators [6, 16, 17]. Although TTX accumulation demonstrates sexual dimorphism in gonads of pufferfish, the molecular mechanism underlying and the relationship between TTX accumulation and steroidogenesis, sex differentiation and reproduction-related genes transcriptome analysis remained unclear.
Transcriptome analysis technologies are important systems-biology methods for the investigation of the organisms' stress response, gene expression, and regulatory pathways changes. Transcriptome refers to the set of all RNA molecules from protein-coding (mRNA) to noncoding RNA, including rRNA, tRNA, lncRNA, pri-miRNA, and others [18, 19] Piwi-interacting RNAs (Piwi/piRNAs) are single-stranded, 26–32 nt length small RNAs that function by interacting with Piwi proteins to form RNA–protein complexes [20, 21]. Animal studies revealed that piRNAs are limited expressed in a few tissues, such as tumors [22,23,24], brain and nervous tissue [23, 25,26,27], whereas generally abundant in gonads of invertebrates [28,29,30,31,32] and vertebrates [23, 33,34,35,36,37,38]. These results were also consistent with teleosts, including zebrafish (Danio rerio) [30, 39], tilapia (Oreochromis niloticus) , Japanese flounder (Paralichthys olivaceus) , tongue soles (Cynoglossus semilaevis) . Previous studies have shown that fertility and normal gonadal morphology in the female progeny of Drosophila requires maternal piRNAs, but not males . piRNAs derived from the W chromosome are expressed more abundantly in the ovary than in the testis of silkworms .
The piRNAs play a role in targeting the transcripts of active TEs (Transcriptional elements) and maintaining the methylation pattern of maternal genomic DNA. It was proved that ovarian and testicular piRNAs are two different classes with different functions, and their expression appears to be regulated by distinct transcription factors in humans . Above all, at least a large number of piRNAs showed sexual dimorphism in the testes and ovaries. In addition, putative piRNAs were found abundant in the gonads of Takifugu rubripes , whereas the function and connection between TTX and piRNAs in the gonads of pufferfish have not been investigated. To gain insights into sexual dimorphism of TTX accumulation in gonads, steroidogenesis, sex differentiation and reproduction-related genes and small non-coding RNAs were chosen to further investigate the molecular mechanisms.
In the present study, we investigated the expression changes of piRNAs and mRNAs in both ovary and testis of Takifugu flavidus on the transcriptomic level after the intramuscular injection with TTX. The differentially expressed piRNAs and genes were identified, function and pathway of these piRNAs and genes were also analyzed. At last, apoptosis detection was performed in the TTX-treated gonads. Through transcriptome analysis, the gene expression pattern's feminization of testis in pufferfish Takifugu flavidus after exogenous TTX administration was observed, and these data will deepen our understanding on the accumulation of TTX sexual dimorphism in Takifugu.
Exogenous TTX administration significantly elevated the tissue TTX level of Takifugu flavidus
After 8 h of TTX administration, the LC–MS/MS (Liquid Chromatography with tandem mass spectrometry) method was performed to detect the concentration of TTX in different groups and tissues, the results showed that the TTX concentration of TG (the TTX-treated group) is significantly higher than CG (the control group) in kidney, cholecyst, skin, liver, heart, and muscle with mixed-sex (Fig. 1A).
To further explore the accumulation of TTX in the ovary and testis of Takifugu flavidus, histological analysis was performed to identify the sex of gonads in 6 months old of Takifugu flavidus. Twelve gonads were dissected from each TG and CG. Five ovaries and seven testes were confirmed in the CG, and six ovaries and six testes were confirmed in the TG. All of these gonads are immature, the ovarian cavity, oogonia and oocytes were observed in the ovary (Fig. 1B-a, b). The spermatogonia without spermatocytes or spermatids were observed in the testis (Fig. 1B-c, d). Immunofluorescence histochemical analysis used an anti-TTX antibody demonstrated that TTX was localized at the somatic cells around the oocytes in the FCG (the female control group) and FTG (the female TTX-treated group) ovary with intense signals (Fig. 1B-e, f). Weak and strong TTX signals were detected in the MCG (the male control group) and MTG (the male TTX-treated group) testis, respectively (Fig. 1B-g, h). Algorithms for unified illumination, focus correction, and tiling and stitching enabled us to analyze large fields of cells and quantify the average cellular staining intensity for different TTX, the statistical analysis showed that the total TTX signals in the TG gonads were significantly higher than the CG gonads after TTX administration (Fig. 1C).
Transcriptome sequencing and novel piRNA identification
Total RNA of gonads (testis or ovary) was isolated and sequenced for piRNA and mRNA analysis from the CG and TG of Takifugu flavidus. Overall, 94.91% of FCG, 94.91% of MCG, 98.09% of FTG, and 94.32% of MTG of the clean reads had scored at the Q30 level in mRNA sequencing (Table 1). 85.06% (FCG), 86.10% (MCG), 86.04% (FTG) and 85.63% (MTG) clean reads were matched onto the Takifugu rubripes reference genome (fTakRub1.2).
piRNA sequencing results showed that the peak value of all groups was mainly concentrated at 28 bp (Figure S1A), with the character of 1U-bias of the molecule (Figure S1B), which was consistent with the characteristic length and structure bias of piRNAs. The numbers of unique/novel piRNAs for the four groups were calculated, there were 32,089, 20,003, 30,020 and 31,691 in the FCG, MCG, FTG and MTG, respectively. The results demonstrated that TTX administration resulted in a decreased number of piRNAs in females, while an increased number of piRNAs in males. The total unique piRNA is 65,535 (Table 1 and Table S1), including known 943 piRNAs, about 1.439% of the total unique piRNAs. 98 female-specific and 90 male-specific piRNAs (Table S1). There were 60,497 unique piRNAs that display less than 10 reads in all groups (Table S1). A total of 79 piRNA clusters were isolated based on the characteristics of the piRNA clusters (Table S2). These piRNA clusters exist and are unevenly distributed on each chromosome (Table S2). Two PIWI protein-coding genes, piwil1 and piwil2, were identified from transcriptome data in Takifugu flavidus. The transcriptome analysis revealed that piwil1 and piwil2 were both expressed in the ovary and testis, and the expression level in the testis is relatively higher than that in the ovary. Among them, the expression level of piwil1 was always higher than piwil2 in both ovaries and testes. After TTX administration, both piwil1 and piwil2 were down-regulated in the testis, whereas piwil1 was up-regulated in the ovary with no significant difference (Table S4). These results indicated that piwil1 and piwil2 were more sensitive to TTX in the testis than the ovary of Takifugu flavidus.
Transcriptome analysis of differentially expressed genes and piRNAs between the FCG and MCG
Overall distribution of differentially expressed mRNAs and piRNAs was observed on the volcano plot. By comparing MCG and FCG, a total of 219 differentially expressed genes were identified with P < 0.05 and |log2 (fold change) |> 1 using the base mean method, and 660 differentially expressed piRNAs were identified with P < 0.05 and |log2 (fold change) |> 1 using the TPM method. 41 genes and 335 piRNAs were differentially expressed in MCG, 178 genes and 315 piRNAs were differentially expressed in FCG (Fig. 2A-B). The female dominant genes such as cyp19a1 (cytochrome P450 19A1-like) and foxl2 (Forkhead Box L2) were identified in the FCG, and the male dominant genes such as LOC101075863/cyp11b2 (cytochrome P450 11B, mitochondrial) and dmrt1 (doublesex and mab-3 related transcription factor 1) were identified in the MCG (Table S4).
KEGG (Kyoto Encyclopaedia of Genes and Genomes) pathway enrichment analysis of FCG and MCG differentially expressed genes showed that steroid hormone biosynthesis, oocyte meiosis and GnRH signaling pathway were both enriched in MCG and FCG. However, the involved genes are different, such as cyp19a1 and hsd17b1 (hydroxysteroid 17-beta dehydrogenase 1) enriched in FCG’s, cyp11b2 and cyp11a1/LOC101063139/ (cholesterol side-chain cleavage enzyme, mitochondrial) enriched in MCG’s steroid hormone biosynthesis. Folate biosynthesis, p53 signaling and some DNA replication/repair were found to be specifically enriched in FCG. Apoptosis and autophagy were found to be specifically enriched in MCG (Fig. 2C-D).
Transcriptome analysis of up and down regulated genes and piRNAs between the FCG and FTG
By comparing FTG and FCG, TTX administration resulted in 23 piRNAs and 223 genes expression up-regulated, meanwhile, 80 piRNAs and 126 genes expression down-regulated in the ovaries (Fig. 3A-B and Table S2, S4).
To investigate the reliability of differentially expressed genes and piRNAs in transcriptomes. Eight relative highly expressed piRNAs (uniq_1163081, uniq_3850196, uniq_556841, uniq_1137493, uniq_1958827, uniq_1013268, uniq_479370 and uniq_1713020) and six genes (amhr2, anti-Mullerian hormone receptor type 2; cacna1c, calcium voltage-gated channel subunit alpha1 C; myc, MYC proto-oncogene bHLH transcription factor a; ppard, peroxisome proliferator activated receptor delta; LOC101078894/adh1, alcohol dehydrogenase 1-like and LOC101069468/lpr5, low-density lipoprotein receptor-related protein 5-like) were screened for further verification by qRT-PCR (Quantitative Real-Time PCR). The piRNAs and genes’ expression changes of trend between the FTG and FCG in the qRT-PCR results were almost consistent with the transcriptome data, these results indicating that the transcriptome data was reliable (Fig. 3C-D and Table S1, S4).
KEGG and GO (Gene Ontology, which was shown in Table S7-8) pathway enrichment analysis of FTG compared FCG down-regulated piRNA target genes showed that these genes mainly enriched in signaling pathways, such as MAPK, cAMP, Rap 1, calcium, thyroid hormone, insulin, and AMPK signaling pathway (Fig. 4A and Table S7-8). KEGG pathway enrichment analysis of FTG compared FCG up-regulated genes showed that these genes mainly enriched in metabolism pathways, such as nitrogen, histidine, retinol and tyrosine metabolism. Furthermore, amhr2, and myc were enriched in the TGFβ signaling pathway, LOC101064673/map2k6 (dual specificity mitogen-activated protein kinase kinase 6), LOC101067856/prkaca (cAMP-dependent protein kinase catalytic subunit alpha), cacna1c were enriched in GnRH signaling pathway (Fig. 4B and Table S5-6).
To further confirm the correlation between down-regulated piRNAs and up-regulated expression genes, KEGG pathway enrichment analysis of FTG compared FCG down-regulated piRNA and up-regulated genes were compared. Up-regulated genes enriched in the oocyte meiosis, progesterone-mediated oocyte maturation, and steroid hormone biosynthesis were selected for mapping with the down-regulated piRNAs in FTG. The results showed that uniq_1958827 predicated target gene cacna1c; uniq_1013268 predicated target gene LOC101078894/adh1; uniq_479370 predicated target gene LOC101069468/lpr5; uniq_1713020 predicated target gene ppard; uniq_671648 predicated target gene amhr2 were identified. cacna1c was involved in the GnRH signaling pathway, adh1 was involved in retinol metabolism; both lpr5 and ppard were involved in the Wnt signaling pathway; and amhr2 was involved in the TGFβ signaling pathway (Fig. 4C and Table S5, S7). These genes/piRNAs' relative expression values in FTG and FCG were shown using RPKM/TPM (Fig. 4D). In addition, the uniq_671648-amhr2 mRNA pairing structure was shown in Fig. 4E.
Transcriptome analysis revealed gene expression feminization of testis after TTX administration
By comparing MTG and MCG, TTX administration resulted in 224 piRNAs and 443 genes expression up-regulated, meanwhile, 286 piRNAs and 2 genes (ddah2, dimethylarginine dimethylaminohydrolase 2 and ing4, an inhibitor of growth family member 4) expression down-regulated in the testes, the female dominant genes such as cyp19a1, gdf9 (growth differentiation factor 9) and foxl2 were identified up-regulated in the MTG (Fig. 5A-B and Table S1, S4).
To investigate the reliability of differentially expressed genes and piRNAs in transcriptomes, nine relative highly expressed piRNAs (uniq_2692448, uniq_904526, uniq_1754588, uniq_1320158, uniq_1010678, uniq_1323470, uniq_2869791, uniq_1678838 and uniq_766456), the up-regulated genes including hsd17b1, ccnb1 (cyclin B1), aurka (aurora kinase A), fbxo43 (F-box protein 43), ccne2/LOC101065671 (cyclin N-terminal domain-containing protein 2-like), pygl (glycogen phosphorylase L), gadd45β/LOC101080011 (growth arrest and DNA damage-inducible protein GADD45 beta) and the down-regulated two genes ing4 and ddah2 in the MTG were screened for further verification by qRT-PCR. The qRT-PCR results confirmed a good correlation with transcriptome sequencing results in piRNA and gene expression level (Fig. 5C-D and Table S1, S4), these results also indicated that the transcriptome data was reliable.
KEGG pathway enrichment analysis of MTG compared MCG down-regulated piRNA target genes showed that these genes also mainly enriched in signaling pathways, such as calcium, oxytocin, adrenergic, Hippo, thyroid hormone, wnt, foxO, estrogen, neurotrophin, GnRH and erbB signaling pathway (Fig. 6A and Table S7-8). KEGG pathway enrichment analysis of MTG compared MCG up regulated genes showed that these genes mainly enriched in ovary differentiated related pathways, such as ovarian steroidogenesis and steroid hormone biosynthesis, oocyte meiosis and progesterone-mediated oocyte maturation; the cell apoptosis-related pathways, such as cell cycle, tight junction, PPAR and p53 signaling pathway (Fig. 6B). Furthermore, ccnb1, ccnb2/LOC101076205 (G2/mitotic-specific cyclin-B2-like), cdc20/LOC101063487 (cell division cycle 20), mad2l1, zfp36l2/LOC101079369 (ZFP36 ring finger protein-like 2), pttg1/LOC101072182 (PTTG1 regulator of sister chromatid separation, securin), aurka, fbxo43 were enriched in oocyte meiosis; cyp19a1, hsd17b1, scarb1 (scavenger receptor class B member 1), gdf9, pla2g4c/LOC101069259 (cytosolic phospholipase A2-gamma) were enriched in ovarian steroidogenesis; ccnb1, ccnb2, mad2l1 and aurka were enriched in progesterone-mediated oocyte maturation; cyp19a1, hsd17b1, hsd17b2/LOC101067779 (very-long-chain 3-oxoacyl-CoA reductase-A) were enriched in Steroid hormone biosynthesis (Fig. 6B and Table S5-6). These results indicated the gene expression feminization in testis after exogenous TTX administration.
To further confirm the correlation between down-regulated piRNA and up-regulated expression genes, KEGG pathway enrichment analysis of MTG compared MCG down-regulated piRNA and up-regulated genes were compared. Up-regulated genes enriched in oocyte meiosis, progesterone-mediated oocyte maturation, steroid hormone biosynthesis, p53 signaling pathway, ferroptosis, necroptosis and cellular senescence were selected for mapping with the down-regulated piRNAs in MTG. The results showed that uniq_1010678 predicated target gene ccne2/LOC101065671; uniq_1138005 predicated target gene ccnb2; uniq_1678838 predicated target gene aurka; uniq_2509649 and uniq_1115741 predicated target genes mad2l1; uniq_1504863 and uniq_2344671 predicated target gene ccnb1; uniq_988037 predicated target gene cdc20; uniq_1323470 predicated target gene fbxo43; uniq_766456 predicated target gene hsd17b1; uniq_554482 predicated target gene cyp19a1; uniq_1713020 predicated target gene gadd45β; uniq_2869791 and uniq_1253913 predicated target gene pygl and uniq_1938531 predicated target gene acyl-CoA synthetase long chain family member 5 (acsl5) were identified (Fig. 6C and Table S3). The genes ccnb1, mad2l1 (mitotic arrest deficient 2 like 1), zfp36l2/LOC101079369, pttg1/ LOC101072182, aurka, fbxo43, ccnb2/ LOC101076205 and cdc20/LOC101063487 were involved in the oocyte meiosis or progesterone-mediated oocyte maturation pathway, hsd17b1 and cyp19a1 were involved in the steroid hormone biosynthesis, the gene gadd45β/LOC10108001 was involved in the p53 signaling pathway, pygl was involved in necroptosis and acsl5 was involved in ferroptosis (Fig. 6C). These genes/piRNAs’ relative expression values in MTG and MCG were shown using RPKM/TPM (Fig. 6D). In addition, the uniq_554482-cyp19a1/ccnb1 RNA pairing structures were shown in Fig. 6E.
Transcriptome analysis of sexual dimorphism distribution mechanism of TTX in gonads
To gain insight into the sexual dimorphism distribution mechanism of TTX in the ovary and testis. Differentially expressed genes (DEGs) and piRNAs (DEpiRNAs) in MTG vs MCG and FTG vs FCG were compared analysis. Venn diagram analysis showed that 401 and 304 DEGs were found to be significantly differentially expressed in MTG vs MCG and FTG vs FCG groups, respectively. And 44 DEGs were common between these two groups (Fig. 7A). There were 212 and 11 DEpiRNAs significantly differentially expressed in MTG vs MCG and FTG vs FCG groups, respectively. And 298 DEpiRNAs were common between these two groups (Fig. 7B).
Common and differentially expressed DEGs/DEpiRNAs among six groups, including FCG and MCG highly expressed genes/piRNAs groups, FTG vs FCG and MTG vs MCG up-regulated genes/piRNA groups, and FTG vs FCG and MTG vs MCG down-regulated genes/piRNA groups, were compared analyzed. As the Venn diagram showed that 47, 2 and 303 differentially expressed genes were found in FCG highly expressed genes, MTG vs MCG down- and up-regulated genes groups, the FCG highly expressed genes group shared 140 common genes with MTG vs MCG up-regulated genes group (Fig. 7C). The same pattern is also seen in the piRNA comparison group (Fig. 7E). In contrast, 39, 126 and 220 differentially expressed genes were found in MCG highly expressed genes, FTG vs FCG down and up-regulated genes groups, the MCG highly expressed genes group only shared 2 common genes with FTG vs FCG up-regulated genes group (Fig. 7D). The same pattern was also seen in the piRNA comparison group (Fig. 7F). These results further indicated that exogenous TTX can cause feminization of testis.
In the KEGG enrichment analysis of differentially expressed DEGs in FTG vs FCG and MTG vs MCG groups (304 and 401 genes), the results showed that compared with FTG vs FCG group, differentially expressed DEGs in MTG vs MCG group were more enriched in metabolism pathway, including nitrogen, vitamin B6, taurine and hypotaurine, starch and sucrose, glyoxylate and dicarboxylate and glycerolipid metabolism pathways (Fig. 7G-H). KEGG pathway enrichment analysis of 44 shared common DEGs of MTG vs MCG and FTG vs FCG, the results showed that they were also mainly enriched in metabolism pathways, including nitrogen, glyoxylate and dicarboxylate, alanine, aspartate and glutamate, drug-other enzymes, pyrimidine and purine metabolism pathways (Fig. 7I). These results suggested that the effect of TTX on gonads was direct. Compared with the ovary, the testis produced more metabolic pathways in response to exogenous TTX, which might be a reason for the sexual dimorphism of TTX distribution in gonads.
Exogenous TTX induces cell apoptosis in testis
The piRNA and mRNA transcriptome analysis showed that cell apoptosis-related genes were up-regulated in testis and related piRNA were down-regulated. To further confirm the damaging effects of TTX on the testis, TUNEL assay (TdT-mediated dUTP-biotin nick end labeling) was used for identifying apoptotic signals or cells in situ. Immunohistochemical analysis showed that positive signals were detected in the somatic cells nuclei of the FTG, FCG ovary and MTG testis. Little to no signal of apoptosis was observed in the MCG testis (Fig. 8A). Statistical analysis of the area of the positive signals, the results showed no significant differences between FTG and FCG, whereas, the positive signals in MTG were significantly higher than that in the MCG (Fig. 8B). These results indicated that exogenous TTX administration induces cell apoptosis in testis.
Sex differentiated Takifugu flavids showed sexual dimorphism on response to exogenous TTX
Processes involving sex determination and reproduction were complicated and regulated by various internal and/or external factors, particularly in pufferfish which showed a sexual dimorphism of TTX accumulation in gonads. The molecular data available was limited for Takifugu to know about the exact internal mechanism or role of TTX in the development of gonads. The accumulation and tissue-specific distribution of TTX in pufferfish, mainly in the genus Takifugu have been widely investigated from the viewpoint of the TTX-resistance VGSCs expression and distribution [2, 10, 47]. Evidence also showed that TTX is usually significantly at a high level during the maturation/spawning period than in the ordinary period [1, 4, 5, 9, 17]. The present study also demonstrated that the TTX concentration level is quite low in 6-month-old Takifugu flavids, and after TTX administration, TG level of TTX was significantly higher than CG in kidney, cholecyst, skin, liver, heart, muscle, ovary and testis, which was reflected by LC–MS/MS and fluorescence immunohistochemistry analysis, these results consistent with the exogenous TTX administration in Takifugu rubripes [48,49,50].
The artificially cultured 6-month-old Takifugu flavids were chosen for the first time to investigate the piRNAs and mRNAs expression changes in separate sex of pufferfish. Firstly, from a seasonal perspective, whole body TTX level was low in this period, accumulation of TTX also demonstrated sex dimorphism, ovary accumulated a higher level of TTX than the testis. Meanwhile, “Endogenous” TTX was also found in the ovary at this stage, TTX immunofluorescence signals were detected indicating the TTX was also accumulated in the ovary before the exogenous TTX treatment. At this time, endogenous TTX levels in the gonads do not reach the high levels compared with sexual maturation period, which can avoid the potentially serious or to be ignored effects of exogenous TTX administration. Secondly, TTX in the developing stage of pufferfish, especially for the gonads marginally been investigated. Also, at this time, Takifugu flavids has been sex differentiated, and it was reflected by clear histological differences, i.e., the ovary has an ovarian cavity and oocytes have initiated meiosis and formed growing oocytes, the testis only available the spermatogonia without initiated meiosis. Thirdly, after exogenous TTX administration, the TTX amount of both the testis and ovary was significantly increased. Meanwhile, the transcriptome data and analysis also showed more piRNAs and genes in the testis were affected by TTX than that in ovaries. Therefore, the response to exogenous TTX in this period is sufficient and worthy of in-depth study.
piRNAs play a role in the regulation of sex related genes’ expression after TTX administration
Previous studies demonstrated that piRNAs are indispensable for generating functional germ cells, gametogenesis, and female infertility [22, 33, 51]. and often bound to members of the piwi protein family to realize its regulation function [22, 24, 27, 28, 52]. Two PIWI proteins, piwil1 and piwil2, were identified in Takifugu flavidus, the results were consistent with another pufferfish Takifugu fasciatus, which also reported having two Piwils that were abundantly expressed in the immature gonads . piwil1 and piwil2 were both expressed in the ovaries and testes, and the expression level in the testis is relatively higher than that in the ovary of Takifugu flavids. After TTX administration, both piwil1 and piwil2 were down-regulated in the testis, whereas piwil1 was up-regulated in the ovary with no significant difference. The present study predicted piRNA uniq_1504863, uniq_1302131 and uniq_2344671 might be bound to these two Piwi proteins, and these results indicated that piwil1 and piwil2 were more sensitive to TTX in the testis in Takifugu flavidus.
Interestingly, up-regulated genes in the MTG were enriched in steroid hormone biosynthesis and ovarian steroidogenesis pathways, and the related female dominant gene cyp19a1 was identified. A piRNA uniq_554482 was also screened in the down-regulated piRNAs in MTG, and the uniq_554482 predicted the target gene including cyp19a1. In addition, another ovarian factor gdf9 and steroidogenesis related gene hsd17b1 were also significantly up regulated after TTX administration. The down-regulated piRNAs uniq_2093102 (targeted gdf9) and uniq_766456 (targeted hsd17b1) were also identified and involved in the upregulated genes in the MTG, showing the impact of TTX on genes’ expression may be through the piRNAs in the testis, and indicated that piRNAs play a role in the regulation of sex related genes after TTX administration.
Ovaries are more tolerant to exogenous TTX than testis
The present study revealed that ovaries are more tolerant to exogenous TTX than testis, this is reflected by a smaller number of genes and piRNAs that were changed after TTX administration and not significantly elevated TUNEL positive signals in the ovary. In addition, the total number of piRNAs is greatly increased in the testis but slightly decreased in the ovary after exogenous TTX administration. These results indicated that TTX has a bigger influence on the testis. There are two possible reasons for this phenomenon. First, the ovary itself has the accumulation of TTX and was more tolerant to exogenous TTX, so there is little change in piRNA quantity. The second is exogenous TTX administration can lead to feminization of testicular gene expressions, such as up-regulated expression of female-dominant genes. We speculate that the increased piRNAs in the testis silence more transposons or genes to resist/reduce the toxicity of TTX.
As previous studies reported that the p53 signaling pathway is involved in the regulation of cell cycle arrest, senescence and apoptosis [54,55,56], necroptosis and ferroptosis pathways were also involved in the regulation of cell death [56,57,58]. Transcriptome analysis in this study indicated that the apoptosis, p53 signaling and necroptosis pathways were found in MTG vs MCG up-regulated genes group in testis, whereas only the necroptosis pathway was found in FTG vs FCG up-regulated genes group in the ovary. These results were consistent with the TUNEL analysis of the present study that exogenous TTX induced more cell apoptosis in the testis than the ovary.
Furthermore, it’s hard to observe the changes in somatic and germ cells from the view of histological or HE results after 8 h of TTX administration. In general, the changes of nucleic acids are always earlier than protein and histological changes, for example, tilapia had formed primitive gonads 5 days after hatching, and there was no difference at the histological level between XX and XY individuals at this period . However, at the molecular level, XX individuals had begun to express Cyp19a1a, while XY individuals did not. In contrast, XY individuals had begun to express Amh (anti-Müllerian hormone) and Dmrt1, while XX individuals did not .
Testis demonstrated feminized gene expression patterns and up-regulated metabolic pathways after TTX administration
The molecular mechanism of sexual dimorphism in TTX accumulation in gonads of pufferfish is largely unknown, the relationship between TTX accumulation and sexual steroidogenesis, sex differentiation related genes also needs to be clarified. In terms of piRNA and mRNA levels, the present study found that exogenous TTX administration would feminize gene expression in the testis, female-dominant genes foxl2, cyp19a1/ cyp19a1a and gdf9 were up regulated in the MTG. The cyp19a1 is the key gene encoding cytochrome P450 aromatase, which is responsible for estrogen production in the ovary [60,61,62], the increase of cyp19a1a expression might lead to an increased level of endogenous estrogen E2 (17beta-estrodiol), excessive estrogen also induces the initiation of cell apoptotic pathways by Stabilizing Schlafen-12 Protein Turnover , leading to a significant increase in cell apoptosis in the testis . Estrogen is also considered a key factor in sex differentiation [64, 65]. It is well documented that foxl2 and cyp19a1 play a vital role in promoting ovarian development and maintaining feminization afterward in teleosts [38, 63, 64, 66, 67]. Growth differentiation factor 9 (gdf9) is a member of the TGFβ superfamily. As an oocyte-derived growth factor, gdf9 plays key role in regulating follicle development [68, 69]. Steroidogenesis is critical for gamete maturation in the teleost, more essentially oogenesis , except cyp19a1, hsd17b1 was found to be up-regulated in the MTG. Moreover, KEGG analysis of differentially expressed genes/piRNAs in MTG also enriched in Wnt signaling and oocyte meiosis pathway. Studies have shown that the components of the Wnt signaling pathway are responsible for reproductive functions, including the regulation of follicular maturation and the production of steroids in the ovaries [71,72,73,74]. The B-type cyclins take important roles in oocyte meiosis , the ccnb1 and ccnb2 are essential for the is critically required for the proliferation of gonocytes , and both were upregulated by the TTX. These up-regulated key female dominant genes and pathways induced feminization of testicular gene expression, which may need to mimic the female endocrine environment to counter the damaging effects of TTX on the testis.
The TGFβ signaling pathway regulates growth, division and proliferation in the ovary [77, 78], one of the TGFβ signaling pathway member amhr2 were up-regulated after TTX administration FTG, the amhr2 was associated with sex determination in Takifugu rubripes, and showed regulation of germ cell proliferation [77,78,79,80]. Combined with the differential analysis base on piRNA sequencing, the up-regulated gene involved in the down-regulated piRNA in the ovary, indicated TTX may impact the gene expression of the ovary through the regulation of piRNA and showed a positive effect.
The present study compared differentially expressed piRNAs and genes in MTG vs MCG and FTG vs FCG, these two groups shared 298 common piRNAs and 44 common genes with the changes. KEGG pathway enrichment analysis of these 44 common genes showed that TTX directly affects piRNAs and mRNAs changes by influencing nucleotide, lipid, glycan, carbohydrate and protein metabolism as well as cell genetic information processes and signals. On the other hand, compared with the ovary, the testis produced more metabolic pathways in response to exogenous TTX. KEGG pathway enrichment analysis of 401 DEGs significantly differentially expressed in MTG vs MCG group, indicating the testis have richer metabolic and/or transport pathways than ovary to against the accumulation of TTX. These results indicated that TTX plays a detrimental role in testicular development, the testes themselves may be repellent to TTX.
Above all, our research revealed that the response of ovary and testis to TTX administration was largely different, ovary are more tolerant to TTX (has fewer piRNA and mRNA expression changes), whereas testes are more senstive to TTX. Exogenous TTX administration mainly resulted in testis gene expression pattern's feminization, up regulation of more metabolic pathways and increased apoptosis. These data will deepen our understanding on the accumulation of TTX sexual dimorphism in Takifugu.
In conclusion, the present research highlights the differentially expressed piRNAs and genes after TTX administration in gonads. The response of ovary and testis to TTX administration was largely different, ovary are more tolerant to TTX, whereas testes are more sensitive to TTX. The differentially expressed piRNAs potential regulated genes and the differentially expressed gene were partly overlapping. The regulatory role of TTX may through piRNA on genes related to reproduction support the notion that piRNAs play crucial roles in reproduction and the TTX novel function in pufferfish reproduction. Exogenous TTX administration mainly resulted in testis gene expression pattern's feminization, up regulation of more metabolic pathways and increased apoptosis in Takifugu rubripes. Further confirmed that excessive TTX is harmful to testicular development. Overall, this research will deepen our understanding on the accumulation of TTX sexual dimorphism in Takifugu.
Materials and methods
Animals and chemicals
The experimental animals, artificial cultured marine pufferfish Takifugu flavidus juveniles (6-month-old, 100.0 ± 10.0 g body weight, 100 individuals) were purchased from Shanghai Aquatics Research Institution. Crystalline of TTX (Supelco, USA) was used in the administration experiment and as a standard sample for the liquid-chromatography-tandem mass spectrometry (LC–MS/MS) analysis. All other chemicals were reagent grade (Sinopharm, China). The present study on Takifugu flavidus was carried out following the relevant guidelines and regulations (following the ARRIVE guidelines ). The protocols were approved by the academic ethics committee of Shanghai Ocean University.
TTX administration and sample preparation
The administration experiment was performed at the laboratory of the Shanghai Ocean University. The TTX was directly dissolved in 0.1% acetic acid to 1 mg/ml, then diluted with 0.7% saline solution to 0.25 μg/μl. The test fish received an intramuscular administration of 0.25 μg TTX/g body weight into the caudal muscle as described previously and were maintained in an 80-L oxygenated plastic tank of aerated artificial seawater at 20 ℃ . After 8 h of administration, twelve fish of control and TTX-treated group were randomly collected, followed ice-cold anesthetized and dissected. The gonads were dissected immediately and separated into two parts, 1/4 gonad was fixed in 4% PFA for histological and immunofluorescence histochemical observation; another 3/4 of gonad were stored at -80 ℃ until RNA extraction. Kidney, cholecyst, skin, liver, heart, and muscle were also dissected and stored at -20 ℃ until RNA exaction.
TTX extraction and quantification
Tissues (kidney, cholecyst, skin, liver, heart, and muscles) were extracted with 0.1 M HCl as reported previously [48, 82]. The extracts from the tissues were defatted with ethoxyethane and centrifuged at 12,000 rpm for 10 min. The resulting supernatants were ultrafiltered (Amicon Ultra, Millipore, USA), then submitted to LC–MS/MS (Liquid Chromatography with tandem mass spectrometry) for TTX determination with some modification [83, 84]. In the LC–MS/MS analysis, chromatography was carried out using a Xevo TQ-XS system (Waters, Milford, USA) with an XBridge Amide column (2.1 mm × 100 mm, particle size 1.7 μm,) and mobile phase comprising 45% 0.1 mM methanoic acid and 55% acetonitrile at a flow rate of 350 μL/min. The eluate was introduced into a MassLynxTM N detector (Waters, Milford, USA) in which the TTX was ionized by positive-mode electrospray ionization with a desolvation temperature of 550 °C, the fragment ion at m/z 162 that results from the dissociation of the parent ion of TTX at m/z 320 was detected by the multiple reaction monitoring modes. These tissue extracts were spiked with TTX standard [83, 85, 86].
Histological observation and TTX immunofluorescence histochemistry
For histological observation and immunofluorescence histochemical analyses, the Takifugu flavidus gonads were dissected and fixed in 4% PFA (Paraformaldehyde) for 12 h at 4℃, then dehydrated, embedded in paraffin, and sectioned at 5 μm. Hematoxylin–eosin staining and immunofluorescence histochemical staining were performed as described previously [87, 88]. The antibody against TTX was purchased from CD (Creative diagnostics, USA) and was diluted at 1:200, Hoechst 33,258 (Sigma, USA) was diluted at 1:1000 staining the nucleic acid. A Nikon ECLIPES 80i light microscope (Nikon, Japan) was used to image the H&E (hematoxylin–eosin staining) stained sections. Leica TCS SP8 X (Leica, German) confocal microscope was used to image fluorescence-stained sections. The TTX signals ratio of gonads was statistical analysis of the three different areas of each sample using ImageJ (NIH, USA) and Graph Prism 8.0 (San Diego, USA).
TdT-mediated dUTP-biotin nick end labeling (TUNEL) technique was applied to evaluate the apoptotic response of TTX administration in the ovary and testis. Reactions were performed on Sects. (5 µm) of 4% Paraformaldehyde (PFA)-fixed (overnight at 4 °C) and paraffin-embedded ovaries and testes. Apoptotic cells with DNA breaks were detected using Colorimetric TUNEL Apoptosis Assay Kits (#C1091, Beyotime Institute of Biotechnology, Shanghai, China). TUNEL assay was then performed according to the instructions by the manufacturer. Images were used a Nikon ECLIPES 80i light microscope (Nikon, Japan). The apoptosis signal ratio of gonads was a statistical analysis of the three different areas of each sample under TUNEL staining by using ImageJ (NIH, USA) and Graph Prism8.0 (San Diego, USA) as described before.
After sex is determined by histological observation, 3 or 4 gonad samples are mixed into one pool for each FTG, MTG, FCG and MCG group. Total RNAs of each group were extracted using RNAiso Plus (Takara, Japan) . The agarose electrophoresis was used to detect the integrity of total RNA, and the Agilent 2100 Bioanalyzer (Agilent Technologies, USA) was used to detect the concentration and purity of total RNA. Then, four cDNA libraries were constructed following the manufacturer’s recommendations and sequenced on the Illumina HiSeq platform, generating 150-bp paired-end reads, a summary of the raw and clean data were shown in Table 1. High-quality clean reads from each sample were mapped onto the assembled transcript sequences using Bowtie software (Version 2.2.2)  with default parameters. The RSEM program  was then used to estimate the expression abundance of the transcripts. The total number of mapped reads for each transcript was determined and then normalized to determine FPKM (Fragments Per Kilobase Million).
For piRNA, four small RNA libraries were constructed using TruSeq Small RNA Donor Prep Kits (Illumina, USA). The small RNA sequencing analysis was performed on the Illumina Hiseq 2500 sequencing platform by the Novegene Company (Beijing, China). The small RNA with 15–41 nt were selected after screening the sequence of Q-score > 20, removing the N base, and filtering out the inappropriate lengths. The clean reads were classified and annotated by comparing them with relevant databases using Bowtie and BLAST . The small RNAs were classified by filtering clean reads using the Rfam database  as well as aligning with the genome of fTakRub1.2 genome (closely related species Takifugu rubripes, fTakRub1.2, E-value < 0.01, Bioproject number PRJEB31988/PRJNA543527, Assembly under accession number GCA_901000725.2 at NCBI, National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov). Sequences without aligned transcripts were matched with the miRBase database by zero-base mismatch [93, 94], and the miRNA in the sequence was filtered as far as possible. Then, the remaining sequences with 26–32 nt unannotated reads analyzed by Piano were predicted as novel piRNAs  based on using the structure and sequence features of transposon-piRNAs interactions, the novel piRNA aligned with piRBase (http://bigdata.ibp.ac.cn/piRBase/) using bowtie by zero-base mismatch , the aligned sequences were considered as known piRNAs. The genome was used to match with the clean sequences for piRNA to determine various small RNAs. Annotated non-coding RNAs included rRNAs, tRNAs, snRNA, miRNA, piRNA, etc. These aligned RNAs were subjected to BLAST against Rfamv.10.1 (http://www.sanger.ac.uk/software/Rfam) , GenBank (http://www.ncbi.nlm.nih.gov/genbank/) databases and miRBase v.21 database (http://www.mirbase.org/) . The known piRNAs were identified by aligning against the piRBase release v2.0 database (http://www.regulatoryrna.org/database/piRNA/) , and the expression patterns of known piRNAs and novel piRNA in four groups were analyzed.
Bioinformatic analysis and data statistics
For mRNA, DESeq software standardized the counts’ number of each group gene and calculated the fold of difference and used NB (negative binomial distribution test) to test the significance of the difference in the number of reads . The expression level of genes can be calculated using FPKM (Fragments Per Kilobase Million). The formula is FPKM = (total exons reads)/ (mapped reads, millions × exon length, kb) . Finally, the differential protein-coding genes were screened according to the different folds and the different significance test results. DEGs were identified with a threshold of P < 0.05 and |log2 (fold change) |> 1 based on the BaseMean method . Due to no biological replicates in this study, the P value was calculated with the DEG algorithm  in the R package including the Audic-Claverie statistic as above . The differential transcripts/genes were classified according to the official annotation of the Gene Ontology (GO)  and Kyoto Encyclopedia of Genes and Genomes (KEGG) [32, 102,103,104] was analyzed.
Differentially expressed piRNAs (DEpiRNAs) were identified with a threshold of P < 0.05 and |log2 (fold change) |> 1. Due to no biological replicates in this study, the P value was calculated with the DEG algorithm  in the R package including the Audic-Claverie statistic . In addition, a quantitative analysis of piRNAs was performed. The expression level of piRNAs is directly reflected by their abundance. The higher the abundance of piRNAs is, the higher the expression level is. In piRNA sequencing analysis, the expression level of piRNAs can be calculated by locating the piRNA sequence of Takifugu flavidus using the transcripts per million (TPM). The formula is TPM = (number of reads aligned for each piRNA) / (sample Total aligned reads) × 106 . Differential expression analysis aims to identify differentially expressed piRNAs among different groups based on the TPM method . After obtaining differential expression, all the expressed piRNA sequences in four groups are compared with the human genome using sRNAmapper software (version 1.0.5) , the parameters are default, the obtained map file, and then use proTRAC software (version 2.4.2) to do piRNA cluster prediction . Then the target genes of differentially expressed piRNAs were predicted using the MiRanda software with a score of 150 or higher and free energy -30 or lower .
The enrichment analysis of the Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was performed as previously reported [32, 101,102,103,104]. The false discovery rate (FDR) was used to adjust the P value and GO or KEGG terms (FDR < 0.01) were considered significantly enriched.
Gonadal piRNA and mRNA expression were also assayed using qRT-PCR. Total RNA was extracted following the above description. cDNAs were synthesized using HiScript® II 1st Strand cDNA Synthesis Kit (Vazyme, China). Hieff UNICON® Universal Blue qPCR SYBR Green Master Mix (Yeason, China) was used for quantitative qRT-PCR assays. The gene names and the primers are listed in Table S9. Gadph was used as an internal control [89, 109, 110]. The relative abundances of mRNA transcripts were evaluated using the formula: R = 2−ΔΔCt .
For piRNA, 1 μg of total RNA above as a template for synthesis of the first-strand cDNA with a Mir-X miRNA First-Strand Synthesis (Takara, Japan) following the manufacturer’s instructions. The expression of the piRNAs used a miScript SYBR Green PCR Kit (Takara, Japan) and followed the manufacturer’s instructions on a CFX96 Real-Time PCR System (Bio-rad, USA). The primers for piRNA qRT-PCR were listed in Table S9. We used U6 as an internal reference gene to control for differences among samples. Each group sample of four or three mixed fishes was run triple times and the relative quantification of piRNA expression was calculated using the 2−ΔΔCt method .
mRNA and piRNA expression levels of each sample were normalized to the levels of gadph and U6 gene in the same samples, respectively. the control group was set as 1 and was omitted in the figures, the mRNA and piRNA level folds to control were shown in figures. This method was reported in previous research .
Except where otherwise indicated, Statistical analyses were performed in the GraphPad Prism 9.0 (San Diego, California USA) environment. Each measurement was repeated at least three times. data are presented as mean ± SEM, and the P values were determined by two-tailed Student’s t-tests. P < 0.05 was considered to be statistically significant.
Availability of data and materials
All sequencing raw data have been deposited in the Sequence Read Archive (SRA) database under Bioproject number PRJNA769942 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA769942). the Takifugu rubripes reference genome fTakRub1.2 under Bioproject number PRJEB31988/PRJNA543527, Assembly under accession number GCA_901000725.2 at NCBI, National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov).
Miyazawa K, Noguchi T. Distribution and origin of tetrodotoxin. J Toxicol-Toxin Rev. 2001;20(1):11–33.
Noguchi T, Arakawa O, Takatani T. TTX accumulation in pufferfish. Comp Biochem Physiol Part D Genomics Proteomics. 2006;1(1):145–52.
Itoi S, Ishizuka K, Mitsuoka R, Takimoto N, Yokoyama N, Detake A, Takayanagi C, Yoshikawa S, Sugita H. Seasonal changes in the tetrodotoxin content of the pufferfish Takifugu niphobles. Toxicon. 2016;114:53–8.
Ikeda K, Emoto Y, Tatsuno R, Wang JJ, Ngy L, Taniyama S, Takatani T, Arakawa O. Maturation-associated changes in toxicity of the pufferfish Takifugu poecilonotus. Toxicon. 2010;55(2–3):289–97.
Itoi S, Yoshikawa S, Tatsuno R, Suzuki M, Asahina K, Yamamoto S, Takanashi S, Takatani T, Arakawa O, Sakakura Y, et al. Difference in the localization of tetrodotoxin between the female and male pufferfish Takifugu niphobles, during spawning. Toxicon. 2012;60(6):1000–4.
Williams BL. Behavioral and chemical ecology of marine organisms with respect to tetrodotoxin. Mar Drugs. 2010;8(3):381–98.
Itoi S, Kozaki A, Komori K, Tsunashima T, Noguchi S, Kawane M, Sugita H. Toxic Takifugu pardalis eggs found in Takifugu niphobles gut: Implications for TTX accumulation in the pufferfish. Toxicon. 2015;108:141–6.
Nagashima Y, Toyoda M, Hasobe M, Shimakura K, Shiomi K. In vitro accumulation of tetrodotoxin in pufferfish liver tissue slices. Toxicon. 2003;41(5):569–74.
Mahmud Y, Okada K, Takatani T, Kawatsu K, Hamano Y, Arakawa O, Noguchi T. Intra-tissue distribution of tetrodotoxin in two marine puffers Takifugu vermicularis and Chelonodon patoca. Toxicon. 2003;41(1):13–8.
Venkatesh B, Lu SQ, Dandona N, See SL, Brenner S, Soong TW. Genetic basis of tetrodotoxin resistance in pufferfishes. Curr Biol. 2005;15(22):2069–72.
Ikeda K, Murakami Y, Emoto Y, Ngy L, Taniyama S, Yagi M, Takatani T, Arakawa O. Transfer profile of intramuscularly administered tetrodotoxin to non-toxic cultured specimens of the pufferfish Takifugu rubripes. Toxicon. 2009;53(1):99–103.
Noguchi T, Arakawa O. Tetrodotoxin-distribution and accumulation in aquatic organisms, and cases of human intoxication. Mar Drugs. 2008;6(2):220–42.
Noguchi T, Arakawa O, Takatani T. Toxicity of pufferfish Takifugu rubripes cultured in netcages at sea or aquaria on land. Comp Biochem Physiol D: Genomics Proteomics. 2006;1(1):153–7.
Feroudj H, Matsumoto T, Kurosu Y, Kaneko G, Ushio H, Suzuki K, Kondo H, Hirono I, Nagashima Y, Akimoto S, et al. DNA microarray analysis on gene candidates possibly related to tetrodotoxin accumulation in pufferfish. Toxicon. 2014;77:68–72.
Okita K, Takatani T, Nakayasu J, Yamazaki H, Sakiyama K, Ikeda K, Arakawa O, Sakakura Y. Comparison of the localization of tetrodotoxin between wild pufferfish Takifugu rubripes juveniles and hatchery-reared juveniles with tetrodotoxin administration. Toxicon. 2013;71:128–33.
Itoi S, Suzuki M, Asahina K, Sawayama E, Nishikubo J, Oyama H, Takei M, Shiibashi N, Takatani T, Arakawa O, et al. Role of maternal tetrodotoxin in survival of larval pufferfish. Toxicon. 2018;148:95–100.
Itoi S, Yoshikawa S, Asahina K, Suzuki M, Ishizuka K, Takimoto N, Mitsuoka R, Yokoyama N, Detake A, Takayanagi C, et al. Larval pufferfish protected by maternal tetrodotoxin. Toxicon. 2014;78:35–40.
Stahl F, Hitzmann B, Mutz K, Landgrebe D, Lübbecke M, Kasper C, Walter J, Scheper T. Transcriptome analysis. Adv Biochem Eng Biotechnol. 2012;127:1–25.
Kalucka J, de Rooij L, Goveia J, Rohlenova K, Dumas SJ, Meta E, Conchinha NV, Taverna F, Teuwen LA, Veys K, et al. Single-Cell Transcriptome Atlas of Murine Endothelial Cells. Cell. 2020;180(4):764-779.e720.
Kim VN, Han J, Siomi MC. Biogenesis of small RNAs in animals. Nat Rev Mol Cell Biol. 2009;10(2):126–39.
Iwasaki YW, Siomi MC, Siomi H. PIWI-Interacting RNA: Its Biogenesis and Functions. Annu Rev Biochem. 2015;84:405–33.
Ozata DM, Gainetdinov I, Zoch A, O’Carroll D, Zamore PD. PIWI-interacting RNAs: small RNAs with big functions. Nat Rev Genet. 2019;20(2):89–108.
Perera BPU, Tsai ZT, Colwell ML, Jones TR, Goodrich JM, Wang K, Sartor MA, Faulk C, Dolinoy DC. Somatic expression of piRNA and associated machinery in the mouse identifies short, tissue-specific piRNA. Epigenetics. 2019;14(5):504–21.
Liu Y, Dou M, Song X, Dong Y, Liu S, Liu H, Tao J, Li W, Yin X, Xu W. The emerging role of the piRNA/piwi complex in cancer. Mol Cancer. 2019;18(1):123.
Qiu W, Guo X, Lin X, Yang Q, Zhang W, Zhang Y, Zuo L, Zhu Y, Li CR, Ma C, et al. Transcriptome-wide piRNA profiling in human brains of Alzheimer’s disease. Neurobiol Aging. 2017;57:170–7.
Kim KW. PIWI Proteins and piRNAs in the Nervous System. Mol Cells. 2019;42(12):828–35.
Lee EJ, Banerjee S, Zhou H, Jammalamadaka A, Arcila M, Manjunath BS, Kosik KS. Identification of piRNAs in the central nervous system. RNA. 2011;17(6):1090–9.
Britton C, Laing R, Devaney E. Small RNAs in parasitic nematodes - forms and functions. Parasitology. 2020;147(8):855–64.
Menor MS, Baek K, Poisson G. Prediction of mature microRNA and piwi-interacting RNA without a genome reference or precursors. Int J Mol Sci. 2015;16(1):1466–81.
Tan CH, Lee TC, Weeraratne SD, Korzh V, Lim TM, Gong Z. Ziwi, the zebrafish homologue of the Drosophila piwi: co-localization with vasa at the embryonic genital ridge and gonad-specific expression in the adults. Gene Expr Patterns. 2002;2(3–4):257–60.
Wei Z, Liu X, Zhang H. Identification and characterization of piRNA-like small RNAs in the gonad of sea urchin (Strongylocentrotus nudus). Mar Biotechnol (NY). 2012;14(4):459–67.
Waiho K, Fazhan H, Zhang Y, Li S, Zhang Y, Zheng H, Ikhwanuddin M, Ma H. Comparative profiling of ovarian and testicular piRNAs in the mud crab Scylla paramamosain. Genomics. 2020;112(1):323–31.
Lau NC, Seto AG, Kim J, Kuramochi-Miyagawa S, Nakano T, Bartel DP, Kingston RE. Characterization of the piRNA complex from rat testes. Science. 2006;313(5785):363–7.
Yang L, Ge Y, Cheng D, Nie Z, Lv Z. Detection of piRNAs in whitespotted bamboo shark liver. Gene. 2016;590(1):51–6.
Gong J, Zhang Q, Wang Q, Ma Y, Du J, Zhang Y, Zhao X. Identification and verification of potential piRNAs from domesticated yak testis. Reproduction. 2018;155(2):117–27.
Li B, He X, Zhao Y, Bai D, Bou G, Zhang X, Su S, Dao L, Liu R, Wang Y, et al. Identification of piRNAs and piRNA clusters in the testes of the Mongolian horse. Sci Rep. 2019;9(1):5022.
He X, Li B, Fu S, Wang B, Qi Y, Da L, Te R, Sun S, Liu Y, Zhang W. Identification of piRNAs in the testes of Sunite and Small-tailed Han sheep. Anim Biotechnol. 2021;32(1):13–20.
Wu S, Guo J, Zhu L, Yang J, Chen S, Yang X. Identification and characterisation of microRNAs and Piwi-interacting RNAs in cockerels’ spermatozoa by Solexa sequencing. Br Poult Sci. 2018;59(4):371–80.
Houwing S, Kamminga LM, Berezikov E, Cronembold D, Girard A, van den Elst H, Filippov DV, Blaser H, Raz E, Moens CB, et al. A role for Piwi and piRNAs in germ cell maintenance and transposon silencing in Zebrafish. Cell. 2007;129(1):69–82.
Zhou Y, Zhong H, Xiao J, Yan J, Luo Y, Gan X, Yu F. Identification and comparative analysis of piRNAs in ovary and testis of Nile tilapia (Oreochromis niloticus). Genes Genomics. 2016;38(6):519–27.
Wang CL, Wang ZP, Wang JQ, Li MY, Chen XW. Identification of candidate piRNAs in the gonads of Paralichthys olivaceus (Japanese flounder). Zool Res. 2016;37(5):301–6.
Zhang B, Zhao N, Jia L, Che J, He X, Liu K, Bao B. Identification and application of piwi-interacting RNAs from seminal plasma exosomes in Cynoglossus semilaevis. BMC Genomics. 2020;21(1):302.
Gonzalez LE, Tang X, Lin H. Maternal Piwi Regulates Primordial Germ Cell Development to Ensure the Fertility of Female Progeny in Drosophila. Genetics 2021.
Hara K, Fujii T, Suzuki Y, Sugano S, Shimada T, Katsuma S, Kawaoka S. Altered expression of testis-specific genes, piRNAs, and transposons in the silkworm ovary masculinized by a W chromosome mutation. BMC Genomics. 2012;13:119.
Yang Q, Li R, Lyu Q, Hou L, Liu Z, Sun Q, Liu M, Shi H, Xu B, Yin M, et al. Single-cell CAS-seq reveals a class of short PIWI-interacting RNAs in human oocytes. Nat Commun. 2019;10(1):3389.
Wongwarangkana C, Fujimori KE, Akiba M, Kinoshita S, Teruya M, Nezuo M, Masatoshi T, Watabe S, Asakawa S. Deep sequencing, profiling and detailed annotation of microRNAs in Takifugu rubripes. BMC Genomics. 2015;16(1):457.
Maruta S, Yamaoka K, Yotsu-Yamashita M. Two critical residues in p-loop regions of puffer fish Na+ channels on TTX sensitivity. Toxicon. 2008;51(3):381–7.
Matsumoto T, Kiriake A, Ishizaki S, Watabe S, Nagashima Y. Biliary excretion of tetrodotoxin in the cultured pufferfish Takifugu rubripes juvenile after intramuscular administration. Toxicon. 2015;93:98–102.
Wang J, Araki T, Tatsuno R, Nina S, Ikeda K, Hamasaki M, Sakakura Y, Takatani T, Arakawa O. Transfer profile of intramuscularly administered tetrodotoxin to artificial hybrid specimens of pufferfish. Takifugu rubripes and Takifugu niphobles Toxicon. 2011;58(6–7):565–9.
Tatsuno R, Gao W, Ibi K, Mine T, Okita K, Nishihara GN, Takatani T, Arakawa O. Profile differences in tetrodotoxin transfer to skin and liver in the pufferfish Takifugu rubripes. Toxicon. 2017;130:73–8.
Tóth KF, Pezic D, Stuwe E, Webster A. The piRNA Pathway Guards the Germline Genome Against Transposable Elements. Adv Exp Med Biol. 2016;886:51–77.
Czech B, Munafò M, Ciabrelli F, Eastwood EL, Fabry MH, Kneuss E, Hannon GJ. piRNA-Guided Genome Defense: From Biogenesis to Silencing. Annu Rev Genet. 2018;52:131–57.
Wen X, Wang D, Li X, Zhao C, Wang T, Qian X, Yin S. Differential expression of two Piwil orthologs during embryonic and gonadal development in pufferfish, Takifugu fasciatus. Comp Biochem Physiol B Biochem Mol Biol. 2018;219–220:44–51.
Stegh AH. Targeting the p53 signaling pathway in cancer therapy - the promises, challenges and perils. Expert Opin Ther Targets. 2012;16(1):67–83.
Zhang H, Zhang X, Li X, Meng WB, Bai ZT, Rui SZ, Wang ZF, Zhou WC, Jin XD. Effect of CCNB1 silencing on cell cycle, senescence, and apoptosis through the p53 signaling pathway in pancreatic cancer. J Cell Physiol. 2018;234(1):619–31.
Liu J, Zhang C, Wang J, Hu W, Feng Z. The Regulation of Ferroptosis by Tumor Suppressor p53 and its Pathway. Int J Mol Sci. 2020;21(21):8387.
Liang C, Zhang X, Yang M, Dong X. Recent Progress in Ferroptosis Inducers for Cancer Therapy. Adv Mater. 2019;31(51):e1904197.
Li J, Cao F, Yin HL, Huang ZJ, Lin ZT, Mao N, Sun B, Wang G. Ferroptosis: past, present and future. Cell Death Dis. 2020;11(2):88.
Wang XG, Orban L. Anti-Müllerian hormone and 11 beta-hydroxylase show reciprocal expression to that of aromatase in the transforming gonad of zebrafish males. Dev Dyn. 2007;236(5):1329–38.
Trant JM, Gavasso S, Ackers J, Chung BC, Place AR. Developmental expression of cytochrome P450 aromatase genes (CYP19a and CYP19b) in zebrafish fry (Danio rerio). J Exp Zool. 2001;290(5):475–83.
Callard GV, Tchoudakova AV, Kishida M, Wood E. Differential tissue distribution, developmental programming, estrogen regulation and promoter characteristics of cyp19 genes in teleost fish. J Steroid Biochem Mol Biol. 2001;79(1–5):305–14.
Hsu HJ, Hsu NC, Hu MC, Chung BC. Steroidogenesis in zebrafish and mouse models. Mol Cell Endocrinol. 2006;248(1–2):160–3.
Li D, Chen J, Ai Y, Gu X, Li L, Che D, Jiang Z, Li L, Chen S, Huang H, et al. Estrogen-Related Hormones Induce Apoptosis by Stabilizing Schlafen-12 Protein Turnover. Mol Cell. 2019;75(6):1103-1116.e1109.
Wang DS, Kobayashi T, Zhou LY, Paul-Prasanth B, Ijiri S, Sakai F, Okubo K, Morohashi K, Nagahama Y. Foxl2 up-regulates aromatase gene transcription in a female-specific manner by binding to the promoter as well as interacting with ad4 binding protein/steroidogenic factor 1. Mol Endocrinol. 2007;21(3):712–25.
Nagahama Y. Gonadal steroid hormones: major regulators of gonadal sex differentiation and gametogenesis in fish. In: International Symposium on the Reproductive Physiology of Fish: 2000; 2000.
Guiguen Y, Fostier A, Piferrer F, Chang CF. Ovarian aromatase and estrogens: a pivotal role for gonadal sex differentiation and sex change in fish. Gen Comp Endocrinol. 2010;165(3):352–66.
Li MH, Yang HH, Li MR, Sun YL, Jiang XL, Xie QP, Wang TR, Shi HJ, Sun LN, Zhou LY, et al. Antagonistic roles of Dmrt1 and Foxl2 in sex differentiation via estrogen production in tilapia as demonstrated by TALENs. Endocrinology. 2013;154(12):4814–25.
Chen W, Liu L, Ge W. Expression analysis of growth differentiation factor 9 (Gdf9/gdf9), anti-müllerian hormone (Amh/amh) and aromatase (Cyp19a1a/cyp19a1a) during gonadal differentiation of the zebrafish. Danio rerio Biol Reprod. 2017;96(2):401–13.
Yu H, Wang Y, Wang M, Liu Y, Cheng J, Zhang Q. Growth differentiation factor 9 (gdf9) and bone morphogenetic protein 15 (bmp15) are potential intraovarian regulators of steroidogenesis in Japanese flounder (Paralichthys olivaceus). Gen Comp Endocrinol. 2020;297:113547.
Rajakumar A, Senthilkumaran B. Steroidogenesis and its regulation in teleost-a review. Fish Physiol Biochem. 2020;46(3):803–18.
Safian D, Bogerd J, Schulz RW. Regulation of spermatogonial development by Fsh: The complementary roles of locally produced Igf and Wnt signaling molecules in adult zebrafish testis. Gen Comp Endocrinol. 2019;284:113244.
Hernandez Gifford JA. The role of WNT signaling in adult ovarian folliculogenesis. Reproduction. 2015;150(4):R137-148.
Kossack ME, Draper BW. Genetic regulation of sex determination and maintenance in zebrafish (Danio rerio). Curr Top Dev Biol. 2019;134:119–49.
Sreenivasan R, Jiang J, Wang X, Bártfai R, Kwan HY, Christoffels A, Orbán L. Gonad differentiation in zebrafish is regulated by the canonical Wnt signaling pathway. Biol Reprod. 2014;90(2):45.
Bouftas N, Wassmann K. Cycling through mammalian meiosis: B-type cyclins in oocytes. Cell Cycle. 2019;18(14):1537–48.
Tang JX, Li J, Cheng JM, Hu B, Sun TC, Li XY, Batool A, Wang ZP, Wang XX, Deng SL, et al. Requirement for CCNB1 in mouse spermatogenesis. Cell Death Dis. 2017;8(10):e3142.
Wu GC, Chiu PC, Lyu YS, Chang CF. The expression of amh and amhr2 is associated with the development of gonadal tissue and sex change in the protandrous black porgy. Acanthopagrus schlegeli Biol Reprod. 2010;83(3):443–53.
Mullen RD, Ontiveros AE, Moses MM, Behringer RR. AMH and AMHR2 mutations: A spectrum of reproductive phenotypes across vertebrate species. Dev Biol. 2019;455(1):1–9.
Pfennig F, Standke A, Gutzeit HO. The role of Amh signaling in teleost fish–Multiple functions not restricted to the gonads. Gen Comp Endocrinol. 2015;223:87–107.
Kamiya T, Kai W, Tasumi S, Oka A, Matsunaga T, Mizuno N, Fujita M, Suetake H, Suzuki S, Hosoya S, 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(7):e1002798.
Percie du Sert N, Hurst V, Ahluwalia A, Alam S, Avey MT, Baker M, Browne WJ, Clark A, Cuthill IC, Dirnagl U, et al. The ARRIVE guidelines 2.0: updated guidelines for reporting animal research. BMJ Open Sci. 2020;4(1):e100115.
Gao W, Kanahara Y, Yamada M, Tatsuno R, Yoshikawa H, Doi H, et al. Contrasting toxin selectivity between the marine pufferfish Takifugu pardalis and the freshwater pufferfish Pao suvattii. Toxins (Basel). 2019;11(8):470.
Yang G, Xu J, Liang S, Ren D, Yan X, Bao B. A novel TTX-producing Aeromonas isolated from the ovary of Takifugu obscurus. Toxicon. 2010;56(3):324–9.
Zhang N, Wei F, Baolong B. Extraction and detection of tetrodotoxin in Takifugu niphobles and its symbiotic bacteria. Journal of shanghai ocean university. 2017;26(06):801–7.
Ma T, Xiaoling G, Guan Z, Baolong B. Analysis of the tetrodotoxin content of tissues in Yongeichthys criniger. J Shanghai Ocean Univ. 2014;23(5):675–9.
Wei F, Ma T, Gong X, Zhang N, Bao B. Identification of tetrodotoxin-producing bacteria from goby Yongeichthys criniger. Toxicon. 2015;104:46–51.
Xie QP, He X, Sui YN, Chen LL, Sun LN, Wang DS. Haploinsufficiency of SF-1 Causes Female to Male Sex Reversal in Nile Tilapia. Oreochromis niloticus Endocrinology. 2016;157(6):2500–14.
Mao J, Gong X, Bao B. Identification of TTX Anisakis pegreffii parasites in Takifugu xanthopterus from the East China Sea. J Shanghai Ocean Univ. 2020;29(4):585–92.
Gao L, Penglee R, Huang Y, Yi X, Wang X, Liu L, et al. CRISPR/Cas9-induced nos2b mutant zebrafish display behavioral abnormalities. Genes Brain Behav. 2021;20(5):e12716.
Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics 2010. Chapter 11:Unit 11.17.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):1–16.
Griffiths-Jones S, Bateman A, Marshall M, Khanna A, Eddy SR. Rfam: an RNA family database. Nucleic Acids Res. 2003;31(1):439–41.
Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Database issue):D140-144.
Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36(Database issue):D154-158.
Wang K, Liang C, Liu J, Xiao H, Huang S, Xu J, Li F. Prediction of piRNAs using transposon interaction and a support vector machine. BMC Bioinformatics. 2014;15(1):419.
Wang J, Zhang P, Lu Y, Li Y, Zheng Y, Kan Y, Chen R, He S. piRBase: a comprehensive database of piRNA sequences. Nucleic Acids Res. 2019;47(D1):D175-d180.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Tino P. Basic properties and information theory of Audic-Claverie statistic for analyzing cDNA arrays. BMC Bioinformatics. 2009;10:310.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.
Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 1999;27(1):29–34.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480-484.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545-d551.
Sun J, Wang S, Li C, Ren Y, Wang J. Novel expression profiles of microRNAs suggest that specific miRNAs regulate gene expression for the sexual maturation of female Schistosoma japonicum after pairing. Parasit Vectors. 2014;7:177.
Rosenkranz D, Han CT, Roovers EF, Zischler H, Ketting RF. Piwi proteins and piRNAs in mammalian oocytes and early embryos: From sample to sequence. Genom Data. 2015;5:309–13.
Rosenkranz D, Zischler H. proTRAC–a software for probabilistic piRNA cluster detection, visualization and analysis. BMC Bioinformatics. 2012;13:5.
John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human microRNA targets PLoS biol. 2004;2(11):e363.
Shao C, Bao B, Xie Z, Chen X, Li B, Jia X, Yao Q, Orti G, Li W, Li X. The genome and transcriptome of Japanese flounder provide insights into flatfish asymmetry. Nat Genet. 2017;49(1):119–24.
Jiang C. Selection of appropriate reference genes for real-time quantitative PCR analysis in Takifugu flavidus. Journal of Fisheries Research. 2020;42(2):105.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.
Luo Y, Li C, Landis AG, Wang G, Stoeckel J, Peatman E. Transcriptomic profiling of differential responses to drought in two freshwater mussel species, the giant floater Pyganodon grandis and the pondhorn Uniomerus tetralasmus. PLoS ONE. 2014;9(2):e89481.
The authors would like to thank Dr. Qing-ping Xie (Institute of Hydrobiology, Zhejiang Academy of Agricultural Sciences) for his critical reading and revision suggestions.
This work was supported by grants from the National Natural Science Foundation of China (31872546, 41176108) and the Key Innovation Projects of Shanghai Municipal Education Commission (14zz145).
Ethics approval and consent to participate
Shanghai Aquatics Research Institute (shanghai, China) provided artificial cultured marine pufferfish Takifugu flavidus juveniles for research. The methods for Takifugu flavidus research was carried out following the relevant guidelines and regulations The present study on Takifugu flavidus was carried out following the relevant guidelines and regulations (following the ARRIVE guidelines). The protocols were approved by the academic ethics committee of Shanghai Ocean University.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Relative expression of all piRNAs identified and DEpiRNAs in gonads.
Prediction of piRNA clusters.
List of preidicted target genes of DEpiRNAs.
Relative expression of mRNAs and DEGs in gonads.
KEGG analysis of DEGs in gonads.
GO analysis of DEGs in gonads
DEpiRNA target gene KEGG analysis.
GO analysis of DEpiRNA target genes.
Primers and for RT-qPCR.
Length range and first base bias of piRNA in the gonad of juvenile Takifugu flavidus.
About this article
Cite this article
He, X., Wu, H., Ye, Y. et al. Transcriptome analysis revealed gene expression feminization of testis after exogenous tetrodotoxin administration in pufferfish Takifugu flavidus. BMC Genomics 23, 553 (2022). https://doi.org/10.1186/s12864-022-08787-z
- Tetrodotoxin (TTX)
- Piwi-interacting RNAs
- Sexual dimorphism
- Takifugu flavidus