Skip to main content

Maize miRNAs and their putative target genes involved in chilling stress response in 5-day old seedlings



In the context of early sowing of maize as a promising adaptation strategy that could significantly reduce the negative effects of climate change, an in-depth understanding of mechanisms underlying plant response to low-temperature stress is demanded. Although microRNAs (miRNAs) have been recognized as key regulators of plant stress response, research on their role in chilling tolerance of maize during early seedling stages is scarce. Therefore, it is of great significance to explore chilling-responsive miRNAs, reveal their expression patterns and associated target genes, as well as to examine the possible functions of the conserved and novel miRNAs. In this study, the role of miRNAs was examined in 5d-old maize seedlings of one tolerant and one sensitive inbred line exposed to chilling (10/8 °C) stress for 6 h and 24 h, by applying high throughput sequencing.


A total of 145 annotated known miRNAs belonging to 30 families and 876 potentially novel miRNAs were identified. Differential expression (DE) analysis between control and stress conditions identified 98 common miRNAs for both genotypes at one time point and eight miRNAs at both time points. Target prediction and enrichment analysis showed that the DE zma-miR396, zma-miR156, zma-miR319, and zma-miR159 miRNAs modulate growth and development. Furthermore, it was found that several other DE miRNAs were involved in abiotic stress response: antioxidative mechanisms (zma-miR398), signal transduction (zma-miR156, zma-miR167, zma-miR169) and regulation of water content (zma-miR164, zma-miR394, zma-miR396). The results underline the zma-miRNAs involvement in the modulation of their target genes expression as an important aspect of the plant’s survival strategy and acclimation to chilling stress conditions.


To our understanding, this is the first study on miRNAs in 5-d old seedlings’ response to chilling stress, providing data on the role of known and novel miRNAs post-transcriptional regulation of expressed genes and contributing a possible platform for further network and functional analysis.

Peer Review reports


The negative consequences of climate change on crop yield, mainly due to increases in average and extreme temperature, changes in rainfall patterns, and elevated CO2 concentrations are already noticeable and the future crop yield projections are being recognized as potentially a major societal concern [1, 2]. By using ensembles of latest-generation crop and climate models, Jägermeyr et al. suggested that climate impacts will emerge earlier than previously projected [3]. At the same time, maize was detected as the most vulnerable crop, because as a C4 crop it has a small capacity to benefit from the elevated CO2. Also, it is grown across a wide range of low latitudes which are recognized as the most endangered regions due to current proximity to crop-limiting temperature thresholds. As maize is the most important global crop in terms of total production, these projections impose great risks to food security.

However, it is expected that cropping system adaptation and/or better-adapted varieties can considerably reduce climate change impacts [4, 5]. In the latter context, early sowing of chilling tolerant early maturity maize hybrids in temperate regions with hot summers coupled with absence/low rainfall could ensure avoiding the negative effects of drought stress during the flowering period and ensure sustainable yield. In maize, low temperature in the range of 5–15 °C is sufficient to compromise growth from the seedling to the mature stage of development. For this reason, shifting the sowing date forward (sowing at sub-optimal temperatures below 15 °C) requires tolerant hybrids because low temperatures can reduce germination, emergence rate, and seedling vitality. Furthermore, low temperatures affect cell survival, cell division, photosynthetic efficiency, and water transport, subsequently leading to a reduction in plant growth and lower productivity [6].

Disruption of normal cell functions caused by chilling stress demands a brisk and comprehensive reprogramming at the molecular level as the result of transcriptional, post-transcriptional, and translational regulation of stress-responsive genes [7]. MicroRNA (miRNA) are a class of small non-coding RNA molecules approximately 18–24 nucleotides in length which regulate gene expression by cleavage or/and translational inhibition of target mRNA [8]. They control complex regulatory networks, are involved in a broad range of biological processes and are shown to be critical regulators of developmental process and response to abiotic/biotic stress [9]. As they regulate various biological processes by targeting multiple genes simultaneously, miRNAs are considered a promising strategy for complex trait improvement, including yield, biomass production, and stress tolerance [10]. The possibility of harnessing the diverse functions of miRNAs to achieve desirable agronomic traits in important crops, particularly abiotic stress response, was also emphasized by Zhang et al., by reviewing their role and the role of their targets in response to low temperature, high temperature, drought, soil salinity, and heavy metals [8].

It has been shown that miRNAs are involved in the regulation of chilling stress in different plant species. Megha et al. summarized the current understanding of miRNA-mediated modulation of the expression of key genes as well as genetic and regulatory pathways, involved in chilling stress responses in plants [11]. The authors accentuated growing evidence in the literature that miRNAs reprogramming of gene expression is a major defense mechanism in plants enabling them to respond to stresses. Among the miRNAs found to be cold responsive in various plant species are miR319, miR394, miR398 and miR156 [11]. In rice, sugarcane, and cassava, miR319 positively regulates cold tolerance by repressing the expression of two TCP genes. miR394 is highly conserved miRNA in both monocots and dicots, and its involvement in cold stress response was found in Arabidopsis. miR398 are known regulators of the antioxidative response in plants, controlling the expression of Cu/Zn superoxide dismutase 1 and 2 (CSD1 and CSD2), shown to affect the tolerance of winter wheat to cold stress. miR156, which are key regulators of squamosa promoter-binding (SPB) like proteins (SPL) expression, were identified as chilling stress-responsive in tomato. Further studies may reveal more detailed functions of miRNAs involved in chilling stress response.

However, despite research done on the role of miRNAs in the chilling response in maize during the V3 stage [12], to our knowledge there are no studies on earlier seedling stages. Considering that early sowing of maize hybrids is a promising adaptation strategy that could significantly reduce the negative effects of climate change, the development of genotypes tolerant to low-temperature stress in the early stages of plant growth and development is of major importance. To achieve this goal, an in-depth understanding of mechanisms underlying plant response to stress is required.

As miRNAs are being recognized as a major defense mechanism in plants enabling them to respond to stress [11], the aim of the experiment was examining the role of miRNAs in the chilling response of young maize seedlings, by applying high throughput sequencing. This study compares the expression of miRNA in the 5-day old seedlings of a chilling-tolerant and a chilling-sensitive maize inbred line, under both optimal and chilling conditions, in order to explore chilling-responsive miRNAs and reveal expression patterns of miRNAs and associated target genes, as well as to examine the possible functions of the conserved and novel miRNAs.


Chilling response assessment of maize inbred lines

Six elite maize inbred lines (L1-L6), parental components of commercial ZP hybrids, developed in Maize Research Institute Zemun Polje, were chosen for the assessment of the chilling stress response. Lines L1 to L3 belong to Lancaster, L4 to Iowa dent, L5 to BSSS/Iowa dent, and L6 to BSSS heterotic groups. With respect to the kernel type, L1 to L4 are dents, and L5 and L6 are semi-dents. According to the breeders' experience in the field regarding chilling stress tolerance, L1 and L6 are considered as highly tolerant genotypes, L2 and L3 as genotypes with low tolerance, while L4 and L5 as intermediately tolerant genotypes.

