The Impact of Exogenous Tetrodotoxin on the piRNAs and mRNAs Proles of Gonads in Puffersh Takifugu Flavidus

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 puffersh. The molecular mechanism underlying sexual dimorphism of TTX accumulation in gonads is complex and unclear. Piwi/piRNA complexes are essential for germline specication, gametogenesis, and gonadal development, they also demonstrate sexual dimorphism in teleosts. Hence, the present study investigated the expression of piRNAs and mRNAs by transcriptomics in cultured puffersh Takifugu avidus after intramuscular administration of exogenous TTX. The results showed 80 piRNAs were down-regulated and 223 genes were up-regulated in the ovary after TTX administration. By contrast, 286 piRNAs were down-regulated and 445 genes were up-regulated after TTX administration in testis. Functional and pathway analyses indicated that the TTX up-regulated genes were enriched in the Wnt, ErbB and GnRH signaling pathways in ovary, while were enriched in the oocyte meiosis, estrogenesis biosynthesis and cell apoptosis-related pathways in testis. Interestingly, these genes were also the potential target genes of TTX downregulated piRNAs. Amhr2 and cyp19a were also identied as sex-related genes, involved in TTX administration. These results showed a certain level of consistency with the enrichment pathway analysis, which indicated that the TTX could affect the expression of sex-related genes, and work as an inhibitor for the testicular meiosis, but as a promoting factor for ovary development through piRNAs in puffersh. In addition, TUNEL staining showed that signicant apoptosis was detected in the TTX treated testis, and the role of the cell apoptotic pathway was further conrmed. Overall, this research will contribute to an understanding of the piwi function in TTX sexual dimorphism accumulation in puffersh and provide a basis for further studies. (amhr2) pairing MiRnada. target


