PceRBase: a database of plant competing endogenous RNA.

Competition for microRNA (miRNA) binding between RNA molecules has emerged as a novel mechanism for the regulation of eukaryotic gene expression. Competing endogenous RNA (ceRNA) can act as decoys for miRNA binding, thereby forming a ceRNA network by regulating the abundance of other RNA transcripts which share the same or similar microRNA response elements. Although this type of RNA cross talk was first described in Arabidopsis, and was subsequently shown to be active in animal models, there is no database collecting potential ceRNA data for plants. We have developed a Plant ceRNA database (PceRBase, http://bis.zju.edu.cn/pcernadb/index.jsp) which contains potential ceRNA target-target, and ceRNA target-mimic pairs from 26 plant species. For example, in Arabidopsis lyrata, 311 candidate ceRNAs are identified which could affect 2646 target-miRNA-target interactions. Predicted pairing structure between miRNAs and their target mRNA transcripts, expression levels of ceRNA pairs and associated GO annotations are also stored in the database. A web interface provides convenient browsing and searching for specific genes of interest. Tools are available for the visualization and enrichment analysis of genes in the ceRNA networks. Moreover, users can use PceRBase to predict novel competing mimic-target and target-target interactions from their own data.

INTRODUCTION

Non-coding RNAs (ncRNAs) are emerging as important regulators of gene expression. Their length has been used to subdivide them into two major classes namely (i) small noncoding RNAs and (ii) long non-coding RNAs (lncRNAs) (1). MicroRNAs are important small noncoding RNAs which play a key role in the regulation of post-transcriptional pathways in development (2). For example, founder-member miRNAs, lin-4 and let-7 regulate the proper timing of stem cell fate decisions in the nematode Caenorhabditis elegans (3–5), while miR156 regulates developmental timing in Arabidopsis by repressing the expression of functionally distinct squamosa promoter binding protein-like (SPL) transcription factors (6).

MiRNAs regulate the stability of mRNA transcripts by pairing to complementary sequences known as miRNA response elements (MREs) (7). In plants, perfect base-pairing between a miRNA and a mRNA transcript leads to cleavage and degradation of that mRNA. Most miRNAs target multiple mRNA transcripts, while individual transcripts may also contain multiple MREs for either the same or different miRNAs (8). The abundance of mRNA transcripts can also be indirectly regulated by the concentration of other mRNA transcripts if there is competition for binding to the same miRNAs. There is thus great potential for miRNAs to be part of an intricate network which regulates gene expression.

The amount of free miRNAs in a cell is dependent on their expression levels, the concentration of their targets, and whether they bind to these targets with or without base mismatches. Imperfect pairing does not lead to mRNA transcript degradation, and thus such RNA molecules can act as miRNA sponges (or decoys) to mop up free miRNAs in the cell. For example, in Arabidopsis, the lncRNA IPS1and mRNA PHO2 compete to bind miR399. miR399 can bind to PHO2 and induce its cleavage, but can also bind to IPS1 with a central three-nucleotide bulge (a hallmark of miRNA target mimics in plants) (9–11). The IPS1 lncRNAs are not cleaved, and instead serve to sequester miR399, thereby preventing it from binding to PHO2 mRNA (9) and allowing production of PHO2 protein under phosphate-replete conditions. IPS1 and PHO2 are thus classical examples of ‘competing endogenous RNAs’ (1,12–14) or target-mimic ceRNAs. As the knowledge of the transcriptome space is increasing, it is becoming evident that a large number of MREs exist in a wide variety of RNA transcripts, including mRNAs, lncRNAs, pseudogenes and transposable elements (9). In this work we focus on identifying ceRNAs, their target miRNAs and the potential regulatory networks they form.

Recently, several databases dedicated to the prediction and curation of ceRNAs have been developed. CeRDB (15) stores information about potential MRE-containing mRNAs. In ceRDB, the ceRNA pairs are listed according to a score based on the number of shared MREs. However, it is evident that not only mRNA transcripts but also many lncRNA transcripts act as ceRNAs (9,14,16). LnCeDB (17) comprises a dataset of human lncRNAs (from GENCODE) that potentially act as ceRNAs. Unlike ceRDB, which mainly contains putative predicted miRNA–mRNA interactions, LceDB provides some AGO-CLIP supported miRNA-mRNA/lncRNA pairing interactions. Both of these databases provide relative expression levels of ceRNAs facilitating user evaluation of the potential ceRNA influence (15,17). StarBase v.2.0 (18) is a comprehensive RNA interaction network, including CLIP-seq verified ceRNA interaction networks.

Although the first ceRNA pair was discovered in Arabidopsis (9), the above ceRNA-related databases are limited to animal species; no plant-specific ceRNA database has been developed to date. In this paper, we present a database of miRNA associated plant competing endogenous RNA interactions (Plant ceRNA database or PceRBase) (Figure 1). PceRBase is designed to provide the plant research community with easy access to a large amount of resources regarding candidate ceRNA pairs in order to build ceRNA networks, and inform future experimental work in this area. In PceRBase, two types of potential RNA interaction between each pair of RNA transcripts are considered for a ceRNA relationship: (i) ‘target-target’, where the common miRNA binds nearly perfectly to both transcripts (7,19); or (ii) ‘target-mimic’, where a bulge exists in the middle of the corresponding miRNA, so that only its two ends can bind to the mimic transcripts (9,10,20,21). The database currently stores predicted ceRNAs from 26 plant species. The biological importance of these candidate ceRNAs can be further evaluated by considering the overlap in their associated GO annotations, and whether they are co-expressed in particular tissues under the same conditions. Furthermore, a web-tool is provided in PceRBase allowing users to predict potential ceRNA pairs from their own sequence data.

Figure 1.

Overview of PceRBase core framework. (A) Detection of miRNA targets. (B) Prediction of ceRNA pairs. (C) Features of PceRBase, which integrates various data to evaluate the predicted ceRNA pairs. C(i) miRNA base pairing to ceRNAs. C(ii) Relative expression levels of ceRNA pairs across RNASeq samples. C(iii) GO annotations of ceRNA pairs. C(iv) GO term enrichment of gene set of interest. C(v) Network formed by ceRNA pairs.

MATERIALS AND METHODS

Data collection

RNA transcript information, including sequence, genome location, protein-coding potential and length were retrieved for 24 plant species from Phytozome10 (22) (http://phytozome.jgi.doe.gov/pz/portal.html) and from TAIR10 (23) (http://www.arabidopsis.org/) and MSU RGCP7 (24) (http://rice.plantbiology.msu.edu/index.shtml) for Arabidopsis thaliana and Orysa sativa respectively. With the exception of A. thaliana and O. sativa, gene ontology annotations for all transcripts were derived from Blast2GO (25). The GO annotation information for Arabidopsis and rice were retrieved from TAIR10 (23) and MSU RGCP7 (24).

Known miRNAs for all of these plant species were retrieved from miRBase (26). For A. thaliana, an additional 6480 lincRNAs identified by Liu et al. (27) were included. Similarly for Zea mays, 2492, 664, 11105, 1703 lncRNAs identified by Boerner et al. (28), Zhang et al. (29), Wang et al. (30) and Li et al. (31) respectively were also included in our database. All the data used in this study are summarized in Supplementary Table S1.

Prediction of miRNA and lncRNA

We predicted additional novel A. thaliana and O. sativa miRNAs present in smRNA-sequencing data sets (Supplementary Table S1), using a module from Mtide (32), which is suitable for plant-specific miRNA prediction.

Novel lncRNAs in A. thaliana, O. sativa, Sorghum bicolor and Z. mays were predicted from RNA-sequencing datasets that were obtained from SRA (http://www.ncbi.nlm.nih.gov/sra, Supplementary Table S1). This prediction work was executed in three steps. Firstly, the reads, filtered to exclude low quality data, were aligned to the respective plant genomes using Tophat (33). Secondly, the transcripts were assembled according to the alignment results using Cufflinks (34), and these assembled transcripts were then merged together using Cuffmerge. The merged transcripts provide a uniform basis for calculating transcript expression in each data set and the transcript class code. The transcripts with class code of ‘u’,’x’ or ‘i’ were selected as lncRNA candidates. Thirdly, CNCI (35) was used to assess the coding potential of the lncRNA candidates. An lncRNA candidate was retained if it satisfied the three following conditions: (i) the length had to be more than 200 nucleotides; (ii) the maximum ORF length had to be no longer than 120 amino acids; (iii) the coding probability had to be noncoding.

MicroRNA target and mimic scanning

The sequences of mRNA transcripts for 26 plant species and lncRNAs for four plant species served as the target database for the following two series of predictions:

to identify targets perfectly bound by miRNA, the targets of all known and novel miRNAs were predicted using TargetFinder (36), using a score parameter less than three.

to identify potential target mimics, the improved Tapir software (37) was used with the following three rules (a) the bulge in the complementary site should be composed only of three sequential nucleotides and located in the middle of the corresponding miRNAs. The middle position was defined as ninth to 10th or 10th to 11th or 11th to 12th. (b) Perfect pairing was required from the second position to the bulge starting position including G/U pairs at the miRNA 5′ end. (c) Except for the central bulge, the total mismatched and G/U pairs within mimic target and miRNA pairing regions should be no more than three and the consecutive mismatches should not exceed two (20,21).

Prediction of competing endogenous transcripts

To identify the ceRNA ‘X’ for a transcript ‘M’ of interest the following procedure is used: (i) identify all miRNAs targeting transcript ‘M’ from the results of TargetFinder, (ii) identify all miRNAs that could bind transcript ‘X’ from the results of TargetFinder (perfect binding) and Tapir (potential mimics), (iii) If one or more miRNA could bind to both ‘X’ and ‘M’, we defined transcript ‘X’ as a ceRNA of transcript ‘M’ and they thus represent a ceRNA pair. If the shared miRNA can bind nearly perfectly to ‘X’, this ceRNA pair is defined as ‘target–target’, while if the shared miRNA can instead bind imperfectly (as defined above) to transcript ‘X’, the ceRNA pair is defined as ‘target-mimic’.

To assess the reliability of the ceRNA pair prediction, we generated three scores (details are provided in Supplementary File S1). Firstly, we have calculated the mutually targeted MRE enrichment (MuTaME) score (38), which evaluates a candidate ceRNA ‘X’ for a transcript ‘M’ of interest. Secondly, to calculate the probability that a transcript is a target ceRNA or mimic ceRNA, a hypergeometric test is executed for each pair separately. Thirdly, the simultaneous expression of the transcript and its ceRNAs in a sample is an important indicator of their potential biological relevance. To evaluate the co-occurrence of transcript and ceRNA pair, we have computed an index of co-expression (ICE) (39) between them.

We considered a pair as co-expressed if the ICE ≥ 0.5. In this database, the ICE parameter was calculated only for the ceRNA pairs for Arabidopsis, rice, sorghum and maize, because sufficient experimental datasets were only available for these species (66, 72, 22 and 21 RNA-Seq datasets respectively).

DATABASE CONTENTS

ceRNA information

Potential significant ceRNA target–target and target–mimic pairs were predicted for the 26 plant species (Supplementary Table S2). From the browse page, users can visualize a list of all ceRNA pairs in a selected plant species. Alternatively, users can get the ceRNA list for a specific transcript by using the ‘Search’ module, as described in the following section.

The ceRNA information page is divided in five main sections: (i) ceRNA overview, (ii) transcript expression profile, (iii) transcript information, (iv) miRNA target site information and (v) GO annotation. The ceRNA overview section contains a summary of the relationship between two transcripts, including the species, the MuTaME score, the P-value, the total number of shared miRNAs and validation type. Validation categories are as follows: computational prediction (only in single plant species), orthologous (computational prediction of ceRNA pair in more than one species) and pubmed relationship (previously described in literatures). In the transcript expression profile section, the graphs illustrate the expression profiles for the two transcripts across multiple RNA-Seq experiments. This section also contains the ICE score and Pearson correlation coefficient of expression. The transcript information section shows detailed information for both transcripts, including their genetic loci and nucleotide sequences. The miRNA target site information section illustrates the potential miRNA and target base-pairing structures. The GO annotation section provides useful information for evaluating the function of ceRNA pair, because these pairs may be expected to have similar ‘Molecular Function’ or to be involved in related ‘Biological Process’.

Search modules

PceRBase provides two query interfaces for retrieving potential ceRNAs: ‘Simple Search’ and ‘BLAST Search’. The ‘Simple Search’ module includes two modes: transcript search and miRNA search. Once users obtain the transcript- or miRNA-related ceRNA pairs, they can filter the ceRNA list by adjusting the ‘Filter opions’, such as P-value and number of shared miRNAs. In addition, users can perform a Blast sequence search to retrieve ceRNA data in the ‘BLAST search’ module, which facilitates for the identification of novel and/or orthologous transcripts and their ceRNAs.

GO functional analysis

To help to understand the potential biological function of candidate ceRNA pairs, we provide the GO annotation for most transcripts in our database. The ceRNA pairs are expected to have similar functions, so users can identify potential novel functions of the second transcript according to the GO term assigned to the other transcript of the pair. Additionally, a combination of the Chi-Square test and Fisher's exact test is used to identify enrichment of GO term in user defined gene sets.

ceRNA regulatory networks

To help the users construct and display a complex miRNA-mediated ceRNA regulatory network of interest, a graphical browser is available by using the option ‘My network’. It has been developed on the basis of Cytoscape Web program (40). Users can choose one or more ceRNA pairs of interest from the ‘browse’ page or from the search result page, then click on the button ‘Add to My Network’ and eventually click on the button ‘Open My Network’. Then a new explorer window will present the network, where the nodes are coloured according to their types.

A web-tool for user driven ceRNA prediction

To enable users to assemble and analyse miRNA-mediated ceRNA networks, we have integrated the workflow we have developed for the identification of ceRNA pairs into the web interface. Though the input form, users can paste the miRNA and transcript sequences in fasta format or upload their miRNA and transcript sequences file in fasta-formatted files. After clicking on ‘GO’, the users will receive both cleaved and mimic ceRNA pairs predictions as a result web page, or in an email, if they have provided their address in the input form.

An example—identification of a potential novel ceRNA in Arabidopsis

SPL proteins are transcription factors that play important roles in plant development, notably the transition from juvenile to adult vegetative development and floral induction (41). The Arabidopsis genome contains 16 SPL genes, 10 of which contain an MRE for the highly similar miRNAs miR156 and miR157, which differ at only two nucleotide positions (42). A functional role for these miRNAs has been demonstrated, with over-expression of miR156 resulting in decreased SPL mRNA levels and delayed flowering in Arabidopsis (43). In Arabidopsis miR156 levels decline as the plant ages, thereby allowing miR156-targeted SPL mRNAs to increase in abundance. While the cause of this decline is not clear, it has been suggested that sugar-mediated repression of miR156 expression may play a role (41). An alternative hypothesis is that age-controlled expression of a ceRNA which can compete for miR156 and miR157 binding could also lead to an increase in SPL mRNA levels thereby allowing these developmental transitions to occur. Interrogation of PceRBase using SPL6 (At1g69170.1) as the search input identified the lncRNA At1NC083510 as a potential ceRNA for miR157 binding to the MRE of this and other SPL mRNAs in Arabidopsis. This lncRNA is predicted to bind to miR157 with a three nucleotide bulge (and could also bind to miR156 with a four nucleotide bulge). Measurement of miR156/7-targeted SPL mRNA levels in At1NC083510 over-expressor and knock-out plants, coupled with analysis of the timing of the vegetative and floral phase transitions in these lines would determine whether this potential ceRNA is indeed functionally relevant.

CONCLUSIONS AND FURTHER DIRECTIONS

Recent reports have described the microRNA-associated interplay among diverse RNA transcripts, and its crucial importance in the development and pathogenesis of eukaryotes (9,19,21,38,44,45). We have developed PceRBase, a collection of predicted plant ceRNA pairs from 26 plant species. This database not only provides convenient browse and search functions, and links to expression and GO enrichment analysis, but also contains a web-tool to predict significant ceRNA pairs from datasets provided by users.

In order to make the database as comprehensive as possible, a large number of smRNA-sequencing and RNA-sequencing datasets have been analyzed, to predict novel miRNAs and lncRNAs. A large number of novel miRNAs were predicted for Arabidopsis (1176) and rice (1822). Furthermore, we have identified a number of novel lncRNAs for Arabidopsis (987), rice (4808), sorghum (3534) and maize (3056) not previously reported. We predict that many more ceRNA pairs will be identified as more sequence information on miRNAs and lncRNAs in different plant species becomes available. For example, although sequence information is included for Citrus clementime mRNA transcripts in PceRBase, only five miRNAs have been identified to date in this species, limiting the power of PceRBase to identify significant C. clementime ceRNAs. Similarly, our analysis of lncRNAs acting as potential ceRNAs is currently limited to 4 plants: Arabidopsis, rice, sorghum and maize.

We will update PceRBase regularly by providing more ceRNA pairs, as well as related miRNA and lncRNA in additional plant species, as more datasets become available, to expand the database. However, a major shortcoming of PceRBase is the lack of plant CLIP-seq (cross-linking immunoprecipitation sequencing) datasets to validate the ceRNA predictions. In future, we plan to collect and integrate experimentally validated ceRNA pairs, to test the predictive value of the PceRBase platform.