Identiﬁcation and application of piwi-interacting RNAs from seminal plasma exosome in Cynoglossus semilaevis

Piwi-interacting RNAs (piRNAs) have been linked to epigenetic and post-transcriptional gene silencing of retrotransposons in germ line cells, particularly in spermatogenesis. Moreover, exosomic piRNAs are promising biomarkers for disease diagnosis and physiological status indication. To determine whether piRNAs from fish germ line cells have similar features, seminal plasma exosomes from half-smooth tongue sole, Cynoglossus semilaevis, were identified, and their small RNAs were sequenced and analysed. We used C. semilaevis because of its commercial value and its sexual dimorphism, particularly the sex reversed ''pseudo-males'' who have a female karyotype, produce sperm, and copulate with normal females to produce viable offspring. We identiﬁed six signature piRNAs as biomarkers in seminal plasma exosomes from males and pseudo-male C. semilaevis. Bioinformatic analysis showed that all six signatures were sex-related, and four were DNA methylation-related and transposition-related piRNAs. Their expression proﬁles were veriﬁed using real-time quantitative PCR. The expression of the signature piRNAs was markedly higher in males than in pseudo-males. The signature piRNAs could be exploited as male-specific biomarkers in this fish. Exosomes are important mediators of vesicle transport, and the piRNAs in exosomes might play an important role in cell communication and signal pathway regulation. These signatures provide an effective tool to explore the regulatory mechanism of sex development in C. semilaevis and may provide guidance for future research on the function of piRNAs in the generative mechanism of sex reversed ''pseudo-males'' in C. semilaevis.


Background
Half smooth tongue soleCynoglossus semilaevis) is, a commercially valuable flatfish that is widely distributed in Chinese coastal waters, and is commonly found in shallow waters on a muddy or sandy bottom 1 . Many previous studies have shown that C. semilaevis employs a heterogametic sex determination system (ZW/ZZ) 2,3 and has significant sexual dimorphism, with a larger female body size and faster growth rate 4,5 . This species is attracting more attention in rReproductive and sexrelated research are attracting more attention to this species, and it could become a tool to study sex determination in fish 6 . Furthermore, this species also exists as pseudo-male fish, both in nature and aquaculture 7 . The high proportion of males in populations of C. semilaevis was partly attributed to the considerable number of pseudo-males. The pseudo-male has the same karyotype as the female fish, butfish but has the physiological characteristics of males 8 . Interestingly, pseudo-male fish are fertile and can pass on its their pseudo-male HYPERLINK "https://fanyi.baidu.com/" \l "en/zh/characteristic" \t "_blank" characteristiccharacteristic to its their offspring. When pseudo-male fish are used as parents, an imbalance in the proportion between the female and male tongue soles will arise 9 . Distinguishing pseudo-male from male fish and inhibiting them from mating with females could maintain the sex balance of sex proportion in C. semilaevis populations, which has great commercial value in aquaculture 10 . AlsoIn addition, it is of great significanceimportant to explore the influencing factors and determining mechanism of pseudo-male occurrence for to enriching the theory ofobtain further details on the sex determination mechanism of fish.

5
Piwi interacting RNAs (piRNAs) are single-stranded, 25-to 33 -nt-long small RNAs that play a function via forming RNA-protein complexes through interactions with piwi proteins 11 . PiRNAs are distinct from microRNAs (miRNAs) in terms of their size (26-31 nt rather than 21-24 nt), lack of sequence conservation, and increased complexity 12 . Previous profiling studies showed that miRNAs are widely expressed in different tissues, while piRNAs are abundant in gametes 13,14 . PiRNAs have been found in the testes and ovaries in mammals 15 , and were detected in both male and female germlines 16,17 . PiRNAs play roles in spermatogenesis in Caenorhabditis elegans, mice, and humans;. As we know,The piRNA pathway relies on the specificity provided by the piRNAs to identify transposon element (TEs TE) targets, while the effector function is provided by the piwi protein. Different piRNAs recognize different TEs target gene sequences to play different regulatory roles 18 . This leads us to think aboutwhether theWe hypothesized that piRNAs are might be differentially expressed in germlines between males and pseudo-males in Cynoglossus semilaevis, whether and piRNA might plays a role in the cross generation inheritance of pseudo-males. If so , how are piRNAs, as a crossgenerational sex-related regulatory factor of cross generation, how tos, ensure its stabley transmissiontransmitted.? So Consequently, we paid our attentions oninvestigated piRNAs in exosomes from the seminal fluid in on C. semilaevis.
Previously, sSmall RNAs, including piRNAs in exosomes from body fluid, have been reported as biomarkers 19-20 previously. Exosomes from mouse and human seminal fluid have been isolated and analysed, representing an excellent system to study piwi-interacting genes and their regulatory network 21 . Although several different kinds of sex molecular markers have been developed in C. semilaevis, such as 6 amplified fragment-length polymorphism (AFLP)AFLP markers 3 , co-dominant microsatellite markers 22 , single nucleotide polymorphisms (SNPs) 23 , and miRNAs markers 24 , identifying more different types of sex markers is also of some significanceimportant for to studying the sex determination mechanism of this fish.
In this the present studywork, we compared the piRNA profiles in exosomes from males and pseudo-males seminal plasma to identify the differentially expressed piRNAs. PiRNAs with significant differential expression in between males and pseudo-males were selected from candidate piRNAs for sex identification and their expression profiles were verified using quantitative real-time reverse transcription PCR (qRT-PCR).The signatures could be used as biomarkers to distinguish pseudomale from male fish in C. semilaevis.,, also,In addition, the interaction and regulatory mechanism between piRNAs and target genes would play an important role in explaining the cross generational genetic mechanism of pseudo-male Cynoglossus semilaevis..

