Conserved and Widespread Expression of piRNA-Like Molecules and PIWI-Like Genes Reveal Dual Functions of Transposon Silencing and Gene Regulation in Pinctada fucata (Mollusca)

PIWI proteins and PIWI-interacting RNAs (piRNAs) suppress transposon activity in animals, thus safeguarding the genome from detrimental insertion mutagenesis. Recent studies revealed additional targets and functions of piRNAs in various animals. piRNAs are ubiquitously expressed in somatic tissues of the pearl oyster Pinctada fucata, however, the role of somatic piRNAs has not well characterized. This study reports the PIWI/piRNA pathway, including piRNA biogenesis and piRNA-mediated transposon silencing, and gene regulation in P. fucata. The biogenesis factors of PIWI, Zucchini, and HEN1, which are ubiquitous in somatic and gonadal tissues, were first identified in P. fucata using transcriptome analysis. Bioinformatics analyses suggested that different populations of piRNAs participate in the ping-pong amplification loop in a tissue-specific manner. In addition, a total of 69 piRNA clusters were identified in the genome of P. fucata based on the expression of piRNAs, which contained 26% transposons and enhanced for DNA/Crypton, LINE/CR1, SINE/Deu, and DNA/Academ. The expression patterns of the piRNAs and piRNA clusters in somatic tissues were not substantially different, but varied significantly between the somatic and gonadal tissues. Furthermore, locked-nucleic-acid modified oligonucleotide (LNA-antagonist) was used to silence single piRNA (piRNA0001) expression in P. fucata. Hundreds of endogenous genes were differentially expressed after piRNA silencing in P. fucata. Target prediction showed that some endogenous genes were targeted by piRNA0001, including twelve upregulated and nine downregulated genes after piRNA0001 silencing. The results indicated that piRNAs from somatic tissues may be related to gene regulation, whereas piRNAs from gonadal tissues are more closely associated to transposon silencing. This study will enhance our understanding of the role of piRNAs in mollusks, transposon silencing, and the regulatory function of the PIWI/piRNA pathway on protein-coding genes outside of germ line cells in P. fucata.

PIWI proteins and PIWI-interacting RNAs (piRNAs) suppress transposon activity in animals, thus safeguarding the genome from detrimental insertion mutagenesis. Recent studies revealed additional targets and functions of piRNAs in various animals. piRNAs are ubiquitously expressed in somatic tissues of the pearl oyster Pinctada fucata, however, the role of somatic piRNAs has not well characterized. This study reports the PIWI/piRNA pathway, including piRNA biogenesis and piRNA-mediated transposon silencing, and gene regulation in P. fucata. The biogenesis factors of PIWI, Zucchini, and HEN1, which are ubiquitous in somatic and gonadal tissues, were first identified in P. fucata using transcriptome analysis. Bioinformatics analyses suggested that different populations of piRNAs participate in the ping-pong amplification loop in a tissue-specific manner. In addition, a total of 69 piRNA clusters were identified in the genome of P. fucata based on the expression of piRNAs, which contained 26% transposons and enhanced for DNA/Crypton, LINE/CR1, SINE/Deu, and DNA/Academ. The expression patterns of the piRNAs and piRNA clusters in somatic tissues were not substantially different, but varied significantly between the somatic and gonadal tissues. Furthermore, locked-nucleic-acid modified oligonucleotide (LNA-antagonist) was used to silence single piRNA (piRNA0001) expression in P. fucata. Hundreds of endogenous genes were differentially expressed after piRNA silencing in P. fucata. Target prediction showed that some endogenous genes were targeted by piRNA0001, including twelve upregulated and nine downregulated genes after piRNA0001 silencing. The results indicated that

