Altered expression profiles of microRNA families during de-etiolation of maize and rice leaves

MicroRNAs (miRNAs) are highly conserved small non-coding RNAs that play important regulatory roles in plants. Although many miRNA families are sequentially and functionally conserved across plant kingdoms (Dezulian et al. in Genome Biol 13, 2005), they still differ in many aspects such as family size, average length, genomic loci etc. (Unver et al. in Int J Plant Genomics, 2009). In this study, we investigated changes of miRNA expression profiles during greening process of etiolated seedlings of Oryza sativa (C3) and Zea mays (C4) to explore conserved and species-specific characteristics of miRNAs between these two species. Futhermore, we predicted 47 and 42 candidate novel miRNAs using parameterized monocot specific miRDeep2 pipeline in maize and rice respectively. Potential targets of miRNAs comprising both mRNA and long non-coding RNA (lncRNA) were examined to clarify potential regulation of photosynthesis. Based on our result, two putative positive Kranz regulators reported by Wang et al. (2010) were predicted as potential targets of miR156. A few photosynthesis related genes such as sulfate adenylytransferase (APS3), chlorophyll a/b binding family protein etc. were suggested to be regulated by miRNAs. However, no C4 shuttle genes were predicted to be direct targets of either known or candidate novel miRNAs. This study provided the comprehensive list of miRNA that showed altered expression during the de-etiolation process and a number of candidate miRNAs that might play regulatory roles in C3 and C4 photosynthesis.

Background miRNAs are about 21-nt long non-coding RNAs that regulate gene expression by binding to complementary sequences within their target mRNAs [3]. Although miR-NAs only occupy about 0.01% of total RNA mass, their average copy numbers per cell were estimated to be higher than mRNAs [4]. The majority of verified miRNA binding targets are transcription factors [5]. Many evidences suggest that miRNAs play important roles during plant growth and development, see recent reviews by Jones-Rhoades et al. [6], Eldem et al. [7], Ameres and Zamore [8]. In addition, miRNAs were also reported to be involved with many other aspects such as siRNA biogenesis [9], signal transduction [10], plant disease [11] and environmental stress responses [12,13].
Plant miRNAs differ from animal miRNAs in many aspects. Firstly, plant miRNA genes are often located in intergenic regions while animal miRNA genes are often located in introns or coding sequences [14]. Secondly, plant miRNA genes are generally monocistronic while animal miRNA genes are more often clustered together [15]. Thirdly, animal miRNAs tend to bind to 3′UTR regions of target transcripts while plant miRNAs follow more strict reverse complementary match with targets thus possessing limited number of targets compared with animal miRNAs [15]. Fourthly, plant and animal miR-NAs also have different lengths and stabilities of precursors and involve different enzymes during the process of biogenesis and regulation [2,3]. Furthermore, monocot and dicot miRNAs also differ in statistical features such as minimal free energies and stabilities of miRNA precursors [16].
Plant miRNAs are quite conserved across species. Through large scale survey of available sequences, expressed sequenced tags, and nonredundant nucleotides, a large number of miRNA families have been identified across different plant species, e.g. 15 conseved miRNA families in 11 plant species, at least 21 miRNA families conserved across monocot and dicots [17]. Zhang et al. [18] reported that they have identified 481 miRNAs that belong to 37 miRNA families across 71 different plant species. Based on miRBase release version 21 [19], 28 conserved miRNA families were identified between maize and rice while both species possess species-specific miRNAs. Plant miRNAs were reported to be involved in the regulation of many different developmental and metabolic processes, including responses of plants to different environmental conditions, such as temperature [12], abiotic stress [20] and salt stress [21,22]. Examination of variations of categories and abundances of miRNAs under different environmental conditions is a common strategy to gain insights regarding the role of miRNAs during plant growth and development.
Light is a major environmental signal influencing many aspects of plant growth and development. Many earlier studies demonstrated that a large number of genes including transcription factors show dramatic changes after light induction [23,24]. Some miRNAs have been implied to be involved in this process, e.g. LONG 14 HYPOCOTYL 5 (HY5), a global regulator of light responsive transcription, was reported to regulate expressions of eight miRNA genes in Arabidopsis [25]. One major plant biological process heavily influenced by light is photosynthesis, which was suggested to be under the regulation of miRNAs [26]. Besides this possible role of miRNAs in regulating C 3 photosynthesis in general, it is of special interest that miRNAs might also be involved in the regulation of C 4 photosynthesis. The cell specific expression of C 4 related proteins or enzymes is under multiple layers of regulation, see recent reviews on this by Hibberd and Covshoff [27] and Williams et al. [28]. Apart from cisregulatory motifs that have been identified to play important role in controlling C 4 specific expression patterns [29], post-transcriptional control, post-translational control and epigenetic control might also contribute to the cell specific accumulations of C 4 enzymes, see a recent review by Wang et al. [30]. Though miRNAs are regarded as highly possible signals controlling the expression of C 4 related genes, so far, experimental evidences suggesting such an association have not been established.
The present study explores the dynamic changes of expression levels of both known and candidate novel miRNAs and their predicted target genes during the greening of etiolated Zea mays (C 4 ) and Oryza sativa (C 3 ) leaves. We examined the possibilities of miRNAs involving with regulation of C 4 photosynthesis. To aid this analysis, we parameterized a monocot specific miR-NAs prediction pipeline. Our results indicated that although the majority of miRNA families are sequentially and functionally conserved between maize and rice, their expression profiles in response to light induction during de-etiolation process are quite different. Potential mRNA targets and mimic lncRNA targets of both known and candidate novel miRNAs were predicted. By checking the overlap of predicted targets with C 4 cycle genes, reported regulators of C 4 trait and photosynthesis related genes, we identified potential miRNAs that might be related to C 4 photosynthesis.