Identification and characterization of Exosomes.
We collected about 20 ml of seminal plasma from 60 male and 60 pseudo-male donors, respectively (Fig. 1). Exosomes were isolated and purified using an SBI ExoQuick-TC kit after the samples were filtered through a 0.45-µm membrane. We used transmission electron microscopy (TEM) and nanoparticle tracking analysis (NTA) to identify exosomes and capture images and data for granularometric analysis of the exosome preparations ( Fig. 2A). The exosome particles had diameters ranging from 30 to 150 nm, which was consistent with the characteristic size range (30-120 nm) of exosomes, and represented almost 70% (67.043% and 69.693%) of all the particles in both samples, with mean sizes of 108.2 ± 0.57 nm and 116.3± 0.32 nm in males and pseudo-males, respectively. The main peaks of particle size in the NTA analysis were 47 and 48 nm, respectively. The concentrations of two samples were 24.9 ± 0.94×10 8 and 23.9 ± 2.64 ×10 8 particles/ml, respectively. ( Fig. 2B and C). All these measurements suggested that the exosome preparations isolated from male and pseudo-male C. semilaevis seminal plasma contain a heterogeneous mixture of exosomes and microvesicles, which was similar to that in previous reports 24 .
We also investigated the presence of three tetraspanins as exosome markers using western blotting, including CD63, CD9, and heat shock protein 90 (HSP90)to confirm the existence of exosomes. Immunoreactive bands corresponding to CD63 and HSP90 were observed, whereas CD9 had no obvious immunoreactive band ( Small RNA sequencing and the nucleotide composition in the exosomes. RNA was isolated from the exosome preparations from male (ZZ♂) and pseudo-male (ZW♂) C. semilaevis and sequenced for small RNA analysis. The reads that aligned reads to the genome of half-smooth tongue sole were employed to determine the length distribution of the two groups: we found that the peak values of both groups were mainly concentrated at 31 bp, which corresponded with the characteristic length of piRNAs (supplementary Supplementary figures Figure 4). The numbers of different types of mature piRNAs for the species were calculated as follows: the number of unique known piRNA aligned reads was 56484, representing about 22.71% of all clean reads in the pseudo-males donor group, while 55324 (26.6%) came from male donors. We constructed pie charts for the classification and annotation of the small RNA reads of each donor group (Fig. 3). In both donor groups, the known piRNAs occupied about 25% of the unique reads. A piRNA annotation program (Piano) was developed using piRNAtransposon interaction information. The piRNA-transposon interaction was predicted using RNAplex. We employed the unaligned sequences filtered from piRBase to carry out novel piRNA prediction. The predicted novel piRNAs are were between 21 and 38 bp and can could be mapped to the genome. In total, 14006 non non-repetitive novel piRNAs were predicted from both pseudo-males (ZW♂) and male (ZZ♂) donors Identification of Signature piRNAs between male and pseudo-male C. semilaevis.
The differential expression profiles of piRNA were investigated between male and pseudo-male C. semilaevisusing the TPM (transcript per million) values. In total, 26135 differentially expressed piRNAs were identified according to the criteria detailed in the methods section (supplementary datasets4dataset 4). Among these piRNAs, 15373 were upregulated and 10762 were downregulated in ZZ♂ compared with ZW♂. Using further screening conditions, we filtered out the novel piRNAs (6622) because of its their unproven existence and narrowed the highly expressed 9 piRNAs down to 87 known piRNAs under the condition of at least a TPM of one group ≥ 150 and a fold-change ZZ♂ / ZW♂ ≥ 100. Then, considering that 87 was too many piRNAs for subsequent analysis,we determined the final filtration conditions bywe adjusting adjusted the fold-change (ZZ♂ / ZW♂) to ≥ 200 and at least a TPM of one group ≥ 400, which identified narrowed the dataset 44 candidate piRNAs ( Supplementary Fig 4) . We predicted the target genes of the 44 candidate signature signature piRNAs that were differentially expressed in males and pseudomales. In the present study, "signature" meant a piRNA marker with significant differential expression as verified by qRT PCR. After piRNA-targetst prediction, we obtained 12145 piRNA-target pairs and 6231 target genes (supplementary datasets 5).
Target prediction for candidate signature piRNAs and GO enrichment and KEGG pathway enrichment analysis.
The 6231 target genes were employed for gene ontology (GO) enrichment and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway enrichment analysis. GO term categories generated from the 6231 genes targeted by the 44 candidate signature piRNAs showed that most of the target genes are involved in plasma membrane, integral component of membrane, extracellular exosome, metal ion, and ATP binding, transcription, and DNA−templates in the cellular component and biological process categories ( Supplementary Fig. 5). Among the target genes of the 44 candidate signature piRNAs, those related to sex differentiation, sex determination, and sex development in subsystems of GO enrichment and KEGG pathway enrichment analysis were identified, which allowed us to reduce the number of candidate signature piRNAs to 37 ( supplementary datasets 6). The 37 candidate 10 signature piRNAs included: 17 piRNAs with high expression in the ZZ♂ group but little expression in ZW♂; and(2) 20 piRNAs with a non-zero expression in the ZZ♂ group, but much higher expression in the ZW♂ group. The expression profiles of these piRNAs are shown in Fig. 4 and supplementary datasets 7... The KEGG pathway enrichment analysis showed that lipidcarbohydrate metabolism and signal transduction were the top two functional categories of the target genes ( Supplementary Fig. ure 6). Meanwhile, we also investigated the target genes related to DNA methylation and transposition, because previous research showed that epigenetic regulation plays multiple crucial roles in the sex reversal of halfsmooth tongue sole, and piRNAs may be involved in transposon silencing. Fiftythree DNA methylation related target genes were identified that were predicted to interact with 15 regulating piRNAs (supplementary datasets 8). Meanwhile, 16 target genes related to transposition, especially heterochromatin formation, were also identified together with their 10 interacting piRNAs (supplementary datasets 9). We carried out Venn diagram analysis among the 37 sex-related, 15 methylationrelated, and 10 transposition-related piRNAs, which identified seven piRNAs in the intersection of all three sets. The result implied that piRNAs might regulate the sex development of C.semilaevis through epigenetic regulation or transposition (Fig. 5).
To investigate the candidate signature piRNAs identified in the present study, we chose 15 candidate signature piRNAs for further verification by qRT-PCR from among the 37 most differentially expressed candidate piRNAs. We used RNA from 10 male and 10 pseudo-male fish to carry out qRT-PCR to quantitatively measure the expression of marker piRNAs. The results of qRT-PCR showed that the expression of 11 six marker piRNAs in 10 male and 10 pseudo-male fish were significantly higher in males than pseudo-males ( Fig. 6 and supplementary datasets 10), which was consistent with the results obtained from the piRNA profiling in the small RNA sequencing analysis. Therefore, these six signature piRNAs (piR-mmu-29271668, piR-mmu-6643660, piR-xtr-979116, piR-mmu-32360528, piR-mmu-72274, and piRmmu-31018127) could be considered as male molecular biomarkers for C. semilaevis.

Discussion
The sex determining mechanisms of fish are complex and diverse.Research using model organisms has revealed that gender determination is influenced by many factors. In C. semilaevis, there are so called "pseudo-male" fish with male sex characteristics and a female karyotype, which generate fertile sperm and mate with normal females to produce viable offspring 25 . We chose half-smooth tongue sole as a model to characterize reproductive regulation differences at the subcellular and molecular level between male and pseudo-male fish. As a result, several signature biomarkers were developed based on small RNAs sequencing. We successfully isolated and captured exosomes derived from seminal plasma in C. semilaevis, from which we identified six piRNAs with significant differential expression for development as biomarkers to distinguish males from pseudo-males in sex identification. piRNAs from exosomes as biomarkers in fish. We identified six piRNAs with significant differential expression as biomarkers for males in sex identification.
Target genes prediction showed that there was a high coincidence between between piRNAsA targets related to sex development and DNA methylation. Earlier research indicated that the piRNA pathway relies on the specificity provided by the piRNA sequence to identify complementary transposable element (TE) targets, while the effector function is provided by the PIWI protein. PIWI silences TE transcription at 13 the chromatin level by directing inhibited histone marker deposition and DNA methylation to the TE copy 31,32 . Whether the sex regulation mechanism of half smooth tongue sole depends on epigenetic regulation or DNA methylation through transcriptional silencing by piRNAs requires further study.
Our signature piRNAs: piR-mmu-6643660, piR-mmu-32360528, piR-mmu-72274, piR-mmu-31018127, piR-mmu-29271668, and piR-xtr-979116, were all highly expressed in male C. semilaevis donors, but showed very low expression in pseudomale fish. The first four signatures were are related to both sex development and methylation through target gene prediction, while the last two are related to sex development only. Currently, the functions of the six piRNAs have are unclear functions in all species, including fish, which illustratesindicating that further studies should be focus on the regulation regulatory mechanism of these piRNA during their interaction with their targets.

Conclusion
In this study, exosomes derived from C. semilaevis seminal plasma were successfully isolated. Small RNAs were sequenced and used to identify signature piRNAs as sexual biomarkers to distinguish male and pseudo-male C. semilaevis. We identified 44 candidate signature piRNAs with extremely significant differential expression profiles. Target genes were predicted and then subjected to GO enrichment and KEGG pathway enrichment analysis. Furthermore, seven piRNAs appeared at the intersection of 37 sex-related, 15 methylation-related, and 10-transposition related piRNAs, which implied that these piRNAs might regulate sex development of C. semilaevis through epigenetic regulation or transposition. Finally, six markers that were verified by qRT-PCR were selected from theas signature 14 miRNAs. This work provides the basis for a method to identify the sex of fish by employing piRNAs derived from exosomes as biomarkers, which might prove to be applicable to other species in the future.

Materials and Methods
Isolation of exosomes from C. semilaevis. C. semilaevis specimens were obtained from Weizhuo Ltd., Hebei province, China. As Transmission electron microscopy.
Exosomes were fixed to formvar-carbon coated 300 mesh copper grids 36 .The absorbed exosomes were then negatively stained with 3% phosphotungstic acid and dried at room temperature for 20 min 37 . Subsequently, the exosomes were observed under a transmission electron microscope (Hitachi, H600IV, Japan) and the images were captured using a digital camera (Sony, Tokyo, Japan) 41-43 .

Nanoparticle tracking analysis.
Quantification and size analysis of the purified exosomes (5 μl, 3 three times) were performed using a NanoSight NS300 instrument (Malvern Instruments, Westborough, MA, USA) 44 . Vesicles were visualized by using light scattering usingunder a light microscope 45 . Measurement of the exosome concentration was performed in this study by calculating particle size on a particle-by-particle basis in a 60-second video recorded at a frame rate of 25 frames/s to provide accuracy and statistics for further analysis 46,47 .The results were subsequently analysed using the NTA 2.3 software (Malvern Instruments).

Western blotting analysis.
Exosomal extracts extracts (each sample contains contained at least 10 μg of total protein at least 10 μg) were separated by using sodium dodecyl sulphate- . The basic reads were converted into sequence data (also called raw data/reads) by base calling. Low-quality reads were filtered 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 15 nt, and those longer than 41 nt were filtered out of the raw data to obtain clean reads.
The expression patterns of known piRNAs in the different donors were analysed. Unannotated reads were analysed using Piano 51 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 identified with the threshold of a p value of p < 0.05. The p value was calculated using the DEG algorithm 52 in the R package with the Audic-Claverie statistic 53 without biological replicates. The targets of the differentially expressed piRNAs were predicted using the MiRanda software 54 for animals, with the following parameters: S ≥ 150; ΔG ≤ −30 kcal/mol and strict 5 seed pairing. Differentially expressed piRNAs were screened out according to the criteria detailed in the methods section section shown in (supplementary datasets4dataset 4).GO enrichment and KEGG pathway enrichment analysis of differentially expressed piRNA-target genes were performed using R based on the hypergeometric distribution.
Real Real-time quantitative PCR to verify signature piRNAs expression.
Exosomal piRNA expression was assayed using real qRT-PCR. Total RNAs (20 ng) from the two donor groups were used for reverse transcription. The RT reactions from total RNA were performed using specific primers Taqman MicroRNA primers (seen shown in the supplementary datasets 11) and a thermal cycler under the following conditions: 30 min at 65 °C, 50 min at 42 °C, and 5 min at 95 °C. The products were stored at −20 °C for later use or immediately processed according to the manufacturer's protocol. Quantitative PCR was performed in 96-well reaction plates using a QuantStudio6 Flex Real-Time PCR System (Thermo Fisher Scientific).
No template controls were used to evaluate background signal. The qPCR program consisted of 95 °C for 5 min, followed by 40 cycles each of denaturation at at 95 °C for 5 s and annealing and extension for 30 s at 60 °C. The expression level of U6 was used as a stable endogenous control for normalization. Each donor sample was run in triplicate and the relative quantification of piRNA expression was calculated using the 2 − ΔΔCt method 55 . Venn analysis among 37 sex-related, 15 methylation-related, and 10 transposition-related pi