piRTarBase: a database of piRNA targeting sites and their roles in gene regulation.

PIWI-interacting RNAs (piRNAs) are a class of small noncoding RNAs that guard animal genomes against mutation by silencing transposons. In addition, recent studies have reported that piRNAs silence various endogenous genes. Tens of thousands of distinct piRNAs made in animals do not pair well to transposons and currently the functions and targets of piRNAs are largely unexplored. piRTarBase provides a user-friendly interface to access both predicted and experimentally identified piRNA targeting sites in Caenorhabditis elegans. The user can input genes of interest and retrieve a list of piRNA targeting sites on the input genes. Alternatively, the user can input a piRNA and retrieve a list of its mRNA targets. Additionally, piRTarBase integrates published mRNA and small RNA sequencing data, which will help users identify biologically relevant targeting events. Importantly, our analyses suggest that the piRNA sites found by both predictive and experimental approaches are more likely to exhibit silencing effects on their targets than each method alone. Taken together, piRTarBase offers an integrative platform that will help users to identify functional piRNA target sites by evaluating various information. piRTarBase is freely available for academic use at http://cosbi6.ee.ncku.edu.tw/piRTarBase/.

INTRODUCTION

Animals have evolved a robust mechanism to silence transposons and thereby guard their genomes against mutation using small noncoding RNAs called PIWI interacting RNAs (piRNAs) (1). piRNAs guide the PIWI subclass of Argonaute protein to search for mRNA targets and regulate their expression (2–5). Surprisingly, tens of thousands of piRNAs produced by animals, from nematodes to mice, do not match transposon mRNA sequences (6–8). Recently, piRNAs have been shown to regulate the expression of endogenous mRNAs. For example, it has been reported in mice that piRNAs produced during the pachytene stage of meiosis initiate the destruction of thousands of mRNAs to prepare spermatocytes for sperm production (9). In flies, piRNAs regulate the localization and degradation of mRNAs in the early embryo (10,11). Additionally, recent studies have also reported the function of piRNAs in the regulation of sex determination, axonal regeneration, and nictation in silkworms and nematodes (12–14). However, predicting transcriptome-wide piRNA targeting sites has been difficult due to limited knowledge of piRNA targeting.

Our group has recently determined the targeting rules that piRNAs use to find mRNA targets in Caenorhabditis elegans (15). We found that piRNAs, like microRNAs (16), possess a second to seventh nucleotide seed sequence that requires nearly perfect complementarity between the piRNA and its target. Though unlike most miRNA targeting events, piRNAs can only tolerate a few mismatches outside of the seed sequence. We have used the piRNA targeting rules to successfully express otherwise frequently silenced transgenes in the C. elegans germline by avoiding piRNA targeting and constructed a tool, pirScan, that allows other researchers to do the same (15,17). Additionally, the cross-linking, ligation, and sequencing of hybrids (CLASH) approach has recently been applied to the PIWI Argonaute in C. elegans, which has revealed endogenous piRNA-mRNA interactions in vivo (18). Both of these methods predict that many endogenous mRNAs can be targeted by piRNAs. However, both methods have been shown to identify many miRNA target sites that do not repress the expression of their targets, suggesting additional approaches are needed to reveal functional piRNA sites (19).

