Homology of miRNAs to Boechera ESTs

The BLASTn search using a reference set of 8433 non-redundant known conserved plant miRNAs against flower-specific sexual and apomictic Boechera EST libraries led to the identification of 282 sexual and 301 apomictic transcripts with high homology to miRNAs of other plant species (Figure 1). Of these, 13 sexual and 16 apomictic transcripts could fold back into stable hairpins containing conserved miRNAs (Table 1 & 2; Additional file 1, Figure S1). Many EST sequences were found that could not fold back into stem loops, although it is unclear whether this was due to the fact that they were not pre-miRNAs or whether this was due to sequencing errors introduced by the 454 FLX system. Predominantly, the less conserved miRNA families (e.g. miR444 to miR869) matched a small number of cDNA sequences which in most cases were found to be truncated precursor sequences in the EST libraries, and thus could not fold into stable stem-loops (Table 1 & 2).

Bioinformatically-identified conserved miRNA families

In all, 44 miRNA families across 67 plant species were found to match at least one Boechera 454 EST read, with A. thaliana being the predominant species (Figure 2). Conserved plant miRNA families in Boechera were identified to a large extent based upon high homology with reported conserved A. thaliana miRNAs (Figure 2). In cases where Boechera and A. thaliana did not share particular miRNA families, a search for conserved miRNA families was performed in other plant species. The predominant miRNA families which shared similarity with the highest number of Boechera 454 reads were miR156, miR157, miR160, miR167 and miR172 (Figure 3). It was observed also that the Boechera miRNAs exhibit a wide variation in the length of pre-miRNA sequences (Table 1 & 2; Figure 4, 5, 6 & 7). 29 families were found to be common between the sexual and apomictic genotypes. Of these, 17 mature miRNAs (miR156, 160, 167, 170, 172, 395, 396, 408, 415, 529, 824, 835, 841, 846, 859, 860 and 865) were similar in sequence, whereas 12 were different in sequence constitution due to nucleotide differences between the two reproductive modes. These included miRNAs miR157, 159, 161,166, 319, 394, 398, 400, 414, 854, 861 and 869 (Table 1 & 2). Pre-miRNA lengths varied from 66 to 233 nucleotides, with most between 66 and 184 nucleotides, a length similar to that of pre-miRNAs in other species. The location of the mature miRNAs in the precursor pre-miRNAs also varied among the miRNA families. In 12 pre-miRNAs, the miRNAs were found in the 3' arm while 18 were in the 5' arm of the stem-loop hairpin structures (Table 1 & 2; Figure 4, 5, 6 & 7).

Evaluation of the pre-miRNAs was also based on A+U content. The miRNA precursors have A+U content ranging from 41.96 to 63.46% (Table 1 & 2; Figure 4, 5, 6 & 7), similar to proportions observed in other plant species [41]. Consistent with general notion, the majority of identified Boechera miRNA precursors contain more A+U nucleotides than G+C [42, 43]. It is also important to note that the formation of a stem-loop structure is not a unique feature of miRNAs, since other RNAs such as mRNA, rRNA, and tRNA can also form similar structures. For this reason, uniform systems for annotating new miRNAs comprising negative minimal fold energy (MFE), adjusted minimal fold energy (AMFE) and the minimal fold energy index (MFEI) have been developed [42–45] and have become generally accepted. Zhang et al [43] indicated that most identified miRNA precursors have an MFEI greater than 0.85, which is much higher than in tRNA (0.64), rRNA (0.59), or mRNA (0.65). However, a number of pre-miRNAs with lower MFEIs have been reported, provided the number of nucleotide substitutions in the particular conserved miRNA compared with other species does not exceed three (Table 1 & 2; [41]).

Microarray analysis of conserved miRNA families

