Computational analysis of drought stress-associated miRNAs and miRNA co-regulation network in Physcomitrella patens.

miRNAs are non-coding small RNAs that involve diverse biological processes. Until now, little is known about their roles in plant drought resistance. Physcomitrella patens is highly tolerant to drought; however, it is not clear about the basic biology of the traits that contribute P. patens this important character. In this work, we discovered 16 drought stress-associated miRNA (DsAmR) families in P. patens through computational analysis. Due to the possible discrepancy of expression periods and tissue distributions between potential DsAmRs and their targeting genes, and the existence of false positive results in computational identification, the prediction results should be examined with further experimental validation. We also constructed an miRNA co-regulation network, and identified two network hubs, miR902a-5p and miR414, which may play important roles in regulating drought-resistance traits. We distributed our results through an online database named ppt-miRBase, which can be accessed at http://bioinfor.cnu.edu.cn/ppt_miRBase/index.php. Our methods in finding DsAmR and miRNA co-regulation network showed a new direction for identifying miRNA functions.

Introduction

Plant growth and productivity are severely affected by drought stress. A large number of genes that respond to drought stress at the transcriptional level have been identified 1, 2. The products of these genes play important roles not only in protecting cells from drought stress (, but also in regulating genes for signal transduction in the drought stress response 4, 5.

miRNAs are 21~24 nt non-coding small RNAs. In plant, mature miRNAs are generated from stem–loop regions of longer RNA precursors mainly by an endoribonuclease III-like enzyme, dicer like-1 (DCL1). The processed and methylated miRNA/miRNA* duplex is then exported to the cytosol. miRNAs that are incorporated into an argonaute protein containing RNA-induced silencing complex can affect the target gene expression 6, 7. miRNA-mediated regulations rely on specific complementary base pairing of miRNAs to target mRNA (, which can degrade the target transcript ( and/or attenuate translation (. Plant miRNAs bind to target mRNAs with near-perfect complementarity (. Based on this observation, several groups have developed computational methods to predict miRNA targets in Arabidopsis thaliana and rice 11, 12, 13.

miRNAs serve critical roles in plant development, including adaptation to biotic and abiotic stresses 6, 14, 15. It is now considered that, developmental processes are primarily regulated by RNA regulatory networks (, and miRNAs are key components of large gene regulatory networks 6, 17. However, our understanding of miRNA regulation network in plant is still poor.

Unlike A. thaliana, the model moss Physcomitrella patens occupies a key phylogenetic position with almost equal distance between green algae and flowering plants 18, 19. It is highly tolerant to drought ( and serves as a valuable source for identifying the mechanism for plant drought resistance. However, currently little is known about the molecular traits underlying the drought resistance. Therefore, identifying gene regulatory processes involved in the drought resistance of P. patens may ultimately help design genetically modified crops with drought tolerance.

In this work, we adopted a gene-oriented strategy to identify drought stress-associated miRNAs (DsAmRs) in P. patens and discovered their new relationships with gene regulations, using bioinformatics techniques. We share our results through an online P. patens miRNA database named ppt-miRBase to promote community development of drought stress response.

Results and Discussion

Potential DsAmRs in P. patens

By implementing the strategy described in Materials and Methods, we identified 6 miRNA families (including 6 miRNAs) that target 8 drought down-regulated genes (Table 1) and 10 miRNA families (including 19 miRNAs) that target 11 drought up-regulated genes (Table 2).

Table 1

miRNAs and their down-regulated gene targets under drought stress

miRNA	Target protein ID	Protein annotation
miR1039-5p	105818	Oxyanion transporter-like protein
miR1054	194139	Glutamate synthase
miR1028a-3p	204491	Chloroplast ribosomal protein L5
miR1050	226242	Chloroplast translation elongation factor EF-G
miR2083-5p	109367	Rubisco small subunit
miR2083-5p	205107	Rubisco small subunit
miR2083-5p	210883	Rubisco small subunit
miR902a-5p	235275	Glycine decarboxylase

Note: Target protein IDs were retrieved from JGI’s annotation (http://genome.jgi-psf.org/Phypa1_1/Phypa1_1.home.html).

Table 2

miRNAs and their up-regulated gene targets under drought stress

miRNA	Target protein ID	Protein annotation
miR1023a/b/c/d/e-5p	163734	No homology
miR1043-5p	176638	“SOUL” haem-binding protein
miR1035	177939	ATP-binding chaperone
miR1059	212537	Ubiquitin associated protein
miR1025	228672	Little protein 1
miR1221-3p	228678	WCOR413 cold acclimation protein
miR414	228695	Sucrose transporter
miR414	228697	Group 3 LEA
miR537a/b/c/d	228711	Organellar nucleoside-diphosphate-sugar epimerase
miR1031a/b	228736	Group 3 LEA
miR477f/g-5p	228741	19S proteasome subunit 9

Note: Target protein IDs were retrieved JGI’s annotation (http://genome.jgi-psf.org/Phypa1_1/Phypa1_1.home.html).

Most down-regulated target genes encode chloroplast proteins, especially photosynthesis-associated proteins (. For example, Ribulose-1,5-bisphosphate carboxylase/owygenase (Rubisco) small subunit is the target of miR2083-5p. Previous studies observed that during the drought stress, the activity of Rubisco decreased (. However, the underlying molecular mechanism was unclear. Our results suggest that miR2083-5p might inhibit the translation of Rubisco small subunit and consequently decrease the activity of Rubisco.

miRNAs may also up-regulate translation during stress 22, 23, 24 and cell cycle (, depending on how miRNA ribonucleoprotein complexes (miRNPs) are associated with other RNA-binding proteins 26, 27. Recently, a new study demonstrated that the miRNA target site may not be cleaved and instead negatively regulates miRNA activity through mimicry (. In Arabidopsis, miR398 targets two Cu/Zn-superoxide dismutase (CSD2 and CSD1) transcripts. During oxidative stress, miR398 is down-regulated and releases its suppression of CSD1 and CSD2 genes. The up-regulation of CSD1 and CSD2 genes also depends on expression levels of miR398 (. Our results raise the possibility that miRNAs might also be activators of certain genes in P. patens.

miRNA co-regulation network in P. patens

Biological systems are driven by complex biomolecular interactions. Multiple miRNAs could target the same gene 30, 31, 32, 33 and their target sites may partially overlap (, suggesting that these miRNAs could combine together to regulate the expression of miRNA target genes (. In this study, we identified 98 genes that were each targeted by more than one miRNA family.

General description

In the miRNA co-regulation network, each node represents a different miRNA family or the same miRNA family member from opposite arms of the miRNA stem–loop structure (for example, miR902-3p and miR902-5p). Two are connected, if they share at least one target gene (. The P. patens miRNA co-regulation network (Figure 1) contained 82 nodes and 98 edges.

Figure 1

miRNA co-regulation network in P. patens. Each pair of connected nodes indicates that they share at least one common target. A red node represents that at least one of its targets is up-regulated under drought stress; whereas a blue node indicates that at least one of its targets is down-regulated under drought stress. The node size is proportional to the number of neighboring nodes. The group of nodes connected with green lines represents that they form a module.

Targets of miR902-5p and miR414

The most intriguing feature of the P. patens miRNA co-regulation network was that miR902-5p and miR414 serve as two network hubs (nodes with a high number of edges connected), suggesting that they may play important roles in regulating drought resistance for P. patens.

In this study, both miR902-5p and miR414 were DsAmRs. miR902a-5p had 86 target genes. Gene Ontology (GO) term enrichment analysis (Table S1) indicated that targets of miR902a-5p are primarily involved in lipid metabolism. Target glycine decarboxylase (GDC) (Table 1) is down-regulated under drought stress. GDC is a key enzyme of photorespiration pathway (, which could protect plants against drought stress (.

miR414 had 145 predicted targets. They mainly involve in transportation of protein or sugar (Table S1). Sucrose transporter and group 3 late embryogenesis abundant (LEA) proteins are up-regulated under drought stress (Table 2). Sugar transporters have been recognized as pivotal targets for regulatory roles of distribution and partitioning of carbon resources in plants (. Sucrose metabolism is extremely responsive to internal and external signals and consequently alters development and stress adjustment (. LEA proteins are accumulated in plants during drought stress and have been considered to be associated with drought tolerance (.

Here we must mention, due to moderately poor precursor hairpin structure, low Northern blot signal, the lack of identified targets, and the lack of conservation in genomes other than Arabidopsis and rice, the status of miR414 sequence as a miRNA has been questioned. Therefore, it is still suspicious whether miR414 is a potential DsAmR in P. patens.

Modules in P. patens miRNA co-regulation network

Many cellular processes are performed not by individual molecules, but by groups of functionally associated molecules, often referred to as functional modules. In the P. patens miRNA co-regulation network, modules appear as groups of densely interconnected nodes. Some modules overlap with each other and form a network of their own.

We identified six modules in the P. patens miRNA co-regulation network (Figure 2). Among them, three modules overlap with each other and form a sub-network (Figure 2A). Target genes of this module are mainly involved in cellular component organization. No enriched GO terms were found in target genes of module B (Figure 2B); targets of module C (Figure 2C) regulated translational initiation; and targets of module D (Figure 2D) played roles mainly in DNA replication.

Figure 2

Modules in P. patens miRNA co-regulation network. Three modules share nodes to form a sub-network (A), their target genes mainly play roles in cellular component organization. No enriched GO terms were detected in target genes of module B. The functions of module C and module D target genes are translational initiation regulation and DNA replication, respectively.

ppt-miRBase database

To share our discoveries with the plant research community, we developed a P. patens microRNA database named ppt-miRBase, which could be accessed online from http://bioinfor.cnu.edu.cn/ppt_miRBase/index.php.

Here we show how to use the query capabilities supported by ppt-miRBase to explore the function of miRNAs included in the database (Figure S1). A user can start at the home page of ppt-miRNA, select a query miRNA, such as “miR902a-5p”, from the pull-down list of miRNA and click the “Query” button (Figure S1A). On the query result page (Figure S1B), the user could browse 86 predicted target genes and 26 co-regulation miRNAs for miR902a-5p, which is listed as a DsAmR. By clicking the linked numbers “86” and “26”, the user can also navigate into the detailed target page (Figure S1C) or the co-regulation miRNA page (Figure S1D). On the detailed target page (Figure S1C), the user can read that miR902a-5p had one drought-responsive target, protein 171132. According to GO annotation, protein 171132 involves in glycine metabolism.

Materials and Methods

Data preparation

P. patens EST sequences were retrieved from NCBI. Of all 381,669 EST sequences, 227,942 are 5′ ESTs, and 153,727 are 3′ ESTs. A total of 275 mature P. patens miRNA sequences were downloaded from miRBase (http://microrna.sanger.ac.uk/sequences/) 42, 43. P. patens genome sequence was retrieved from the Physcomitrella genome sequence assembly, version 1 (JGI Genome Portal; http://genome.jgi-psf.org/Phypa1_1/Phypa1_1.home.html). The genome encodes 35,938 proteins. The DNA sequences of drought stress responsive genes were downloaded from JGI according to Cuming’s annotations (. There are 82 down-regulated genes and 125 up-regulated genes (Tables S2 and S3).

Prediction of miRNA target

We modified the criteria that were developed for plant miRNA target prediction 8, 44, 45 to predict P. patens miRNAs targets. We mapped all 381,669 EST sequences onto P. patens genome using BLAT ( and clustered EST sequences into groups according to their genome locations before assembling EST sequences of each group into contigs using Phrap program (.

An initial pool of 275 P. patens miRNA candidate targets was created in these contigs using PatScan (, with at most four unpaired nucleotides. Using Perl scripts, we aligned all 275 P. patens mature miRNA sequences against the candidate targets. This method leads to identification of all potential target sites of a single miRNA on the same target sequence, including false positive target sites. Consequently, a filter was applied to these targets, and the parameters were used according to Schwab et al. (: no mismatch at positions 10 and 11, no more than one mismatch at positions 2-12, and no more than two consecutive mismatches downstream of position 13. Finally, a mispair scoring system was applied to these filtered targets according to Wu et al. (: mismatches and single-nucleotide bulges were each scored as 1, and G:U pairs were each scored as 0.5. Candidate targets with a mispair score less than 3 were selected as putative miRNA targets.

Discovery of DsAmRs

We adapt a target-guided strategy to discover DsAmRs in P. patens. Since many drought stress responsive genes had been identified, our starting point is the drought stress-induced genes instead of related miRNAs.

Microarray analysis by Cuming et al. revealed 207 differentially expressed genes in P. patens under drought stress (. After predicting targets of all 275 P. patens miRNAs, we mapped all the predicted target genes onto the P. patens genome. If a predicted target gene is located within the same genomic region, in which a known drought stress-induced gene co-exists, we trace the predicted target gene back to its miRNA, and consider its targeting miRNA as a potential DsAmR (Figure S2).

Construction of miRNA co-regulation network

Two miRNAs are co-regulation miRNAs if they share at least one target gene (. We constructed the miRNA co-regulation network using Cytoscape (. In the network, nodes represent different miRNA families or the same miRNA family member from opposite arms of the miRNA stem–loop structure, and links indicate the sharing of at least one target gene by two nodes. Network modules in the network were identified using CFinder (.

GO term enrichment analysis

To understand the function of an individual miRNA, we performed GO term enrichment analysis on target genes using topGO (. We implemented a weighted counting method and Fisher’s exact test to assess the significant of over-representation of GO terms. The enriched GO term (P<10−5) represented the function of the miRNA (Table S1).

Construction of ppt-miRBase

The ppt-miRBase was implemented on a HP Core 32 Linux computer using MySQL as the database management system. Apache web server and PHP code were implemented to generate data-driven dynamic web contents.

Authors’ contributions

PW and JW conceived the idea. PW, JW, YZ, JX, JF, WZ, SX and GJ carried out the experiments. PW and JYC supervised the work and wrote the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors have declared that no competing interests exist.