Deep sequencing of small RNAs
In this study, about 85.8 and 87.5 million reads were sequenced from maize and rice small RNA libraries, which yielded about 24.2 and 19.2 million collapsed identical reads (termed tags) for maize and rice, respectively. After quality filtration, 3′ adapters were trimmed for each sample. The length distribution of the trimmed reads and tags showed that the most abundant length regions are between 20-nt and 24-nt ( Fig. 1). No significant difference in length distribution between maize and rice trimmed reads was detected. Contaminant small RNA reads were removed by BLAST against rRNA, tRNA and snRNA databases (see "Methods" section for details). For each sample, numbers of reads and tags that were aligned with other small RNAs were listed in Additional file 1 and plotted in Additional file 2. In general, no significant variations between maize and rice were found in categories or amounts of identified rRNAs, tRNAs and snR-NAs. Retained reads were then mapped against genomes for further analysis.

Monocot specific miRDeep2
Plant miRNAs differ from animal miRNAs in many aspects such as biogenesis [2], precursor length [31], target binding patterns [3]. Besides, monocot and dicot plants also differ in statistical features [16]. We therefore parameterized a widely used tool miRDeep2 [32] to be monocot specific in this study. Several parameters and scoring functions were replaced in the original scripts. The modified pipeline could be used in the same way as the original one and it could be accessed via https:// github.com/Rossifumi/monocot_specific_miRDeep2.git.

