Antiquity of microRNAs and their targets in land plants

Abstract

MicroRNAs (miRNAs) affect the morphology of flowering plants by the posttranscriptional regulation of genes involved in critical developmental events. Understanding the spatial and temporal dynamics of miRNA activity during development is therefore central for understanding miRNA functions. We describe a microarray suitable for detection of plant miRNAs. Profiling of Arabidopsis thaliana miRNAs during normal development extends previous expression analyses, highlighting differential expression of miRNA families within specific organs and tissue types. Comparison of our miRNA expression data with existing mRNA microarray data provided a global intersection of plant miRNA and mRNA expression profiles and revealed that tissues in which a given miRNA is highly expressed are unlikely to also show high expression of the corresponding targets. Expression profiling was also used in a phylogenetic survey to test the depth of plant miRNA conservation. Of the 23 families of miRNAs tested, expression of 11 was detected in a gymnosperm and eight in a fern, directly demonstrating that many plant miRNAs have remained essentially unchanged since before the emergence of flowering plants. We also describe an empirical strategy for detecting miRNA target genes from unsequenced transcriptomes and show that targets in nonflowering plants as deeply branching as ferns and mosses are homologous to the targets in Arabidopsis. Therefore, several individual miRNA regulatory circuits have ancient origins and have remained intact throughout the evolution and diversification of plants.

Figure 1.

The miRNA Array Is Semiquantitative and Highly Reproducible. (A) Mean array values for each miRNA detected (average of two hybridizations) are plotted versus the number of times the corresponding RNAs were cloned and sequenced from wild-type, mixed stage C. elegans (Lim et al., 2003). The linear regression line shown has a correlation coefficient of 0.623; this regression excludes two outliers, miR-2 and miR5-2, which are noted as triangles. (B) Technical replicate array hybridizations. Log₂-transformed array values from two hybridizations of independently labeled C. elegans wild-type, mixed stage RNA from the same total RNA sample. (C) Technical replicate array hybridizations. Log₂-transformed array values from two hybridizations of independently labeled C. elegans glp-4(bn2) RNA from the same total RNA sample. (D) Mean array values from Arabidopsis rosette leaves (diamonds), inflorescences (triangles), and seedlings (squares) are plotted versus the number of times the corresponding RNAs were cloned and sequenced from the corresponding Col-0 tissues. The linear regression line shown has a correlation coefficient of 0.593. (E) Biological replicate array hybridizations. Log₂-transformed array values from two hybridizations of independently labeled Arabidopsis inflorescence RNA from two RNA samples derived from different crops of plants grown under identical conditions. (F) Biological replicate array hybridizations. Log₂-transformed array values from two hybridizations of independently labeled Arabidopsis rosette leaf RNA from two RNA samples derived from different crops of plants grown under identical conditions. Lines in (A) to (F) represent linear regressions with the given correlation coefficients (r²) and slopes.

Figure 2.

A Global Expression Profile of Arabidopsis Small RNAs. (A) Log₂-transformed array values for Arabidopsis spots that were above the detection threshold (see Methods) in both replicates for at least one organ were hierarchically clustered both by gene and by hybridization sample. Log₂-transformed values are displayed as a color gradient from gray (low values) through white (intermediate values) to bright red (high values), with black indicating not detected (N.D.) (B) Differentially expressed small RNAs in Arabidopsis tissues. Arabidopsis small RNAs whose array values in at least one organ were significantly different than those in other organs were F3-3_B01-5, miR157, miR172, miR156, miR396, miR398, miR160, and miR163 (single-factor analysis of variance with Bonferroni-Holm corrected P-values of < 0.05). The other differentially expressed small RNAs with lower confidence were miR167 (P = 0.056), miR169 (P = 0.069), miR394 (P = 0.064), miR158 (P = 0.073), siR480(+) (P = 0.118), and miR171 (P = 0.117). To highlight the relative differences between tissues, the expression values for each miRNA were normalized on a per-gene basis to their median level of expression and log₂-transformed. Yellow indicates high relative expression, gray indicates moderate relative expression, and blue indicates low relative expression. N.D., not detected.

Figure 3.