Introduction
Tetrodotoxin (TTX), one of the most potent neurotoxins which is speci c to voltage-gated sodium channels (VGSCs) on excitable membranes of muscle and nerve tissues, was long believed to occur exclusively in TTX-bearing puffer sh, mainly the genus Takifugu [1,2]. TTX is produced primarily by marine bacteria, and accumulated extremely high in puffer sh through the food chain with tissue-speci c and seasonal change [3][4][5]. It is suggested that TTX may function as a chemical defense against predators and as pheromone during spawing [6,7].
The accumulation of TTX in puffer sh has been studied extensively [1,2,[8][9][10][11]. Interestingly, the ovaries of puffer sh accumulate large amounts of tetrodotoxin, whereas the testes accumulate few or none, especially the ovarian toxicity was extremely high during the maturation period, indicating that the accumulated TTX showed sex difference in the toxicity of tissues [3][4][5]. Previous studies had found that tissue-speci c distribution and the amount of TTX in the mature puffer sh Takifugu niphobles were sexdependent, female gonads and male liver showed the highest concentrations of the toxin followed by male skin [12]. The larvae of Takifugu puffer sh were protected by maternal TTX which had accumulated in the eggs during their development in the ovary [13,14]. In addition, studies have shown that TTX as a toxin or pheromone accumulated in the ovary protects puffer sh larvae from predators [6, 13,14].
Although TTX was accumulation demonstrates sexual dimorphism gonads of puffer sh, the molecular mechanism underlying and the relationship between TTX accumulation and steroidogenesis/reproduction is complex and unclear.
To gain insights into sexual dimorphism of TTX accumulation in gonads, steroidogenesis/reproductionrelated genes and small non-coding RNAs were chosen to further investigate the molecular mechanisms. Piwi-interacting RNAs (piRNAs) are single-stranded, 26-32 nt long small RNAs that play a role function via the formation of RNA-protein complexes through interactions with piwi proteins [15,16]. In animals, previous studies revealed that piRNAs are limited expressed in a few tissues, such as tumors [17][18][19], brain and nervous tissue[18, [20][21][22], whereas generally abundant in gonads of invertebrates [23][24][25][26][27] and vertebrates [18,[28][29][30][31][32][33], including zebra sh (Danio rerio) [25,34], tilapia (Oreochromis niloticus) [35], Japanese ounder (Paralichthys olivaceus) [36], tongue soles (Cynoglossus semilaevis) [37] and other teleosts. Previous studies showed that maternal Piwi /piRNA is required for the fertility and normal gonad morphology of female, but not male, progeny in Drosophila [38]. piRNAs derived from the W chromosome are expressed more abundantly in the ovary than in the testis of silk worm [39]. The piRNAs play a role in targeting the transcripts of active TEs and maintaining the methylation pattern of maternal genomic DNA, and proving that ovaries and testes piRNAs are two different classes of piRNAs with different functions and that their expression appears to be regulated by distinct transcription factors in human [40]. Above all, at least a large number of piRNA show sexual dimorphism in the testis and ovaries. In addition, putative piRNAs were found abundant in the gonads of Takifugu rubripes [41], but the function and connection between tetrodotoxin and piRNAs in the gonads of puffer sh have not been investigated.
In the present study, we investigated the expression changes of piRNAs and mRNAs in both ovary and testis of Takifugu avidus on the transcriptomic level after the intramuscular injection with TTX. The up and down-regulated piRNAs and genes were identi ed, function and pathway of these piRNAs and genes were also analyzed. At last, apoptosis detection was performed in the TTX treated gonads. These data would provide us some clues to understand the speci c effects of TTX accumulated and function in the gonads of puffer sh.

Results
Exogenous TTX accumulated in both ovary and testis of arti cially breeding puffer sh A TTX treatment was performed in 6 months old of Takifugu avidus, after TTX administration, the TTX concentration of the TTX treated group is signi cantly higher than the control group in kidney, cholecyst, skin, liver, heart and muscle tissue by LC-MS/MS (Liquid Chromatography with tandem mass spectrometry) methods ( Figure 1A).
To further research the TTX accumulation effect in gonads of Takifugu avidus, histological analysis was performed to identify the sex and gonadal features in 6 months old of Takifugu avidus. Twelve gonads were dissected from each of the control and TTX treated groups. Five ovaries and seven testes were con rmed in the control group, six ovaries and six testes in the TTX-treated group through the gonadal sections histological observation. All of these gonads are immature, the ovarian cavity is obvious, oogonia and oocytes were observed in the ovary ( Figure 1B-a, b). In the testis, the spermatogonia were observed but without spermatocytes or spermatids ( Figure 1B-c, d). Immunohistochemical observations with anti-TTX antibody suggested that TTX was localized at the somatic cells around the oocytes in the control ovary with intense signals ( Figure 1B-e, f). However, just a few TTX signals were detected in the control testis ( Figure 1B-g, h). After TTX administration, the total TTX signals in the TTXtreated gonads were signi cantly higher than the control gonads ( Figure 1B-i, j, k, l &1C). piRNA and mRNA sequencing and novel piRNA identi cation Total RNA of gonads (testis or ovary) was isolated and sequenced for piRNA and mRNA analysis from the TTX-treated and control groups of Takifugu avidus. Overall, 94.91% of con_O (control ovary group), 94.91% of con_T (control testis group), 98.09% of TTX_O (TTX-treated ovary group) and 94.32% of TTX_T (TTX-treated testis group) of the clean reads had scored at the Q30 level in mRNA sequencing (Table 1). Of all the clean reads, 85.06% (con_O), 86.10% (con_T), 86.04% (TTX_O) and 85.63% (TTX_T) were matched onto the fugu reference genome sequence (fTakRub1.2).
For piRNA sequencing, the peak value of all groups was mainly concentrated at 28 bp (Figure s1-A) with the character of 1U-bias of the molecule (Figure s1-B), which corresponded with the characteristic length and structure bias of piRNAs. The numbers of unique piRNAs for the species were calculated as follows: 507,004 in the control ovary group, 530,323 came from the TTX-treated testis group; 346,106 came from the control testis group while 472,764 in the TTX-treated ovary group, the total unique piRNA is 1,247,865 (Table 1 and Additional le 1). In contrast, more than 288,864 unique piRNAs display less than 10 reads in all groups. These results indicate that expression varied signi cantly among different piRNA families. Multiple piRNAs with at least 10 piRNAs located on the same chromosome and no more than 5000 nt apart are de ned as a piRNA cluster. A total of 5,457 piRNA Clusters were isolated. These piRNA clusters exist on each chromosome and are unevenly distributed on each chromosome (Figure s2 The region consistent with the gene sequence may be the generation region of the piRNA clusters, and the gene sequence that is reversed complementary to the piRNA cluster may be the target site of the piRNA clusters. TTX affecting the piRNA/mRNAs pro les in the ovary In the piRNA and mRNA transcriptomes, the differential expression in the ovaries of the control group and the TTX treatment group were analyzed. TTX administration resulted in 23 piRNAs and 223 genes were up-regulated, meanwhile, 80 piRNAs and 126 genes were down-regulated in the ovaries (Figure 2 To investigate the reliability of differential expression data in the piRNA and mRNA transcriptomes. Eight 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, alcohol dehydrogenase 1-like and LOC101069468, low-density lipoprotein receptor-related protein 5-like) were screened for further veri cation by qRT-PCR. The piRNA and mRNA expression change between the TTX-treated group and control group in the qRT-PCR results are consistent with the transcriptome data indicating that the transcriptome data is reliable (Figure 2-C, D and Additional le 3).
Studies have been suggested that Piwi/piRNA complexes are important in gene silencing, cell differentiation and gonadal development regulation in animals. Similar to the method of miRNA and other small RNAs to predict target genes, we obtained differential piRNAs between TTX-treated and control ovary for target gene prediction and enrich the target genes with KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene ontology) analysis. KEGG and GO analysis were also used for the prediction target gene of TTX down-regulated piRNA in ovaries. As revealed in Figure 3 TTX up-regulated pathways in the ovary were compared with the piRNA transcriptome data, the oocyte meiosis, progesterone-mediated oocyte maturation, and steroid hormone biosynthesis in TTX-upregulation genes were selected for mapping with the TTX-down-regulation piRNA. The results showed that uniq_1958827 predicated target gene including calcium voltage-gated channel subunit alpha1(cacna1c); uniq_1013268 predicated target gene including alcohol dehydrogenase 1-like (LOC101078894); uniq_479370 predicated target gene including low-density lipoprotein receptor-related protein 5 (LOC101069468); uniq_1713020 predicated target gene including peroxisome proliferatoractivated receptor delta (ppard). uniq_671648 predicated target gene including anti-Mullerian hormone receptor type 2 (amhr2). The gene cacna1c was involved in the GnRH signaling pathway, the gene LOC101078894 was involved in Retinol metabolism; both LOC101069468 and ppard were involved in the Wnt signaling pathway; The gene amhr2 was involved in the TGFβ signaling pathway ( TTX affecting the piRNA/mRNAs pro les in testis In the piRNA and RNA transcriptome, the differential expression of the control group and the TTX treatment group testes were analyzed. In the testis, TTX up-regulated the expression of 224 piRNAs and down-regulated the expression of 286 piRNAs ( Figure 4-A and Additional le 5). Only the ddah2 (dimethylarginine dimethylaminohydrolase 2) and ing4 (inhibitor of growth family member 4) in the testis were signi cantly down-regulated. However, 445 genes were up-regulated ( Figure 4-B and Additional le 5).
To investigate the reliability of differential expression data in the piRNA and RNA transcriptomes, nine piRNAs (uniq_2692448, uniq_904526, uniq_1754588, uniq_1320158, uniq_1010678, uniq_1323470, uniq_2869791, uniq_1678838 and uniq_766456) in the testis were selected for further veri cation by qRT-PCR. For the differential expressed genes, the hsd17b1(hydroxysteroid 17-beta dehydrogenase 1), ccnb1(cyclin B1), aurka (aurora kinase A), fbxo43 (F-box protein 43), LOC101065671(cyclin N-terminal domain-containing protein 2-like), pygl(glycogen phosphorylase L), LOC101080011(growth arrest and DNA damage-inducible protein GADD45 beta) and the TTX down-regulated two gene ing4 (inhibitor of growth family member 4) and ddah2 (dimethylarginine dimethylaminohydrolase 2) were selected. The qRT-qPCR results con rmed the good correlation between the piRNA expression level and piRNA sequencing. Besides, the expression changes between the TTX-treated group and control group in the qPCR results are consistent with the transcriptome data ( Figure 4-C, D and Additional le 3), showed the piRNA/mRNA transcriptome data were reliable.
KEGG and GO analysis were used for prediction potential target genes of TTX down-regulated piRNAs in testis. Results showed as in Figure 5-A and Additional le 4. Focal adhesion, Wnt signaling pathway, Calcium signaling pathway, Apoptosis, Autophagy, GnRH signaling pathway, Thyroid hormone signaling pathway, Estrogen signaling pathway and oocyte meiosis were enriched. KEGG and GO analysis showed that the enrichment of TTX up-regulated gene in the testis are mainly in cell apoptosis-related pathways, such as cell cycle, tight junction, PPAR signaling pathway, p53 signaling pathway, Ferroptosis and Necroptosis. Besides, oocyte reproduction-related pathways were also enriched, such as oocyte meiosis, progesterone-mediated oocyte maturation and Steroid hormone biosynthesis ( Figure 5-B and Additional le 4). Compare the differential expression piRNAs potential target gene enrichment with the differential expression gene enrichment. Both piRNA and mRNA showed that after TTX injection, the gene expression changes in the testis develop toward the "ovarian" pattern and rises the apoptosis pathway (p53, Ferroptosis and Necroptosis).
Exogenous tetrodotoxin induced the cell apoptosis of testis As we have known that puffer sh testis usually accumulates few tetrodotoxins, the excessive exogenous tetrodotoxin may be harmful to testis. Besides, the piRNA and mRNA transcriptomes data showed that cell apoptosis-related genes were up-regulated in testis and related piRNA were down-regulated. These results indicated that abundant exogenous tetrodotoxin will induce the apoptosis of testis. TUNEL assay (Terminal deoxynucleotidyl Transferase Biotin-dUTP Nick End Labeling) was used for identifying apoptotic cells in situ. As shown in Figure 6-A, most nuclei were stained as a discernible brown in the TTX-treated testis but few signals were found in the control testis. The TTX-treated and control ovary were both found the signals at the somatic cells, very few signals were found in the oocytes. In addition, the area of the positive signals was statistics and analysis, the results showed that no signi cant difference of cell apoptosis between the TTX-treated and control ovary, but the signals of cell apoptosis in TTX-treated testis is signi cantly higher than that in the control testis ( Figure 6-B), and the cell apoptosis induced by exogenous tetrodotoxin administration.

Discussion
In this study, the exogenous tetrodotoxin was administration to the arti cially cultured Takifugu avids. TTX level elevation was observed at the testis and other tissues after eight hours. The piRNA and transcriptome sequence results showed that TTX down-regulated piRNAs and up-regulated genes in testis were biased toward the ovarian expression pattern, it is re ected by piRNAs and genes enriched in Wnt signaling pathway, Estrogen biosynthesis signaling pathway and oocyte meiosis pathway. Additionally, the cell apoptosis pathway was also enriched and con rmed by the TUNEL assay in testis.
To our knowledge, the present study is rst time investigated the piRNA and transcriptome expression changes in Takifugu avids' separate sex. As the ovary can accumulate abundant TTX while the testis has just been detected little TTX in the control group, we speculate a high level of tetrodotoxin will hamper the development of testis. Processes involving sex determination and reproduction were complicated and regulated by various internal and/or external factors, particularly in puffer sh which show sexual dimorphism of TTX accumulation in gonads. With limited molecular data available for Takifugu to know about the exact internal mechanism or role of tetrodotoxin in the development of gonads. The accumulation and tissue-speci c distribution of tetrodotoxin in puffer sh, mainly in the genus Takifugu has been widely investigated from the viewpoint of the TTX-resistance VGSCs expression and distribution [2,10,42]. Evidence also showed that tetrodotoxin is usually signi cantly at a high level during the maturation/spawning period than in the ordinary period, exogenous tetrodotoxin administration also showed consistent results [1,4,5,9,14]. But the tetrodotoxin in the developing stage of puffer sh, especially for the gonads marginally been investigated. The present study was performed to elucidate the possible TTX impact on the ovary and testis of juvenile puffer sh development stage through piRNA and mRNA transcriptome.
The piRNA and mRNA pro le change of both ovary and testis were analyzed from the view of transcriptomics, total of 1,247,865 unique piRNAs was identi ed from four groups, but unique piRNAs from every group is different. The most interesting is the total number of uniq piRNA is increase after exogenous TTX administration in testis, but decrease after exogenous TTX administration in the ovary. it is also showed that the number of DEpiRNA (differentially expressed piRNA) which is affected by TTX in the testis (510) is more than that in the ovary (103), this indicated TTX has a bigger in uence on the testis. The increased piRNAs in testis silence more transposons or genes to resist/reduce the toxicity of TTX. In addition, TUNEL analysis showed that exogenous tetrodotoxin induced the cell apoptosis of testis, but there is no signi cant effect on the apoptosis of the ovary.
Despite fewer genes were down-regulated by the TTX in the testis, many genes involved in the oocyte meiosis, reproduction and apoptosis-related gene/pathway were remarkably up-regulated in the testis, showed the "ovarian" developmental tendency and cell apoptosis tendency. These results of KEGG pathway analysis in differential expression gene/piRNA in TTX treated testis revealed gametogenesis and reproduction-related pathways, including the Wnt signaling pathway, MAPK signaling pathway, oocyte meiosis pathway. Especially, the up-regulated oocyte meiosis pathway was found in the TTXtreated testis. The B-type cyclins take important roles in oocyte meiosis [43]. Cyclin B1(ccnb1) and Cyclin B2 like (ccnb2L) are essential for the is critically required for the proliferation of gonocytes [44], both of them were upregulated by the TTX. Combined with the predication of differential expression piRNA target gene, piRNA uniq_1504863, uniq_1302131, uniq_2344671 may target them. Another interesting upregulation in the TTX-treated testis pathway was steroid hormone biosynthesis. Steroidogenesis is critical for gamete maturation in teleosts, more essentially oogenesis [45], the cyp19a is the key gene encoding cytochrome P450 aromatase, which is responsible for estrogen production in the ovary [46-48], differential expression piRNA uniq_554482 predicted target also including cyp19a. In addition, related steroidogenesis gene (nr5a1), growth factor (LOC101078012 or gsdf), germ cell-related marker (piwil1; nanos2 or LOC101063734) and meiosis related gene (scyp3, LOC101063989 or cyp26a1) were nonsigni cant down-regulated by TTX. TTX also up-regulated the apoptosis-related pathway and genes in the testis, the p53 signaling pathway, Necroptosis and Ferroptosis pathway. The p53 signaling pathway has been widely accepted that it is involved in the regulation of cell cycle arrest, senescence and apoptosis [49][50][51]. Both Necroptosis and Ferroptosis were pathways discovered form of regulated cell death [51][52][53]. Interestingly, these down-regulated piRNA-predicted target mRNAs were involved in the upregulated genes in the testis, showed the impact of TTX on gene expression may through the role of piRNA in the testis.
Studies have shown that the components of the Wnt signaling pathway affect reproductive functions, including the regulation of follicular maturation and the production of steroids in the ovaries. Both GnRH signaling pathways and Retinol metabolism are related to meiosis, indicating that a certain increase TTX in the ovary contributes to ovarian meiosis and development. Similarly, the Wnt signaling pathway, GnRH signaling pathway, TGFβ signaling pathway and Retinol metabolism related genes were up-regulated after TTX administration. The Wnt and GnRH signaling pathways are required normal fertility and reproductive cycles [54][55][56][57]. Above all, TTX accumulation in gonads of puffer sh showed sexual dimorphism. TTX may work as a promoting factor to puffer sh ovary development and ovarian meiosis but as an inhibitor to puffer sh testis development, too much TTX will induce the apoptosis of testis. The metabolism of TTX in males and the molecular mechanism of testicular resistance to TTX accumulation still require further investigation.

Conclusion
In conclusion, the present research highlights the differentially expressed piRNAs and genes after TTX administration in gonads. 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 tetrodotoxin novel function in puffer sh reproduction. A certain level of consistency with the enrichment pathway analysis, which indicated that the TTX could affect the expression of sex-related genes, and work as an inhibitor for the testicular meiosis, but as a promoting factor for ovary development through piRNAs in puffer sh. TTX administration also induces signi cant apoptosis in testis, further con rmed that excessive TTX is harmful for testicular development. Overall, this research will contribute to an understanding of the piwi function in TTX sexual dimorphism accumulation in puffer sh and provide a basis for further studies.

Animals and chemical
Cultured marine puffer sh Takifugu. avidus Juveniles (6-month-old, 100.0 ± 10.0 g body weight, 100 individuals) purchased from Shanghai Aquatics Research Institution were used as the experimental sh. Crystalline TTX from (Supelco, USA) was used in the administration experiment and as a standard for the liquid-chromatography-tandem mass spectrometry (LC-MS/MS) analysis. All other chemicals were of reagent grade. All other chemicals were of reagent grade (Sinopharm, China). The methods for Takifugu. avidus research was carried out following the relevant guidelines and regulations (following the ARRIVE guidelines[69]). The protocols were approved by the academic ethics committee of Shanghai Ocean University.

TTX treatment and sample preparation
The administration experiment was performed at the laboratory of the shanghai ocean university. The tetrodotoxin was directly dissolved in 0.1% acetic acid to 1 mg/ml, then diluted with 0.7% saline solution to 0.25 ug/ul. The test sh received an intramuscular administration of 0.25 ug TTX/g body weight into the caudal muscle as described previously and were maintained in an 80-L oxygenated plastic tank of aerated arti cial seawater at 20 ℃ [70]. At 8 hours after administration, each group of 12 sh were randomly collected, ice-cold anesthetized, and dissected. The gonads of each were immediately separated into two parts, 1/4 gonad was stored in 4% PFA for histological observation; 3/4 gonads were stored at -80 ℃ until RNA extraction. TTX extraction and quanti cation TTX in the tissue samples from puffer sh was extracted with 0.1% acetic acid as reported previously [71,72]. The TTX assay from tissues was determined by the LC-MS/MS analysis according to the method used previously [71,73,74].

Histological observation, TTX immuno uorescence
For histological observation and immuno uorescence analyses, the gonads of Takifugu avidus were dissected, xed in 4% PFA (Paraformaldehyde) for 12 hours at 4℃, dehydrated, embedded in para n, and sectioned at 5 µm. Hematoxylin-eosin staining and immuno uorescence staining were performed as described previously [75,76]. The antibody against TTX (tetrodotoxin) was purchase from CD (Creative diagnostics, USA) and was diluted at 1:200 to use, Hoechst 33258 (Sigma, USA; 1:1000) was used to stain the nucleic acid. The TTX signals ratio of gonads was statistics the three different areas of one sample under TTX immuno uorescence by using ImageJ (NIH, USA) and Graph Prism8.0 (San Diego, USA).

TUNEL analysis
The sections of gonads were used for apoptosis analysis by TUNEL staining, the terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) assay was conducted by a colorimetric TUNEL Apoptosis Assay Kit (Beyotime, China) according to the manufacturer's instructions. Brie y, The sections of gonads were depara nized twice in dimethyl benzene and then graded in a series of ethanol, after repair with proteinase K (2 mg/mL) and then moved the slides into the 0.3% H 2 O 2 in methanol for 20 min at RT (room temperature). Washed 3 times with 1×PBS (Phosphate Buffered Saline). That was followed by incubation with TdT enzyme solution for 60 min at 37°C. The reaction was terminated by incubation in a stop/wash buffer for 10 min at 37°C. Then used the Streptavidin-HRP (Streptavidin-Horseradish Peroxidase) for coloration. Incubated the sections in the Streptavidin-HRP working solution for 30 min at RT, following with the color reagent incubation for 3 min, and then restaining by hematoxylin for 3 min, hydrated and nally sealed and documented as digital images on a Nikon Eclipse 80i microscope (Nikon, Japan). The apoptosis signal ratio of gonads was statistics the three different areas of one sample under TUNEL staining by using ImageJ (NIH, USA) and Graph Prism8.0 (San Diego, USA) as described before.

piRNA/RNA library construction and sequencing
For piRNA and RNA transcriptome analysis, including library construction, sequencing, and bioinformatic analysis were performed by Novegene Company (Beijing, China). Brie y, total RNAs from every group (four sh gonads mixed to one sample) were extracted using RNAiso Plus (Takara, Japan). The quality of RNA was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and agarose gel electrophoresis.
For piRNA, Small RNA libraries were constructed and sequenced using TruSeq Small RNA Donor Prep Kits (Illumina, USA) by the Novegene Company (Beijing, China), sequenced using the basic reads were converted into sequence data (also called raw data/reads) by base calling. Low-quality reads were ltered out, and the reads with 5′ primer contaminants and poly (A) regions were removed. Reads without a 3′ adapter and insert tag, reads shorter than 24 nt, and longer than 32 nt were ltered out of the raw data to obtain clean reads.
For piRNA primary analysis, the length distribution of the clean sequences in the reference genome was determined (closely related species Takifugu rubripes, fTakRub1.2, PRJNA543527). Non-coding RNAs of 26-32 nt were annotated as piRNAs. The known piRNAs were identi ed by aligning against the piRNABank (http://pirnabank.ibab.ac.in), and the remaining unannotated reads were characterized by piRNA predictor [77]. The novel piRNAs reads were analyzed using Piano to predict novel piRNAs (http://ento.njau.edu.cn/Piano.html).
Previous studies showed that by using the structure and sequence characteristics of transposon-piRNA interactions, Piano demonstrated excellent predictive performance for piRNAs. Differentially expressed piRNAs were identi ed with the threshold of |log2 fold change|>1 and pvalue<0.05. The p-value was calculated using the DEG algorithm [78]. In the R package with the Audic-Claverie statistic without biological replicates. The targets of differentially expressed piRNAs were predicted using the MiRanda software [79]for animals, with the following parameters: S ≥ 150; ΔG ≤ −30 kcal/mol and strict 5′ seed pairing. Differentially expressed piRNA-target genes were classi ed according to the o cial annotation of the Gene Ontology (GO) Consortium, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was performed using phyper, a function of R [80,81]. False discovery rate (FDR) was used to determine the threshold of the P-value and GO or KEGG terms (FDR≤0.01) were considered signi cantly enriched. Validation of DEG results of piRNA/mRNA Gonadal piRNA and mRNA expression were assayed using real qRT-PCR. For mRNA, total RNA was extracted using RNAiso Plus (TaKaRa, Japan) following the protocol provided by the manufacturer. 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 Additional le 6. gadph was used as an internal control [83,84].
For piRNA, we validated the differential expression of four randomly-selected piRNAs (uniq_2692448, uniq_1288101, uniq_1694274 and uniq_1744355), four RNA-seq DEG related pathway piRNA (uniq_1958827, uniq_1013268, uniq_479370, uniq_1713020) in the ovary and four randomly-selected piRNAs (uniq_2692448, uniq_904526, uniq_1754588 and uniq_1320158) ve RNA-seq DEG related pathway piRNA (uniq_1010678, uniq_1678838, uniq_1323470, uniq_766456 and uniq_2869791) in the testis and with qRT-PCR. we used 1µg total RNA above as a template for the synthesis of the rst-strand cDNA with a Mir-X miRNA First-Strand Synthesis (Takara, Japan) following the manufacturer's instructions. we measured the expression of the above piRNAs on a CFX96 Real-Time PCR System (Bio-rad, USA) with a miScript SYBR Green PCR Kit (Takara, Japan), following the manufacturer's instructions, the primers for piRNA qPCR were listed in Additonal le 6. We used U6 as an internal reference gene to control for differences among samples. Each sample was run in triple times and the relative quanti cation of piRNA expression was calculated using the 2 −ΔΔCt method [ respectively; i, j, the ovary of the treated group; k, l, the testis of the treated group. OG, oogonia; OC, the oocyte; SG, spermatogonia; The green signals are TTX and the blue signals are Hoechst33258 (staining for nucleic acid), the white dotted-rectangle indicate the magni ed area, black arrows showed the germ cells, the white triangle arrows are indicated the TTX signals, the scale bar is 100, 200 and 50 μm. C, TTX signals of gonads in juvenile Takifugu avidus after TTX administration. Data were shown as mean±SEM. *, indicated that ρ<0.05; **, indicated that ρ<0.001; ***, indicated that ρ<0.0001.

Figure 2
The Differential expression of piRNAs and mRNA in the TTX-treated ovary and control ovary. In piRNA/mRNA sequencing, the red dot represents up-regulated gene/piRNA, the green dot is downregulated gene/piRNA, the gray dot is no signi cant difference gene/piRNA, DEG analysis threshold isρvalue < 0.05 && | log2 fold change | > 1; Expression levels of piRNA/mRNA measured with qRT-PCR were depicted with column charts, and the red line represents the log2 fold changes/fold changes in sequencing. Error bars represent one standard deviation of three different biological replicates. the internal reference gene is U6 for piRNA and Gapdh for mRNA. set as **ρ < 0.01, ***ρ < 0.001 and ****ρ < 0.0001 vs controls. Statical signi cance is reported for each piRNA or RNA. DEpiRNA, differential expression piRNA; DEG, differential expression mRNA. A Volcano diagram of differential piRNA analysis in the ovary after TTX treatment; B, Differential expression of genes in the ovary of TTX treatment group and control group; C, Relationship between relative expression levels of selected four piRNA in the TTXtreated ovary and control ovary validated by qRT-PCR and log2fold changes derived using piRNA sequencing, control ovary as a calibrator; D, RT-qPCR and mRNA sequencing of selected three genes between ttx treated ovary and control, control ovary as a calibrator.  The Differential expression of piRNAs and mRNA in the TTX-treated and control testis. Expression levels of piRNA or RNA measured with qRT-PCR were depicted with column charts, and the dotted line