Piwi-interacting RNA (piRNA) expression patterns in pearl oyster (Pinctada fucata) somatic tissues

Piwi-interacting RNAs (piRNAs) belong to a recently discovered class of small non-coding RNAs whose best-understood function is repressing transposable element activity. Most piRNA studies have been conducted on model organisms and little is known about piRNA expression and function in mollusks. We performed high-throughput sequencing of small RNAs extracted from the mantle, adductor muscle, gill, and ovary tissues of the pearl oyster, Pinctada fucata. RNA species with sequences of approximately 30 nt were widely expressed in all tissues. Uridine at the 5′ terminal and protection from β-elimination at the 3′ terminal suggested that these were putative piRNAs. A total of 18.0 million putative piRNAs were assigned to 2.8 million unique piRNAs, and 35,848 piRNA clusters were identified. Mapping to the reference genome showed that 25% of the unique piRNAs mapped to multiple tandem loci on the scaffold. Expression patterns of the piRNA clusters were similar within the somatic tissues, but differed significantly between the somatic and gonadal tissues. These findings suggest that in pearl oysters piRNAs have important and novel functions beyond those in the germ line.

and closely related organisms. In particular, it will aid understanding of the physiological roles of small RNAs in mollusks, particularly in biomineralization.

Results
Sequencing of putative piRNAs from somatic and gonadal tissues of pearl oysters. Eight small-RNA libraries were sequenced using the Ion Proton system. We obtained 50.31 million raw reads. After removing the adapters, low-quality reads, reads with unknown nucleotides, and reads outside the 15-35 nucleotides range, 42.7 million reads ( Table 1, Step 1) remained for length-distribution analysis (Fig. 1). Two obvious peaks (at 21-23 nt and 29-31 nt) were observed in the distribution of the lengths of small RNAs from somatic tissues, and one peak at 29-31 nt was found in small RNAs from the gonadal tissue (Fig. 1A). After removing known RNAs, including miRNAs, rRNAs, tRNAs, snRNAs, and snoRNAs, 31.5 million clean reads (Step 2) remained for analysis. The size profile of small RNAs after removed all the known RNAs were shown in Supplementary Fig. S1, which have a low percentage of 21-23 nt RNAs. Then 22.6 million reads were mapped to the pearl oyster reference genome (Step 3). Finally, 18.0 million reads were considered valid putative piRNA reads, ranging from 26 to 31 nt in length (Step 4). Consistent with the signature of piRNAs observed in previous studies, 59.31% unique putative piRNA reads were uniquely mapped to the reference genome and the remaining unique reads were mapped to multiple loci of the reference genome (Fig. 1B). Of the unique putative piRNAs, 70.1-87.8% began with a U in each sequencing library (Fig. 1C). This was also consistent with previously documented features of piRNAs. Therefore, we deduced that the 18.0 million reads were mostly derived from piRNAs, and we therefore used them for our subsequent analyses. The final dataset contained 2.8 million unique sequences. For simplicity, in the following sections we refer to the putative piRNA reads as "piRNAs", and unique putative piRNA sequences as "unique piRNAs".
Further examination of the piRNAs revealed that a piRNA (piRNA0001) species with sequence 5′-UACUUUAACAUGGCACAGAUAUAAUGACCU-3′, which was mapped on scaffold1502.1 with multiple tandem manners, had the highest number of reads per million (RPM) of the piRNAs expressed in the somatic tissues, and that piRNA 5′-ACAGAAAUCUCGGAUCCAUCUAUAGACAGG-3′ was the most strongly expressed in the gonadal tissues. The top 20 most-highly expressed piRNAs (7.7% of reads) in the two types of tissues are shown in Supplementary Table S1. piRNAs were highly expressed in both somatic and gonadal tissues (Fig. 2). Approximately 18,179-33,301 unique piRNAs had an RPM > 5. The gonadal tissues had more unique piRNA species with low RPMs than the somatic tissues did.
Putative piRNA sequence mapping with the reference genome. It was thought that, in the organisms in which this has been studied, mature piRNAs originate from the processing of longer RNA precursors 27 . We mapped the putative piRNAs to the pearl oyster reference genome to identify piRNA gene clusters. The numbers of piRNA clusters identified based on different merge distances and criteria for minimum coverage are shown in Supplementary Fig. S2. We used moderate criteria (<1 kb merge distance and a minimum coverage of six reads) to identify piRNA clusters from 18.0 million piRNAs. A total of 35,848 piRNA clusters were defined, and the locations and reads per kilobase per million (RPKM) of each piRNA gene clusters are shown in Supplementary Table S2. Clusters are named in order in each scaffold. The sizes of the clusters ranged from 30 bp (single unique piRNA) to 60.9 kb, and 49.1% of the piRNA clusters were longer than 1 kb. piRNAs from 22,629, 24,933, 27,789, and 30,671 clusters were observed in the mantle, adductor muscle, gill, and gonad tissues, respectively. Putative piRNAs from 14,125 clusters were commonly expressed in all four tissues, whereas 2,896 clusters were exclusively observed in the gonadal tissues ( Supplementary Fig. S3). A principal component analysis (PCA) of piRNA-cluster expression patterns from all eight libraries clustered those of the somatic tissues together and differentiated the somatic and gonadal tissues ( Supplementary Fig. S4).
Mapping the piRNAs to the reference genome revealed two different mapping patterns: bidirectional and unidirectional (Fig. 3). We found that 40.69% of the unique piRNAs were mapped to the reference genome at multiple loci, and 25% of the unique piRNAs were mapped with multiple and tandem patterns on the same scaffold. Dot-matrix analysis of the piRNA gene clusters revealed that piRNAs were often encoded in a tandem manner (Fig. 4).
Putative piRNA clusters expressed differently in somatic and gonadal tissues. Analysis of the expression of piRNA clusters aids understanding of their expression patterns in the somatic and gonadal tissues of the pearl oyster. PCA analysis indicated that piRNA clusters from the somatic tissues have similar expression patterns but differ from gonadal tissues. To investigate the difference between piRNA cluster patterns in the somatic and gonadal tissues, we conducted a differential expression analysis using edgeR, with the following criteria: FDR < 0.01, and an absolute log2Ratio value of >1. The results showed no difference in the expression patterns of piRNA clusters among the somatic tissues. However, 127, 30, and 377 piRNA clusters were significantly differently expressed in the mantle, adductor muscle, and gill tissues than in the gonadal tissues, respectively (Supplementary Table S3).
β-elimination reaction confirms putative piRNA. piRNAs are modified by methylation at the 3′ terminal for protection from periodate oxidation. The β-elimination reaction identifies methylation at the 2′-OH of the 3′ terminal of small RNAs. The products of the β-elimination reaction of miR-279a showed increased migration in a 15% polyacrylamide gel, while the mobility of piRNA0001 was unaffected by the β-elimination reaction (Fig. 5), indicating that there is an important difference between miRNAs and putative piRNAs in terms of their molecular structure at the 2′-OH of the 3′ terminal.  Detection of putative piRNA in somatic tissues. When RNA extracted from gill, adductor muscle, mantle, intestine, and abdominal foot tissues was used for piRNA detection by northern blotting and in situ hybridization, piRNA0001 was widely expressed in all of these tissues, particularly the gill tissue (Fig. 6A). The in situ hybridization of small RNAs also revealed the expression of piRNAs in mantle tissue (Fig. 6B).