The miRNAs identified from cDNA sequencing of floral tissues, using the bioinformatics described above, were further verified using a proprietary microarray analysis with LC Sciences, in order to validate their expression in sexual and apomictic Boechera flower tissues. The LC Sciences proprietary miRNA microarray chip that was used was designed by spotting all known plant miRNAs that were available in the miRBase Release 14 (total 1117 unique mature miRNAs) and the Plant miRNA Database, PMRD (total 5690 unique mature miRNAs). Subsequently hybridization was performed as described in Methods using isolated enriched Boechera small RNAs to confirm expressed conserved miRNAs. As expected, most (n = 50) mature miRNAs representing 22 miRNA families were identified to be conserved mainly compared to A. thaliana. The microarray assay confirmed 15 conserved families identified with the bioinformatics techniques. It is also noteworthy that 7 and 29 other miRNA families were respectively detected separately by the microarray and bioinformatics approaches (Additional file 2, Figure S2).

Boechera-specific miRNA nucleotide substitutions (NSs) enhance pre-miRNA stem-loop stability

The stability of a secondary structure is quantified as the amount of free energy released or used by forming base pairs. The more negative the free energy of a structure, the more likely is formation of that structure and its stability, because more stored energy is released, and this principle is used to predict the secondary structure of a particular sequence [46, 47]. Out of the 30 stable Boechera pre-miRNA stem-loop structures obtained, 19 contain miRNAs with nucleotide substitutions (NSs) when compared with corresponding Arabidopsis or other plant miRNAs. The frequency of A, C and G substitutions were similar between sexual and apomictic mature miRNAs, while U appeared to show a higher rate of substitution in the apomictic mature miRNAs (Figure 8). Considering that a single nucleotide change in the sequence of a target site can affect miRNA regulation, NS could conceivably be under selection pressure to enhance the conformation and thermodynamic stability of the pre-miRNA stem-loop structure. We thus examined whether these Boechera-specific nucleotide changes had any effect on the structure and thermodynamic stability of their corresponding pre-miRNAs.

In all, this nucleotide substitution-stability phenomenon was most common in our analyses of both apomictic (8 out of 9 miRNA families) and sexual Boechera (5 out of 9 miRNA families; Table 3), in addition to Arabidopsis (9 out of 11 miRNA families; Table 4). Naturally occurring miRNA NSs thus appear to confer optimal thermodynamic stability on pre-miRNA stem-loop structures in Boechera, and is consistent with similar analyses in other plants. For example, a similar comparison of the ΔG of predicted secondary structures of two variants of barley miR1137 precursor with a C and a G in the 13th position showed differences in stability between the variants [49]. Interestingly, Thakur et al [50] reported that species background may also be correlated with the calculation of both the minimum free energy and miRNA hairpin stability, although this difference appeared to be manifested at the level of mono- and dicots. Thus, at least with respect to the comparisons between closely related Boechera and Arabidopsis used here, our data imply that natural selection has guided sequence variation in these regulatory elements.

In one case pre-miRNA stability was also manifested on the intraspecific level, comparing sexual and apomictic Boechera. In the family miR394, the pre-miRNA of the sexual Boechera species has the same miRNA sequence as in Arabidopsis, however that of the apomictic species shows one C to A NS change at position seven (Figure 5A). The pre-miRNA stability was examined by introducing the apomictic NS into the sexual sequence at the same position and ΔGs compared. As expected there was a decrease in the negative ΔG by 6.5 kcal/mol in the "new" sexual pre-miRNA with the introduced apomictic NSs, suggesting that the sexual pre-miRNA is perhaps at its optimal thermodynamic stable state. This final evidence is consistent with trans-acting regulatory differences between sexual and apomictic ovules, the result of sequence variation in regulatory factors in the sexual (homozygous) versus apomictic (hybrid) genomes, as suggested by Sharbel et al [40].

Conserved Boechera miRNAs target many transcription factors (TFs)

The EST libraries from which the Boechera miRNAs were mined were flower-specific [39, 40], and expectedly, a number of identified TF-targeting miRNAs have been associated with flower development in other species. For example, miR156 and miR157, the homologues of the squamosa-promoter binding proteins and whose function is well conserved across plant species [43], were identified in both apomictic and sexual Boechera (Table 5; Additional file 3, Table S1). In Arabidopsis these TF regulatory miRNAs have been reported to regulate the Antirrhinum floral meristem identity squamosa promoter binding protein-like (SPL) genes [56]. Other TF regulatory miRNA families which have regulatory roles during flower development in other species were also identified (Table 5; Additional file 3, Table S1), including miR156, miR159, miR164 and miR172, which have been implicated in the control of LFY expression, floral organ identity, and flowering time [27, 57, 58]. miR172 has furthermore been reported to regulate stem cell fate, and defines the inner boundary of the APETALA3 and PISTILLATA expression domains in Arabidopsis floral meristems [38].