Prediction of novel miRNAs in rice and maize
In total, 288 known maize miRNAs and 658 known rice miRNAs were detected in our samples. By applying the modified miRDeep2, we predicted 111 and 141 potential candidate novel miRNAs for maize and rice, respectively. Genomic loci of the precursors of these potential candidate novel miRNAs were listed in Additional file 3. Conservations of mature sequences of the predicted candidate novel miRNAs were checked across six plant species, i.e., Arabidopsis thaliana, O. sativa, Populus trichocarpa, Triticum aestivum, Vitis vinifera, Z. mays. When setting the similarity cutoff to be 90%, 6 maize and 4 rice potential miRNAs showed hits with at least one known miRNA mature sequence in miRBase (Additional file 4). Those predicted miRNAs were further filtered with expression cutoffs. A potential candidate novel miRNA was considered to be candidate novel miRNA and retained for further study only when its total TPQ (transcripts per quarter million) value in all samples together exceeded 10 (Additional file 5). Length, mature sequence, genomic loci of precursor sequence and TPQ value across samples of 47 maize candidate novel miR-NAs and 42 rice candidate novel miRNAs were identified (Additional file 6).
Prior to this study, a number of genome-wide predictions of novel miRNAs in maize and rice have been conducted. For example, Wang et al. [33] predicted 167 novel miRNAs in imbibed maize seeds. Jiao et al. [34] reported 66 novel miRNAs in mixed tissues, embryo and endosperm of maize. Kang et al. [35] identified 54 novel miRNAs in developing seeds and growing leaves of maize. Sunkar et al. [36] reported 63 novel miRNAs in rice seedlings and seedlings exposed to drought and salt stresses. Here we compared our predictions with these previous studies. When setting the similarity cutoff to be 90%, 14 out of 47 maize candidate novel miRNAs were able to match with previous predictions (Additional file 7). However, no similarity was identified between 42 rice candidate novel miRNAs and the predictions from Sunkar et al. [36] at the same similarity cutoff. When lowering the similarity cutoff, we were able to identify a few matched pairs (Additional file 7). However, similarity cutoff lower than 90% is not recommanded because it would allow too many mismatches. Thus matched miRNA pairs might be quite distinguished from each other.
The length distributions of known and candidate novel miRNAs ( Fig. 2) show that although the length of the most abundant mature sequence for both maize and rice miRNAs is 21-nt, rice miRNAs however are usually slightly longer than maize miRNAs. Taking into account that the majority of plant miRNAs regulate their target mRNAs by reverse complementary match [37], the increased mature sequence length of miRNAs in rice indicate that rice miRNAs might have higher specificity in terms of target binding and thus affecting regulation.

miRNA families and their altered expression profiles during de-etiolation
Many plant miRNA families possess multiple family members within the same species [38]. Different members from the same family are highly conserved in terms of mature sequence within or across species. However, the precursors of those miRNAs might differ at not only precursor sequence level but also their genomic loci [18]. In our study, 25 multi-member miRNA families and 3 solo-member miRNA families were detected in maize. Their homolog miRNAs were also identified in rice. Apart from these 28 conserved miRNA families between maize and rice, 30 more multi-member and 133 more solo-member miRNA families were detected in rice as well. Thus, the number of members in each miRNA family was larger in maize than rice (about 11 members per miRNA family for maize, 2 members per miRNA family for rice). Compared to rice, the majority of maize miR-NAs were conserved across species while rice possess more species-specific miRNAs, according to miRBase release 21 [19].
However, for the 28 conserved miRNA families between maize and rice, more miRNA members were detected in rice than maize (Fig. 3). Although members of plant miRNA families rarely form clusters [39], in contrast to animal miRNAs which usually cluster together [15], the conservation of genomic loci of miRNA families were also detected in plants. This is reflected in the presence of at least one member of 22 out of 28 miRNA families residing within maize and rice syntenic gene blocks (Fig. 3). More interestingly, miRNA families that form clusters in maize, such as miR2118, miR395 and miR159, tend to form clusters in rice as well, indicating that those miRNA families already existed in the last common ancestor of maize and rice.
The majority of members of the same miRNA families usually share similar mature sequences and similar regulatory roles, although a few exceptions of members of the same miRNA family with different functions were detected as well [35]. Due to sequence similarity of mature miRNAs, it is difficult to precisely detect the expression profiles of different members of the same miRNA families. To maximun the reliability of expression levels, family-wise expression profiles of 28 conserved miRNA families in maize and rice were provided in Additional file 8. Comparing the fold changes of expression values of etiolated leaves to control samples showed that although the majority of miRNA families had relatively slight changes during de-etiolation, some of them showed rather dramatic expression changes (Fig. 4).
miRNAs families with different expression patterns between maize and rice were further examined. Among miRNA families with considerable expression levels, some of them were not only differentially expressed during de-etiolation process between maize and rice but also differentially expressed between control sample and 24-h illuminated sample, i.e., their expression levels did not converge to normal condition after 24-h illumination. To be more specific, miR167 family, repressors of auxin responsive factor 8 (ARF8) [40], was expressed higher in rice than maize during de-etiolation process (Additional file 8). Besides, after 24-h illumination, the expression level of miR167 was still lower than control sample in maize, but significantly higher than control sample in rice (Fig. 4). Similarly, expression levels of miR159 family, regulators of myeloblastosis (MYB) genes [41], and miR166 family, regulators of Homeodomain-Leucine zipper (HD-ZIP) transcription factor families [42], were also lower than control sample during de-etiolation process in maize, but higher than control sample in rice (Additional file 8). On the contrary, miR528 family, regulators of cupredoxin domains involving oxidative stress responses [43], was expressed higher in maize than rice during deetiolation process (Additional file 8). Similarly, miR827 family, verified regulators of SPX (named after SYG1/ PHO81/XPR1)-major facilitator superfamily (SPX-MFS) genes [44], was expressed during de-etiolation in maize but not in rice (Additional file 8). Different from miRNA families discussed above, expression levels of miR164 family, reported regulators of the CUP-SHAPED COTY-LEDON (CUC) genes [45], and miR390 family, regulators Color codes stand for different miRNA families conserved between maize and rice. Grey lines linked a pair of syntenic gene blocks between these two species. Color lines indicated that the pair of corresponding miRNA genes resides between these two gene blocks of lateral root development by cleaving TRANS-ACTING SIRNA3 (TAS3) precursor RNA [46], converged to normal condition after 24-h illumination although different expression patterns between maize and rice were detected during de-etiolation process (Additional file 8, Fig. 4).