Comparison of miRNA and Target mRNA Expression Levels. (A) The median relative level of miRNA expression plotted versus the median relative level of target mRNA expression within seven Arabidopsis organs. (B) Analysis as in (A), except the median expression levels of miRNA targets were replaced with the median expression levels of the closest paralogs of each target that do not have recognized miRNA complementary sites. Paralogous nontargets were chosen as the top nontarget Arabidopsis BLASTp hits at an e-value of < 1e-10 to a target that were called present in greater than half of the array experiments (four maximum per miRNA). (C) Analysis as in (A), except the median expression levels of miRNA targets were replaced with median expression levels of nonparalogous mRNAs. These control RNAs were randomly selected members of the WRKY and MADS box families of transcription factors. (D) Plots showing the distributions of linear correlation coefficients between relative miRNA expression and expression of target (n = 43), paralogous nontarget (n = 33), and control (n = 35) mRNAs within seven Arabidopsis organs. Central lines indicate median values, boxes bound the 25th to 75th percentiles, and the bars indicate the 10th and 90th percentiles. The median value of the targets was significantly lower than the medians of both the paralogous nontargets and the controls (P = 0.0038 and P < 0.0001, respectively; Mann-Whitney U-test). The median value of the paralogous nontargets was also significantly less then the median of the controls (P = 0.0079; Mann-Whitney U-test).

Figure 4.

Detection of Ancient miRNAs. (A) Duplicate array hybridizations were used to probe RNA samples from various specimens representative of major clades of land plants. RNAs detected in the species indicated at right for each row are shaded gray, whereas those not detected are left unshaded. Arabidopsis is represented by two rows, one for the composite of all organs (all), and the other for the analysis of rosette leaves only (leaves). Names of the major groups represented by each species are shown at the far right. An abbreviated cladogram displaying the evolutionary relationships of the sampled species is shown at left. Major evolutionary innovations are marked on the cladogram by perpendicular hash marks. Detection of three endogenous siRNAs is indicated. Unclassified small RNAs, of which none were detected outside of Arabidopsis, were omitted from this figure. (B) RNA gel blot using probes specific for Arabidopsis miR390, miR160, and the U6 small nuclear RNA. The blot was sequentially probed and stripped. Markers (M) are ³³P-labeled RNAs with sizes indicated at left. Initials are abbreviations of the species listed above.

Figure 5.

Identification of miRNA Targets in Nonflowering Land Plants. (A) Scheme of strategy for empirical discovery of miRNA targets. Black rectangles represent adapter sequences added before or during cDNA synthesis, and gray rectangles represent miRNA complementary sites. (B) Ethidium bromide–stained gel showing 5′-RACE products from candidate miRNA targets. Pixel values were inverted for ease of viewing. Bands in lanes 1 to 8 are identified in (C). Numbers at bottom indicate the predicted sizes of cleavage products (in kilobases) assuming termination at position 10 relative to the alignment with the Arabidopsis miRNA. Cloning and sequencing showed that products with higher molecular weights in lanes 4 and 6 were PCR artifacts resulting from amplification of nontarget cDNAs. Marker lanes (M) contain DNA standards, with sizes (kilobases) noted at the left. (C) Mapping 5′ ends of cleavage fragments of miRNA targets. In each duplex, mRNA sequence of the target is shown at top, aligned with the Arabidopsis miRNA sequence on the bottom. Base pairs are indicated by solid lines and G:U wobbles by circles. Predominant positions of 5′ target ends are indicated. Fractions refer to the number of independently cloned RACE products whose 5′ end terminated at the indicated position (numerator) over the total number of sequenced clones matching the target gene (denominator). In most cases, no 5′-RACE products were obtained upstream of the miRNA complementary site, precluding determination of the target sequence that pairs with the 3′ of the miRNA. Nucleotides in red indicate the 5′ ends of fern miR171 and fern miR172 as determined by PCR and sequencing.

Figure 6.

Moss Small RNAs. The lengths and 5′ residues of 214 nonredundant small RNAs from the moss P. juniperinum. RNA fragments that were obvious tRNA or rRNA fragments were eliminated before this analysis, and sequences cloned multiple times are only represented once. Sequences are listed in Supplemental Table 7 online.