Functional and evolutionary analyses of the miR156 and miR529 families in land plants
© Morea et al. 2016
Received: 17 September 2015
Accepted: 18 January 2016
Published: 3 February 2016
The Erratum to this article has been published in BMC Plant Biology 2016 16:153
MicroRNAs (miRNAs) are important regulatory elements of gene expression. Similarly to coding genes, miRNA genes follow a birth and death pattern of evolution likely reflecting functional relevance and divergence. For instance, miRNA529 is evolutionarily related to miRNA156 (a highly conserved miRNA in land plants), but it is lost in Arabidopsis thaliana. Interestingly, both miRNAs target sequences overlap in some members of the SQUAMOSA promoter-binding protein like (SPL) family, raising important questions regarding the diversification of the miR156/miR529-associated regulatory network in land plants.
In this study, through phylogenic reconstruction of miR156/529 target sequences from several taxonomic groups, we have found that specific eudicot SPLs, despite miRNA529 loss, retained the corresponding target site. Detailed molecular evolutionary analyses of miR156/miR529-target sequence showed that loss of miR529 in core eudicots, such as Arabidopsis, is correlated with a more relaxed selection of the miRNA529 specific target element, while miRNA156-specific target sequence is under stronger selection, indicating that these two target sites might be under distinct evolutionary constraints. Importantly, over-expression in Arabidopsis of MIR529 precursor from a monocot, but not from a basal eudicot, demonstrates specific miR529 regulation of AtSPL9 and AtSPL15 genes, which contain conserved responsive elements for both miR156 and miR529.
Our results suggest loss of functionality of MIR529 genes in the evolutionary history of eudicots and show that the miR529-responsive element present in some eudicot SPLs is still functional. Our data support the notion that particular miRNA156 family members might have compensated for the loss of miR529 regulation in eudicot species, which concomitantly may have favored diversification of eudicot SPLs.
MicroRNAs (miRNAs) are small RNAs important to transcriptional and post-transcriptional regulation in animals, plants, and viruses. MiRNAs bind complementarily to their target mRNA sequences, leading to translational repression, RNA degradation, or RNA cleavage . Most plant miRNAs are encoded by gene families, and mature miRNAs often have multiple target genes with similar complementary motifs in their mRNAs. The almost perfect complementarity between miRNAs and targets facilitates computational prediction and could be due to their evolutionary origins. One well-accepted model for MIR gene evolution is the inverted duplication of target gene sequences in plant genomes. These duplicated regions become MIR genes over time through sequence variations, which permit the generation of hairpin-like transcripts to become substrates for DICER-LIKE enzymes . Another interesting model suggests that MIR genes can randomly originate from various inverted repeats throughout the genome, independently of target gene duplications [3–6]. For instance, recent evidences indicate that MIR genes may be generated as a result of transposon activity [7, 8].
It has been recently shown that evolutionary patterns of miRNA genes (including lineage-specific gain or loss) can potentially influence the evolution of their target genes . Moreover, synonymous codons in target genes have been found near conserved miRNA target sites in at least four plant genomes, indicating selection constraint on synonymous codons for efficient miRNA binding and proper miRNA function . Several miRNAs are conserved throughout large evolutionary distances from embryophyta to core rosids. However, some miRNAs appear to be species or lineage specific [11–13]. For instance, miR156 is conserved in all Angiosperms studied thus far. Interestingly, miR156 is correlated at the nucleotide level with miR529, sharing 14–16 nt . Curiously, although both miRNAs have a common ancestor in embryophytes, miR529 seems to have been lost in some taxonomic groups, including core eudicots such as Arabidopsis thaliana [12–14]. Both miRNAs share target genes, which are members of the SQUAMOSA promoter-binding protein like (SPL) family . SPL genes are plant-specific transcription factors defined by a highly conserved region of 76 amino acids called SBP domain , and by their crucial and widespread functions in development [17–23]. SPL genes can be roughly separated into two major groups–long and short–the latter containing responsive elements for miR156 and, in some species, for miR156 and miR529. Interestingly, sites for miR156 reside in coding as well as untranslated regions of target sequences, whereas miR529 binding sites are chiefly located in coding regions and overlap with miR156 sites [24, 25].
In plant lineages containing both miRNA genes, differential expression of miR156 and miR529 in vegetative and reproductive organs/tissues might have favored lineage-specific retention of miR156/529 sequence variants in their genomes . Another possibility is that the combinatory action of miR156 and miR529 leads to the regulation of distinct targets in specific lineages, such as in Physcomitrella patens . Recently, the evolutionary differences between these two miRNA families in monocot species have been investigated , but the consequences of miR529 loss in core eudicots such as Arabidopsis are not yet clear. A broader evolutionary analysis of MIR156 and MIR529 genes and their targets including eudicot species should offer valuable insight into this issue.
In this study, we extended our knowledge regarding evolutionary and functional divergences between miR529 and miR156 regulation on their conserved targets. Through phylogenic reconstruction of target sequences, we confirmed that specific eudicot SPL family members retained miR529-target sites, independently of the presence of MIR529 genes in their genomes. Loss of miR529 function in eudicots is correlated with a more relaxed selection of the miRNA529-specific target element, while miRNA156-specific target sequence is under stronger selection. Additionally, we showed that A. thaliana plants overexpressing MIR529 precursor from a monocot, but not from a basal eudicot, display similar phenotypes as the spl9;spl15 mutant due to the specific down-regulation of these miR156/529-targeted SPLs. Based upon functional and evolutionary analyses, we proposed that the loss of MIR529 genes might have favored diversification of SPLs in eudicot species. It is also possible that new miR156 family member(s) have replaced miR529 functions in eudicots.
Results and discussion
Sequence and phylogenetic analyses reveal that a miR529-responsive element in SPL genes is broadly conserved across land plants
To further explore the evolutionary relationship of miR156/529 common targets, phylogenetic inference of SBP-box genes containing the miR156/529-responsive element was estimated using maximum-likelihood (ML) and Bayesian inference methods. The percentage of pairwise identity of well-aligned sequence blocks (see Methods) was 73.3 %, and the substitution saturation test reported that the alignment was not saturated (data not shown). We observed two groups of SPLs in our consensus tree (Fig. 1b). Group I included known miR156/529 SPL targets in bryophyte, whereas group II contained various monocot and core eudicot SPLs harboring conserved binding sites for miR156/529. This analysis indicated that SPLs containing miR156/529 target sites have a common origin in land plants (Fig. 1b). A. thaliana SPL9 and SPL15 are closely related and most likely form a pair of paralogous genes [29, 30]. Accordingly, both SPL9 and SPL15 as well as their orthologs retained the miR156/529-responsive element (Fig. 1b and c).
It has been proposed for monocot species and P. patens that SPLs containing miR156/529 sites evolved conservatively with a slow rate when compared with SPLs harboring only the miR156-responsive element . To further elucidate the evolutionary fates of eudicot SPLs containing the miR156/529-responsive element, we analyzed two blocks in SPL sequences: “SBP domain” block, which contains nucleotides of the SBP domain  plus few nucleotides upstream/downstream, and “target site” block, which contains nucleotides that comprise both miR156/529-responsive elements (see Methods). For the “SBP domain” block, we estimated nonsynonymous (Ka) and synonymous substitution (Ks) ratios (Ka/Ks) of representative SPLs. We chose Arabidopsis as a representative eudicot due to the fact that, even after extensive sequencing efforts, precursors or canonical mature sequences of miR529 have not been found in either A. thaliana or its closest relative A. lyrata . Pairwise alignments of best-aligned blocks (249 to 2061 nt) among SPL genes from A. thaliana, monocots, and P. patens indicated that the “SBP domain” block is under purifying selection for all comparisons, reflecting functional constraints, which is expected since SBP is likely a DNA binding domain (Additional file 2; ).
Overexpression of a basal eudicot microRNA529 precursor in A. thaliana phenocopies miR156 overexpressor
Although miR156 is highly conserved, being present in all plant species assessed thus far, miR529 seems to be restricted to particular taxonomic groups . To get a better view of the MIR156/529 gene evolution, the phylogenetic relationship of these miRNAs was accessed using the maximum-likelihood (ML) approach. For phylogenetic analyses, we included MIR156 and MIR529 precursors from Physcomitrella patens, monocot species (Oryza sativa, Zea mays, Brachypodium distachyon, and Sorghum bicolor), a basal eudicot (Aquilegia coerulea), and precursors of MIR156 genes from Arabidopsis thaliana. A consensus ML tree was generated in which two general, broad groups were readily observed (Additional file 3A). Group I comprised MIR156 precursors from different species, while group II contained MIR529 precursors from monocots and moss. Monocot MIR529 precursors were grouped into a distinct subset of group II (Additional file 3A), suggesting that evolutionary divergence occurred in a common ancestor of land plants, which led to the split between MIR529 genes of moss and flowering plants.
It had been proposed earlier that a key feature of miRNA evolution is that, once evolved, families and family members are rarely lost . However, not all miRNAs are equally conserved and it has been recently shown that several miRNA losses occurred in families that evolved prior to the origin of spermatophytes . Our data suggest that miR156 and miR529 families experienced dynamic duplications and losses across embryophytes, through which clade- or species-specific miRNA gene subgroups have arisen or were eliminated. For instance, A. thaliana has at least 10 MIR156 loci and 10 miR156-targeted SPLs, whereas rice has at least 12 MIR156 loci, two MIR529 loci, and eight miR156-targeted and four miR156/miR529-targeted SPLs .
Interestingly, the predicted MIR529 precursor from the basal eudicot A. coerulea  was grouped into group I, with A. thaliana MIR156h and Aquilegia MIR156a and b precursors, indicating a common origin of these miRNAs (Additional file 3A). Moreover, Aquilegia MIR529 seems to be highly conserved with Arabidopsis MIR156h at both nucleotide and secondary structure levels (Additional file 3B). These observations raised the question of whether this MIR precursor of Aquilegia defined as MIR529 is indeed a MIR156 homolog. To test this hypothesis, we investigated whether loci flanking Aquilegia pre-miR529 are localized into syntenic blocks when comparing with monocot species. We firstly searched for such syntenic groups among distinct monocot species, including Z. mays, O. sativa, B. distachyon, and S. bicolor (Additional file 4A). The colinearity of genes around pre-miR529 locus is relatively conserved, therefore defining conserved syntenic block in monocots. Conversely, we did not detect any conservation of syntenic block including pre-miR529 from Aquilegia, although some orthologous genes were found to be present in this species (Additional file 4B). Moreover, there is no detectable synteny between MIR156h-associated regions from eudicots (A. thaliana and V. vinifera) and the MIR529-associated region from basal eudicot A. coerulea (Additional file 4B, C). This complete loss of synteny, along with our phylogenetic and sequence analyses, suggests that the basal eudicot Aquilegia lost a bona-fide MIR529 gene, perhaps after the regulatory role of miR529 was perturbed due to mutations in the miR529 sequence as recently proposed .
In line with the observed phenotypes and expression analyses, RACE analysis of SPL15 cleavage sites demonstrated that transcripts are chiefly targeted by miR156 in 35S::AqcMIR529 plants (Fig. 3c). Based on these observations and given that AqcMIR529 precursor is more conserved with MIR156-like precursors (Additional file 3), we propose at least two non-exclusive possibilities: (1) the predicted AqcMIR529 precursor is more likely paralogous to AqcMIR156a and b genes and would produce a miRNA156-like small RNA; (2) AqcMIR529 precursor might have accumulated mutations in the miR529 sequence as recently proposed , leading to a loss of miR529 biogenesis and/or function. As other basal eudicots seem to accumulate miR529-like small RNAs , we cannot rule out the possibility that miR529-like small RNAs still accumulate in specific A. coerulea tissues, since there is no available information regarding large-scale identification of small RNA populations in this basal eudicot.
MiR156 and miR529 show overlapping expression patterns during rice vegetative development
OsmiR156 is dynamically expressed during rice leaf development, and a gradual increase of OsmiR156 expression might be essential for regulating the temporal expression of target genes, including OsSPL14 . We evaluated the expression pattern of OsmiR529b in similar developed leaves in seedling (L3-L5) and tillering stages. In contrast with the observed OsmiR156 expression patterns , OsmiR529b was ubiquitously expressed in all leaf developmental stages (Fig. 4b), suggesting that this miRNA has a minor or negligible contribution for temporal control of the expression of SPL genes during rice leaf maturation. Nevertheless, it is also possible that miR529 function as a dampening miRNA to establish the correct balance of SPL targets during temporal leaf development in monocots.
Given that OsmiR529b and OsmiR156a-j transcripts accumulated in the vegetative apex (Fig. 4a), we decided to evaluate their spatial expression patterns in the shoot apical meristem (SAM) via in situ hybridization using specific probes as described (; see Methods). Both miRNAs are expressed in incipient (P0) and developing leaf primordia, but not in the meristem proper (i.e., the peripheral and central zones) (Fig. 4c). Such expression pattern strengthened the data from Arabidopsis, which showed that miRNA regulation of SPLs occurs mainly in leaf primordia and that SPL activity may nonautonomously inhibit initiation of new leaves at the SAM, perhaps via auxin signaling pathways . Our in situ data showed that OsmiR529a-b and OsmiR156a-j have overlapping spatial expression patterns in leaf primordia, which suggest that these miRNAs can redundantly regulate or collaborate to fine-tune regulation of target expression in these organs. However, based on reported higher levels of OsmiR156 expression compared with OsmiR529 expression , the contribution for SPL expression modulation is unlikely to be comparable for both miRNAs, mainly during early stages of rice vegetative development.
Modified 5′-RACE procedure can be used to access cleavage products of miRNA targets as well as the processing of miRNA precursors [40, 41]. Parallel analysis of RNA end (PARE) signatures that are derived from rice degradome and that only mapped to the pre-miRNAs can give additional evidences of roles of miR156 and miR529 in rice development. We therefore collected PARE signatures of OsMIR156a-l and OsMIR529a-b precursors from publicly available resources (see Methods). Based on available data, it seems that rice MIR156 and MIR529 precursors are differentially processed, which may lead to differential miRNA accumulation across rice tissues/organs (Fig. 4 and Additional file 5). Even within the OsMIR156 family, specific members are differentially processed. For example, OsMIR156k and –l showed fewer PARE signatures when compared with the remaining MIR156 precursors. Likewise, OsMIR529a and b precursors have a smaller amount of signatures when compared with most OsMIR156 precursors (Additional file 5). It is possible that differential precursor processing also account for the distinct functions of OsmiR156 and OsmiR529 in development. It would be interesting in the future to address the question of whether miR529 regulation of SPLs is crucial for rice vegetative development. In contrast to miR156, the effects of over-expressing or down-regulating miR529 have yet to be examined in transgenic/mutant monocot plants.
Conserved miR529-responsive element is functional in AtSPL9 and AtSPL15 genes
Phenotypic evaluation of 35S::OsMIR529b lines in comparison with Col-0, spl9;spl15, and 35S::MIR156a plants under LD growing conditions
Rosette leaves (DAG)
6.8 ± 0.5c
7.7 ± 0.7c
8.8 ± 0.9c
13.3 ± 0.8c
14.3 ± 0.9c
16.5 ± 1.2d
23.6 ± 1.6c
19.4 ± 1.8a
15.0 ± 0.9a
19.5 ± 1.6a
35.6 ± 4.6a
52.9 ± 2.6a
45.0 ± 1.8a
54.0 ± 3.3a
10.3 ± 0.7b
11.0 ± 1.2b
13.3 ± 1.8b
21.1 ± 2.3b
29.4 ± 2.2b
20.0 ± 0.7bc
25.8 ± 1.9b
10.3 ± 0.5b
12.3 ± 1.1b
13.0 ± 0.6b
20.3 ± 1.1b
26.8 ± 2.0b
19.4 ± 0.7c
24.8 ± 1.5bc
10.0 ± 1.1b
11.4 ± 1.0b
12.6 ± 1.1b
20.3 ± 2.7b
27.5 ± 1.6b
18.9 ± 0.9c
23.4 ± 1.3c
9.8 ± 0.7b
11.0 ± 1.3b
13.2 ± 0.9b
19.6 ± 1.7b
26.5 ± 1.9b
20.9 ± 1.2b
26.8 ± 2.2b
In addition to estimating the number of leaves formed before the appearance of the first flowers, we also determined the time that 35S::OsMIR529b lines needed to bolt as well as to anthesis (Table 1). On average, transgenic plants overexpressing OsmiR529b showed a slight delay in the transition to flowering (3.3 days) when compared with Col-0 (wild type). In line with the observations of Schwarz and co-workers , we also observed that the double mutant spl9;spl15 showed an intermediate behavior between Col-0 and miR156 overexpressor (Table 1). Importantly, spl9;spl15 plants did not differ statistically from 35S::OsMIR529b lines in terms of transition to flowering and leaf development (Table 1). These data reinforced the observation that OsMIR529b overexpressors display similar vegetative and reproductive phenotypes as spl9;spl15 mutant, likely due to the low levels of SPL9 and SPL15 transcripts in both genotypes.
MiR156 targets, besides SPL9 and SPL15, exclusively other eight SPL family members  and these were shown to be down-regulated in AtMIR156b-overexpressing plants . In comparison with spl9;spl15 double mutant and 35S::OsMIR529b lines, 35S::MIR156a line displays more severe, aberrant vegetative and reproductive phenotypes (Table 1; ), which is likely due to the fact that additional miR156-targeted SPL genes act redundantly to regulate leaf initiation and phase change . Conversely, as in the spl9;spl15 double mutant, only AtSPL9 and AtSPL15 genes may be repressed in 35S::OsMIR529b lines, rendering them a less aberrant phenotype (Fig. 5a). To confirm this hypothesis, we evaluated the expression patterns of several SPL family members in leaf tissues of 35S::OsMIR529b lines, spl9;spl15 mutant, 35S::MIR156a line, and Col-0. We also evaluated the presence of transcripts from OsMIR529b precursor and the accumulation of its respective mature miRNA (Fig. 5b). As expected, several SPL genes were down-regulated in the 35S::MIR156a plants (Fig. 5b). However, in line with our phylogenetic (Fig. 1b) and phenotypic data (Table 1), only transcript levels of AtSPL9 and AtSPL15 showed a severe reduction in 35S::OsMIR529b plants. Moreover, levels of miR156 transcripts were similar in Col-0 and 35S::OsMIR529b plants (Fig. 5b), suggesting that the accentuated reduction in SPL9 and SPL15 expression in OsMIR529b overexpressors is most likely due to the accumulation of rice miR529b transcripts. Together, these results substantiated the fact that monocot MIR529 precursor is correctly processed in a core eudicot to generate mature miR529 transcripts and down-regulate specific SPL genes.
To confirm the post-transcriptional regulation of AtSPL genes by miR529b, we mapped cleavage sites in SPL15 transcripts in 35S::OsMIR529b plants employing the modified 5′-RACE approach . RACE analyses showed that cleavage sites occurred between the base 10 and 11 of the miR529b (Fig. 5c), similarly to what has been described for monocots . These results robustly showed that SPL15 is mainly regulated by OsmiR529b in Arabidopsis 35S::OsMIR529b plants, demonstrating a conserved function of miR529 in post-transcriptionally regulating specific SPL family members. Importantly, our data imply that the miR529-responsive element is conserved and functional in Arabidopsis SPL9 and SPL15 genes, likely due to the selective constraint on the amino acid or RNA secondary structure of the region surrounding miR156/529-responsive element.
We have shown that, although MIR529 genes have been lost in Arabidopsis and perhaps in all eudicot species, particular SPL genes in these species retained the miR529-responsive element, possibly due to the maintenance of synonymous codons for efficient miR156 binding and proper function . More specifically, A. thaliana SPL9 and SPL15 genes retained a functional miR529-responsive element, even in the absence of a miR529-generated locus. Similarly to monocot SPLs, eudicot SPL genes containing the miR156/529-responsive element appear to be under evolutionary constraints distinct from those containing only the miR156-responsive element. Such tendency would be indicative of target evolution constrained by miRNA-mediated regulation.
It is possible that the interplay between miR156 and miR529 regulation of specific SPLs be important to fine-tune flower architecture development in monocots, particularly in grasses [26, 37]. As Arabidopsis does not have miR529, perhaps particular miR156 family members (such as miR156h.2, which is preferentially expressed in flower tissues; ) functionally replace miR529. It is conceivable that other core rosids and/or closely related species of A. thaliana share similar miR156/miR529/SPL evolutionary history, though such confirmation requires future studies.
Plant material and growth conditions
Arabidopsis thaliana plants (ecotype Columbia-0 or Col-0) were grown at 21 °C (day)/19 °C (night) under long-day conditions (16 h light/8 h dark). Transgenic plants 35S::MIR156a and the double mutant spl9-1;spl15-2 were described . Transgenic plants were confirmed by PCR genotyping.
For transgenic Arabidopsis plants, the binary constructs 35S::OsMIR529b and 35S::AqcMIR529 were delivered into Agrobacterium tumefaciens GV3101 (pMP90) by the electroporation method. Transgenic plants were generated by the floral dipping method  and screened with 50 mg/mL kanamycin on half-strength MS plates. At least six independent kanamycin-resistant lines were selected for transgene integration by PCR and subsequently examined for transgene expression levels (data not shown). Further analyses were performed with selected lines in the T3 generation.
Rice seeds (Oryza sativa ssp japonica) were germinated on soil, and plants were grown under greenhouse conditions.
Oligonucleotide primers for all constructs are given in the Additional file 6. A 1000-bp fragment encompassing the OsMIR529b precursor was amplified from genomic DNA of O. sativa. The PCR product was subcloned into pGEM (promega) and sequenced. The confirmed OsMIR529b precursor was digested with BamHI and SacI restriction enzymes and subsequently cloned into the binary vector pBI121 behind the CaMV35S promoter. For 35S::AqcMIR529 construction, a 125-bp fragment encompassing the annotated AqcMIR529 precursor  was amplified from genomic DNA of A. coerulea, sequenced, and further cloned into the plant binary destination vector pK7WG2 (Gateway System; ) behind the CaMV35S promoter.
RNA extraction and stem–loop pulsed reverse transcriptase (RT)-PCR
Total RNA from Arabidopsis (leaf tissues) and rice (vegetative apices, leaf, and panicle tissues) was extracted using Trizol reagent (Life Technologies, USA) according to manufacturer’s instructions and subsequently treated with DNAse I (Life Technologies, USA). For miRNA and mRNA detection, DNAse I-treated RNA (2.0 μg) was reverse-transcribed to generate the first-strand cDNA, according to Varkonyi-Gasic et al. . Oligo(dT) primer was also added to the reaction for detecting target mRNAs and internal controls. cDNA dilutions were used for PCR reactions as follows: 1.0 μL cDNA, 1.5 mM Magnesium Sulfate, 0.25 mM each dNTP, 10 pmol each primer, and 1 U of Taq DNA Polymerase (Promega, USA). The reactions were done under the following conditions: 94 °C for two minutes and appropriate cycle numbers of 94 °C for 20 s, 60 °C for 30 s, and 72 °C for 45 s. All reactions were repeated at least twice with two biological samples. Primer sequences are described in Additional file 6.
Analysis of 5′-RACE
Five micrograms of total RNA from rosette leaves of Arabidopsis plants (Col-0, 35S::OsMIR529b and 35S::AqcMIR529) was ligated to a RNA adapter, in a reaction mixture containing 0.5 U/μL of T4 RNA Ligase, 4 U/μL RNAse inhibitor, and 1 mM ATP. The subsequent steps were performed according to the manufacturer’s guide of the GeneRacer kit (Invitrogen). The first PCR was done using the following AtSPL15 specific primer: 5′-AGCCATTGTAACCTTATCGGAGAATGAG. The PCR reaction was subsequently used as a template for a semi-NESTED PCR with an internal AtSPL15-specific primer (5′-TCATCGAGTCGAAACCAGAAGAT). After amplification, 5′-RACE products were gel-purified and cloned, and at least eight independent clones were randomly chosen and sequenced.
The number of rosette leaves was measured during several developmental stages (15, 20, 25, and 30 days after germination or DAG). Flowering and bolting time as well as the number of juvenile leaves were estimated as described. Data were subjected to statistical analyses by using the program ASSISTAT version 7.6 beta (t-student P <0.01).
In situ hybridization
Non-radioactive in situ hybridization was done as described . Oryza sativa vegetative apices were collected from seedlings 25 days after germination. Locked nucleic acid probes with sequences complementary to miR56 and miR529 as described  and negative control Scramble-miR (5′-GTGTAACACGTCTATACGCCCA) were synthesized by Exiqon (USA) and digoxigenin-labeled with the DIG Oligonucleotide 3′-end Labeling kit (Roche Applied Science, USA). Ten picomoles of each probe was used for each slide. Hybridization and washing steps were performed at 55 °C.
Phylogenetic and sequence analysis
Sequences of MIR156 and MIR529 precursors (pre-miRNAs) were retrieved from miRBase v.21 (http://www.mirbase.org/) whereas SPL sequences were retrieved from PHYTOZOME v. 9.1 (http://www.phytozome.net/) and TAIR (https://www.arabidopsis.org/). For miRNA precursors, retrieved sequences were aligned using ClustalW  using default values. Phylogenetic analyses were performed in MEGA v. 5.05  using default values. Phylogenetic inference was done using maximum-likelihood (ML) method with Bootstrap analysis (1000 trees).
The DNA sequences of SPLs were aligned using Muscle algorithm  and the well-aligned blocks were selected using Gblocks server (http://molevol.cmima.csic.es/castresana/Gblocks_server.html) by the most stringent option. Multiple sequence alignment is depicted in Additional file 7. The alignment was submitted to the estimation of proportion of invariant sites and substitution saturation test using the algorithm of Xia test implemented in DAMBE5 software . The option for the best-fit evolutionary model was performed using Akaike information criterion implemented in jModelTest . The phylogenetic reconstruction was determined by ML and Bayesian inference methods, using PhyML v3.0  and Beast v1.8.0 , respectively, the latter being implemented in CIPRES Science Gateway (https://www.phylo.org/). The approximate likelihood ratio test or aLRT  was used for ML analysis. The posterior probability estimates were calculated for Bayesian inference. The software Tracer was applied to determine the burn-in (using the log likelihood scores) in Bayesian method generation and the TreeAnnotator  to summarize the data after the exclusion of the trees that appeared outside the convergence area. The proportion of invariable sites and gamma distribution (number of categories = 4) was estimated and random local clock model for Bayesian analysis was also used.
PARE signatures mapping to OsMIR156 and OsMIR529 precursors and RNA-seq and sRNA signatures were retrieved from Rice Next-Gen sequence DBs (http://mpss.udel.edu/). Sequence abundance was estimated by normalizing all samples to TP10M (transcript per 10 million reads).
Nonsynonymous (Ka) and synonymous (Ks) substitution calculation
The selective pressure analysis (Ka/Ks) was performed using best-aligned blocks for nucleotide sequences encompassing the SBP domain plus few nucleotides upstream/downstream (named “SBP domain”). Eudicot, monocot, and bryophyte SPL sequences were carefully selected based on the phylogenetic tree depicted in Fig. 1b and phylogenetic analyses reported previously . Codons of each DNA sequence for each edited alignment were selected using an in-house Python script. Values of Ka/Ks were estimated by comparing sequences among and within eudicot, monocot, and bryophyte groups through the software KaKs_Calculator using the model selection (MS) method . Selected “target site” sequences were aligned using the Muscle algorithm, and Logos were generated using Geneious tools (http://www.geneious.com).
Based on coordinates of neighbor genes of sites of pre-miR529 in O. sativa (OsMIR529a and OsMIR529b), A. coerulea, and pre-miR156h from V. vinifera and A. thaliana, the conservation of syntenic blocks among and within monocot and eudicot species was searched in Genomicus Plants v.16.03 . For syntenic mapping of Aquilegia coerulea, we firstly used coordinates of pre-miR529 sites (scaffold_4:4,760,784..4,810,783) from Phytozome database to map flanking genes around aqc-MIR529 locus. Orthologous genes for those loci in selected eudicots and monocots were queried in Genomicus Plants v.16.03 . Phytozome database was also used to search for homologs in A. coerulea of pre-miR529 and pre-miR156h neighbor genes from O. sativa, V. vinifera, and A. thaliana, respectively.
Availability of supporting data
All published datasets referred to in the manuscript are cited in the reference list. All the supporting data are included as additional files.
AGI identifiers for Arabidopsis thaliana genes are as follows: SPL2, At5g43270; SPL3, At2g33810; SPL4, At1g53160; SPL5, At3g15270; SPL6, At1g69170; SPL9, At2g42200; SPL10, At1g27370; SPL13, At5g50570; SPL15, At3g57920; Actin-2, At3g18780.
We thank Dr. Scott Poethig for 35S::AtMIR156a seeds; Dr. Peter Huijser for spl9;spl15 seeds; Dr. Renato Vicentini for initial bioinformatic analyses and helpful discussions; and Dr. Luiz Del Bem for initial phylogenetic analyses. This work was supported by the State of Sao Paulo Research Foundation, FAPESP, Brazil (grants no. 07/58289-5 and 12/51146-2). EGOM was a recipient of a fellowship from Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil). GFFS (from Centro de Energia Nuclear na Agricultura –CENA/USP) and EMS were recipients of a fellowship from the State of Sao Paulo Research Foundation, FAPESP, Brazil.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Axtell MJ, Westholm JO, Lai EC. Vive la differénce: biogenesis and evolution of microRNAs in plants and animals. Genome Biol. 2011;2:221.View ArticleGoogle Scholar
- Allen E, Xie Z, Gustafson AM, Sung GH, Spataford JW, Carrington JC. Evolution of microRNAs genes by inverted duplication of target gene sequences in Arabidopsis thaliana. Nat Genet. 2004;36:1282–90.View ArticlePubMedGoogle Scholar
- Lu C, Tej SS, Luo S, Haudenschild CD, Meyers BC, Green PJ. Elucidation of the small RNA component of the transcriptome. Science. 2005;309:1567–9.View ArticlePubMedGoogle Scholar
- Kasschau KD, Fahlgren N, Chapman EJ, Sullivan CM, Cumbie JS, Givan SA, et al. Genome-wide profiling and analysis of Arabidopsis siRNAs. PLoS Biol. 2007;l5:e57.View ArticleGoogle Scholar
- Felippes FF, Schneeberger K, Dezulian T, Huson DH, Weigel D. Evolution of Arabidopsis thaliana microRNAs from random sequences. RNA. 2008;14:2455–9.View ArticlePubMedGoogle Scholar
- Piriyapongsa J, Jordan IK. Dual coding of siRNAs and miRNAs by plant transposable elements. RNA. 2008;14:814–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Rajagopalan R, Vaucheret H, Trejo J, Bartel DP. A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana. Genes Dev. 2006;20:3407–25.View ArticlePubMedPubMed CentralGoogle Scholar
- Creasey KM, Zhai J, Borges F, Van Ex F, Regulski M, Meyers BC, et al. miRNAs trigger widespread epigenetically activated siRNAs from transposons in Arabidopsis. Nature. 2014;508:411–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhao M, Meyers BC, Cai C, Xu W, Ma J. Evolutionary patterns and coevolutionary consequences of MIRNA genes and MicroRNA targets triggered by multiple mechanisms of genomic duplications in soybean. Plant Cell. 2015;27:546–62.View ArticlePubMedPubMed CentralGoogle Scholar
- Gu W, Wang X, Zhai C, Xie X, Zhou T. Selection on synonymous sites for increased accessibility around miRNA binding sites in plants. Mol Biol Evol. 2012;29(10):3037–44.View ArticlePubMedGoogle Scholar
- Axtell MJ, Bartel DP. Antiquity of microRNAs and their targets in land plants. Plant Cell. 2005;17:1658–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Cuperus JT, Fahlgren N, Carrington C. Evolution and functional diversification of MIRNA genes. Plant Cell. 2011;23:431–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Montes RA, de Fátima Rosas-Cárdenas F, De Paoli E, Accerbi M, Rymarquis LA, Mahalingam G, et al. Sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs. Nat Commun. 2014;5:3722.View ArticleGoogle Scholar
- Ortiz-Morea FA, Vicentini R, Silva GF, Silva EM, Carrer H, Rodrigues AP, et al. Global analysis of the sugarcane microtranscriptome reveals a nique composition of small RNA associated with axillary bud outgrowth. J Exp Bot. 2013;64:2307–20.View ArticlePubMedPubMed CentralGoogle Scholar
- Rhoades MW, Reinhart BJ, Lim LP, Burge CB, Barte B, Barte DP. Prediction of plant microRNA targets. Cell. 2002;110:513–20.View ArticlePubMedGoogle Scholar
- Klein J, Saedler H, Huijser P. A new family of DNA binding proteins includes putative transcriptional regulators of the Antirrhinum majus floral meristem identity gene SQUAMOSA. Mol Gen Genet. 1996;250:7–16.PubMedGoogle Scholar
- Schwab R, Palatnik JF, Riester M, Schommer C, Schmid ME, Weigel D. Specific effects of MicroRNAs on the plant transcriptome. Dev Cell. 2005;8:517–27.View ArticlePubMedGoogle Scholar
- Chuck G, Cigan A, Saeteurn K, Hake S. The heterochronic maize mutant Corngrass1 results from overexpression of a tandem microRNA. Nat Genet. 2007;39:544–9.View ArticlePubMedGoogle Scholar
- Wang JW, Czech B, Weigel D. miR156-regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana. Cell. 2009;138:738–49.View ArticlePubMedGoogle Scholar
- Nodine MD, Bartel DP. MicroRNAs prevent precocious gene expression and enable pattern formation during plant embryogenesis. Genes Dev. 2010;24:2678–92.View ArticlePubMedPubMed CentralGoogle Scholar
- Yu N, Cai WJ, Wang S, Shan CM, Wang LJ, Chen XY. Temporal control of trichome distribution by microRNA 156-targeted SPL genes in Arabidopsis thaliana. Plant Cell. 2010;22:2322–35.View ArticlePubMedPubMed CentralGoogle Scholar
- Gou JY, Felippes FF, Liu CJ, Weigel D, Wang JW. Negative regulation of anthocyanin biosynthesis in Arabidopsis by a miR156-targeted SPL transcription factor. Plant Cell. 2011;23:1512–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Silva GFF, Silva EM, Azevedo MS, Guivin MAC, Ramiro DA, Figueiredo CR, et al. microRNA156-targeted SPL/SBP box transcription factors regulate tomato ovary and fruit development. Plant J. 2014;78:604–18.View ArticlePubMedGoogle Scholar
- Ling LZ, Zhang SD. Exploring the evolutionary differences of SBP-box genes targeted by miR156 and miR529 in plants. Genetica. 2012;140:317–24.View ArticlePubMedGoogle Scholar
- Zhang SD, Ling LZ, Zhang QF, Xu JD, Cheng L. Evolutionary comparison of two combinatorial regulators of SBP-Box genes, MiR156 and MiR529, in plants. PLoS One. 2015;10:e0124621.View ArticlePubMedPubMed CentralGoogle Scholar
- Jeong DH, Park S, Zhai J, Gurazada SG, De Paoli E, Meyers BC, et al. Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell. 2011;23(12):4185–207.View ArticlePubMedPubMed CentralGoogle Scholar
- Arif MA, Fattash I, Ma Z, Cho SH, Beike AK, Reski R, et al. DICER-LIKE3 activity in Physcomitrella patens DICER-LIKE4 mutants causes severe developmental dysfunction and sterility. Mol Plant. 2012;5:1281–94.View ArticlePubMedPubMed CentralGoogle Scholar
- Dai X, Zhao PX. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011;39(Web Server issue):W155–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Schwarz S, Grande AV, Bujdoso N, Saedler H, Huijser P. The microRNA regulated SBP-box genes SPL9 and SPL15 control shoot maturation in Arabidopsis. Plant Mol Biol. 2008;67((1-2)):183–95. doi:10.1007/s11103-008-9310-z.View ArticlePubMedPubMed CentralGoogle Scholar
- Preston JC, Hileman LC. Functional evolution in the plant SQUAMOSA-PROMOTER BINDING-LIKE (SPL) gene family. Front Plant Sci. 2013;4(80):1–13.Google Scholar
- Fahlgren N, Jogdeo S, Kasschau KD, Sullivan CM, Chapman EJ, Laubinger S, et al. MicroRNA gene evolution in Arabidopsis lyrata and Arabidopsis thaliana. Plant Cell. 2010;22(4):1074–89.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen K, Rajewsky N. The evolution of gene regulation by transcription factors and microRNAs. Nat Rev Genet. 2007;8:93–103.View ArticlePubMedGoogle Scholar
- Taylor RS, Tarver JE, Hiscock SJ, Donoghue PCJ. Evolutionary history of plant microRNAs. Trends Plant Sci. 2014;19:175–82.View ArticlePubMedGoogle Scholar
- Puzey JR, Kramer EM. Identification of conserved Aquilegia corulea microRNAs and their targets. Gene. 2009;448:46–56.View ArticlePubMedGoogle Scholar
- Odell JT, Nagy F, Chua N-H. Identification of DNA sequences required for activity of the cauliflower mosaic virus-35S promoter. Nature. 1985;313:810–2.View ArticlePubMedGoogle Scholar
- Barakat A, Wall K, Leebens-Mack J, Wang YJ, Carlson JE, Depamphilis CW. Large-scale identification of microRNAs from a basal eudicot (Eschscholzia californica) and conservation in flowering plants. Plant J. 2007;51:991–1003.View ArticlePubMedGoogle Scholar
- Jeong DH, Szchmidt SA, Rymarquis LA, Park S, Ganssmann M, German MA, et al. Parallel analysis of RNA ends enhances global investigation of microRNAs and target RNAs of Brachypodium distachyon. Genome Biol. 2013;24:R145.View ArticleGoogle Scholar
- Xie K, Shen J, Hou X, Yao J, Li X, Xiao J, et al. Gradual increase of miR156 regulates temporal expressin changes of numerous genes during leaf development in rice. Plant Physiol. 2012;158:1382–94.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang JW, Schwab R, Czech B, Mica E, Weigel D. Dual effects of miR156-targeted SPL genes and CYP78A5/KLUH on plastochron length and organ size in Arabidopsis thaliana. Plant Cell. 2008;20:1231–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol. 2008;18:758–62.View ArticlePubMedPubMed CentralGoogle Scholar
- Meng Y, Gou L, Chen D, Wu P, Chen M. High-throughput degradome sequencing can be used to gain insights into microRNA precursor metabolism. J Exp Bot. 2010;61:3833–7.View ArticlePubMedGoogle Scholar
- Jones-Rhoades M, Bartel D. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14:787–99.View ArticlePubMedGoogle Scholar
- Poethig RS. Vegetative phase change and shoot maturation in plants. Curr Top Dev Biol. 2013;105:125–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Chuck G, Whipple C, Jackson D, Hake S. The maize SBP-box transcription factor encoded by tasselsheath4 regulates bract development and the establishment of meristem boundaries. Development. 2010;137:1243–50.View ArticlePubMedGoogle Scholar
- Jeong DH, Thatcher SR, Brown RS, Zhai J, Park S, Rymarquis LA, et al. Comprehensive investigation of microRNAs enhanced by analysis of sequence variants, expression patterns, ARGONAUTE loading, and target cleavage. Plant Physiol. 2013;162:1225–45.View ArticlePubMedPubMed CentralGoogle Scholar
- Clough SJ, Bent AF. Floral dip: a simplified method for Agrobacterium-mediated transformation of Arabidopsis thaliana. Plant J. 1998;16:735–43.View ArticlePubMedGoogle Scholar
- Karimi M, Inzé D, Depicker A. GATEWAY vectors for Agrobacterium-mediated plant transformation. Trends Plant Sci. 2002;7:193–5.View ArticlePubMedGoogle Scholar
- Varkonyi-Gasic E, Wu R, Wood M, Walton EF, Hellens RP. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:1–12.View ArticleGoogle Scholar
- Javelle M, Timmermans MCP. In situ localization of small RNAs in plants by using LNA probes. Nat Protoc. 2012;7:533–41.View ArticlePubMedGoogle Scholar
- Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23:2947–8.View ArticlePubMedGoogle Scholar
- Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28(10):2731–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Xia X. DAMBE5: a comprehensive software package for data analysis in molecular biology and evolution. Mol Biol Evol. 2013;30:1720–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.View ArticlePubMedPubMed CentralGoogle Scholar
- Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systs Biol. 2010;59(3):307–21.View ArticleGoogle Scholar
- Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Anisimova M, Gascuel O. Approximate likelihood-ratio test for branches: a fast, accurate, and powerful alternative. Systs Biol. 2006;55:539–52.View ArticleGoogle Scholar
- Zhang Z, Li JL, Zhao X-Q, Wang J, Wong GK-S, Yu J. KaKs calculator: calculating Ka and Ks through model selection and model averaging. Genomics Proteomics Bioinfomatics. 2006;4:259–63.View ArticleGoogle Scholar
- Louis A, Nguyen NTT, Muffato M, Crollius HR. Genomicus update 2015: KaryoView and MatrixView provide a genome-wide perspective to multispecies comparative genomics. Nucleic Acids Res. 2015;43:D682–9. doi:10.1093/nar/gku1112.View ArticlePubMedGoogle Scholar
- Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003;31:3406–15.View ArticlePubMedPubMed CentralGoogle Scholar