mRNA targets and mimic lncRNA targets of miRNAs
The majority of plant miRNAs bind to its target mRNA through perfect reverse complementary matching [37]. A number of other factors also influence target binding, such as target accessibility, minimal free energy and secondary structure of miRNA-target pair. As a result, all these factors need to be considered when their binding targets are predicted [47]. In this study, we identified 47 and 42 candidate novel miRNAs for maize and rice respectively (Additional file 6). The mRNA targets of both known miRNAs and candidate novel miRNAs were predicted by psRNATarget [48] and listed in Additional file 9. We compared the obtained correlation coefficients between miRNAs and their binding targets. More rice miRNA-target pairs showed negative correlation coefficients than maize (Additional file 10). We further calculated the correlation coefficients between the expression patterns of miRNAs and their binding targets given by PTMED, a plant miRNA-target database based on microarray datasets [49]. Similar results were obtained, i.e., rice showed more negative correlations than maize.
Some of the correlation coefficients between expression values of miRNAs and their predicted targets were not negative. One possible explanation is that miR-NAs can bind to not only mRNAs but also lncRNAs as target mimicry and thus affecting the miRNA-target regulation [50]. Given that a few previous studies have predicted genome-wide lncRNAs for maize [51,52], here we tested whether those previously reported lncRNAs can serve as mimic targets of maize miRNAs (Additional file 11). In addition, by applying an in-house pipeline (see "Methods" section for details), we also identified 25 conserved lncRNAs precursors (13 from maize and 12 from rice genome). Two of the 12 predicted rice lncRNAs can potentially act as mimic target for two rice miRNAs (Additional file 12).

The binding of miRNA to C 4 photosynthesis related genes
To examine if miRNAs could be potentially involved with the regulation of photosynthesis, we examined the potential miRNAs that target photosynthesis related genes. We developed the list of photosynthesis related genes from a number of sources: (a) genes classified into photosynthesis category in MapMan annotation in maize and rice; (b) genes differentially expressed in bundle sheath cell and mesophyll cell in maize [53][54][55]; (c) rice orthologs of cell specific accumulated maize genes. miRNAs that predicted to have binding targets belonging to the list of photosynthesis related genes were listed in Table 1.
Our anlaysis did not find any C 4 cycle gene being a direct target of miRNA. However, a number of C 4 photosynthesis related genes were predicted to be targets of miRNAs in both maize and rice. MiR395 has been reported to be crucial for sulfate homeostasis during growth and development in Arabidopsis [56]. Our data showed that miR395 regulates sulfate adenylytransferase (APS3) in both maize (GRMZM2G051270, GRMZM2G149952) and rice (LOC_Os03g53230), suggesting that this regulatory mechanism is conserved between maize and rice. miR1432 was predicted to target calmodulin-binding protein (CaMBP) and EF-hand protein in rice [36]. Our data indicate that it might also target a chlorophyll a/b binding protein (GRMZM2G131489) in maize (Table 1). Chlorophyll a/b binding protein (LOC_Os03g39610) in rice was also predicted to be target of a candidate novel miRNA os_03_07 in rice. No sequence similarity was detected between miR1432 (CUCAGGAGAGAUGACACCGAC) and os_03_07 (AAUGACUUACAUUGUGGAACGGAG) mature sequences. This suggestes the mechanisms of chlorophyll a/b binding protein regulation might have been evolved through different routes between maize and rice.
Beta-CA5 (GRMZM2G145101) was predicted to be a target of miR528 in maize. Notably this is not the beta-CA gene (GRMZM2G121878) that was recruited into C 4 cycle in maize based on cell specific and leaf gradient data [53,57]. miR528 was reported to regulate seed development in rice and response to drought/waterlogging stress in maize and rice [43,58]. Thus the regulation of miR528 to beta-CA is more likely a response to environmental stress, in this case, a sudden exposure to light during de-etiolation process.
In addition, we also compared target gene lists of miRNAs with previously identified transcription factors, which might serve as putative positive or negative regulators of the development of Kranz anatomy [55]. Analysis showed that two target genes of miR156 (GRMZM2G097275 and GRMZM2G126018) were reported to be putative positive Kranz regulators [47].

