Identification and comparative analysis of drought-associated microRNAs in two cowpea genotypes

Identification of miRNAs in cowpea

In order to study the role of miRNAs in drought tolerance, we grew cowpea plants (CB46 and IT93K503-1) in green house under well-watered and drought stress conditions. Drought stress was applied to 30-day-old plants. After 10 to 15 days of moderate drought stress (ψ_w = -1.5 MPa), the two genotypes showed apparent differences in drought tolerance. While IT93K503-1 plants continued to grow relatively well, CB46 plants displayed severe drought stress symptoms such as chlorotic leaves (Figure 1).

We constructed four small RNA libraries (2 genotypes × 2 growth conditions) and obtained on average 6.9 million (range: 6.5 - 7.3 million) clean small RNA reads from each library (deep-sequencing data have been deposited into the NCBI/GEO database with accession number GSE26402). The average number of unique reads per library is 4.3 million (range: 3.9 - 4.6 million). Using the procedure and criteria described in the materials and methods section, we mapped unique small RNA reads to a cowpea EST assembly, BAC end sequences and methylation filtration sequences, GSS sequences in dbGSS, and a draft cowpea genome assembly, and predicted 14, 78, 6, and 125 miRNA precursors, respectively. These four sets of putative miRNA precursors were then compared with each other to remove redundancy, and we obtained 157 candidate miRNA genes (for detailed information, see Additional file 1). Based on similarity of mature miRNA sequences, these miRNA genes were clustered into 89 families. Whereas 27 families (93 miRNAs) have match to miRNAs from other plants in the miRBase (release 16) [37], 62 families (64 miRNAs) appear to be cowpea-specific. Using a cowpea EST assembly, we have also identified putative target protein-coding genes for 112 (71%) miRNAs.

Genotype-specific expression of miRNAs

Because small RNA libraries were sequenced to great depth, counts of mature miRNAs can be used to evaluate their relative expression levels in different genotypes and growth conditions. We first applied Principal Component Analysis (PCA) to the log2 normalized counts (transcripts per ten million, TPTM) of 91 unique mature miRNAs that had combined expression of at least 50 TPTM in four libraries. As shown in Figure 2, the first two components account for over 93% of variation in the data set, with the first component accounting for 63%. The first component (PC1) separates two samples of one genotype from two samples of the other genotype, indicating genotype is the main factor that determines miRNA expression levels. Indeed, nine miRNAs account for 75% of variation in PC1 and they show clear genotype-specific expressions (Table 1, for predicted hairpin structures and mapping of small RNA reads to precursors, see Additional files 2 and 3). Whereas two miRNAs (vun_cand014 and vun_cand054) are predominantly expressed in CB46, the other seven miRNAs are exclusively or predominantly expressed in IT93K503-1 plants. The expression pattern of IT93K503-1 specific miRNA, vun_cand058, was confirmed with northern blot assay (Figure 3b).

Because perfect matches were required when small RNA reads were mapped to cowpea sequences for miRNA prediction, genotype-specific expression of miRNAs could be caused by inter-genotype single nucleotide polymorphisms (SNPs) in mature miRNAs. To address this possibility, we re-mapped clean small RNA reads from each library to the precursors of nine miRNAs in Table 1, allowing up to one mismatch. The normalized counts of these mature miRNAs were essentially unchanged (data not shown). Therefore, genotype-specific expression of these miRNAs was genuine and was not an artifact of the reads mapping process.

Drought-associated miRNAs

To identify drought-associated miRNAs, we tested for differential expression of miRNAs in drought-stressed and corresponding control samples in each genotype using the statistical method developed by Audic and Claverie [38]. We used the following criteria to identify drought-associated miRNAs: (1) adjusted p-value was less than 0.01 in at least one of the two comparisons; (2) normalized counts (TPTM) was at least 100 in one of the four libraries; (3) log2 ratio of normalized counts between drought and control libraries was greater than 1 or less than -1 in one of the two genotypes. For differential expression analysis, we considered only unique mature miRNAs as they are the active form of the miRNA and in some cases, identical mature miRNA can be generated from two or more homologous miRNA genes. We found 44 drought-associated unique mature miRNAs that belong to 28 families (Additional file 4). The direction of statistically significant change was the same in both genotypes for all 44 miRNAs, indicating that miRNA gene expression in IT93K503-1 and CB46 had similar overall response to drought stress. Whereas thirty of 44 miRNAs were upregulated in the drought-stressed condition, fourteen were downregulated in one or both genotypes.