Discussions
In the present study, we examined the expression patterns of small RNAs in the mantle (Ma), adductor muscle (Ad), gill (Gi), and ovary (Go) tissues of pearl oysters. Our first experiment examined the expression patterns of miRNA species in pearl oysters. Surprisingly, we found that RNA species of approximately 30 nt were more abundant than those of approximately 22 nt, which is the length expected for miRNA and/or endogenous siRNA. The 30 nt RNA species were rich in uridine (U) residue at the 5′-end. In addition, β-elimination reactions confirmed the presence of 2′-O-methylation at the 3′ terminal of a representative 30-nt RNA. Since these are typical features of piRNA species, we regarded the 30 nt RNA species as putative piRNA in P. fucata. After removing reads that were of poor quality or not within the target size range, two obvious peaks in the distribution of lengths of small RNAs were observed in somatic tissues, and dominant 29-31 nt small RNAs observed in gonadal tissues. Four steps of sequence processing identified 18.0 million putative piRNAs and 2.8 million unique piRNAs. Abundant putative piRNAs were identified and characterized in pearl oyster somatic and gonadal tissues. Northern blotting and in situ hybridization revealed widespread expression of these RNA species in pearl oyster somatic tissues. Mapping of the putative piRNAs to the reference genome of the pearl oyster revealed that they were encoded in multiple loci, with tandem patterns on the genome. Also, their expression patterns were notably different in somatic and gonadal tissues. These results provide useful information for further research on piRNA in pearl oysters, and on the functions and molecular regulatory mechanisms of small RNAs in mollusks.  The main characteristics of most of the putative piRNAs-that they started with uridine at the 5′ terminal and were within the range 26-31 nt-are consistent with piRNAs in other animals [28][29][30] . A common feature of piRNAs in all systems analyzed thus far is that they derive from repetitive transposon rich regions of the genome 20,31 . In the piRNA biogenesis pathway, long piRNA precursors are transcribed from specific genomic loci known as piRNA clusters, cleaved and modified in the cytoplasm, and then transported into the nucleus in a complex with Piwi protein 32 . Due to the limited database of transposable elements in the pearl oyster, we performed a piRNA cluster analysis: we merged piRNA mapping loci within a distance of <1 kb and with a minimum coverage of six reads. A total of 35,848 piRNA clusters were identified on the reference genome, varying in length from 30 bp (single unique piRNA) to 60.9 kb, similar to piRNA clusters that have been identified in vertebrates and other invertebrates 29,33 . Twenty-five percent of the unique piRNAs mapped to multiple tandem loci on the scaffold, and comparisons between the dot-matrix analyses of piRNA clusters and piRNA genomic loci showed that most piRNAs were located on repeats on the scaffold. Previous report indicated that 54-65% of the piRNAs were mainly generated from repeat regions 22 . The result indicated that the biogenesis of piRNA in P. fucata is essentially common with those of other reported species. Our low repeat-mapping percentage may be the result of the imperfect pearl oyster reference genome, with 29,306 scaffolds 4 . Another important biogenesis region of piRNAs is the 3′-UTR regions 34,35 , but the relationship between piRNAs and the exonic regions of genes was beyond the scope of this study.
To gain better understanding of the putative piRNA expression pattern in the somatic and gonadal tissues, the RPKM of each piRNA cluster was used to compare piRNA expression among tissues. We found that 39.4% of piRNA clusters were commonly expressed in all tissues, and that 49.4% occurred in all somatic tissues. The PCA analysis of piRNA cluster density clustered the somatic tissues' libraries together, indicating that overall, the piRNA clusters expressed in somatic tissues were not substantially different. Moreover, piRNA clusters were expressed differently in the mantle, adductor muscle, and gill tissues, respectively, from the gonadal tissue. Thus expression patterns in the pearl oyster are different in the somatic and gonadal tissues.
The piRNAs tested from Drosophila, zebrafish, and mouse are resistant to periodate treatment, indicating the occurrence of a modified 2′-OH at the 3′ terminal 36 . This 2′-O-methylation at the 3′ terminal is recognized as a defining feature of piRNAs in all animals studied to date [37][38][39] . Consistent with these findings, a β-elimination reaction confirmed the presence of a 2′-OH modification at the piRNA 3′ terminal. Thus, it is likely that piRNAs from pearl oysters carry a 2′-O-methylation at their 3′ terminal and display similar characteristic features to the piRNAs of mammals and flies 20,40 . Lastly, northern blotting and in situ hybridization further verified the expression of piRNAs in the somatic tissues of pearl oysters.
In this study, millions of putative piRNAs revealed a fascinating and unanticipated dimension of pearl oyster biology. The piRNA pathway is commonly thought to be germline specific, even though the somatic function of Piwi proteins was documented when they were first discovered 41,42 . However, new studies are currently re-exploring the Piwi-piRNA pathway in the somatic cells of a diverse range of organisms, and notable insight has recently been gained into the somatic function of the pathway in the ovarian somatic cells of Drosophila 21 .
In addition, groundbreaking work in the lower eukaryotes has demonstrated a conserved function for Piwis and piRNAs in the somatic tissues, particularly in stem cells 32,43,44 . Piwi genes are expressed in the planarian totipotent stem cells that repopulated all somatic and germline lineages in planarians 45 . Furthermore, recent studies suggest that piRNAs are also associated with neuronal Piwi proteins and contribute to neuronal development and function 46,47 . Thus, it is evident that the diversity of functions of piRNAs is greater than previously thought. In the present study, putative piRNAs partially generated from repeats shared characteristic features with piRNAs identified in other organisms, and were strongly expressed in somatic and gonadal tissues, although with different expression patterns. However, the roles of piRNAs in the pearl oyster have remained unknown. They may play important roles beyond those of the germline in mollusks. Further studies are needed to determine whether piR-NAs play a critical role in epigenetic silencing through DNA/chromatin methylation, similar to their germline homologs in mollusks 48 .