Experimental validation of predicted miRNA-target pairs
To validate our prediction of novel miRNAs and their regulation of target genes, we constructed short tandem targe mimic (STTM) vector sequences for four miRNAtarget gene pairs, i.e., osa-miR395b-Os03g53230, osa-miR444a-Os05g47560, osa-miR6251-Os03g532230, Os_03_07-Os03g39610, and verified expression levels of both miRNAs and target genes among different transgenic lines in O. sativa japonica (Fig. 5). STTM sequence and detection primers were illustrated in Fig. 6.
Among four pairs we chose, two of them showed relative significant correlation between miRNA and predicted target gene with p value equals 0.05 and 0.06 indicating they might be ture miRNA-target pairs (Fig. 5a, b). miR444a, regulator of STN7, chloroplast serine-threonine protein kinase, showed negative correlation with its predicted target Os05g47560, whose orthologous gene in maize showed enriched expression in mesophyll cells in leaves ( Fig. 5b; Table 1). While miR395b, regulator of sulfate adenylytransferase, showed positive correlation with its predicted target Os03g53230, whose maize ortholog showed enriched espression in bundle sheeth cell in leaves. Although the majority of known miRNAs repress gene expression by complementary mapping to their targets, miRNAs that positvely regulate target gene expression have been reported in other species as well [59,60]. Our finding suggested this mechanism might also exist in rice.
Os_03_07, a novel miRNA predicted by our monocot specific pipeline, showed a negative correlation with its predicted target Os03g39610, a chlorophyll a/b binding protein (Fig. 5d). Other miRNA-target gene pairs predicted in this study as showed in Table 1 is still under experimental verfication.

The de-etiolation system as a model to study regulation of C 4 photosynthesis
We chose the greening process of the etiolated leaves for this study for the following seasons. First, light not only provides energy for photosynthesis, but also function as an environmental signal during development of photosynthesis [61]. Many photosynthetic genes showed altered expression during light induction [62] together with changes of many other genes, which makes it possible to use data collected during the de-etiolation process to predict miRNA-gene regulation pairs. During the de-etiolation process, we sampled 7 time points, i.e., 0, 0.5, 1, 3, 6, 12 and 24 h, following some previous studies [63][64][65]. Time interval was shorter during the beginning phase of de-etiolation process to better capture the initial responses of de-etiolated leaves to light illumination.

Substantial level of conservation exists between maize and rice miRNAs
miRNAs serve as important regulators during plant growth and development [7]. There is substantial level of conservation in plant miRNAs. The existence of one of the most conserved miRNA family, miR166, was estimated to date back to at least 400 million years old [66]. More families of miRNAs were conserved in plants than animals [17]. This study showed a number of aspects of miRNAs are conserved between maize and rice. Firstly, maize and rice share 28 conserved miRNA families according to miRBase release 21 [19]. The majority of these families share not only similar mature sequences but also similar biological functions [58,67,68]. Many of them, i.e., miR156, miR160, miR164, miR166, miR167, miR171, miR172 and miR390, were suggested to play highly evolutionary conserved roles across plant species [54]. Secondly, in both maize and rice, the majority of plant miRNA targets are genes coding transcription factors rather than protein coding genes [5]. Thirdly, at least one member of 22 out of 28 conserved miRNA families were identified within maize and rice syntenic gene blocks, i.e., conserved blocks of order of genes within two sets of chromosomes that are being compared (Fig. 3), which suggested that they were inherited from the same ancestor.
For some miRNA families, their members form clusters together on the chromosome (Fig. 3). Though in both maize and rice, the majority of members of the same miRNA family scattered across different chromosomes, which is in contrast to animals where members of the same miRNA family tend to cluster together [69]. The fact that the same miRNA families of which their members form clusters on maize chromosomes also form clusters on rice chromosomes indicated the common evolutionary origins of these miRNA families between maize and rice (Fig. 3). In contrast, the scattered distributions of miRNA family members indicated that these members might have experienced fragment rearrangements possibly due to transposon activity during evolution.