Among 44 drought-associated miRNAs, the expression of 22 miRNAs (17 families) in drought condition changed at least two-fold compared to the control in both genotypes (Additional file 4). Some of these miRNA families have been found to be associated with drought stress in previous studies, including miR156 and miR166 [39], miR159 [40], miR167 [24], miR169 [22], miR171 [25,41], miR395 [40], miR396 [24,39], and miR482 [29]. Most of the predicted targets encode transcription factors (Additional file 4). Other miRNA families, miR162, miR164, miR319, miR403, miR828, and four cowpea-specific miRNA (vun_cand001, vun_cand010, vun_cand041, and vun_cand057) were found to be associated with drought stress for the first time (northern blot confirmation of vun_cand001 was shown in Figure 3b).

We also found that 12 miRNAs showed at least two-fold change only in IT93K503-1 (Table 2), and 10 only in CB46 (Table 3). Although statistical tests indicated that some of these miRNAs (e.g. miR1515 which was validated by northern blot as shown in Figure 3b) were up- or down-regulated under drought stress in both genotypes without having two-fold change, 11 miRNAs were clearly regulated in only one genotype. Whereas miR160a, miR160b, miR171e, vun_cand015, vun_cand033, and vun_cand048 were significantly regulated by drought stress in IT93K503-1 plants only, miR171b, miR171d, miR2111b, miR390b, and miR393 were regulated only in CB46.

Plant materials

CB46 and IT93K503-1 plants were grown in a greenhouse at the University of California Riverside campus in Spring 2009. The temperature was 35°C during the day and 18°C at night with no artificial control of day length. Four seeds were germinated in 2 gallon-pots filled with steam-sterilized UC Riverside soil mix UCMIX-3 and thinned to two plants per pot two weeks after planting. Three replicate pots per treatment were arranged in a completely randomized block design. When plants were 30 days old, corresponding to late vegetative stage, deficit irrigation treatments were applied by withholding watering on the stressed pots while controlled pots were water daily to soil capacity. Third leaf water potential was monitored using a pressure chamber (Cornallis, PMS instruments, USA) [55] as the indicator of the stress level. Fresh leaves (second from apex) of three replicates were sampled and frozen in liquid nitrogen from control plants (well watered, ψ_w = -0.5 MPa) and moderately stressed plants (ψ_w = -1.5 MPa) for RNA extraction.

Small RNA library construction and sequencing

Total RNA was extracted with the TRIzol reagent (Invitrogen) according to the manufacturer's instructions. Small RNA libraries were constructed from cowpea leaves using the procedure used by Sunkar and Zhu [20] with minor modifications [56]. Briefly, for each treatment/genotype group, equal amount of total RNA was pooled from three replicates to generate ~700 μg of RNA. Pooled total RNA was resolved in a 15% denaturing polyacrylamide gel and the 20-30 nt small RNA fraction was extracted and eluted. A preadenylated adaptor (linker 1, IDT) was ligated to the 3' end of small RNAs with the use of T4 RNA ligase. Ligation products were then gel purified and subsequently ligated to an RNA adaptor at the 5'end. After ligation and purification, the products were used as template for RT-PCR. After synthesis and purification, the PCR products were quantified and sequenced using an Illumina Genome Analyzer.

miRNA identification