A number of well-defined TF targeting miRNAs were also identified in Boechera. For example, miR160 and miR167 (Table 5; Additional file 3, Table S1) are associated with post-transcriptional regulation of the A. thaliana auxin response transcription factor (ARF) family genes [26, 59]. miR319 is known to regulate the expression of TCP transcription factor genes whose down-regulations cause abnormalities in leaf development [16]. Vierstra [60] showed that miR394 regulates the messages of F-box proteins, which in turn target specific proteins for proteolysis by making them substrates for ubiquitination by SCF E3 ubiquitin ligases. Growth Regulating Factor genes, the targets of the miR396 family, are putative transcription factors that regulate cell expansion in leaf and cotyledon [61]. Argonaute, one of the important proteins in the regulation of miRNA biogenesis, is a target of miR403 whereas miR408 regulates a copper ion binding protein. The miR414 family regulates a number of other genes including the transcription factors, transducin family protein/WD-40 repeat family protein and peptidyl-prolyl cis-trans isomerase cyclophilin-type family protein.

Expression patterns of transcription factor (TF) targets and apomixis in Boechera

The switch from sexual to apomictic seed production is hypothesized to involve global regulatory changes during ovule development which are induced by hybridization and/or polyploidy [9, 62], both common characteristics of apomictic plants and parthenogenetic animals. Using data from a previously-published SuperSAGE analysis [39, 40], the ovule expression patterns of putative target TFs for the miRNAs identified here were compared between sexual and apomictic Boechera across four ovule developmental stages. Of the 17 TFs identified as potential miRNA targets, expression data for 6 were found in the SuperSAGE libraries, including: the squamosa promoter binding protein like SPL6, SPL11 and SPL15, Myb domain protein 120 (MYB120), RAP2.7, TOE1 RAP2.7 (RELATED TO AP2.7) DNA binding and TCP10 (TCP family transcription factor 10), which are targets of the miRNA families miR156/157, miR159, miR172 and miR319 respectively.

It is noteworthy that, whereas the other genes showed no significant differential expression levels between sexual and apomictic species, SPL11 was found to be significantly (p ≤ 0.05) up-regulated at the stage two of ovule development in apomictic species in all libraries studied (Figure 9). SPL11 also showed low level expression in all the other apomictic ovule stages and at only stage two of the sexual ovules. Using six apomictic and five sexual genotypes of Boechera, the differential expression of SPL11 at ovule stage two of floral development was further validated using quantitative Real Time-PCR. With the exception of a single sexual B. divaricarpa from Mule Ranch, Montana, all apomictic accessions clearly showed relatively higher expression of SPL11 than the sexuals (Figure 10), result which is consistent with the expression pattern observed with the SPL11 SuperSAGE tag (Figure 9). The single sexual outlier (Figure 10) for SLP11 implies that the expression pattern of this TF may not be a key factor associated with apomixis expression, but rather is associated with DNA sequence variation in regulatory factors in the hybrid B. divaricarpa. Alternatively, population-level variation for TF expression could be associated with the penetrance of the apomictic phenotype, which has been shown to be genotype-specific in Boechera[63].

Flower-specific Boechera 454 cDNA libraries used

Floral cDNA libraries used in this study are those previously reported by Sharbel et al. [39, 40]. These libraries were sequenced from pooled flower stages 1-12 [64] of three diploid sexual plants (Accessions ES910-2-2 K, 105.6-1 K and B07261) and three apomictic plants (Accessions 67.5-K, 300.6.1-1 K and 218.2-2 K).

Conserved miRNA reference set for bioinformatics procedures