Diverse features between rice and maize miRNAs
Maize and rice miRNAs show many differences. Firstly, family sizes of maize miRNAs are much larger than those of rice miRNA families. On average, a maize miRNA family possesses about 11 members while a typical rice  [19]). The majority of identified maize miRNAs families are conserved across species while rice possess more species-specific miRNA families (miRBase release 21; Griffiths-Jones et al. [19]). Secondly, the average length of mature miRNA sequences is also 0.7-nt longer in rice than maize (Fig. 2). Rice possesses more 24-nt long miRNAs than maize (Fig. 2). Chen et al. [70] also reported a dominant portion of small RNAs of rice fall into 21-nt to 24-nt region, especially 24-nt. Montes et al. [71] reported that the majority of plant 24-nt sRNAs are heterochromatic siRNAs. Thus, one possible explanation might be some 24-nt siRNAs are mis-identified as miRNAs in rice as previous studies reported that 24-nt siRNAs are quite active in gene regulation in rice [72][73][74]. Given that the majority of plant miRNAs bind to targets by complementary reverse matching [37], much more targets were predicted in maize than in rice since maize miRNAs are usually shorter (Additional file 9). Thirdly, the correlations between miRNA expression profiles and target gene expression profiles are better in rice than maize (Additional file 10) suggesting a more sophisticated miRNA regulation network in rice. For those miRNAs that showed similar expression patterns between maize and rice, i.e., miR156, miR166, miR168, miR172, miR2275 and miR528, GO enrichment analysis of their predicted targets was applied (Additional file 13). Our analysis showed that their targets enriched in much less number of biological processes in maize compared to those in rice (Additional file 13), which also suggests that miRNAs might play more sophisticated regulatory roles in rice compared to maize.
It is worth mentioning that in the current miRBase, there are more rice miRNAs in miRBase than maize miR-NAs, likely due to the larger number of publications on rice compared to maize. Interestingly, though a numbers of genome-wide identification of novel miRNAs have been conducted in maize [33][34][35]75], the number of newly identified maize miRNAs did not dramatically increase the number of miRNAs. Thus, though maize genome (~2500 MB) is much larger than rice genome (~390 MB), maize might not necessarily possess more miRNAs than rice.

Some miRNA families conserved between maize and rice showed drastically different expression patterns during de-etiolation
Although the majority of conserved miRNA families across species show similar functions [58,67], their transcriptional responses during de-etiolation can be quite different. Here we used one example to illustrate this point. Earlier study showed that HY5, a global regulator of light responsive transcription factor, regulates expressions of eight miRNA genes, i.e., miR156d, miR172b, miR402, miR408, miR775, miR858, miR869 and miR1888, in Arabidopsis [25]. Among these 8 miRNAs, miR156, miR172 and miR408 are conserved miRNA families between maize and rice. Though miR156, regulators of flowering time, phase changing modulation, later embryonic maturation and root development [76], together with miR168, regulators of stress responses and signal transductions in plant development [67], showed very similar expression patterns between maize and rice during deetiolation (Additional file 14). miR172, regulator of seed development and phase change in shoot [67], showed drastic difference in responses during de-etiolation, i.e., it was barely detectable in maize while constantly expressed in rice during de-etiolation (Additional file 8).
Dramatically different expresison patterns have also been detected in other miRNA families. miR156, miR160, miR164, miR166, miR167, miR171, miR172, and miR390, had been earlier reported to play evolutionarily conserved roles in plant development [54]. Though all these miRNA families share functions related to flowering control, embryonic maturation, root development, leaf primordial and development etc. [58,67,68], most of them showed different expression patterns upon exposure to light (Additional file 14) except (a) miR156, miR166, miR172, which showed almost identical expression curves, and (b) miR171 and miR390, which showed shifted expression patterns.