Methods
Ethics. This study was conducted in strict accordance with the recommendations in the University of Tokyo's Guide for the Care and Use of Laboratory Animals. All efforts were made to minimize animal suffering.
Experimental samples and tissue collections. Two female pearl oysters (approximately two years of age, with body length and body weight of 13 ± 2 cm and 95 ± 10 g, respectively) were collected from a pearl oyster culture population and cultured at the Mikimoto Pearl Research Institution, Mie Prefecture, Japan. The pearl oysters were put in ice water for anesthetization and then dissected for tissue collection. Mantle, adductor muscle, gill and ovary tissues were collected and stored in 2 ml tubes for later separation of the RNA (Thermo Fisher, Japan). The samples were stored overnight at 4 °C, and then preserved at −80 °C until use.
Small-RNA extraction, library construction, and sequencing. Eight total RNA samples were extracted from the pearl oyster somatic and gonadal tissues. The subsequent small RNA extraction conformed to the protocol for the mirVana miRNA Isolation Kit (Life Technologies, America). Briefly, the sequencing libraries were constructed as per the Ion Total RNA-seq Kit v2 small-RNA library construction protocol (Life Technologies, America). Adapter ligation, synthesis of cDNA by reverse transcription, purification of the small-RNA fraction, and amplification by PCR were then performed to construct the cDNA libraries. Agilent 2200 Bioanalyzer (Agilent Technologies, Germany) was used for cDNA concentration checking. Using these cDNA libraries, a template for loading onto the analysis chip was prepared according to the protocol of the Ion PI Template OT 2200 Kit v3 (Life Technologies, America). Sequencing was then performed according to the protocol for the Ion Proton system (Ion Proton Sequencer and Ion Proton Torrent Server, Life Technologies, America).
Putatvie piRNA sequence processing. The raw sequencing data were saved as FASTQ files. The PRINSEQ program was used to exclude data with an average quality score of ≤20, and reads with lengths outside the 15-35-nt range were removed from FASTQ files ( Number of a unique piRNA reads Total number of piRNA reads from given library 6 was used to compare the expression of unique piRNAs between libraries. The density of unique piRNAs in each library was investigated. Normalized putative piRNA cluster identification. piRNA clusters were identified from putative piR-NAs using the following procedure: Firstly, all piRNA reads were mapped to the reference genome using "bwa mem". The piRNA genomic mapping loci were merged, and we defined a region with a distance of <1 kb between any piRNAs as a merged region (bedtools merge -d 1000). Total reads were then re-mapped to the merged region to calculate the coverage (bedtools coverage -a -b). Here, we defined merged regions with a minimum coverage of six reads as a piRNA cluster, and counted the normalized number of piRNAs in each cluster using the measure reads per kilobase per million mapped reads (RPKM) 50 . The normalized piRNA enrichment in each cluster was defined as follows: × × Number of reads mapped to a piRNA cluster 10 Total number of piRNA reads from given library piRNA cluster length in kb 6 Our definition of RPKM was slightly different from the conventional one, since piRNA precursor transcript information was not available, and our RPKM was normalized by genomic DNA length rather than transcript length. Finally, Dotter 51 and the Integrative Genomics Viewer (IGV) tracks tool 52 , respectively, were used to visualize the repeats and the mapping density on the scaffold. Differential expression of putative piRNA clusters in the somatic and gonadal tissues. The mapped reads were normalized by cluster length according to RPKM for each piRNA cluster in the somatic and gonadal tissues, which facilitated comparisons of piRNA levels between tissues. The edgeR package 53 was used to identify differently expressed piRNA clusters between two tissues. Clusters that were differently expressed in two tissues were identified using the following filter delimitations: FDR (false discovery rate) <0.01; and absolute value of log2Ratio > 1.  54,55 . The β-elimination reaction was used to identify methylation at the 3′ terminal of small RNAs. The detailed process of the β-elimination reaction has been described previously 56 . Briefly, an NaIO 4 solution was used to cleave the vicinal hydroxyls to a dialdehyde, and then the dialdehyde was treated with borate buffer at pH 9.5 (β-elimination reaction) in the dark at room temperature for 30 min. This resulted in a 3′-monophosphate of the miRNA that was shorter by 1 nt, which had a higher migration rate in gel electrophoresis than the miRNA did prior to the treatment 57 . Although the β-elimination reaction also resulted in methylation at the 2′-OH of the 3′ terminal of the piRNAs, these methylated piRNAs had the same migration rate in gel electrophoresis. The β-elimination reaction was used to investigate the modification of the 3′ terminal of small RNAs. We confirmed the most strongly expressed piRNA (piRNA0001: 5′-UACUUUAACAUGGCACAGAUAUAAUGACCU-3′) displayed 3′ terminal methylation, and we used miR-279a as a control. Four total RNAs from P. fucata gill tissues were used for small RNA β-elimination reaction analysis. To increase binding specificity and affinity to the targets, locked nucleic acid (LNA) probes were designed and used to detect small RNAs by northern blotting. Both terminals of the probes were labeled with digoxigenin. The process of northern blotting was conducted as previously described 58,59 , and the northern blotting probes are shown in Table 2.
Detection of putative piRNA in somatic tissues. The most strongly expressed putative piRNA (piRNA0001) in somatic tissues was selected for detection analysis in pearl oyster. Six pearl oysters (three of each sex) were collected from Mie Prefecture, Japan. After anesthesia with ice water, five somatic tissues (gill, mantle, adductor muscle, intestine, and abdominal foot) were dissected for sampling. The process of total RNA extraction was as described above. When sampling, part of the mantle tissue was stored in 4% paraformaldehyde at 4 °C for 12 hours, which was then replaced with 10%, 15%, and 20% sucrose for long-term preservation. We conducted the in situ hybridization of small RNAs using the IsHyb In Situ Hybridization Kit protocol (Biochain, America).
Accession codes. All sequencing data are deposited in the DNA Data Bank of Japan (DDBJ) 491 database under accession DRA006953.