A total of 9274 previously reported non-redundant 21-24 nucleotides long miRNAs (including their precursor sequences) collected from 121 plant species were obtained from the Plant MicroRNAs Database (PMRD as of February 8, 2011; [65]). These miRNAs were defined as a reference set of miRNA sequences for the identification of potentially conserved miRNAs in Boechera. To avoid redundant miRNAs, duplicated miRNAs shared between different species within the database were removed. In all, 8433 non-redundant miRNAs were obtained, and these were used as query sequences for a BLASTn search against all original 454 sequence reads from the apomictic and sexual Boechera libraries.

Identification of conserved miRNAs

The bioinformatics approaches used for identification of conserved miRNAs in Boechera species are outlined in Figure 1. The length of the EST sequences used to search for conserved miRNAs ranged between 51 and 478 nucleotides, with about 80% of them around 200 nucleotides long. In order to exclude all ESTs having exact matches to tRNA or rRNA sequences from further BLASTn searches, the sexual and apomictic EST libraries were first queried against ribosomal RNAs database from Rfam (http://www.sanger.ac.uk/Software/Rfam/) and the Arabidopsis transfer RNAs database (http://lowelab.ucsc.edu/GtRNAdb/Athal/). Rather than using the miRNA precursors for BLASTn searches against our databases, the analysis was based mainly on the mature miRNA sequences considering that only mature miRNAs are highly conserved in plants [42, 43]. The following BLASTn parameters which gave the highest and most reliable number of hits (blastall -p blastn -m 8 -e 1 -W 7 -r 1 -q -1 -i) were used. All resulting EST sequences with an alignment length of 20-24 nucleotides, three or fewer mismatches and no gaps compared to previously identified plant miRNAs were selected and compared with each other to eliminate redundancies. The obtained non-redundant sequences were then used for the prediction of secondary structures and screening for miRNA precursor sequences. The secondary structures of pre-miRNAs were generated using the Mfold 3.2 software, which is based on Zuker folding algorithm principles [48].

The secondary structure of candidate pre-miRNA sequences were analysed and scored for their potential to form miRNA precursors. A stem-loop was selected as a candidate miRNA precursor if it satisfied most of the following generally accepted criteria: (1) the mature miRNA is 20-24nt with a maximum of three mismatches compared with the corresponding known miRNA in other plant species; (2) the miRNA precursor (pre-miRNA) sequence folds into a stable hairpin structure such that one arm of the hairpin contains the mature miRNA sequence; (3) the predicted secondary structure of the pre-miRNA has lower minimal free energy (MFE ≤-10 kcal/mol) and minimal free energy index (MFEI) than other types of RNA (e.g. tRNA, rRNA); (4) the predicted mature miRNA has an A+U content of 40-70%; and (5) no loop or gap in the mature miRNA sequences [41].

Microarray validation of conserved plant miRNAs

The bioinformatically-identified miRNAs in floral tissues were further verified using a proprietary microarray analysis with LC Sciences, USA. The microarray assay was performed using 4 to 8 μg total RNA sample from pooled flower tissues of sexual and apomictic genotypes. The total RNA was size fractionated using a YM-100 Microcon centrifugal filter (Millipore) and the isolated small RNAs (< 300 nt) were 3'-extended with a poly(A) tail using poly(A) polymerase. An oligonucleotide tag was then ligated to the poly(A) tail for later fluorescent dye staining, and two different tags were used for two RNA samples in dual-sample experiments. Hybridization was performed overnight on a μParaflo microfluidic chip [spotted with all known plant mature miRNAs that were available in miRBase Release 14 (total 1117 unique mature miRNAs) and the Plant miRNA Database (total 5690 unique mature miRNAs)] using a micro-circulation pump (Atactic Technologies; [66]). On the microfluidic chip, each detection probe consisted of a chemically modified nucleotide coding segment complementary to a target miRNA, and a spacer segment of polyethylene glycol to extend the coding segment away from the substrate. The detection probes were made by in situ synthesis using PGR (photogenerated reagent) chemistry. The hybridization melting temperatures were balanced by chemical modifications of the detection probes, and hybridization was performed using 100 μL 6 × SSPE buffer (0.90 M NaCl, 60 mM Na₂HPO₄, 6 mM EDTA, pH 6.8) containing 25% formamide at 34°C. After RNA hybridization, tag-conjugating Cy3 and Cy5 dyes were circulated through the microfluidic chip for dye staining. The fluorescence data were collected on an Axon GenePix 4000B Microarray Scanner, and then analysed by first subtracting the background followed by normalization of the signals using a LOWESS filter (Locally-weighted Regression; [67]). A detectable miRNA on the array was identified if its signal intensity was higher than 3×(background standard deviation) and spot CV < 0.5, and p < 0.01 for the difference between Cy3 and Cy5 signals (LC Sciences).

Prediction of Boechera gene targets of miRNA families

A BLASTn search (blastall -p blastn -m 8 -e 1 -W 7 -r 1 -q -1 -i) was employed to detect complementarity between the validated miRNAs and predicted target ESTs in sexual and apomictic Boechera (Additional file 3, Table S1). Putative miRNA targets were identified based on the total numbers of mismatched nucleotides between miRNAs and the alignment structures of potential targets. To identify potential regulatory targets, a BLAST search (blastall -p blastn -m 8 -e 1 -W 7 -r 1 -q -1 -i) was performed using the validated (from LC Sciences) conserved Boechera miRNAs against the A. thaliana protein-coding nucleotide databases (TAIR9 cDNA) using the miRU web server [68] from the Arabidopsis Information Resource (TAIR). The total number of allowed mismatches at complementary sites between miRNA sequences and potential mRNA targets in Arabidopsis were limited to a maximum of three, and no gaps were allowed at complementary sites. Finally, the Boechera homologues of potential targeted genes in Arabidopsis were chosen using a BLAST search (blastall -p blastn -m 8 -e 1 -W 7 -r 1 -q -1 -i) based on the degree of similarity of protein-coding mRNAs between A. thaliana and Boechera.

Expression analysis of Transcription factor (TF) targets using SuperSAGE tags

Finally, a comparative gene expression analysis of TF targets from 11 miRNA families was carried out. First, a BLASTn search using TF genes from Arabidopsis against the assembled Boechera EST database was performed (blastall -p blastn -m 8 -e 1 -W 7 -r 1 -q -1 -i). For each Arabidopsis TF, homologous Boechera TFs with an alignment having a bit score ≥100 were selected. Next, 100% sequence matches between the Boechera TFs and expression tags from 8 ovule-specific Boechera SuperSAGE libraries [40] were found using a BLASTn search. Finally, the expression patterns of the selected Boechera TFs corresponding to the obtained SuperSAGE tags were compared across four different ovule developmental stages between a sexual and apomictic Boechera genotype [40].

Quantitative RT-PCR validation of differential SPL11 expression in ovule stage two of Boechera flowers

Six accessions of apomictic and five of sexual Boechera were selected for the validation the differential expression of SPL11 (Additional file 4, Table S2). From these accessions, stage two ovules were micro-dissected, RNAs isolated and cDNAs prepared as described in Sharbel et al. [40]. The forward primer 5'-CAAAGTGCCCAAAAGTTACCGTGAGT-3' and reverse primer 5'-ACGCCTCGCATTATGATGAGAAAGA-3' with amplicon size of 137 nucleotides long were used for the qRT-PCR. Primers were designed avoiding intronic regions (to ensure the elimination of likely DNA contamination in samples) using the following parameters: temperature; 60°C, 20% < CG content < 80%, and PCR product size < 150 bp. For the real-time PCR reactions, the SYBR Green PCR Master Mix (Applied Biosystems) was used. qRT-PCR amplifications were performed in a 7900 HT Fast RT-PCR system (Applied Biosystems) with the following temperature profile for SYBRgreen assays: initial denaturation at 90°C for 10 min, followed by 40 cycles of 95°C for 15 s, and 60°C for 1 min. The Ct, defined as the PCR cycle at which a statistically significant increase of reporter fluorescence is first detected, was used as a measure for the starting copy numbers of the target gene. The mean expression level and standard deviation for each set of three technical replicates for each cDNA was calculated. Relative quantitation and normalization of the amplified targets were performed by the comparative ΔΔCt method in reference to the expression levels of the housekeeping gene ubiquitin [69].