Potential relevance of miRNAs and regulation of C 4 photosynthesis
Given that the majority of miRNA targets are transcription factors genes that might further regulate gene expressions of other genes [5], it is rational to suggest that miRNAs, as upstream regulators, have the potential to function as master switches of some biological processes. Previous studies have reported that some miR-NAs can move between different cells and cause cell specific accumulations [77], which raises the possibility that miRNAs might be related to cell specific expression of C 4 enzymes. The large family size of plant miRNA families and relative frequent duplication events in plant genomes can potentially provide the conditions necessary for the evolution of neo-functions of some redundant miRNAs. Thus, we explored the predicted targets of both known and candidate novel miRNAs to investigate if miRNAs might be potentially associated with regulation of C 4 photosynthesis.
Based on our results, no C 4 cycle genes were identified to be direct target of known or predicted miRNAs (Additional file 9). However, a few photosynthesis related genes were predicted as targets of miRNAs (Table 1). APS3 genes, i.e. GRMZM2G051270 and GRMZM2G149952 for maize and LOC_Os03g53230 for rice, were predicted to be binding targets of miR395 (Table 1). These two maize genes showed cell specific accumulations in maize according to Li et al. [45]. The similar regulation patterns of miR395 to APS3 genes between maize and rice indicated that this regulatory mechanism might have existed long before the speciation of maize and rice, i.e. although APS3 genes showed cell specific accumulation patterns in maize, which should not be direclty linked to C 4 photosynthesis. Beta-CA5 (GRMZM2G145101) was also predicted to be a binding target of miR528 in maize (Table 1). However, it is not the beta-CA gene (GRMZM2G121878) that was recruited into C 4 photosynthesis in maize based on cell specific and leaf gradient data [53,57]. Thus, these genes predicted to be targets of miRNAs are considered being peripheral to C 4 photosynthesis. In fact, out of the total number of predicted miRNA binding targets, the photosynthesis related genes only represented a small fraction of them (Table 1). This indicates that miRNAs might have played a relatively minor role in the regulation of photosynthesis.
We further explored whether miRNAs might be involved in regulation of C 4 photosynthesis by regulating C 4 related transcription factors. We compared our target gene lists with putative positive or negative regulators of the development of Kranz anatomy reported by Wang et al. [47]. Our results showed that two target genes of miR156 (GRMZM2G097275 and GRMZM2G126018) were reported to be putative positive Kranz regulators (Additional file 15). Considering that miR156 is a highly conserved miRNA across plant species [76], we conducted GO enrichment anlaysis of all the predicted 24 and 21 target genes for maize and rice miR156 respectively. Analysis showed that the function of target gene lists is enriched with nucleus related functions in both species (Additional file 16). The functional significance of this miRNA and the corresponding binding targets in C 4 Kranz anatomy formatoin needs to be tested next (Additional file 17).

Conclusions
In this study, we examined the transcriptomic changes of miRNAs and their potential targets at 7 time points during the greening process of de-etiolated leaves in both maize and rice. No biological replicates were available in this study; instead we used 7 time points during the deetiolation process to facilitate identification of miRNAtarget gene pairs. We have parameterized a monocot specific prediction pipeline and identified 47 and 42 candidate novel miRNAs in maize and rice respectively. Further analysis showed that miRNA families are quite conserved between maize and rice in terms of mature sequence and genomic loci. However, there are also some different features between these miRNAs in two species such as the average length of mature sequence, family size and changes of expression profiles during de-etiolation. Our study further indicated that rice might possess a more complex miRNA regulation network than maize. Analysis of potential targets of both known and candidate novel miRNAs did not find any miRNA that direclty binds to C 4 cycle genes. However indirect associations between miRNAs and C 4 trait development related TFs were identified. Some of these identified miRNA and their binding targets need to be experimentally tested to examine their role during C 4 photosynthesis development.