The experiment was done in three replicates of 10 plants per inbred line for each analysis. Seeds were germinated for five days in a growth chamber (MLR-352H-PE Climate Chamber, PHC Europe B.V., The Netherlands), after they were sterilized in 10% sodium hypochlorite (commercial bleach). Germination was performed in the dark at 25/20 °C (12/12 h) and relative humidity of 75%. After reaching the desired growth stage, the 5-day old seedlings were subjected to chilling stress for 24 h at 10/8 °C with a 12/12 h light/dark photoperiod and light intensity up to 700 µmol/m2/s. After recording the survival rate (SR), calculated as the percentage of plants surviving the stress treatment, a subset of ten seedlings per genotype was immediately sampled and used for seedling fresh weight (FW5-d), radicle (Lrad) and coleoptile (Lcol) length measurements. The rest of the seedlings were sown in pots containing a mixture of soil and sand (3:1), moved back to the growth chamber, and grown under the optimal temperatures of 25/20 °C, with the same photoperiod, relative humidity, and light intensity as described, for another seven days to recover.

Recovered plants were harvested and used for measuring root (RFW) and shoot (SFW) fresh weight, as well as root (RL) and shoot (SL) length. Dry root (RDW) and shoot (SDW) weights were determined after a 24-h drying period at 110° C in a drying oven. Control plants of all inbred lines were grown under the optimal conditions until the same developmental stage as the treated ones and sampled at the identical time points.

Chilling treatment of LT and LS lines

Based on the results of the chilling stress response assessment and their statistical significance, two maize inbred lines of contrasting tolerance to chilling were selected for further research. The selection underwent two levels of assessment – following the 24-h stress period and after seven days of recovery. The chosen inbred lines were marked as LT, shown to be the most tolerant during the assessment, and LS, chosen as the most sensitive of the six lines.

The experimental setup for the low-temperature treatment of LT and LS was the same as explained above – after five days of germination under optimal conditions in the dark, the 5-day old seedlings were exposed to the previously described chilling conditions for 24 h. Sampling for further analyses was done 6 h and 24 h after the start of the treatment. Control plants were grown under optimal conditions in the same time period and sampled at identical time points.

Total RNA isolation

Total RNA was extracted from 30 maize seedlings per inbred line (LT and LS) after each time point (6 h and 24 h). The sampled tissue was ground in liquid nitrogen using a pestle and mortar and stored at -80° C until further use. Approximately 100 mg of the frozen tissue was used for the total RNA extraction (GeneJet™ RNA Purification kit, Thermo Scientific, MA, USA). Isolated RNA was further treated with DNase I (Ambion® DNA-free™ kit, Invitrogen, MA, USA), and total RNA concentrations and quality were determined by several methods. Preliminary quantitation was performed with the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), while agarose gel electrophoresis was utilized to check for RNA degradation and potential contamination. Finally, RNA integrity and quantitation were carried out with the 2100 Bioanalyzer and RNA Nano 6000 Assay Kit (Agilent®, CA, USA). All samples with RIN above seven were chosen for the downstream analysis.

NGS library preparation and sequencing

Eight samples (four per genotype) were used for the library preparation and Illumina sequencing, which were carried out at Novogene Bioinformatics Technology Co., Ltd., in Beijing, China. The LT samples encompassed 6 h and 24 h time points, control (C) and treated (T) plants—LT-C-6, LT-C-24, LT-T-6, and LT-T-24, respectively. Accordingly, the LS samples included LS-C-6, LS-C-24, LS-T-6, and LS-T-24. Small RNA library construction was carried out using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB®, USA), adding index codes to attribute sequences to each sample and following the manufacturer’s recommendations. Briefly, small RNA libraries were generated from 3 μg of total RNA by adapter ligation (first the 3', followed by the 5' adapter), first strand cDNA synthesis, PCR amplification, and purification through 8% polyacrylamide gel electrophoresis. Library quality was assessed on the Bioanalyzer 2100 system using DNA High Sensitivity Chips (Agilent®, CA, USA). The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina®, CA, USA) according to the manufacturer’s instructions. After the cluster generation, the single-end 50 bp RNA sequencing was performed on the Novaseq 6000 platform (Illumina®, CA, USA).

sRNAseq bioinformatics data analysis