Only small RNA reads that passed the Illumina quality control and contained clear adaptor sequences were considered good reads for further processing. After adaptor sequence was trimmed, clean small RNA reads of 18nt or more were combined into unique sequences. Reads that match known plant repeats, rRNAs, tRNAs, snRNAs, and snoRNAs were removed. Unique small RNA reads were mapped to four genomic sequence resources with SOAP2 [57]: cowpea EST assembly available in HarvEST:Cowpea [58] (http://harvest.ucr.edu, version 1.17, 18,745 sequences, but we excluded those appear to be protein-coding genes), a combination of 260,642 cowpea gene-space random shotgun sequences [59] and 30,527 BAC end sequences (obtained from M.-C. Luo, UC Davis, http://phymap.ucdavis.edu:8080/cowpea), 54,123 cowpea Genome Survey Sequences (GSS) from dbGSS of GenBank http://www.ncbi.nlm.nih.gov/dbGSS/, and a draft cowpea genome assembly from 63× coverage Illumina pair-ended reads (296,868 contigs with total length of ~186 MB, available at http://www.harvest-blast.org). Perfect match was required.

We used the updated annotation criteria for plant miRNAs [60] and built an in-house pipeline for miRNA prediction. Unique reads with a redundancy of at least 10 copies are used as anchor sequences. With one end anchored at 10 bp from the mapped position, DNA segments of 100 - 300 bp that cover each anchor sequence were sampled with 20 bp as step size. Secondary structure of each segment was predicted with UNAFold [61]. We then examined the structures and only those met the following criteria were considered genuine miRNA candidates: (1) free energy is lower than or equal to -35 kcal/mol; (2) number of mismatches between putative miRNA and miRNA* is 4 or less; (3) number of asymmetrical bulges in the stem region is not greater than 1 and the size of each asymmetrical bulge is 2 or less; (4) strand bias - small RNA reads that map to the positive strand of the hairpin DNA segment account for at least 80% of all mapped reads; (5) precise cleavage - reads that map to the miRNA and miRNA* regions (defined as miRNA or miRNA* plus 2nt on 5' and 3' ends) account for at least 75% of all reads that map to the precursor. If two or more candidate hairpins were predicted from the same region, we compared these hairpins and chose a hairpin that has highest putative mature miRNA expression, lowest free energy, or shortest length.

In order to classify miRNAs into families, all predicted mature miRNAs were compared with themselves using the ssearch35 program in the FASTA package (version 3.5) [62]. Using a single-linkage algorithm, mature miRNAs with up to two mismatches were included in same clusters. Mature miRNAs were then compared with the mature miRNAs in the miRBase (Release 16) [37] using ssearch35. If a member in a cowpea miRNA cluster had a match (allowing up to two mismatches) in the miRBase, the family number of the known miRNA was assigned to the cluster, otherwise the cluster was annotated as a new family.

miRNA Target prediction

Mature miRNA sequences were used as query to search the cowpea EST assembly for potential target sites using miRanda [63]. The alignments between miRNAs and potential targets were extracted from the miRanda output and scored using a position-dependent, mispair penalty system [64-66]. Briefly, miRNA-target duplexes were divided into two regions: a core region that includes positions 2-13 from the 5' end of the miRNA, and a general region that contains other positions. In the general region, a penalty score of 1 was given to a mismatch or a single-nucleotide bulge or gap, and 0.5 to a G:U pair. Scores were doubled in the core region. A match was considered positive if the alignment between miRNA and target meets two conditions: (1) the penalty score is 4 or less; (2) total number of bulges and gaps is less than 2.

Principal component analysis

Counts of each mature miRNA were first normalized to transcripts per ten million (TPTM) according to the total number of clean small RNA reads in each of the four libraries. MiRNAs with combined expression of at least 50 TPTM were chosen for principal component analysis (PCA). We used the log2 values of miRNA normalized counts to build an expression matrix and used the princomp function in MATLAB (MathWorks Inc., Natick, MA) for PCA.

Statistical test for differential expression of miRNAs

Because deep sequencing of small RNAs provides a random sampling of mature miRNAs in the original small RNA pools, counts of miRNAs can be modeled by a Poisson distribution. We applied an established method [38,67] to calculate the p-value for differential expression of miRNAs between a drought-stressed sample and a control sample. The first step was to calculate a conditional probability using the formula:

Where N₁ is total number of clean reads in the control library, N₂ is total number of clean reads in the drought-stressed library, x is number of a mature miRNA in the control library, and y is number of the same mature miRNA in the drought-stressed library. A two-tailed p-value for differential expression was then calculated as p = 2q, where q was the accumulated probability:

Due to the x↔y symmetry of p(y|x), if q was greater than 0.5, p-value could be calculated as p = 2*(1-q). Bonferroni method was used to adjust p-values for multiple comparisons.

Northern blot analysis

~40 μg of total RNA were resolved in 15% denaturing polyacrylamide gels and transferred to neutral nylon membranes (Hybond NX). The RNA was transferred and fixed to the membranes by chemical cross-linking [68] and then hybridized to probes complementary to mature miRNA sequences at 38°C, overnight. After hybridization, the blots were washed twice, 5 minutes each at 38°C with washing solution (2X SSC, 0.1% SDS) and exposed to X-ray film to reveal the signals. Results obtained in Northern blot assays were verified in three replicated samples.