INTRODUCTION
Animal species express three types of endogenous silencinginstigating small RNAs: microRNAs (miRNAs), endogenous siRNAs (endo-siRNAs), and PIWI-interacting RNAs (piRNAs), based on their biogenesis mechanism and type of Argonaute-binding partners (Kim et al., 2009). miRNAs and endo-siRNAs, which are usually 20-23 nucleotides (nt) in length, are generated from double-stranded precursors by Dicer (Peters and Meister, 2007;Kim et al., 2009), whereas piRNAs are generated from single-stranded precursors independent of RNase III enzymes (Houwing et al., 2007), which are necessary for miRNA and endo-siRNA biogenesis. piRNAs are associated with PIWI subfamily members of the Argonaute family of proteins, whereas miRNAs and endo-siRNAs are associated with AGO subfamily members. The large Argonaute superfamily is characterized by a conserved PAZ domain, which is a single-stranded nucleic acid binding motif (Cerutti et al., 2000), and the PIWI domain, which implements RNase H slicing activity (Liu et al., 2004;Song et al., 2004;Tolia and Joshua-Tor, 2007).
The mechanisms underlying piRNA biogenesis and function remain largely unknown, mainly because the process has few similarities with the miRNA and endo-siRNA pathways. However, significant progress has recently been made, particularly in the area of piRNA biogenesis (Izumi et al., 2020). A comprehensive computational analysis of piRNA populations generated two models for piRNA biogenesis in various animals: the primary biogenesis pathway and the amplification loop or ping-pong cycle (Ishizu et al., 2012;Ross et al., 2014). In the primary biogenesis pathway, long piRNA precursors are transcribed from specific genomic loci called piRNA clusters, cleaved, and modified by intricate factors in the cytoplasm before being transported into the nucleus by specific transcription factor complexes (Izumi and Tomari, 2014;Ross et al., 2014). Primary piRNAs undergo an amplification process to induce high piRNA expression, known as the amplification loop or ping-pong cycle (Ishizu et al., 2012). The Zucchini (Zuc) endonuclease potentially forms the piRNA 5 end, whereas the 3 end is 2 -O-methylated by HEN1/Pimet, associated with PIWI proteins (Saito et al., 2007;Ipsaro et al., 2012). In addition, an uncharacterized 3 -5 exonuclease such as Trimmer was observed to trim the 3 end of piRNAs in silkworms (Kawaoka et al., 2011;Izumi et al., 2016). Hsp83/Shu may play a role in the PIWI loading step, and Hsp90/FKBP6 plays a role in secondary piRNA production (Ishizu et al., 2012). Despite the recent characterization of several factors in the piRNA biogenesis pathway, it remains poorly understood. PIWI proteins and their associated piRNAs suppress transposon activity in various animals, thereby safeguarding the genome from detrimental insertion mutagenesis (Aravin et al., 2001(Aravin et al., , 2007Carmell et al., 2007;Yang and Xi, 2017). However, many animals produce piRNAs that do not match transposon sequences determined in Caenorhabditis elegans  and mouse genome (Aravin et al., 2007), suggesting additional piRNA targets and functions. Nucleic small RNA-mediated silencing has gained significant attention in recent years, and numerous advances have been made in this field. Recent findings describing the role of transposon transcriptional regulation downstream of piRNAs have provided a paradigm for studying transcriptional regulation by small RNAs in animals. The use of C. elegans, Drosophila melanogaster, and mice as simple model systems, can offer important new insights into this process (Weick and Miska, 2014;Zhang et al., 2018). For example, Fas3 was postulated to be regulated by piRNAs generated from the 3 untranslated regions (3 UTR) of the traffic jam transcript (Saito et al., 2009), whereas traffic jam protein levels were elevated in PIWI mutants, indicating a cisregulatory mechanism for piRNA action in D. melanogaster (Robine et al., 2009).
Mollusks are one of the largest groups of marine animals, of which Pinctada fucata is well studied due to its economic potential for pearl production and a model organism to investigate the fascinating biology of mollusks. Recent studies have revealed that PIWI proteins are ubiquitously expressed in mollusks (Rajasethupathy et al., 2012;Ma et al., 2017;Jehn et al., 2018). Moreover, abundant piRNA-sized small RNAs were observed in the oysters Crassostrea gigas Zhou et al., 2014), Pinctada martensii (Jiao et al., 2014), and the scallop Chlamys farreri . Our previous study characterized the ubiquitous expression of piRNAs in somatic and gonadal tissues (Huang et al., 2019a). However, piRNA biogenesis and their functions in P. fucata remain largely unknown. This study aimed to demonstrate the biogenesis and function of PIWI/piRNA, and investigate its transposon silencing and endogenous gene regulation in P. fucata. Several key piRNA biogenesis factors (PIWI, Zucchini, and HEN1) have been identified in P. fucata, and the primary and secondary biogenesis pathways were determined by comprehensive computational analysis. To understand piRNA-regulated endogenous gene expression, we examined the effect of silencing of one of the most abundant piRNAs (piRNA0001) on potential target genes, such as miRNAs (predicted by a bioinformatics algorithm), using an LNA-antagonist.

MATERIALS AND METHODS
Transcriptomic Analysis of piRNA Biogenesis Factors in P. fucata Six (female: male = 1: 1) 1-year-old and six (female: male = 1: 13-year-old pearl oysters were collected from the Mikimoto Pearl Research Institution Base, Mie Prefecture, Japan in May 2018. Adductor muscle (Ad), gill (Gi), mantle (Ma), and gonad (Go) tissues were collected and individually transferred to 2-mL tubes containing RNA later (QIAGEN, Maryland, United States) separately. The samples were stored overnight at 4 • C and preserved at −80 • C until further analysis. Total RNA was isolated from each sample using the RNeasy Mini Kit (QIAGEN), according to the manufacturer's instructions. After assessing RNA quality and quantity using the Agilent 2200 TapeStation (Agilent Technologies, Waldbronn, Germany), three isolated total RNA samples were mixed at equivalent concentrations to construct an RNA sequencing library for each tissue sample at different ages. A total of 2 µg of mixed total RNA was used for library construction according to the manufacturer's protocol and subjected to paired-end sequencing on the Illumina HiSeq 4,000 platform. High-quality reads were required for de novo assembly analysis. Before assembly, raw reads were trimmed by removing adapter sequences and low-quality reads. Clean reads were assembled into unigenes using the Trinity software (Grabherr et al., 2011), and transcriptome annotation was performed using the Trinotate pipeline 1 (Bryant et al., 2017). The edgeR package in R software was used to identify and draw the significant differentially expressed genes (DEGs) between two different libraries with the following threshold: P-value < 0.05, and folds > 2 ( Robinson et al., 2010).
The consensus sequences of piRNA biogenesis factors were assembled using CLC Genomics Workbench version 8.0.1 (QIAGEN) and confirmed by published homolog sequences using BLASTp. Open reading frames (ORFs) were identified using the ORF finder available at NCBI. 2 The structural features of the deduced protein sequences were analyzed via SMART webtool 3 (Letunic and Bork, 2017). The predicted molecular weights and isoelectric points were determined using ExPASy (Gasteiger et al., 2005). All nucleotide and deduced amino acid sequences were electronically edited using BioEdit (Hall, 1999). Sequence similarity of the deduced proteins was performed using BLAST. 4 Homologous sequences of piRNA biogenesis factors were downloaded from NCBI Protein databases and the alignment of amino acid sequences was performed using MEGA5.0 (Tamura et al., 2014). The neighbor-joining method from MEGA 5.0 was used for building phylogenic trees with 1,000 bootstrap replications.
piRNA Processing in P. fucata Eight small RNA libraries from the adductor muscle, gill, ovary, and mantle tissues of P. fucata were sequenced using the Ion Proton system in our laboratory (the sequencing raw data were stored in the DNA Data Bank of Japan (DDBJ) under accession number DRA006953). Reads that did not produce a match to known non-coding RNAs were considered putative piRNAs and were merged for piRNA cluster prediction using proTRAC (v2.4.2) with default settings (Jehn et al., 2018). piRNA-sized molecules, approximately 30 nt in length, were widely expressed in all tissues and were used for ping-pong amplification analysis in P. fucata. Small RNA reads were pooled from the same tissues to perform piRNA annotation by the unitas (v1.5.3), which was run with the option -pp (Gebert et al., 2017). To compare ping-pong signatures and the number of sequence reads within different tissue types, PPmeter (v0.4) was used to quantify and compare the extent of ongoing ping-pong amplification (Jehn et al., 2018). Pseudo-replicates were generated by repeated bootstrapping (default = 100) of a fixed number of sequence reads (default = 1,000,000) from a set of original small RNA sequence datasets. Thereafter, the ping-pong signature of each pseudoreplicate was then calculated, and the number of sequence reads that participated in the ping-pong amplification loop was counted. The resulting parameter (ping-pong reads per million bootstrapped reads, ppr-mbr) were subsequently used to quantify and compare ping-pong activity in different tissue datasets. Both tools were operated using default parameter settings.

Transposon Element Annotation
We performed a de novo prediction of repetitive elements in the genome of P. fucata using RepeatMasker (v4.0.7) (Chen, 2004) based on a de novo repeat library, constructed using RepeatModeler (v1.0.11) (Flynn et al., 2020). An analysis was also conducted with the entire repeat data set and repeats localized in predicted piRNA clusters in P. fucata. The enrichment of transposon categories was compared between the genome and piRNA clusters. The relative expression patterns of piRNA clusters were calculated as mapped reads normalized for piRNA cluster length according to transcripts per million reads (TPM) (Li and Dewey, 2011). ggplot2 and pheatmap packages in R were used to perform principal component analysis (PCA) and heat maps of piRNA and piRNA cluster expression performance in the somatic and gonadal tissues. A t-test was used to compare the transposon percentages among different piRNA clusters between the somatic and gonadal tissues.
piRNA0001 Silencing in P. fucata and Quantification of piRNA Expression A single piRNA with the sequence 5 -UACUUUAACAUG GCACAGAUAUAAUGACCU-3 (piRNA0001) showed the highest expression in P. fucata somatic tissues, which was not annotated in any piRNA clusters. To explore its function, an LNA-antagonist was used to silence its expression in P. fucata. Eight P. fucata oysters (approximately 3 years old) were randomly divided into two groups: the LNA-antagonist group (LNA) and the control group (Con). Oysters from the LNA group were injected with 100 µL of 5 mg mL −1 LNA-antagonist probe solution into their body cavity. In contrast, those in the Con group were injected with isotonic PBS solution. The sequence of the high-affinity LNA-antagonist was 5 -AggTcaTtaTatCtgTgcCa-3 (LNA in capitals) (Elmén et al., 2008). 2 weeks after injection, all individuals were subjected to tissue collection. Somatic tissues, including the adductor muscle, gill, and mantle, and gonadal tissues were collected and stored separately in 2-mL tubes containing RNA later (QIAGEN) overnight at 4 • C, and preserved at -80 • C until further use. Stem-loop RT-PCR was used for piRNA0001 relative expression analysis using U6 sRNA as a reference gene. This process was performed as described in our previous study (Huang et al., 2019b).

Transcriptomic Analysis of Gene Expression Profiles After piRNA0001 Silencing in P. fucata
Twenty-four libraries of somatic tissues, including adductor muscle, gill, and mantle, with four replicates for each tissue type from the LNA and Con groups, were constructed for RNA sequencing in 2017. Library construction, sequencing, de novo assembly, and functional annotation were performed as described above. Mapped fragments were normalized for RNA length according to the TPM method (Li and Dewey, 2011), facilitating the comparison of transcript levels between libraries. The analysis of DEGs between the Con and LNA groups was the same as in the previous method using edgeR.
Prediction of piRNA0001 Targeting Sites in P. fucata piRNAs possess second to eighth nucleotide seed sequence that requires near-perfect complementarity between the piRNA and target genes, such as miRNAs (Gou et al., 2014). Additional basepairing outside of the seed is also important for piRNA targeting. piRNAs can only tolerate a few mismatches outside the seed region (Shen et al., 2018;Zhang et al., 2018). In this study, 3 UTR sequences of DEGs were used for piRNA target site predictions. No piRNA target rules have been explored in mollusks, therefore, we used miRNA-like pairing rules to predict potential target sites (Wang et al., 2018). These were predicted between piRNA and mRNA 3 UTR using miRanda tools 5 (Enright et al., 2003) with an alignment score > 150 and energy < -10 kcal mo −1 .

Gene Expression Analysis by RT-PCR
To validate PIWI gene expression in P. fucata, the early development stage was achieved by artificial fertilization, and the somatic and gonadal tissues were dissected from 1-to 3-yearsold oysters. The corresponding primers were designed using Primer5.0, based on the complete sequence with the following parameters: primer length between 18 and 22 bp, melting temperatures (T m ) between 55 and 65 • C, product size between 80 and 150 bp, and avoidance of primer dimers or hairpin structure (Supplementary Table 1). β-actin was used as reference gene for PIWI relative expression calculations. First-stand cDNA was synthesized using the PrimerScript RT Master Mix reagent Kit (TaKaRa, Shiga, Japan), according to the manufacturer's instructions. To determine the expression patterns of PIWI genes in various tissues at different developmental stages, all samples were analyzed using the Applied Biosystem 7,300 Fast Real-Time PCR System (Life Technologies, California, United States). RT-PCR analysis was performed in a 96-well plate, in which each well contained 20 µL of reaction mixture consisting of 10 µL SYBR Premix Ex Taq II (TaKaRa), 0.4 µL ROX Reference Dye (50×) (TaKaRa), 0.8 µL of each primer (10 µM), 2 µL of complementary DNA (cDNA) template and 6 µL of sterilized double-distilled water. RT-PCR conditions were as follows: pre-denaturation at 95 • C for 30 s, followed by 40 cycles of amplification at 95 • C for 5 s and 60 • C for 31 s, a dissociation stage with one cycle at 95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s after amplification. Each sample was tested in triplicate. The average value per gene was calculated from three replicates.
Nine DEGs, predicted to be targeted by piRNA0001, were also selected for relative expression analysis by RT-PCR as described above, with total RNA consistent with RNA_Seq2 libraries. The resulting PCR products contained the predicted interaction sites on the mRNA 3 UTR. The relative expression of each mRNA quantified by RT-PCR was presented as n = 4, from three replicates. Ct values were represented by the mean values of three independent replicates. The relative expression level of each gene was calculated using the 2(-Ct method (Pfaffl, 2011). A t-test was used to compare the differences in the relative expression of LNA/Con between RNA-seq and RT-PCR.

RESULTS
piRNA Biogenesis Pathway in P. fucata Following clean-up and quality checks, sixteen transcriptomic libraries with 191.11 million reads, with a total length of 28.67 Gb, were used for de novo assembly (Supplementary Table 2). A total of 324,779 transcripts were assembled, with a mean length of 564 bp (Table 1, RNA_Seq1). Complete sequences of PIWI (Piwil1 and Piwil2), Zucchini (Zuc), and HEN1 were assembled based on the annotation results of transcripts, and their respective lengths are listed in Table 2. PIWI possessed the conserved PAZ and PIWI domains as homologs in other organisms, whereas the Zuc protein contained the PLDc domain (Figure 1). The conserved identities of the PAZ and PIWI domains were 51.37 and 60.54%, respectively. Piwil1 and Piwil2 clustered into two categories of the PIWI subfamily (Supplementary Figure 1).
The somatic and gonadal tissues of 1-and 3-years-old individuals were used for transcriptome analysis to study the expression profiles of the piRNA biogenesis factors. Gene expression was calculated as transcripts per million (TPM) using transcriptomic analysis. Piwil1 and Piwil2 were widely expressed in the P. fucata somatic and gonadal tissues, with significantly higher expression in the gonadal tissues at 1-and3years-old stages (Figures 2A,B). Piwil1 expression levels were significantly higher than those in Piwil2. Zuc was highly expressed in the gonadal tissues, particularly in3-years-old females (Figure 2C), and HEN1 was expressed in all examined tissues ( Figure 2D). Furthermore, transcriptome and early development stage samples were used for gene expression verification. RT-PCR RNA_Seq1 represents the transcriptomic analysis of piRNA biogenesis factors in P. fucata, RNA_Seq2 represents the transcriptomic analysis of gene expression profiles after piRNA0001 silencing in P. fucata.
analysis validated Piwil1 and Piwil2 expression levels in larvae, juveniles, and adult P. fucata, which were high in fertilized eggs (0 hpf), significantly decreased in embryos (0.5 hpf), and weak from the blastula (4 hpf) to spat (30 dpf) stages ( Supplementary  Figures 2A,B). Piwil1 and Piwil2 were ubiquitously expressed in the somatic and gonadal tissues and were highly expressed in female P. fucata, both at the 1-and 3-years-old stages. These piRNA biogenesis factors are highly expressed in the gonadal tissues of other animals; however, they are weakly but ubiquitously expressed in somatic tissues, indicating their contribution to piRNA biogenesis in P. fucata somatic tissues. P. fucata encodes two ubiquitously expressed PIWI proteins, therefore, we further analyzed the participation of distinct piRNA populations in the ping-pong amplification loop. We employed a bioinformatics approach, under the premise that Piwil1 and Piwil2 bind to piRNAs with different length profiles, similar to their corresponding mouse homologs Miwi and Mili, which particularly bind to 29/30 nt and 26/27 nt long piRNAs, respectively (Vourekas et al., 2012). We analyzed pairs of P. fucata sequenced reads with a 10-bp 5 overlap (ping-pong pairs), which is the typical sequence length of each ping-pong partner ( Figure 3A). A significant 10 nt overlap between sense and antisense strands was detected in the small RNAs derived from P. fucata, which is consistent with piRNAs generated by the ping-pong amplification mechanism (Supplementary Figure 3). In somatic tissues, ping-pong pairs combine 29/30 nt long piRNAs, suggesting Piwil1-Piwil1-dependent homotypic pingpong amplification. In the gonadal tissues, most ping-pong pairs combine piRNAs predominantly 29/30 nt and 25/26 nt long, suggesting both Piwil1-Piwil2-dependent heterotypic and Piwil1-Piwil1-dependent homotypic ping-pong amplification, as shown in Figure 3B. Furthermore, 29/30 nt long piRNAs that presumably bound to Piwil1 were heavily biased for a 5 uridine (U), whereas 25/26 nt long piRNAs that presumably bound to Piwil2 showed a stronger bias for an adenine (A) at position 10.
piRNA-Mediated Transposon Silencing in P. fucata Small RNA sequencing was performed to identify putative piRNAs in P. fucata (Huang et al., 2019a). Reads that did not produce a match to known non-coding RNAs were considered as putative piRNAs and were merged for piRNA processing.
A total of 69 piRNA clusters (total size 0.58 Mb, 0.056% of the genome) were identified in the genome of P. fucata ( Figure 4A and Supplementary Table 3). The P. fucata genome contains 539.3 Mb repetitive sequences accounting for 52% of the genome; however, 23.3% of repetitive sequences are classified as unknown ( Figure 4B). In line with the TE suppressive role of piRNAs, the identified piRNA clusters also showed a lower enrichment for transposon sequences compared to the whole genome situation, which contained 15% known and 11.06% unknown transposons ( Figure 4C). The composition of transposons in piRNA clusters does not at all reflect to the transposon landscape of the whole genome. Therefore, piRNA clusters were enriched for DNA/Crypton, LINE/CR1, SINE/Deu, and DNA/Academ showing up to 28-fold enrichment in piRNA clusters ( Figure 4D). In addition, a total of 43,266 unigenes (13.32% of the total) from RNA-Seq1 were annotated with transposon consensus (Supplementary Figure 4). However, thousands of distinct piRNAs do not map to the transposons in P. fucata.
It is well known that piRNAs suppress transposons in germ line (Ishizu et al., 2012;Yang and Xi, 2017). The expression performance of putative piRNAs was compared between the somatic and gonadal tissues. The PCA of putative piRNA expression patterns from all examined libraries clustered the somatic tissues together and differentiated the somatic and gonadal tissues ( Figure 5A). All putative piRNA reads were re-mapped with piRNA clusters to calculate the relative expression performance using the TPM method. The expression performance of these piRNA clusters was also differentiated the somatic and gonadal tissues; however, there was small difference within different somatic tissues ( Figure 5B). In addition, two different categories of piRNA clusters were clustered in the somatic and gonadal tissues ( Figure 5B). The piRNA clusters of Class 1 were highly expressed in the gonadal tissues, whereas a piRNA clusters of Class 2 were highly expressed in the somatic tissues. The transposon element percentage of each piRNA cluster was compared between Class 1 and Class 2. The piRNA clusters of Class 1 had an average transposon percentage of 32.16%, which was significantly higher than that of the piRNA clusters of Class 2 (1.85×, P < 0.01) ( Figure 5C).

piRNA-Mediated Gene Regulation in P. fucata
LNA nucleotides can be hybridized with DNA or RNA residues in the oligonucleotide to form stable heteroduplexes between the LNA-antagonist and small RNA to silence its expression (Elmén et al., 2008). A modified LNA-antagonist was used to silence the expression of piRNA0001 in P. fucata. All samples were dissected for tissue collection 2 weeks after the injection. Following RNA extraction, stem-loop RT-PCR was used for the relative quantitative analysis of piRNAs. piRNA0001 expression levels in LNA-introduced pearl oysters were decreased in the adductor muscle (0.04×, p < 0.01), gill (0.19×, p < 0.01), gonad (0.28×, p < 0.01), and mantle tissues (0.62×, p < 0.05) compared to those of the control group ( Figure 6A). Moreover, the piRNA0001 precursor detected by transcriptomic analysis of gene expression  The PAZ domain contributes to specific and productive incorporation of miRNAs, endo-siRNAs, and piRNAs into the RNAi pathway; the PIWI domain can be inferred to cleave single-stranded RNA, such as mRNA should be guided by single-stranded small RNAs; PIWI domains place the guide uniquely in the proper position in Argonaute-RNA complexes. PLDc produces phosphatidic acid from phosphatidylcholine, which may be essential for the formation of certain types of transport vesicles or involved in constitutive vesicular transport via signal transduction pathways.
profiles showed no significant differences in the expression between the LNA and Con groups (Supplementary Figure 5). Twenty-four libraries from the Con and LNA groups, including adductor muscle, gill, gonad, and mantle tissues were constructed for RNA sequencing, with three replicates. Following data cleaning and quality checks, RNA sequencing produced 265.71 million clean reads (Supplementary Table 4), that is, a total length of 26.57 Gb, further assembled into 213,323 transcripts by Trinity (Table 1, RNA_Seq2). The Trinotate pipeline annotated 48.63% of these transcripts, based on five public databases (NT, Swiss-Prot, Pfam, GO, and KEGG), and identified 31. 06, 63.37, 80.60, 60.88, and 38.14% transcripts in the NT, Swiss-Prot, Pfam, GO, and KEGG databases, respectively. In addition, most transcripts were blasted with P. fucata and other mollusk species (Supplementary Figure 6).
Analysis of DEGs in somatic tissues between Con and LNA groups was performed to understand further the molecular events involved in the functions of piRNA0001 in P. fucata. We selected the genes from at least one group with a mean TPM value > 5 for analysis using edgeR (P-value < 0.05 and fold > 2). A total of 2,904 transcripts showed significant differential expression between the Con and LNA groups (Supplementary Table 5). Pairwise comparison of the adductor muscle revealed 456 differentially expressed transcripts, including 220 upregulated and 236 downregulated genes after LNAantagonist treatment (in comparison with the Con group). In addition, analysis of 900 transcripts revealed 440 upregulated and 460 downregulated DEGs in the gill tissue and 374 upregulated and 325 downregulated DEGs in the mantle tissues, following LNA-antagonist treatment ( Figure 6B). Among these, 16 and 15 DEGs were typically upregulated and downregulated, respectively, in the three examined LNA-antagonist-treated tissue types.
All DEGs were used for piRNA0001 target prediction, with genes containing target sites in the 3 UTR considered as potential target genes. Four, nine, and eleven DEGs were predicted to be targeted by piRNA0001 in the adductor muscle, gill, and mantle tissues, respectively, with an average alignment score of 162.17 and average energy value of −14.00 kcal mol −1 (Figure 7). Target prediction was anticipated for FHOD3, SRS10, VINC, and ZNF622 in the adductor muscle, for ART2, EF1A, EMC7, ESI1L, NAC2, PA2HB, SYWC, TM87A, and ZNF622 in the gill tissue, and for CAH14, CBF, CDC42, FRRS1, GOLI4, KAT3, NBR1, PA2HB, TYR1, TYR2, and ZNF622 in the mantle tissue. The predicted piRNA-mRNA interaction is shown in Supplementary Table 6, including two sites on both the ART2 and CBF 3 UTR. No polyA signal was observed in the assembly transcripts of ART2, EF1A, ESI1L, FHOD3, PA2HB, SRS10, and TM87A, and predicted interaction sites on CAH14, NBR1, SYWC, and VINC were located after the polyA signal (Supplementary Sequences). The second piRNA:mRNA interaction site on ART2, NAC2, and VINC, were located far from the stop codon (Supplementary Table 6). Alignment pairing from the second to eighth nucleotides, considered to be the piRNA seed region, showed perfect complementarity. Furthermore, base-pairing outside the seed region is also important for piRNA target prediction. Following piRNA0001 silencing in P. fucata, FHOD3, SRS10, and ZNF622 were upregulated, whereas VINC were downregulated in the adductor muscle tissue. In the gill tissue, ART2, EF1A, EMC7, ESI1L, SYWC, TM87A, and ZNF622 were upregulated, whereas NAC2 and PA2HB were downregulated. In the mantle tissue, CBF, CDC42, NBR1, and ZNF622 were upregulated, while CAH14, FRRS1, GOLI4, KAT3, PA2HB, YTR1, and TYR2 were downregulated. ZNF622 was upregulated in all the examined tissues. Both upregulated and downregulated genes were predicted to be targeted by piRNA0001, demonstrating the diverse ways of piRNA-mediated gene regulation in P. fucata.

Validation of Gene Expression by RT-PCR
Total RNA extracted from LNA and Con P. fucata somatic tissues was analyzed using RT-PCR, to determine the authenticity of mRNA expression calculated by RNA sequencing. Inconsistencies were detected in the adductor muscle VINC expression obtained from RNA sequencing and RT-PCR analysis (Figure 8), possibly caused by ineffective RT-PCR primers or by sequencing errors that occurred during the complex process, with the predicted products located after the polyA signal on the 3 UTR, far from the stop codon (Supplementary Sequences: VINC). However, the remaining eight predicted target genes demonstrated consistent levels of relative expression levels through RT-PCR, compared with the results obtained from RNA sequencing. Among these genes, ZNF622 was upregulated in all examined somatic tissues, following piRNA0001 silencing in P. fucata. In contrast, PA2HB and TYR1 were downregulated in the gill and mantle tissues, based on RNA sequencing and RT-PCR data.

DISCUSSION
Analysis of the PIWI/piRNA pathway representative of several animals revealed an extensive diversity of lineage-specific adaptations, challenging the universal validity of data obtained from model organisms. PIWI proteins are characterized by two protein domains, namely the PAZ domain, an RNA-binding motif that binds the 3 end of short RNAs, and the PIWI domain, which is structurally similar to the RNaseH catalytic domain (Parker and Barford, 2006). Numerous PIWI homologs have been identified in various organisms, such as four homologous Hiwi, Hili, Hiwi2, and Hiwi3 have been identified in Homo sapiens, and play a crucial role in human cancer and male germline cell development (Qiao et al., 2002;Sasaki et al., 2003;Liu et al., 2006); homologous Miwi, Mili, and Miwi2 were identified in Mus musculus, and their knockdown led to male sterility Aravin et al., 2006;Lau et al., 2006); and two homolog, Ziwi and Zili, were identified in Danio rerio, and were both were crucial for germ cell differentiation and meiosis (Houwing et al., 2007(Houwing et al., , 2008. These findings indicate that PIWI may play an important role in germline cell development in organisms. PIWI was first discovered in Drosophila, where it functions in germline cell maintenance and self-renewal (Cox et al., 1998). PIWI/piRNA research in germline cells has rapidly advanced. PIWI mutations lead to a profound infertility phenotype in mice Carmell et al., 2007). In addition to PIWI function in germline cells, recent studies have demonstrated the expression and function of PIWI in the soma, from low eukaryotes to mammals. PIWI expression has been detected in human cancer cells (Ross et al., 2014;Krishnan et al., 2016) and mammalian somatic tissues (Yan et al., 2011). Although somatic PIWI was previously detected in Drosophila fat bodies and ovarian somatic tissues (Malone et al., 2009;Jones et al., 2016), further evidence of somatic PIWI/piRNA expression was observed in 16 out of 20 surveyed arthropod species (Lewis et al., 2018). In mollusks, PIWI proteins are widely expressed in Lymnaea stagnalis and C. gigas somatic tissues, and homologous PIWIs were also detected in 11 mollusk species, suggesting the occurrence of somatic PIWI/piRNA expression in early bilaterian ancestors (Jehn et al., 2018). In this study, Piwil1 and Piwil2, assembled by RNA sequencing in P. fucata, were ubiquitously expressed in all of the examined somatic tissues and gonadal tissues. These results provide further evidence of somatic PIWI/piRNA expression as an ancestral bilaterian trait (Lewis et al., 2018).
piRNA biogenesis in metazoa involves synthesizing primary piRNA from a piRNA cluster, assisted by several factors in the cytoplasm and nucleus. Thereafter, the long, singlestranded piRNA precursor was exported from the nucleus to the cytoplasm, where it was shortened to piRNA-like small RNAs by an undetermined endonuclease (Ishizu et al., 2012;Ross et al., 2014). Recent studies have indicated that Zuc may be an endonuclease that forms the 5 end of piRNAs in Drosophila and mice (Ipsaro et al., 2012;Inoue et al., 2017;Izumi et al., 2020). The 3 end was 2-O-methylated by HEN1 (Saito et al., 2007), whereas an uncharacterized 3 -5 exonuclease has been shown to trim the 3 end of piRNAs (Kawaoka et al., 2011;Izumi et al., 2016). In this study, several crucial piRNA biogenesis factors, including PIWI, Zuc, and HEN1, were assembled by RNA sequencing. These were ubiquitously expressed in all examined tissues, particularly in gonads, which previously showed a higher percentage of piRNA expression than somatic tissues (Huang et al., 2019a).
The ping-pong amplification loop is responsible for posttranscriptional silencing of transposable elements (Jehn et al., 2018). In Drosophila and mice, this process normally involves two PIWI proteins (heterotypic ping-pong), one loaded with antisense piRNAs targeting piRNA cluster transcripts, which contain transposon sequences in an antisense orientation (Aravin et al., 2008). The homotypic Aub:Aub ping-pong process occurs in Drosophila  and wild-type prenatal mouse testes (Miwi2:Miwi2 and Mili:Mili) (Aravin et al., 2008). In P. fucata, we also determined the amplification system for a secondary piRNA biogenesis pathway using a comprehensive computational analysis. A Piwil1-Piwil1-dependent homotypic ping-pong amplification loop was observed in the somatic and gonadal tissues, with a Piwil1-Piwil2 dependent heterotypic ping-pong amplification loop, simultaneously occurring in the gonadal tissues; however, we cannot determine whether the binding preferences of PIWI proteins have changed in P. fucata. Both PIWI proteins may bind to the entire range of piRNAs. However, based on the presence of piRNA populations with length profiles and their representation in ping-pong pairs, together with the differences in their number of 1U and 10A reads, we believe the above explanation is a reasonable and parsimonious interpretation of the data while acknowledging the possibility of others.
The expression of piRNA biogenesis factors in somatic tissues implies the potential existence of somatic piRNAs. Recent studies have revealed somatic piRNAs from sponges to humans (Yan et al., 2011;Ross et al., 2014). A functional non-gonadal somatic piRNA pathway in Drosophila fat bodies affects normal metabolism and overall organismal health (Jones et al., 2016). Somatic piRNAs are considered an ancestral trait of arthropods, which predominantly target transposable elements, suggesting that the piRNA pathway was active in the soma of the last common ancestor of arthropods to maintain mobile genetic elements in check (Lewis et al., 2018). Abundant putative piNRAs were observed in the somatic and gonadal tissues of P. fucata in our previous study (Huang et al., 2019a). In this study, we analyzed the expression of piRNAs and piRNA clusters in the somatic and gonadal tissues of P. fucata. The PCA analysis of piRNAs and heatmap of expression patterns of piRNA clusters clustered the somatic tissues together, therefore, the piRNA and piRNA clusters expressed in somatic tissues were not substantially different, but varied significantly between the somatic and gonadal tissues. The transposon annotation results of piRNA clusters also indicated that the gonadal piRNAs are more likely related to transposon silencing. These results imply FIGURE 6 | piRNA0001 silencing in P. fucata. (A) piRNA0001 expression profiles in P. fucata treated with LNA-antagonist. U6 was used as the reference gene. The relative expression of piRNA0001 was normalized by piRNA0001 expression following LNA-antagonist treatment. Differences were statistically analyzed between Con and LNA groups using one-way analysis of variance (ANOVA). Significant differences are marked with * (p < 0.05) and ** (p < 0.01). (B) Number of differentially expressed genes (DEGs) between LNA and Con groups below the following threshold: P-value < 0.05 and folds > 2.
FIGURE 7 | Target predictions of DEGs in P. fucata. TPM shows the relative expression level of target genes analyzed by RNA-seq. Tissue: adductor muscle (Ad), gill (Gi), and mantle (Ma). Large difference in standard deviations of relative expression level of genes analyzed by RNA-seq might cause errors during the complex sequencing protocols or differences among individual samples. Tissue marked with an asterisk (*) represents differentially expressed gene under the following threshold: P-value < 0.05 and folds > 2 in the present tissue.
FIGURE 8 | Validation of potential piRNA0001 predicted target genes in P. fucata. t-test was used to compare the differences in relative expressions of LP/LC for the eight DEGs found common between RNA-seq and qPCR. Gene marked with an asterisk (*) represents significant differences in the gene expression trend observed between RNA-seq and RT-PCR analysis using a t-test (P < 0.05). different functions of piRNAs between the somatic and gonadal tissues in P. fucata.
To explore the function of gene regulation of piRNAs in P. fucata, an LNA-antagonist was used to silence single piRNA (piRNA0001) expression. LNA-antagonists have been used for mouse and non-human primate small RNA silencing without affecting health (Elmén et al., 2008). In this study, we used an LNA-antagonist to silence piRNA0001expression, which was the most highly expressed piRNA in P. fucata somatic tissues, followed by a stem-loop RT-PCR for piRNA0001 quantification analysis, as in previous studies (Hong et al., 2016;Wang et al., 2018;Zhang et al., 2018). Stem-loop RT-PCR is a powerful and reliable tool for quantitatively analyzing and monitoring dynamic changes in piRNAs, specifically at the cellular and tissue levels in invertebrates. In this study, piRNA0001 was effectively downregulated in somatic tissues by a specific LNA-antagonist.
Recent studies have shown that piRNAs may regulate endogenous mRNA expression in Drosophila and mice (Gou et al., 2014;Wu et al., 2018). Thousands of genes were differentially expressed after piRNA0001 silencing in P. fucata, which might be caused by piRNA0001 silencing, or changes in other biological processes, cellular components or molecular functions caused by LNA injection. Although targeting rules have been determined in Drosophila Zhang et al., 2018), the available tools cannot be used for piRNA targeting site prediction in other species. An RNA-RNA interacting prediction tool, such as miRanda (Gou et al., 2014), was alternatively used to predict potential targeting sites between piRNAs and endogenous genes. In this study, we used miRanda to predict potential piRNA-mRNA interaction sites on DEGs in P. fucata. piRNA target recognition was strikingly similar to that of miRNA pairing in the seed sequence (positions 2-8 relative to the 5 end of the piRNA), which is the primary determinant of target recognition. Therefore, seed pairing is not sufficient, and additional basepairing outside of the seed sequence is also important for piRNA target sites.
Among the predicted piRNA0001 target genes in P. fucata, the majority were related to metabolism and biological processes. ZNF622 was up-regulated in all examined tissues following piRNA0001 silencing, and was predicted to be targeted by piRNA0001 with an alignment score of 154.00 and energy value of −10.64 kcalmol −1 . Zinc finger proteins (ZNFs) are one of the most abundant groups of proteins, with a wide range of molecular functions, such as transcriptional regulation, ubiquitin-mediated protein degradation, signal transduction, actin targeting, DNA repair, cell migration, and numerous other processes (Cassandri et al., 2017). A single piRNA (fem piRNA) from a sex chromosome downregulates Masc mRNA, which encodes a C3Htype ZNF that induces masculinization and plays an important role in sex determination in silkworms (Kiuchi et al., 2014). CAH14, NBR1, SYWC, and VINC were predicted targets of piRNA0001, whereas interaction sites were located after the polyA signal on the 3 UTR, suggesting that these interaction sites may be false positive results. This may be due to the insufficient assembly result of RNA sequencing, with the polyA signal not appearing in the assembling gene sequences, as observed in ART2, EF1A, ESI1L, FHOD3, PA2HB, SRS10, and TM87A. Single and double target sites were both observed in these target genes, indicating that piRNA may regulate endogenous genes at multiple sites on the 3 UTR. Except for upregulated genes predicted to be targeted by piRNA0001, nine down-regulated genes may also have been targeted by piRNA0001, extending the uncertainty surrounding piRNAs in P. fucata. The lowest piRNA silencing efficiency was observed in the mantle tissue, possibly explaining the majority of downregulated genes. Although, we cannot finalize the target rules or genes of piRNAs, we have provided new insights into the gene regulatory function of piRNAs in P. fucata. The PIWI/piRNA might have a great diversity of functions in P. fucata, including but not limited to gene regulation. Further studies are needed to fully understand somatic piRNAs in mollusks.
In summary, this study reported a PIWI/piRNA pathway, including piRNA biogenesis and piRNA-mediated transposon silencing and gene regulation, in the pearl oyster P. fucata (Mollusca). These findings have successfully contributed to our understanding of the role of piRNAs in mollusks and will also help to gain further insights into the PIWI/piRNA pathway function outside of germline cells.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ddbj.nig. ac.jp/, DRA008674 and https://www.ddbj.nig.ac.jp/, DRA007432.

AUTHOR CONTRIBUTIONS
SA and SK designed the study. FO, KM, and KN provided the experimental animal samples. MA and SW performed