The quality of the raw reads was checked using FastQC [13], and then processed through fastp [14] to remove reads containing adapter and poly-N sequences and reads with low quality (Phred score > 30, N% < 10%). Sequencing error rates, Phred and Q-scores, as well as the GC content, were calculated. Additionally, the clean reads of high quality were filtered based on their length and only 18–24 nt reads were used for the downstream analyses. Small RNA tags were mapped onto Zea mays B73 NAM reference genome version 5.0 ( using Bowtie 0.12.9 [15] without any mismatches to analyze their expression and distribution on the reference.

Known, conserved miRNAs were identified through the alignment of the mapped sRNA tags to the registered miRNAs in the miRBase database (Release 20, [16] using the miRdeep2 software [17]. Sequences with perfect matches were regarded as conserved miRNAs. These detected known miRNAs were also compared to those of different organisms in the miRbase. Novel miRNAs were predicted from the miRNA precursors by applying software tools miREvo [18] and miRdeep2. Both are based on a modified miRDeep algorithm that uses unannotated sRNA tags to predict potentially novel miRNAs, by exploring the secondary, miRNA hairpin structure, the Dicer cleavage site, and the minimum free energy of the unannotated sRNAs. The algorithm also has stringent score cut-offs that significantly limit the possibility of false positives. Tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, were removed by mapping to RepeatMasker ( and the Rfam database [19]. Target genes of the conserved and potentially novel miRNAs have been predicted using the psRobot 1.2 tool [20]. No mismatches were allowed at miRNA positions 2–17 to increase the reliability of the prediction.

Gene Ontology (GO) enrichment analysis of the potential target genes was performed using the GOseq 2.12 Wallenius-based non-central hyper-geometric distribution [21], which could adjust for gene length bias. KOBAS software 3.0 [22] was used to test the statistical enrichment of the target gene candidates in KEGG pathways [23]. Additionally, potential miRNA secondary structures, counts, and first position bias were obtained through miRdeep2.

Differential expression analysis between the control and treatment samples was performed using the DEGseq R package 1.2.2 [24], after previously executed normalization using the TPM algorithm [25]. The p-value was adjusted according to the Bayesian interpretation [26], and the adjusted value, or the q-value was further used. The threshold for significant differential expression was set as q-value < 0.01 and log2 fold change ≥ 1 or ≤ -1.

Validation of known and novel miRNAs using qRT-PCR

The sequencing results were validated using qRT-PCR to analyze the expression of the selected known and novel miRNAs, as well as their target genes. Five miRNAs were chosen based on their quantity across all eight libraries and differential expression levels: novel_904, zma-miR11970-3p, zma-miR159c-5p, zma-miR166a-5p and zma-miR396c. Chosen corresponding targets genes for these miRNAs were heat stress transcription factor A-2e (hsfA-2e), mitochondrial ADP/ATP carrier protein 1 (aac1), Zea mays Homeobox-leucine zipper protein, ROC7 (roc7), chloroplastic acyl carrier protein 2 (acp2) and growth-regulating factors 5 and 8 (grf5/8), respectively. Additionally, the expression of several known Zea mays cold induced genes (ZmCOI) was analyzed using qRT-PCR, to further confirm the chilling response assessment of the two chosen maize inbred lines.

Total RNA was extracted and treated with DNase I, as previously described. Two μg of total RNA were used in cDNA synthesis using the Revert Aid First Strand cDNA synthesis kit with RNase inhibitor (Thermo Scientific™). Real-time PCR analysis was carried out using cyclophilin (cyp) as the internal reference gene [27] and three biological replicates for each sample. PCR reactions were performed on a StepOnePlus™ Real-Time PCR System (Applied Biosystems™). The reaction mixture (10 μl) contained 1 × HOT FIREPol® EvaGreen® qPCR Mix Plus (ROX) (Solis BioDyne™), 0.2 μM of each primer (forward and reverse) and 1 μl of template cDNA diluted two times. The thermal cycling conditions included the initial denaturation (95° C for 10 min), followed by 40 cycles of denaturation (95° C for 15 s), primer annealing, and extension (appropriate temperature for 60 s). The primers used are listed in the Additional File 1. Primers for miRNA validation were designed using miRprimer [28], while the target gene primers were designed using Primer 3 (v 0.4.0) online software ( and checked in NCBI Primer-BLAST tool ( Relative gene expression was calculated according to Livak and Schmittgen [29] using efficiency correction as in Pfaffl [30].

Statistical analysis

All statistical analyses were done using the stats package in R 4.1.2 [31]. The results obtained from the chilling response assessment of maize inbred lines were checked for normality (Shapiro–Wilk and Kolmogorov–Smirnov test) and statistically processed by applying t-test with a significance level at p < 0.05. t-test was also carried out for mean comparison with a significance level at p < 0.05 for qRT-PCR validation results. Graphs in Figs. 1, 2 and 7 were designed in Microsoft Excel (Windows 10). Cluster analysis, based on the expression patterns of 101 known miRNAs identified in the eight libraries, was performed using normalized read count values and Euclidian distance method, to form both the K-means and hierarchical clusters.

Fig. 1
figure 1

Chilling response assessment of 5-day old seedlings after 24 h stress period. L1 to L6 – six maize inbred lines with different level of chilling tolerance. Control is shown in blue, treatment in yellow. The significance of the difference between the control and treatment of each parameter was determined by t-test and is shown as *** (p < 0.001), ** (p < 0.01), * (p < 0.05), and NS (statistically not significant at p < 0.05)

Fig. 2
figure 2

Chilling response assessment after seven days recovery period. L1 to L6 – six maize inbred lines with different levels of chilling tolerance. The control is shown in blue, and the treatment in yellow. The significance of the difference between the control and treatment of each parameter was determined by t-test and is shown as *** (p < 0.001), ** (p < 0.01), * (p < 0.05), and NS (statistically not significant at p < 0.05)


Chilling response assessment of maize inbred lines

Chilling stress tolerance of the six inbred lines underwent two levels of assessment – following the 24-h stress period (Fig. 1) and after seven days of recovery (Fig. 2). It can be seen from Fig. 1 that lines L1 and L6 displayed the highest levels of tolerance regarding Lrad, Lcol and FW, as all negative effects of the stress were statistically insignificant. However, although insignificant, the percentage of decrease was lower for L6 compared to L1 – 10% vs. 59% for Lrad, 13.5% vs. 18.4% for Lcol, and 3.7% vs. 9% for FW. Besides, line L6 had higher seed vigor than L1, and it is well known that genotypes with high vigor can cope more efficiently with various stress factors since they evolved complex resistance mechanisms [32]. Moreover, SR of L1 and L6 showed a great difference – L1 showed a 20% decrease, while no change was detected for L6 under the stress. On the other hand, lines L3 and L4 showed the highest level of susceptibility, as the values of all analyzed traits were significantly decreased (p < 0.01 or p < 0.001). A high decrease of SR (20%) was also detected. The other two lines (L2 and L5) displayed intermediate tolerance with statistically significant changes only for Lrad, but also with 27% and 10% decrease in SR, respectively.

Considering parameters measured after seven days of recovery, L1 displayed better recovery compared to L6 (Fig. 2). The difference between control and treated plants was insignificant for all root and shoot traits measured in L1, except for RL for which a statistically significant increase (p < 0.01) was detected. However, a statistically significant increase (p < 0.01) in DRW was detected in L6. Still, L6 was chosen as the tolerant inbred line for further studies due to its better performance at the 5-day stage, considering that shifting sowing data backwards requires genotypes tolerant to low temperatures during germination. The L3 line was chosen as the most susceptible line, as it had the poorest performance in both 5-day and recovery stages. From this point onward, L3 will be labeled as LS and L6 as LT.

Identification of known and potentially novel miRNAs

On average, 21.3 million clean reads were generated from the sequencing the eight sRNA libraries and 13.9 million were 18–30 nt in length. Size distribution of the reads (Fig. 3) was similar across all eight libraries and the majority of reads were 20–24 nt long (on average 67.9%), with the most abundant being 24 nt long reads (on average 29.97% of all 18–30 nt reads). Most of the clean reads (81.39%) were successfully mapped sRNA reads, and 56.9% were unique sRNA reads. For annotation purposes, the sRNAs were grouped into several classes, shown for the total and unique reads (Additional File 2).

Fig. 3
figure 3

Size distribution of detected sRNAs. sRNAs are sorted into groups of 18–19 nt, 20–21 nt, 22–23 nt, 24–25 nt, 26–27 nt, 28–29 nt and 30–31 nt

Considering the miRNAs, a total of 145 annotated known miRNAs, belonging to 30 families, were identified (Fig. 4, Additional File 3). miRNA families represented by most members were miR166 and miR171 with 12 members, closely followed by miR156, miR164, and miR169 with 11 members. Only 68 miRNAs were present in all eight libraries, most of them belonging to miR156, miR166 (seven members), and miR159 (six members). Even fewer miRNAs were specific to the low-temperature treatment, and none of them were common for both genotypes. For example, zma-miR164e-5p, zma-miR167j-3p, zma-miR169c-5p, zma-miR171c-3p, and zma-miR393b-3p were expressed only during the treatment in LT, but not Ls. Certain miRNAs were expressed only under the optimal conditions, and they were also specific to the genotype.

Fig. 4
figure 4

Distribution of miRNAs among genotypes and treatments. (A) Known miRNAs in LT under low-temperature treatments and control; (B) Known miRNAs in LS under low-temperature treatments and control; (C) Novel miRNAs in LT under low-temperature treatments and control; (D) Novel miRNAs in LS under low-temperature treatments and control. The total number of miRNAs in each library is shown in italic outside of the diagram

Additionally, detected known miRNAs were compared to those of different organisms. No matches were found with viruses, bacteria, and animals, but members of 25 miRNA families were expressed in the sequenced libraries of other plant species. miRNA families with the most matching miRNAs were miR171, miR169, and miR395, with over ten matching miRNAs across 35 different plant species. Most common miRNAs were found in Oryza sativa (103 matches), Glycine max (97 matches), Populus trichocarpa (92 matches), Malus domestica (91 matches), and Sorghum bicolor (90 matches). Additionally, matching miRNAs found across most species were MIR156a (40 plant species), MIR166a (37 plant species), MIR156b (36 plant species). Interestingly, matches were also found among other vascular plants, such as the lycophyte Selaginella moellendorffii (12 miRNAs across 6 miRNA families) and bryophyte Physcomitrella patens (43 miRNAs across 8 miRNA families).

Unknown sRNAs that could be mapped to the reference sequence were identified as potentially novel miRNAs, and the ones with high confidence included 876 sequences (Fig. 4, Additional File 4). Of the 876 novel miRNAs, 124 were specific to the tolerant genotype and 137 were found only in the sensitive one. In LT, 324 novel miRNA sequences were expressed in all four libraries, under the control and treatment conditions, while in LS the same can be said for 386 miRNA sequences. Only 150 novel miRNA sequences were expressed in all eight sequenced libraries. Out of the 189 novel sequences detected only in the treatment libraries (LT-T-6, LT-T-24, LS-T-6, and LS-T-24), 27 were common for both genotypes, while 71 were found only in LT and 91 only in LS. Analysis of first position nucleotide bias in these miRNA sequences showed that the first nucleotide tended to be adenine (A), while other position biases were not as pronounced. A significant part of these novel miRNAs showed diversified expression levels, and while some showed great abundance (over 1000 reads in at least one of the libraries), most were represented with less than 100 reads. This expression level is lower than that of conserved miRNAs.

Target prediction of the known and novel miRNAs

A total of 1822 genes were identified as targets for both known and novel miRNAs. Out of 145 known miRNAs, 138 were shown to be able to target 678 genes (Additional File 5), while among the novel miRNAs 1241 potential targets were identified for 440 novel miRNAs (Additional File 6). Nearly 70 target genes were regulated by multiple known and novel miRNAs according to prediction.

GO enrichment analysis showed that most target genes were included in the molecular function of binding (over 57%), particularly protein binding (≈18%), as well as regulation of other biological processes, such as localization and transport (Additional File 7). KEGG analysis on the other hand revealed that most of them were active in different metabolic pathways (≈30%) and biosynthesis of secondary metabolites (≈17%) (Additional File 8). Further gene annotation analysis showed that a significant number of these genes were involved in the plants’ stress response. Many target genes encoded different transcription factors important for the stress response, including heat stress transcription factors (HSF), growth regulating factors (GRF), ethylene responsive factors (ERF), as well as many others such as WRKY, HOX, MYB, bZIP, and auxin-responsive (ARF) transcription factors. Additionally, other genes included in the stress response and signal transduction were identified as targets for these miRNAs, such as genes encoding cysteine-rich protein kinases (RLKs), serine/threonine protein kinases (SAPKs), calcium dependent ion transporters and universal stress proteins (USPs).

miRNA expression patterns under chilling stress in maize 5-d old seedlings

Differential expression analysis was performed on 933 miRNAs (119 known and 814 novel miRNAs) that had expression levels higher than 5 TPM. Out of the 933 miRNAs, 649 were expressed differentially (DE) between the control and treatment conditions – 89 known and 560 novel miRNAs (Additional File 9).

In LT, there were 417 miRNAs differentially expressed between the control and treatment conditions: 159 after 6 h treatment and 331 after 24 h (Fig. 5A). There were 73 miRNAs common for both treatment durations with various expression patterns (Fig. 5B). On the other hand, among the 436 DE miRNAs in LS, 355 were differentially expressed after 6 h treatment and only 146 after 24 h (Fig. 5C). Additionally, the ratio of DE miRNAs with the same expression patterns after both treatment durations was higher in LS compared to LT, 42 in LS (Fig. 5D) and 33 in LT (Fig. 5B).

Fig. 5
figure 5

Distribution of differentially expressed miRNAs across the treatments and their expression patterns. (A) Differentially expressed known and novel miRNAs after 6 h and 24 h of treatment in LT. (B) Expression patterns of the 73 miRNAs that are common for both treatment durations in LT. (C) Differentially expressed known and novel miRNAs after 6 h and 24 h of treatment in LS. (D) Expression patterns of the 65 miRNAs that are common for both treatment durations in LS

A total of 55 DE miRNAs were detected in both genotypes after 6 h of treatment. Twenty-eight followed the same expression pattern – 18 up-regulated and 10 down-regulated. On the other hand, 12 were up-regulated in LT and down-regulated in LS, while 15 were up-regulated in LS but not in LT. Considering 24 h treatment, 51 common DE miRNAs were detected. Here, most miRNAs showed the same expression patterns in both LT and LS—12 up-regulated and 19 down-regulated. Only eight miRNAs were expressed in both genotypes and at both time points showing different expression profiles (Fig. 6A).

Fig. 6
figure 6

Differential expression of common miRNAs and cluster analysis. (A) Differential expression of the eight common miRNAs (novel_102, novel_445, novel_625, novel_800, zma-miR159c-5p, zma-miR166a-5p, zma-miR398a-3p, zma-miR398b-5p), between the control and treated samples in LT and LS after 6 h and 24 h shown as the fold change. (B) Cluster analysis, based on the expression patterns of 101 known miRNAs in the eight libraries

Cluster analysis, based on the expression patterns of 101 known miRNAs, revealed that those expression patterns were genotype-specific, more than treatment-specific. Three clusters were observed – all the samples belonging to LS were classified into a single cluster, while LT-C-6 was separated from the rest of LT into a separate cluster. Additionally, the known miRNAs were divided into four clusters based on their expression patterns across all eight libraries. (Fig. 6B).

Validation of known miRNAs, novel miRNAs and target genes using qRT-PCR

Sequencing reliability was confirmed through qRT-PCR analysis of the miRNA expression levels. Five miRNAs (novel_904, zma-miR11970-3p, zma-miR159c-5p, zma-miR166a-5p and zma-miR396c) were selected for the analysis based on their quantity across all eight libraries and differential expression levels. In general, the expression patterns obtained through RT-PCR matched those obtained through NGS sequencing, with small inconsistencies. Additionally, the expression of certain target genes (hsfA-2e, aac1, roc7, acp2, grf5, grf8) was evaluated in the same way and opposite expression profiles between the target genes and their corresponding miRNAs was observed. The results of the validation are shown in Fig. 7.

Fig. 7
figure 7

Validation of known and novel miRNAs, as well as their target genes using qRT-PCR. qRT-PCR validation was performed on selected miRNA (novel_904, zma-miR11970-3p, zma-miR159c-5p, zma-miR166a-5p and zma-miR396c), and their targets genes (hsfA-2e, aac1, roc7, acp2, grf5, grf8). The expression patterns obtained from the next-generation sequencing (sRNAseq) are shown as the log2 fold changes between the control and treated samples in LT and LS after 6 h and 24 h. qRT-PCR expression patterns are shown as 2.–ΔΔCt values [30, 31] obtained from the ΔΔCt values from control and treated samples in LT and LS after 6 h and 24 h. t-test was performed for mean comparison (p < 0.05)

Additionally, the expression levels of several known ZmCOI genes were analyzed using qRT-PCR, at the same time points (Additional file 10). These included: fatty acid desaturases 2 and 7 (FAD2, FAD7), dehydration-responsive element binding protein 2A (ZmDREB2A), calcium-dependent protein kinase 1 (ZmCPK1) and zinc finger protein AN13 (ZmAN13). The expression of the ZmCOI genes was analyzed at the same time points. Statistically significant differences were found between the control and treatment, as well as between the genotypes, at either of the time points for all five genes.


Chilling stress tolerance in maize seedlings is important for implementing early sowing as a strategy for avoiding the negative consequences of climate change on maize yield and biomass. In recent years, many studies have shown that miRNAs play important roles in the response to abiotic stress in plants. Still, miRNA studies in maize seedling tolerance to chilling stress are scarce, and understanding their mode of action would be a valuable reference for future molecular studies and improvement of breeding strategies. This is the first study designed to investigate the role of 5-d old maize seedlings׳ chilling-responsive miRNAs, their expression patterns, associated target genes and possible functions.

miRNA in maize compared to other species

Identified 5d-old maize miRNAs (zma-miRNA) were compared to those of other species. Size distribution analysis showed that most abundant miRNAs were 24 nt long. This is consistent with the results from other research focused on zma-miRNA identification under abiotic [33] and biotic [34] stress conditions, as well as during germination [35].

By comparing the identified miRNAs miRNAs from other species, it was confirmed that the most conserved ones belonged to miR156, mir166, miR169, miR395 and mir396 families as they were found in a broad range of plant species, from bryophytes to eudicots [36]. This concurrence shows that many miRNAs belong to evolutionary conserved regulatory modules that play important roles in plant development [37]. Expectedly, maize shares most of the common miRNAs with other monocots such as rice and sorghum. But, common miRNAs were also found with soybean and tree species (such as Malus domestica and Populus trichocarpa). This suggests a possibility that many miRNAs had independent origins, unrelated to their phylogenetic tree during evolution, meaning that they could be homoplasious rather than homologous.

Roles of zma-mirRNA in the response to chilling stress

Out of the 933 analyzed miRNAs, 98 were differentially expressed in both genotypes at one of the time points – 55 after six hours and 51 after one day of exposure to the treatment conditions. Only eight were differentially expressed in both genotypes at both time points (novel_102, novel_445, novel_625, novel_800, zma-miR159c-5p, zma-miR166a-5p, zma-miR398a-3p, zma-miR398b-5p). Over 300 genes were found to be potential targets of these 98 miRNAs.

GO enrichment analysis revealed that target genes were involved in various biological processes, mainly in different aspects of growth and development (zma-miR396, zma-miR156, zma-miR319, and zma-miR159), and in the direct response to stress factors. Of these, some were involved in the antioxidative response (zma-miR398), signal transduction (zma-miR156, zma-miR167, zma-miR169) and regulation of water content (zma-miR164, zma-miR394, zma-miR396). The rest of the targeted genes had roles in cell wall formation, photosynthesis, protein metabolism, and nutrient assimilation.

zma-miRNAs might regulate growth and organ development during chilling stress

Various miRNAs target genes that encode different components of regulatory networks involved in plant growth, root and shoot development, tissue morphogenesis, and other developmental processes. miRNA families known to be involved in this process are miR396, miR156, mir166, miR319, and miR159.

The role of miR396 is in controlling the growth of multiple tissues and organs in a variety of species by modulating growth regulating (GRFs) and GRF-interacting factors (GIFs) [38]. The miR396–GRF/GIF module fulfills multiple roles in plant development. They include the promotion of shoot and shoot lateral organ growth, particularly leaf development and leaf cell proliferation, as well as overall cell proliferation and expansion. Considering miR156 family members, they are key regulators of squamosa promoter-binding (SPB) like proteins (SPL) expression, controlling the timing of vegetative phase change, leaf initiation rate, shoot branching, and lateral root development [39, 40]. The role of miR319 is important for maintaining cell division functions in the leaf meristem, through expression regulation of several transcription factors (TF) such as gibberellin-and-abscisic-acid-regulated MYB (GAMYB), and teosinte branched1/cycloidea/PCF (TCP) [41]. Besides miR319, miR159 also participates in GAMYB regulation by suppressing its expression during development, ensuring normal growth during the vegetative stage [42]. It was also shown that miR159 was accumulated in response to different abiotic stresses in various crop species and that this increase seems to be regulated by abscisic acid (ABA). Our results of the target prediction analysis coincide with the miRNA-target gene connections found in the literature, giving further support to these regulatory networks. Furthermore, target prediction showed a potential role for zma-miR528a and zma-miR408a in root development regulation, with miR528a possibly controlling bHLH TFs necessary for root hair development and miR408a being involved in Casparian strip lignification and root cell differentiation.

Differential expression analysis showed that the expression of the identified miRNAs was affected by low temperatures, with many being up-regulated in both genotypes after 6 h (zma-miR156d-5p, zma-miR156k-5p) or 24 h (zma-miR319a-3p, zma-miR396c). Other plant species also showed matching expression profiles under abiotic stresses. It was shown that miR396 are up-regulated by various stress conditions, such as drought and UV-B irradiation in Arabidopsis [43]. Similarly, miR156 expression was increased in many plant species during different abiotic stress treatments, including heat [44], high salinity [45], and drought [46].

The described expression profiles indicate the potential role of the analyzed miRNAs in inhibiting seedling growth and development under chilling conditions, suggesting the up-regulation of the specific miRNAs and the consequential decrease in the expression of their target genes is an important aspect of the plant’s survival strategy.

zma-miRNAs involved in antioxidative response regulation

The miR398 are known regulators of the antioxidative response in plants, controlling the expression of Cu/Zn superoxide dismutase 1 and 2 (CSD1 and CSD2) which are responsible for reactive oxygen species (ROS) scavenging [47]. Expression of miR398 is regulated by ABA, as well as by WRKY TFs which are known to be involved in abiotic stress responses, particularly in chilling [48]. Two miR398 family members (zma-miR398a and zma-miR398b) showed a significant difference in their expression levels, at both time points and in both genotypes. In LS, the expression of zma-miR398b was increased at both time points and of zma-miR398a after 24 h. Contrary to LS, the initial up-regulation of both miRNAs after 6 h was followed by down-regulation after 24 h in LT. Down-regulation of miR398 is related to abiotic stress tolerance and is followed by the increase in CSD1 and CSD2 expression in different plant species under many stress treatments [49, 50]. This is concurrent with the DE analysis performed for LT and LS. To our knowledge, this is the first report of miR398 being differentially expressed during chilling stress in maize.

Additionally, according to target prediction, several potentially novel miRNAs could also be involved in antioxidative response. For example, novel_16, which is possibly involved in regulating 2-alkenal reductase necessary for alleviating oxidative stress in maize [51], was up-regulated in both genotypes after 24 h. On the other hand, novel_447 which could be involved in regulating the synthesis of ascorbic acid by targeting L-gulonolactone oxidase 2 and in hydrogen peroxide removal by targeting L-ascorbate peroxidase 2, was down-regulated after 6 h in both genotypes. Both enzymes are involved in the abiotic stress response and alleviating the negative effects of ROS [52], indicating that the down-regulation of this miRNA could be related to increased resistance to oxidative damage during chilling.

zma-miRNAs included in stress signal transduction

Target prediction analysis showed that several DE zma-miRNAs (zma-miR156d-5p, novel_447, zma-miR528a-5p, and novel_342) have a role in Ca2+ signal transduction, known to play a critical role in chilling stress response in plants. zma-miR156d-5p and novel_447 are involved in the regulation of calmodulin-binding transcription activator 3 (CAMTA3), which perceives an increase in calcium level as a response to chilling and freezing stress [53]. CAMTA3, together with CAMTA1, CAMTA2, CAMTA4 and CAMTA5, is required for the chilling-induced expression of TFs of chilling regulated (COR) genes DREB1B/CBF1, DREB1C/CBF2 and GOLS3 [54]. Büyük et al. predicted that among others, miR156 potentially regulates CAMTA factors in various plant species [55]. Several studies noted the potential role of miR156 in regulating several CAMTA during periods of abiotic stress, but not in maize under low temperature conditions. Additionally, Kansal et al. showed that calcium cytosolic levels, through CAMTA4 and CAMTA6, play a role in regulating the expression of os-miR156a and os-miR167h under drought conditions in rice and this could indicate a possibility that a similar relation also occurs in maize [56].

Zhao et al. suggested that the role of miR156 during abiotic stress could be in balancing the vegetative development and chilling or freezing response [57]. Additionally, it was shown that miR156 and their SPL-targets contribute to the activation of ABA signaling pathways, pointing also to their involvement in the chilling stress response [58]. Expression of miR156 after exposure to low temperatures varies among species. Accordingly, down-regulation of miR156 is related to an increase in tolerance in rice [59], but overexpression of miRna156 improved chilling tolerance in tobacco [60]. Herein, zma-miR156d-5p was up-regulated in maize 5d old seedlings of both genotypes after 6 h.

Another two identified miRNAs involved in stress signal transduction, zma-miR528a-5p and novel_342, play a role in regulating calcium-dependent protein kinase 27 (CDPK27). It is well known that CDPKs are major regulators of stress responses in plants, directly regulating target proteins. A study in tomato plants showed that CDPK27, among other CDPKs, can induce crosstalk between ROS and other signaling molecules, leading to the activation of ABA. Silencing of CDPK27 significantly decreased the ROS accumulation, as well as mitogen-activated protein kinases (MPK1/2) further downstream. This suggests that CDPK27 may play a crucial role in adaptation to chilling stress [61]. miR528 have been shown to have a role in various abiotic stresses, such as heat in wheat [62], excess nitrate in maize [63], and drought in rice [61]. Aydinoglu showed that miR528 are up-regulated under chilling conditions in the maize leaf meristem in the later developmental stages [64]. On the contrary, our results show that in the 5-day old seedling stage, their down-regulation is linked to cold tolerance—differential expression showed different expression patterns of zma-miR528a-5p and novel_342 after 24 h in the analyzed genotypes – significant down-regulation in LT, but up-regulation in LS.

Regulation of water loss by zma-mirRNA during chilling stress treatment

Chilling stress is known to invoke similar responses in plants as droughtsince the capacity of water absorption and uptake decreases with temperature decline. Both low temperatures and drought cause an imbalance in the water status and in the short term lead to an upsurge of ABA that promotes stomatal closure and in the long term lead to a decrease in leaf area and stomatal density [65].

Lower transpiration and higher water uptake efficiency (WUE) were associated with reduced stomatal density in the leaf abaxial epidermis, as well as with higher expression of negative stomatal development regulator STOMATAL DENSITY AND DISTRIBUTION 1 (SDD1) in Arabidopsis [66]. Target prediction revealed that several members of zma-miR164 family, as well as zma-miR394a-5p, are involved in regulating SDD1 expression in maize 5-d old seedlings. While zma-miR394a-5p was up-regulated in both genotypes, zma-miR164f-3p was significantly down-regulated after 24 h in LT. The role of miR164 in drought stress and water uptake regulation is well known, but mostly regarding their involvement in stress response through NAC TF modulation [67]. However, their possible role in regulating SDD1 has not been previously described. Herein, it was shown that zma-miR164f could have a role in developing chilling and other abiotic stress tolerance in maize.

The nuclear factor Y (NF-Y) complex is induced by abiotic stressors [68]. Yang et al. showed that NF-Y affects drought tolerance through modulation of ABA signaling pathway and further downstream the modulation of stomata movement, osmolyte accumulation, and ROS metabolism [69]. miR169 respond to and regulate NF-Y in many plant species under different abiotic stress factors [68]. This has also been shown in maize [70], but not under chilling stress. Herein, zma-miR169i-5p was down-regulated in both genotypes after 6 h. Evidence on the role of miR169 in establishing abiotic stress tolerance is conflicting. In some cases, the overexpression of miR169 resulted in increased tolerance [71], but in others, it was vice versa [68]. This suggests that the roles of miR169 in regulating abiotic stress response are more complex and require further research.

Additionally, according to target prediction, zma-miR396c was involved in regulating the protein import into the nucleus in response to chilling by modulating the expression of KPNB1, a homolog of human importin β. KPNB1 acts as a negative effector of drought tolerance in Arabidopsis and its inactivation leads to increased stomatal closure in response to ABA and water loss reduction [72]. While KPNB1 was up-regulated in both genotypes after 24 h, its expression level was higher in LT. The possible role of miR396c in regulating nucleus import during low temperature stress has not been previously described, but their role in drought and other abiotic stress responses is known [73]. As with miR164, there are conflicting results when it comes to the effect of miR396c on abiotic stress tolerance – in some cases, its overexpression led to the increase in tolerance [74], but in the others, tolerance was acquired by down-regulating miR396 [75]. This suggests that the expression profiles and regulatory machinery behind them could vary on the plant species. Based on their role in regulating leaf development and leaf cell proliferation [76], the role of miR396 in regulating water loss could be in decreasing leaf surface and decreasing transpiration.

Based on the information discussed above, the potential regulatory network of the miRNAs (miR156, miR159, miR164, miR166, miR169, miR319, miR394, miR396, miR398, and miR528) and their target genes involved in the maize 5d-old seedlings response to chilling is presented in Fig. 8.

Fig. 8
figure 8

Potential regulatory network of the DE miRNAs and their target genes involved in the maize 5d-old seedlings’ chilling response. Potential regulatory network includes the miRNA (miR156, miR159, miR164, miR166, miR169, miR319, miR394, miR396, miR398, and miR528) and their target genes (CDPK, CAMTA, SPL, ACP2, MYB, GRF, KPNB1, NF1, SDD1, NAC, and CSD) involved in signal transduction, growth and development, waterloss regulation and antioxidative response. Downregulation of a miRNA and/or target gene is labeled green, while the upregulation is labeled orange


Early sowing of maize is projected to reduce the negative effects of climate change and the development of genotypes tolerant to low-temperature stress in early stages of plant growth and development can be considered to be of vital importance. miRNAs are regarded as a promising strategy for complex trait improvement, including yield, biomass production, and stress tolerance. This study presents a first report on the study of miRNAs and their targets involved in chilling stress tolerance of 5d-old maize seedlings. In order to identify miRNAs and understand their potential functions, eight small RNA libraries of a tolerant and sensitive inbred line were constructed and sequenced. Thirty-three known and 65 novel zma-miRNAs, differentially expressed in response to chilling stress, were identified. Target prediction, GO, and KEGG-based enrichment analysis showed that many of the identified miRNAs affect genes involved in regulating plant growth and organ development, stress signal transduction, antioxidative response, and water loss regulation under stress. This research provides a significant basis for further network analysis, with other coding and non-coding RNAs, as well as functional analyses. Furthermore, the findings provide important information for the identification of chilling-responsive genes which could be beneficial for molecular breeding of chilling-tolerant maize.

Availability of data and materials

Sequence data generated during the current study are available in the The European Nucleotide Archive (ENA) ( The data are submitted under the project name “Maize sRNAs involved in cold stress response in 5-day old seedlings” (Accession number: PRJEB68228), with the sequences of each individual library submitted under accession numbers ERR12245651—ERR12245658.


  1. Zhao C, Liu B, Piao S, Wang X, Lobell DB, Huang Y, et al. Temperature increase reduces global yields of major crops in four independent estimates. Proc Natl Acad Sci U S A. 2017;114:9326–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Pörtner H, Roberts DC, Adams H, Adler C, Aldunce P. Climate Change 2022: Impacts, Adaptation and Vulnerability: Summary for Policymakers. Cambridge University: Cambridge, UK. 2022. Accessed 10 Dec 2023.

  3. Jägermeyr J, Müller C, Ruane AC, Elliott J, Balkovic J, Castillo O, et al. Climate impacts on global agriculture emerge earlier in new generation of climate and crop models. Nat Food. 2021;2:873–85.

    Article  PubMed  Google Scholar 

  4. Asseng S, Martre P, Maiorano A, Rötter RP, O’Leary GJ, Fitzgerald GJ, et al. Climate change impact and adaptation for wheat protein. Glob Change Biol. 2019;25:155–73.

    Article  Google Scholar 

  5. Rising J, Devineni N. Crop switching reduces agricultural losses from climate change in the United States by half under RCP 8.5. Nat Commun. 2020;11:4991.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Beck EH, Fettig S, Knake C, Hartig K, Bhattarai T. Specific and unspecific responses of plants to cold and drought stress. J Biosci. 2007;32:501–10.

    Article  CAS  PubMed  Google Scholar 

  7. Jeknić Z, Pillman KA, Dhillon T, Skinner JS, Veisz O, Cuesta-Marcos A, et al. Hv-CBF2A overexpression in barley accelerates COR gene transcript accumulation and acquisition of freezing tolerance during cold acclimation. Plant Mol Biol. 2014;84:67–82.

    Article  PubMed  Google Scholar 

  8. Zhang F, Yang J, Zhang N, Wu J, Si H. Roles of microRNAs in abiotic stress response and characteristics regulation of plant. Front Plant Sci. 2022;13: 919243.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Millar AA. The Function of miRNAs in Plants. Plants. 2020;9:198.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Sun G. MicroRNAs and their diverse functions in plants. Plant Mol Biol. 2012;80:17–36.

    Article  CAS  PubMed  Google Scholar 

  11. Megha S, Basu U, Kav NNV. Regulation of low temperature stress in plants by microRNAs. Plant, Cell Environ. 2018;41:1–15.

    Article  CAS  PubMed  Google Scholar 

  12. Li SP, Dong HX, Yang G, Wu Y, Su SZ, Shan XH, et al. Identification of microRNAs involved in chilling response of maize by high-throughput sequencing. Biologia plant. 2016;60:251–60.

    Article  CAS  Google Scholar 

  13. Andrews S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010.

  14. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47:D155–62.

    Article  CAS  PubMed  Google Scholar 

  17. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.

    Article  PubMed  Google Scholar 

  18. Wen M, Shen Y, Shi S, Tang T. miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics. 2012;13:140.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kalvari I, Nawrocki EP, Ontiveros-Palacios N, Argasinska J, Lamkiewicz K, Marz M, et al. Rfam 14: expanded coverage of metagenomic, viral and microRNA families. Nucleic Acids Res. 2021;49:D192–200.

    Article  CAS  PubMed  Google Scholar 

  20. Wu H-J, Ma Y-K, Chen T, Wang M, Wang X-J. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40:W22–8.

    Article  CAS  PubMed Central  Google Scholar 

  21. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

    Article  CAS  PubMed  Google Scholar 

  23. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2007;36:D480–4.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26:136–8.

    Article  PubMed  Google Scholar 

  25. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-Throughput Sequencing of Arabidopsis microRNAs: Evidence for Frequent Birth and Death of MIRNA Genes. PLoS ONE. 2007;2: e219.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Statist. 2003;31:2013–35.

    Article  Google Scholar 

  27. Lin Y, Zhang C, Lan H, Gao S, Liu H, Liu J, et al. Validation of Potential Reference Genes for qPCR in Maize across Abiotic Stresses, Hormone Treatments, and Tissue Types. PLoS ONE. 2014;9: e95445.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Busk PK. A tool for design of primers for microRNA-specific quantitative RT-qPCR. BMC Bioinformatics. 2014;15:29.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Livak KJ, Schmittgen TD. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  30. Pfaffl MW. Quantification strategies in real-time PCR. In: AZ of quantitative PCR. ed. Bustin, S.A. La Jolla, CA, USA: International University Line (IUL). p. 87–112.

  31. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. 2022.

  32. Li H, Yue H, Xie J, Bu J, Li L, Xin X, et al. Transcriptomic profiling of the high-vigour maize (Zea mays L.) hybrid variety response to cold and drought stresses during seed germination. Sci Rep. 2021;11:19345.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Seeve CM, Sunkar R, Zheng Y, Liu L, Liu Z, McMullen M, et al. Water-deficit responsive microRNAs in the primary root growth zone of maize. BMC Plant Biol. 2019;19:447.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Wang S, Wang X, Chen J. Identification of miRNAs Involved in Maize-Induced Systemic Resistance Primed by Trichoderma harzianum T28 against Cochliobolus heterostrophus. JoF. 2023;9:278.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Liu J, Guo X, Zhai T, Shu A, Zhao L, Liu Z, et al. Genome-wide identification and characterization of microRNAs responding to ABA and GA in maize embryos during seed germination. Plant Biol J. 2020;22:949–57.

    Article  CAS  Google Scholar 

  36. Gramzow L, Theißen G. Plant miRNA Conservation and Evolution. In: De Folter S, editor. Plant MicroRNAs. Springer, New York: New York, NY; 2019. p. 41–50.

    Chapter  Google Scholar 

  37. D’Ario M, Griffiths-Jones S, Kim M. Small RNAs: Big Impact on Plant Development. Trends Plant Sci. 2017;22:1056–68.

    Article  PubMed  Google Scholar 

  38. Liebsch D, Palatnik JF. MicroRNA miR396, GRF transcription factors and GIF co-regulators: a conserved plant growth regulatory module with potential for breeding and biotechnology. Curr Opin Plant Biol. 2020;53:31–42.

    Article  CAS  PubMed  Google Scholar 

  39. Wu G, Park MY, Conway SR, Wang J-W, Weigel D, Poethig RS. The Sequential Action of miR156 and miR172 Regulates Developmental Timing in Arabidopsis. Cell. 2009;138:750–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Gao R, Wang Y, Gruber MY, Hannoufa A. miR156/SPL10 Modulates Lateral Root Development, Branching and Leaf Morphology in Arabidopsis by Silencing AGAMOUS-LIKE 79. Front Plant Sci. 2018;8:2226.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Aydinoglu F, Lucas SJ. Identification and expression profiles of putative leaf growth related microRNAs in maize (Zea mays L.) hybrid ADA313. Gene. 2019;690:57–67.

    Article  CAS  PubMed  Google Scholar 

  42. Millar AA, Lohe A, Wong G. Biology and Function of miR159 in Plants. Plants. 2019;8:255.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Rodriguez RE, Ercoli MF, Debernardi JM, Palatnik JF. Growth-Regulating Factors, A Transcription Factor Family Regulating More than Just Plant Growth. In: Plant Transcription Factors. Elsevier; 2016. p. 269–80.

  44. Zhao J, He Q, Chen G, Wang L, Jin B. Regulation of Non-coding RNAs in Heat Stress Responses of Plants. Front Plant Sci. 2016;7:1213.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Ma Y, Xue H, Zhang F, Jiang Q, Yang S, Yue P, et al. The miR156/SPL module regulates apple salt stress tolerance by activating MdWRKY100 expression. Plant Biotechnol J. 2021;19:311–23.

    Article  CAS  PubMed  Google Scholar 

  46. Jerome Jeyakumar JM, Ali A, Wang W-M, Thiruvengadam M. Characterizing the Role of the miR156-SPL Network in Plant Development and Stress Response. Plants. 2020;9:1206.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Feng K, Yu J, Cheng Y, Ruan M, Wang R, Ye Q, et al. The SOD Gene Family in Tomato: Identification, Phylogenetic Relationships, and Expression Patterns. Front Plant Sci. 2016;7:1279.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Jiang J, Ma S, Ye N, Jiang M, Cao J, Zhang J. WRKY transcription factors in plant responses to stresses. JIPB. 2017;59:86–101.

    Article  CAS  PubMed  Google Scholar 

  49. Park S-Y, Grabau E. Bypassing miRNA-mediated gene regulation under drought stress: alternative splicing affects CSD1 gene expression. Plant Mol Biol. 2017;95:243–52.

    Article  CAS  PubMed  Google Scholar 

  50. Lu Q, Guo F, Xu Q, Cang J. LncRNA improves cold resistance of winter wheat by interacting with miR398. Functional Plant Biol. 2020;47:544.

    Article  CAS  Google Scholar 

  51. Wang Y, Zhao Y, Wang S, Liu J, Wang X, Han Y, et al. Up-regulated 2-alkenal reductase expression improves low-nitrogen tolerance in maize by alleviating oxidative stress. Plant, Cell Environ. 2021;44:559–73.

    Article  CAS  PubMed  Google Scholar 

  52. Thabet SG, Alomari DZ, Börner A, Brinch-Pedersen H, Alqudah AM. Elucidating the genetic architecture controlling antioxidant status and ionic balance in barley under salt stress. Plant Mol Biol. 2022;110:287–300.

    Article  CAS  PubMed  Google Scholar 

  53. Aslam M, Fakher B, Ashraf MA, Cheng Y, Wang B, Qin Y. Plant Low-Temperature Stress: Signaling and Response. Agronomy. 2022;12:702.

    Article  CAS  Google Scholar 

  54. Kidokoro S, Yoneda K, Takasaki H, Takahashi F, Shinozaki K, Yamaguchi-Shinozaki K. Different Cold-Signaling Pathways Function in the Responses to Rapid and Gradual Decreases in Temperature. Plant Cell. 2017;29:760–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Büyük İ, İlhan E, Şener D, Özsoy AU, Aras S. Genome-wide identification of CAMTA gene family members in Phaseolus vulgaris L. and their expression profiling during salt stress. Mol Biol Rep. 2019;46:2721–32.

    Article  PubMed  Google Scholar 

  56. Kansal S, Panwar V, Mutum RD, Raghuvanshi S. Investigations on Regulation of MicroRNAs in Rice Reveal [Ca2+]cyt Signal Transduction Regulated MicroRNAs. Front Plant Sci. 2021;12: 720009.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Cui N, Sun X, Sun M, Jia B, Duanmu H, Lv D, et al. Overexpression of OsmiR156k leads to reduced tolerance to cold stress in rice (Oryza Sativa). Mol Breeding. 2015;35:214.

    Article  Google Scholar 

  58. Yang M, Zhao Y, Shi S, Du X, Gu J, Xiao K. Wheat nuclear factor Y (NF-Y) B subfamily gene TaNF-YB3;l confers critical drought tolerance through modulation of the ABA-associated signaling pathway. Plant Cell Tiss Organ Cult. 2017;128:97–111.

    Article  CAS  Google Scholar 

  59. Zhao J, Shi M, Yu J, Guo C. SPL9 mediates freezing tolerance by directly regulating the expression of CBF2 in Arabidopsis thaliana. BMC Plant Biol. 2022;22:59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Dong H, Yan S, Jing Y, Yang R, Zhang Y, Zhou Y, et al. MIR156-Targeted SPL9 Is Phosphorylated by SnRK2s and Interacts With ABI5 to Enhance ABA Responses in Arabidopsis. Front Plant Sci. 2021;12: 708573.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Lv X, Li H, Chen X, Xiang X, Guo Z, Yu J, et al. The role of calcium-dependent protein kinase in hydrogen peroxide, nitric oxide and ABA-dependent cold acclimation. J Exp Bot. 2018;69:4127–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Ravichandran S, Ragupathy R, Edwards T, Domaratzki M, Cloutier S. MicroRNA-guided regulation of heat stress response in wheat. BMC Genomics. 2019;20:488.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Sun Q, Liu X, Yang J, Liu W, Du Q, Wang H, et al. MicroRNA528 Affects Lodging Resistance of Maize by Regulating Lignin Biosynthesis under Nitrogen-Luxury Conditions. Mol Plant. 2018;11:806–14.

    Article  CAS  PubMed  Google Scholar 

  64. Aydinoglu F. Elucidating the regulatory roles of microRNAs in maize (Zea mays L.) leaf growth response to chilling stress. Planta. 2020;251:38.

    Article  CAS  Google Scholar 

  65. Agurla S, Gahir S, Munemasa S, Murata Y, Raghavendra AS. Mechanism of Stomatal Closure in Plants Exposed to Drought and Cold Stress. In: Iwaya-Inoue M, Sakurai M, Uemura M, editors. Survival Strategies in Extreme Cold and Desiccation. Singapore: Springer Singapore; 2018. p. 215–32.

    Chapter  Google Scholar 

  66. Morales-Navarro S, Pérez-Díaz R, Ortega A, De Marcos A, Mena M, Fenoll C, et al. Overexpression of a SDD1-Like Gene From Wild Tomato Decreases Stomatal Density and Enhances Dehydration Avoidance in Arabidopsis and Cultivated Tomato. Front Plant Sci. 2018;9:940.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65:2119–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Ni Z, Hu Z, Jiang Q, Zhang H. GmNFYA3, a target gene of miR169, is a positive regulator of plant tolerance to drought stress. Plant Mol Biol. 2013;82:113–29.

    Article  CAS  PubMed  Google Scholar 

  69. Yang Y, Zhang X, Su Y, Zou J, Wang Z, Xu L, et al. miRNA alteration is an important mechanism in sugarcane response to low-temperature environment. BMC Genomics. 2017;18:833.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Luan M, Xu M, Lu Y, Zhang L, Fan Y, Wang L. Expression of zma-miR169 miRNAs and their target ZmNF-YA genes in response to abiotic stress in maize leaves. Gene. 2015;555:178–85.

    Article  CAS  PubMed  Google Scholar 

  71. Yu Y, Ni Z, Wang Y, Wan H, Hu Z, Jiang Q, et al. Overexpression of soybean miR169c confers increased drought stress sensitivity in transgenic Arabidopsis thaliana. Plant Sci. 2019;285:68–78.

    Article  CAS  PubMed  Google Scholar 

  72. Luo Y, Wang Z, Ji H, Fang H, Wang S, Tian L, et al. An A rabidopsis homolog of importin β1 is required for ABA response and drought tolerance. Plant J. 2013;75:377–89.

    Article  CAS  PubMed  Google Scholar 

  73. Akdogan G, Tufekci ED, Uranbey S, Unver T. miRNA-based drought regulation in wheat. Funct Integr Genomics. 2016;16:221–33.

    Article  CAS  PubMed  Google Scholar 

  74. Liu D, Song Y, Chen Z, Yu D. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol Plant. 2009;136:223–36.

    Article  CAS  PubMed  Google Scholar 

  75. Fracasso A, Vallino M, Staropoli A, Vinale F, Amaducci S, Carra A. Increased water use efficiency in miR396-downregulated tomato plants. Plant Sci. 2021;303: 110729.

    Article  CAS  PubMed  Google Scholar 

  76. Li S, Tian Y, Wu K, Ye Y, Yu J, Zhang J, et al. Modulating plant growth–metabolism coordination for sustainable agriculture. Nature. 2018;560:595–600.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


Not applicable.


This work was supported by Maize Research Institute Zemun Polje (internal research fund) and Ministry of Science, Technological Development and Innovation under Grant No. 451–03-47/2023–01/200040.

Author information

Authors and Affiliations



MB, DIM and AN conceived and designed the research. MB analysed the data and wrote the main manuscript text. DIM, ND and AN reviewed the manuscript. MB, DIM and AN edited the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Dragana Ignjatović Micić.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Božić, M., Ignjatović Micić, D., Delić, N. et al. Maize miRNAs and their putative target genes involved in chilling stress response in 5-day old seedlings. BMC Genomics 25, 479 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: