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¹ (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.² The structural features of the deduced protein sequences were analyzed via SMART webtool³ (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.⁴ 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.

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 pseudo-replicate 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.

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.

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).

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.

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 base-pairing 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⁵ (Enright et al., 2003) with an alignment score > 150 and energy < –10 kcal mo^–1.

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-years-old 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.

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- and3-years-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 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 ping-pong 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.

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).

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 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 LNA-antagonist 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.

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.