Plant materials, RNA isolations, miRNA and mRNA sequencing
Zea mays L. ecotype B73 and O. sativa japonica seeds were sown and cultured in soil in dark condition at 28/22 °C and 60% humidity for 1 week. The 7-day-old etiolated seedlings were exposed to continuous light (~200 μmol/m 2 /s at the surface of the sampled leaves) and illuminated for 24 h. Seeds for control experiments were sown and cultured in soil for 1 week with a 16-h (7:00-23:00) light/8-h night cycle. Leaf sections of about 2 cm length were harvested from the third leaf at about top one-third of the position from tip to base. Samples were harvested before illumination (termed 0 h, at 9:00 a.m.) and six other time points after illumination, i.e., 0.5, 1, 3, 6, 12 and 24 h. These samples were immediately frozen in liquid nitrogen and stored at −80 °C until use. Total RNA was extracted with TRIzol ® protocol. The RNA integrity was evaluated by agarose gel electrophoresis and the concentration was checked using a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific

Modification of miRDeep2 pipeline and prediction of novel miRNAs
A well-accepted miRNAs prediction package, miRDeep2 [32], was parameterized to be monocot specific. The altered parameters and scoring functions were described in our previous study [16]. The modified version of miRDeep2 is accessible via our website (http://www. picb.ac.cn/PSB/a/DOWNLOAD/). References miRNAs sequences required by miRDeep2 were downloaded from miRBase release 21, June 2014 (http://www.mirbase.org) [78].
We named candidate novel miRNAs by the combination of short forms of species names, chromosome number of precursor and order of candidate novel miRNAs discovered on that chromosome in this paper.

miRNAs expression profiles, family characterization and genomic syntenies
Normalized TPQ (transcripts per quarter million) values were generated for both known and candidate novel miR-NAs. Expression levels of miRNA families were determined by the sum of TPQ values of all family members. Differentially expressed miRNAs were detect by R package edgeR [80]. Expression curves were 3rd polynomial regressed to better illustrate time-dependent expression of miRNAs. This function has been used earlier to describe time-series dataset analysis [81,82]. Maize and rice orthologous gene pairs were identified by BLASTX with maximum e-value as e-10. DAGchainer [83] was used to identify collinear regions amongst those orthologous genes. Two regions are considered to be collinear when at least five orthologous genes were found in the same order with no more than ten genes inserted between each neighbors. Family members and genomic syntenies of miRNAs were plotted by CIRCOS version 0.62 [84]. Two miRNA genes were considered syntenic if they belong to the same miRNA family and located within syntenic collinear regions. miRNA targets prediction, function analysis and lncRNA prediction miRNA targets for known and candidate novel miRNAs were predicted by psRNATarget [48] (http://plantgrn. noble.org/psRNATarget/) using default parameters. GO enrichment analysis of the identified miRNA targets was performed using plantGSEA [85] (http://structuralbiology.cau.edu.cn/PlantGSEA/). Conserved lncRNAs were predicted by the following steps: (1) Identifying ortholog gene pairs that are syntenically conserved across Z. mays, O. sativa and Setaria italica by BBH-LS with default parameters [86]; (2) aligning the intergenic regions between ortholog gene pair with BLAST setting e-value cutoff to be 1e−6 [87]; (3) checking if any reads were mapped to this region by Bowtie2 version 2.2.6 [88]; (4) checking the secondary structure of expressed intergenic regions by RNAz 2.0 [89] with default parameters to identify precursor sequences of potential lncRNAs.

Experimental validation of miRNA-target gene pairs
Short tandem target mimic modules (STTM) for each miRNA was annealed from a pair of complementary long primers which contained the spacer region and two same miRNA binding sites in, adjunct HindIII and PstI sites were at the sides [90]. The miRNA binding sites include a CTA trinucleotide bulge corresponding to positions 10 and 11 from the 5′ end of the mature miRNAs, so that the STTM can effectively bind but will not be cleaved by the miRNAs. This module was then inserted between the 2X35S promoter and rbcS terminator of the binary vector pHB (12 kb) through the HindIII and PstI sites. pHB uses the hygromycin gene and bar gene as selection markers. Single colonies were picked up for plasmid isolation with the detection primers and the constructs were verified by DNA sequencing. Transgenic plants were generated by Agrobacterium tumefaciens-mediated floral dip transformation. For each miRNA, about 10 transgened O. sativa japonica lines were obtained. RT-PCR was then perfomed for each individual plant to access the expression levels of both miRNAs and their corresponding target genes.