Our purpose in constructing piRTarBase (http://cosbi6.ee.ncku.edu.tw/piRTarBase/) is to synthesize the recent developments in piRNA targeting site identification from predictions and experiments into a user-friendly interface that will allow researchers to explore the piRNA targeting sites and their regulatory effects on endogenous genes (Figure 1). piRTarBase displays various information about these piRNA sites (18), allowing users to search for and browse desired endogenous genes for piRNA targeting sites, as well as to identify the predicted mRNA targets of desired piRNAs. Additionally, piRTarBase contains expression data from wildtype and PIWI mutant animals, allowing users to assess the effects that piRNA targeting has on a given gene. Importantly, our analysis suggests that piRNA targeting sites that are both predicted by our targeting rules and identified by CLASH, referred to as common sites, are significantly better at predicting mRNAs that are regulated by PIWI in vivo than either method alone. While numerous miRNA targeting site databases are available, there is currently only one database for piRNA targeting sites, and it only inventories piRNA sequences and the limited published targets (20). By contrast, piRTarBase integrates the results derived from the piRNA targeting rules, PIWI CLASH data, and expression data of mRNAs and small RNAs, allowing researchers of different fields to identify candidates of functional piRNA targeting sites.

Figure 1.

Workflow of piRTarBase.

DATABASE CONTENT

Database entries that pertain to piRNA target sites were collected from transcriptome-wide prediction of piRNA target sites using pirScan (17) and published PIWI PRG-1 CLASH data (18). Database entries corresponding to mRNA and small RNA expression profiles in WT and PIWI mutant animals were collected from published data (6,18,21–23). We will continue to update piRTarBase with new datasets as they become available.

pirScan predicted piRNA targeting sites

pirScan predicts transcriptome-wide piRNA targeting sites of 17849 C. elegans piRNAs (15364 type 1 and 2485 type 2 piRNAs) based on stringent and relaxed rules, resulting in 571204 and 1420256 piRNA targeting sites, respectively (17). The stringent rules are based on an in vivo reporter assay in C. elegans that reveals which sites, when mutated, are sufficient to avoid silencing (15). Therefore, predicted sites based on stringent rules are more likely to contribute to gene silencing. The piRNA targeting score is based on the data obtained in the same reporter assay (15,17). Alternatively, the relaxed rules are based on genome-wide monitoring of secondary siRNA (22G RNA) expression changes surrounding presumptive piRNA targeting sites (15). In C. elegans, when PIWI-piRNA complexes recognize a target gene, they induce the production of 22G RNAs by RNA-dependent RNA polymerases around the targeted site (22,24). The piRNA rules derived from this approach are slightly less stringent than those derived from the piRNA reporter assay and are therefore defined as the relaxed rules.

CLASH predicted piRNA targeting sites

The cross-linking, ligation, and sequencing of hybrids (CLASH) approach is a powerful method that can be used to identify RNA-RNA interactions in vivo (25). Using the CLASH method culminates in the sequencing of chimeric molecules representing a ligation event between the small RNA and a fragment of its mRNA target (26). Recently, CLASH was performed on the C. elegans PIWI Argonaute PRG-1, revealing in vivo interactions between piRNAs and their mRNA targets (18). However, the paper does not provide a list of sites identified by CLASH. Here we reanalyzed their published CLASH data to provide transcriptome-wide piRNA targeting sites. As CLASH data are known to be noisy and chimera reads can result from random RNA ligations, we restricted our analysis to chimeric reads that appear at least 5 times in the ligated libraries. We analyzed CLASH data by first finding a piRNA whose full mature sequence perfectly matches a portion of a chimera. To identify the mRNA target of that piRNA, we mapped the remainder of the chimera sequences that were over 15 nucleotides in length to the C. elegans transcriptome. As steps of CLASH experiments require RNase treatment, the chimeras may experience RNA degradation; therefore, to identify the best possible pairing between the piRNA and the target mRNA, we extended the length of the mapped mRNA region by 21 nucleotides from both the 5′ and 3′ ends of the read and predicted the piRNA targeting site as the interaction with the highest piRNA targeting score, similar to the approach taken in the previous study (25). piRTarBase displays the targeted mRNA region of each chimera, while the predicted pairing between the mRNA and the mapped piRNA is found in a separate table as described below. In total, piRTarBase inventories 10116 CLASH-identified piRNA target sites.

mRNA and small RNA sequencing data

In C. elegans piRNAs trigger silencing of their targets using amplified siRNAs known as 22G RNAs. Therefore, in prg-1 PIWI mutants, which lose all piRNAs, piRNA targets have been shown to have increased levels of mRNAs and reduced levels of 22G RNAs (22,24). In order to allow users to better evaluate how piRNA targeting affects a particular gene, piRTarBase collects mRNA expression data from wildtype (N2) and PIWI mutant (prg-1) animals. Normalized mRNA expression data is displayed for each gene in wildtype and in PIWI mutant animals. piRTarBase also displays the mRNA fold change (averaged between replicates) in PIWI mutant vs. wildtype animals for each gene. Therefore, genes with mRNA fold changes greater than 1 are more likely to be regulated by piRNA targeting events. For example, it has recently been shown that xol-1, a master regulatory gene controlling dosage compensation and sex determination in C. elegans, is regulated by piRNAs (27). Looking at the piRTarBase entry for xol-1, a user can see that the mRNA fold change in PIWI mutants is over 2, meaning that xol-1 is more than twice as highly expressed in PIWI mutants. This confirms that xol-1 is likely to be regulated by piRNAs.

piRTarBase also collects small RNA expression data from wildtype and PIWI mutant animals. ‘22G RNA seq data’ displayed in piRTarBase reflects the accumulation of 22G RNAs on a given mRNA species. This data is also available in wildtype and PIWI mutant animals. The calculated fold change compares the accumulation in 22G RNAs in PIWI mutants to that in wildtype animals. Because successful piRNA targeting events lead to the production of 22G RNAs (22,24), a reduction in 22G RNAs along a given gene in a PIWI mutant would be expected if that gene was being regulated by piRNAs. Therefore, a fold change less than 1 would suggest possible piRNA regulation of that gene. Taking xol-1 as an example again, it can be seen that the 22G RNA fold change is 0.0215, which shows that 22G-RNAs are dramatically reduced in PIWI mutants and suggests that xol-1 is being regulated by piRNAs.

Common piRNA targeting sites

pirScan's prediction algorithm is based on a robust reporter assay that allowed our group to optimize transgene sequences for germline expression in C.elegans (15). While the sites predicted using pirScan robustly identifies sites that contribute to silencing in vivo, it also considers piRNAs that do not actually interact with their predicted mRNA targets. For example, it may be that some piRNAs are expressed at low levels, or a piRNA and its presumptive target are not expressed at the same developmental time, or that mRNAs are protected for piRNA interaction by anti-silencing signals such as periodic An/Tn repeats (PATCs) or CSR-1 Argonaute interaction (15,18,28). Conversely, the targets predicted based on CLASH data may identify targets that interact with mRNA targets but do not actually contribute to silencing, as shown in the case of miRNAs (19). Therefore, we examined whether the piRNA targeting sites identified by both methods (904 sites and 2273 sites using stringent and relaxed rules, respectively), referred to as common piRNA target sites, are more likely to represent targeting events that are sufficient for regulating mRNA targets. Indeed, for reported targets xol-1 and Y40B10A.2a, removing only a single common piRNA targeting site displayed in piRTarBase (21ur-4863 and 21ur-8264 targeting sites, respectively) in each of these genes was sufficient to increase these genes’ expression levels (18). Furthermore, overall mRNA expression is significantly higher in PIWI mutants for genes whose targets were predicted using both pirScan and CLASH, but not for genes identified using each method separately (Figure 2A). Additionally, when we compare the 22G RNA fold changes in genes whose targets were predicted by pirScan, CLASH, or by both independently, we see that overall 22G RNAs are reduced significantly in PIWI mutants for target genes identified using both methods, but not for target genes identified by the individual methods (Figure 2B). Taken together, these observations suggest that the common piRNA target sites provided by piRTarBase are more likely to have regulatory effect on their targets in vivo than sites that are identified only by prediction or experimental approach alone.

Figure 2.

22G RNA and mRNA fold expression changes in germline expressed genes with piRNA targeting sites predicted using different methods. Only genes found to be expressed in the C.elegans germline according to (30) were considered. Distributions are displayed without outliers for clarity. Outliers were defined as points less than Q1–1.5*IQR or greater than Q3+1.5*IQR. Statistics were performed using the Mann–Whitney U test. Comparisons were made between each distribution and the full germline expression gene set. *** indicates P < 1e–9, ** indicates P < 1e–3, * indicates P < 0.05, n.s. indicates not significant. (A) mRNA fold expression changes, from (21), in genes with predicted piRNA targeting sites. (B) 22G RNA fold expression changes, from (22), in genes with predicted piRNA targeting sites. Note that prg-1 mutant animals have a reduced number of germ cells compared to wildtype, which likely contributes to the reduction of overall levels of 22G RNAs and mRNAs among all germline expressed genes.

Database statistics and data analysis

piRTarBase curates target sites that correspond to the WS230 C. elegans genome assembly and the WS230 C. briggsae assembly. The number of target sites compiled in piRTarBase are displayed in Table 1. Expression data was taken from (6,18,21–23,27).

Table 1.

Database statistics

Method	Number of target sites	Number of genes	Number of transcripts
Stringent rules	571 204	19 781	32 825
Relaxed rules	1 420 256	20 097	33 288
CLASH	10 116	7880	16 076
Stringent rules & CLASH	904	903	1898
Relaxed rules & CLASH	2273	2272	4876

The web interface of piRTarBase was constructed using the Python language with the Django MTV framework. The figures were generated by D3.js, which is a JavaScript library for manipulating documents based on data.

Comparison with other piRNA databases

Currently only one database (piRBase) of piRNA targeting sites is established. This piRNA database focuses on curating a list of piRNAs and their targets (20). While such databases successfully condense the existing literature into a more easily searchable format, they lack the predictive power of piRTarBase. piRTarBase uniquely provides candidates of piRNA target sites on endogenous transcripts using pirScan's prediction algorithm and PIWI CLASH data. This allows piRTarBase to collect piRNA targeting sites genome-wide at a much higher confidence level than existing databases. In addition, piRTarBase offers a better integrative platform that allows users to examine various information about piRNA targeting sites, including expression level of mRNA targets, piRNA and 22G-RNA. In addition, users can easily examine the position of mismatches, pairing score, abundance of chimera reads/piRNAs for each piRNA targeting site with sortable tables. While piRTarBase does not inventory the limited piRNA target sites reported in mice or flies, we plan to expand piRTarBase to other organisms once the corresponding piRNA targeting rules and PIWI CLASH data are available. Taken together, we believe that piRTarBase offers a user-friendly interface for researchers of diverse fields to explore the potential regulatory function of piRNAs.

WEB INTERFACE

Searching by gene

piRTarBase supports querying either by gene or by piRNA (Figure 3, 1. Input). The user can select the desired species (currently C.elegans or C.briggsae) and the desired stringency of the targeting rules. For gene/transcript querying, clicking search following input of the gene/transcript name results in a table displaying the queried terms, a table displaying mRNA and small RNA expression data, and a table displaying isoform-specific piRNA targeting site information (Figure 3, 2a. Targeting piRNAs). In the mRNA expression data, expression levels for the input gene are displayed in reads per kilobase million (RPKM) for wildtype (N2) animals and PIWI (prg-1) mutant animals, with two biological replicates for each strain. This information is summarized for the user in a mRNA fold change field, which calculates the average mRNA fold change in PIWI mutants relative to wildtype animals. The user can use this field to determine whether loss of PIWI results in an elevated level of expression (fold change > 1) or a reduced level of expression (fold change < 1). Additionally, 22G RNA sequencing data is available in reads per million (RPM) for wildtype and PIWI mutant animals. The 22G RNA fold change is displayed in a field, which will allow the user to determine whether loss of PIWI results in a reduction of 22G RNAs mapping to that gene (fold change < 1) or an accumulation of 22G RNAs mapping to that gene (fold change >1). As of the publication of this manuscript, piRTarBase curates all four publicly available 22G RNA sequencing datasets in C. elegans (6,18,22,23) and the single available 22G RNA sequencing dataset in C. briggsae (27). The last table on the page displays the piRNA targeting sites for each mRNA isoform of the queried gene. Targeting sites are grouped into three categories: (i) Predicted piRNA target sites: targeting sites predicted by pirScan's algorithm, (ii) CLASH identified piRNA target sites: targeting sites predicted by CLASH data (iii) Common piRNA target sites: targeting sites that were predicted by both pirScan and CLASH data. Additional information about each identified targeting site can be found by clicking ‘Show target sites’ (Figure 3, 3a. Target Site Data). On this page, piRNA targeting sites are first displayed as a graphic, with pirScan predicted target sites in red, CLASH identified target sites in blue, and common piRNA target sites in yellow. Tables corresponding to each type of identified targeting site follow. Hovering over any displayed targeting site gives the user additional information about that site, including which piRNA is being displayed, the expression level of that piRNA, the number of hybrids/chimeras sequenced (if the piRNA was predicted by CLASH), the region targeted in the mRNA identified in the chimera and the predicted pairing region, the piRNA targeting score, and the mismatches and pairing between the piRNA and the targeted site. This information is additionally displayed in the tables that follow the graphic. By clicking the sorting icons, the order of piRNA sites displayed can be sorted by the values in various columns, such as piRNA abundance or piRNA targeting score.

Figure 3.

An overview of piRTarBase's user interface.

Searching by piRNA

For querying by piRNA, the user can input a piRNA name and click search. The average piRNA expression data displayed in this case corresponds to the expression of that piRNA as determined by PIWI immunoprecipitation (PRG-1 IP) small RNA sequencing studies (Figure 3, 2b. Targeted mRNAs). All mRNA targets are displayed in a single table, with each row containing all of the information outlined above for querying by gene, including mRNA and small RNA expression fold changes and numbers of target sites predicted for that particular piRNA on a given gene. Users can click the sorting icons to display the mRNA targets in ascending or descending order of the values in various columns, such as mRNA or 22G RNA fold change. To find more information about the piRNA–mRNA information displayed in each row, the user can click ‘Show target sites’ to open an analogous page as described for querying by gene, except only targeting sites involving the queried piRNA are displayed (Figure 3, 3b. Target Site Data).

Browsing by genes or piRNAs

All transcripts or piRNAs can be displayed in a single sortable table. Browsing by transcript will allow users to find the mRNA and 22G RNA fold change information and the piRNA targeting site information, and to access the target site specific pages for each transcript. Browsing by piRNA will allow users to find expression information for each piRNA and its corresponding mRNA targets, and to access the target site specific pages for each target of the queried piRNA. Rows can be sorted by any field, allowing users to identify which mRNA transcripts contain the most target sites, or whose expression is most affected by loss of PIWI for example.

Downloading piRTarBase content

All data contained in piRTarBase is available to download for advanced analyses.

SUMMARY AND PROSPECTIVE WORKS

piRTarbase provides the first integrative database allowing users to access both predicted and experimentally identified piRNA targeting sites and to explore their regulatory function. Importantly, our analyses demonstrate that the common piRNA sites identified by both bioinformatic and experimental approaches are more likely to have regulatory effects on their targets than those identified by the individual approaches. The user can further examine the change in mRNA and 22G-RNA expression of each specific target gene to assess the likelihood that piRNAs might regulate that target. Together, such an integrative platform uniquely allows users to sort through the hundreds of thousands of possible piRNA target sites and focus hypotheses on those sites that are more likely to be functional in vivo.

Currently, C. elegans is the only species whose piRNA targeting rules have been elucidated and the only species where PIWI CLASH has been performed. As prediction algorithms or CLASH data become available in other organisms, we will expand piRTarBase to include piRNA target sites for those organisms as well. For example, since the C. elegans piRNA targeting rules are likely conserved in C. briggsae due to the high similarity between their PIWIs, we have expanded the piRTarBase to include the piRNA target sites in C. briggsae.

We believe that piRTarBase provides a user friendly and informative interface for researchers of diverse backgrounds to investigate whether genes of interest are regulated by piRNAs or vice versa. For example, recent quantitative trait locus (QTL) analyses in C. elegans have mapped the nictation behavior-controlling and other gene regulatory loci to regions inside the piRNA clusters on chromosome IV (13,29). piRTarBase will allow researchers to identify candidates for the mRNA targets of the piRNAs located within these QTLs that may be contributing to the phenotype. Similarly, the authors of a recent study that linked piRNA regulation to axonal growth in C. elegans (12) will benefit from the ability to search through genes critical for neuronal regeneration and ask whether those genes may be targeted by specific piRNAs that underlie the regeneration defects. Such selective analyses were not possible prior to piRTarBase.