Identification and Characterization of Heat-Responsive miRNAs and their Regulatory Network in Maize

doi:10.21203/rs.3.rs-44831/v1

Download PDF

Research article

Identification and Characterization of Heat-Responsive miRNAs and their Regulatory Network in Maize

https://doi.org/10.21203/rs.3.rs-44831/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 03 Nov, 2021

Read the published version in Plant Growth Regulation →

Version 1

posted

You are reading this latest preprint version

Background: MicroRNAs (miRNAs) are a class of small non-coding RNAs, which have been demonstrated to play essential roles in plant growth and development, and in responses to abiotic stress. Heat stress is one of the most serious stresses that affecting crop yield and quality, however, the related regulatory mechanisms of miRNAs remains poorly understanding in maize.

Results: In this study, a total of 340 miRNAs, including 215 known and 125 novel members, were identified from maize seedlings under heat stress (MH) and control conditions (MC) using high-throughput sequencing approach. The 215 known miRNAs can be further divided into 40 different families, and 21-nt miRNAs were found to be most abundant among the known miRNAs. Thirty-five miRNAs, including 26 known and 9 novel members, were significantly different expressed between MC and MH libraries. Furthermore, 174 target genes were predicted to be cleaved by 115 miRNAs using degradome sequencing. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed for these targets to explore the biological function and pathways involved. Based on the relationships of miRNAs, target genes and the enriched results, a regulatory network was constructed for the miRNAs and their respective target genes, and 16 significantly differently expressed miRNAs (DEMs) were involved in the network.

Conclusions: The results revealed novel insights into the roles of miRNAs in heat stress response and provided a useful foundation for understanding the regulatory mechanisms of heat-responsive miRNAs in maize.

Plant Molecular Biology and Genetics

Plant Physiology and Morphology

Maize

miRNA

heat stress

target gene

regulatory network

Maize (Zea mays L.) is one of the most important cereal crops in the world, and has become the model plant in functional genomics and genetic researches [1]. However, its yield and quality are frequently affected by various abiotic stresses, such as drought, heat, and high salinity. During the process of evolution, plants have developed complex and elaborate regulatory mechanisms to cope with various stresses via the changes at morphological, physiological and molecular levels [2, 3]. Studies indicated that a series of genes can be regulated under stress conditions, and can be divided into structural genes and regulatory genes [4]. For example, many stress-responsive transcription factors (TFs) have been reported, which can regulate plant’s stress responses through regulating the expression of a serious of downstream genes [5–10]. Although a number of stress-responsive genes have been identified in plants, the regulatory mechanism and network of stress response still need to be further explored.

MiRNAs are a group of ~ 18–24 nucleotide (nt) non-coding endogenous small RNAs, which are important negative regulators of gene expression at the post-transcriptional level [11, 12]. MiRNA genes are transcribed by RNA polymerase II (Pol II) to produce primary miRNAs (pri-miRNAs), and then the pri-miRNAs are processed into precursor-miRNAs (pre-miRNAs) which can be further processed to generate miRNA/miRNA* duplexes [12–14]. Finally, the duplexes are transported from the nucleus into the cytoplasm, and mature miRNAs are loaded onto argonaute (AGO) protein to form RNA-induced silencing complex (RISC) that guide target mRNA regulation by either cleaving mRNAs or repressing the translation [12, 15, 16]. Since the first miRNA lin4 was discovered in Caenorhabditis elegans (C. elegans) [17], more and more miRNAs were identified in eukaryotes with the development of high-throughput sequencing technology. In addition, degradome sequencing has been widely used to predict the target mRNAs of miRNA, and results also indicated that miRNAs can regulate one or more target genes, thus identification of target genes can contribute to the understanding of the biological roles of miRNAs [18, 19].

To date, a large number of miRNAs have been reported in plants. Increasing evidences indicated that miRNAs have crucial roles in diverse plant developmental stages and stress responses. For example, Arabidopsis miR169 has been shown to target NF-YA2 and NF-YA10 genes to regulate primary root growth [20]. In tomato, overexpressing miR167 can down-regulate the expression of auxin response factors (ARFs) 6 and 8, which result in floral development defects and female sterility [21]. Especially, numerous stress-responsive miRNAs have been identified in various plant species by high-throughput sequencing. For example, expression patterns of drought-responsive miRNAs have been revealed between different rice cultivars [22]. In tomato, Cao et al. (2014) reported hundreds of chilling-responsive miRNAs, and 49 miRNAs were shown significantly differential expression [23]. In wheat, 12 miRNAs have been identified to be responsive to heat stress [24]. Importantly, functional analysis of some stress-responsive miRNAs has also been reported [25–31].

In recent years, heat stress has become one of the most seriously factors that affecting maize yield and quality. Our previously studies demonstrated that maize inbred line CM1 exhibits high tolerance to high temperature [32]. Although the expression patterns of mRNA genes were investigated by RNA sequencing, the roles of miRNAs in heat response in maize inbred line CM1 were still unknown. In this study, expression patterns of miRNAs was investigated between control condition and heat treatment using high-throughput sequencing, and 35 DEMs were identified. The target genes of miRNAs were further predicted using degradome sequencing, and regulatory network of heat-responsive miRNAs and targets was also constructed. Our result will provide valuable information to explore the miRNA-mediated regulatory mechanisms in plant heat response, which is helpful for breeding heat-resistant maize varieties in the future.

High-throughput sequencing and data analysis of small RNAs

In our previously study [32], we reported that maize inbred line CM1 exhibits high tolerance to high temperature stress. To identify miRNAs involved in the response to heat stress, six small RNA libraries were constructed from the leaves of maize inbred line CM1 under heat treatment (MH1, MH2 and MH3) and normal conditions (MC1, MC2 and MC3). These libraries were sequenced using an Illumina Hiseq2500 platform, and the generated data were shown in Table 1. In total, an average of 17,568,632 and 21,192,599 raw reads of the three biological replicates were obtained from MC and MH libraries, respectively. After filtering low quality reads, adapter sequences, and junk reads, 5,522,709 − 19,120,848 valid reads were remained from each sample with 18–25 nt in length, and subsequently used to identify known and novel miRNAs. The length distribution of the small RNAs (sRNAs) was similar between MC and MH libraries (Fig. 1a). We found that the majority of length distribution of the sRNAs was 21 and 22 nt. In addition, the 24 nt sRNAs were also accounted for a large percentage. We noted that the percentage of 18–21 nt sRNAs in MH library was higher than that in MC library, while the opposite result was found for the range of 22–25 nt sRNAs.

Table 1

Statistics of small RNA sequences from the six libraries.
Type	MC1		MC2		MC3		MH1		MH2		MH3
Type	Total	%	Total	%	Total	%	Total	%	Total	%	Total	%
Raw reads	20,636,944	100.00	13,188,120	100.00	18,880,831	100.00	13,303,420	100.00	11,923,314	100.00	38,351,063	100.00
Valid reads	9,315,467	45.14	5,522,709	41.88	7,125,136	37.74	7,693,944	57.83	6,893,172	57.81	19,120,848	49.86
3ADT&length filter	5,507,005	26.69	2,657,816	20.15	4,303,780	22.79	3,314,282	24.91	2,675,581	22.44	12,435,168	32.42
Junk reads	29,539	0.14	29,112	0.22	34,396	0.18	3,4715	0.26	24,114	0.20	62,523	0.16
Rfam	2,521,197	12.22	1,866,557	14.15	2,249,898	11.92	1,202,084	9.04	1,449,839	12.16	4,165,178	10.86
mRNA	3,746,950	18.16	3,647,915	27.66	5,856,919	31.02	1,223,256	9.20	1,090,629	9.15	3,462,254	9.03
Repeats	20,030	0.10	15,058	0.11	17,083	0.09	22,941	0.17	16,194	0.14	52,363	0.14
rRNA	1,840,815	8.92	1,044,648	7.92	1,184,966	6.28	1,009,375	7.59	1,247,375	10.46	3,155,023	8.23
tRNA	569,939	2.76	733,926	5.57	969,077	5.13	130,267	0.98	135,382	1.14	824,895	2.15
snoRNA	9,604	0.05	13,219	0.10	16,023	0.08	4,435	0.03	5,842	0.05	26,990	0.07
snRNA	11,500	0.06	18,398	0.14	21,375	0.11	2,904	0.02	2,510	0.02	9,273	0.02
Other Rfam RNA	89,339	0.43	56,366	0.43	58,457	0.31	55,103	0.41	58,730	0.49	148,997	0.39

Identification of known and novel miRNAs

To identify known miRNAs, the valid reads were used to compare with known plant miRNA precursors or mature miRNAs by alignment against miRBase (v22.0). Totally, 215 known miRNAs were identified from the MC and MH libraries (Table S1). The length size of the known miRNAs was ranged from 18 to 25 nt, and we found that the 21 nt miRNAs were the most abundant, accounting for 55.81% followed by the 22 nt miRNAs (18.60%) (Fig 1b). These known miRNAs were divided into 40 different miRNA families, and the number of each miRNA family has a large range. The miR166 and miR169 families contained the largest number with 22 members, followed by miR156 and miR171 families, which have 21 and 20 miRNAs, respectively. Among the 40 miRNA families, 14 families contained only one member, such as miR11969, miR395 and miR6199 families (Fig. 2). After identification of known miRNAs from the valid data, the remaining unmapped sequences were further aligned to the maize genome to screen potential novel miRNAs, and the RNAfold software was used to predict the secondary structures of the sequences to formulate a hairpin structure or not. Finally, a total of 125 members were identified as the novel miRNAs. The length distribution of the novel miRNAs ranges from 18 to 24 nt, and the 24 nt miRNA was most abundant (64.00%), followed by the 21 nt miRNAs (20.80%) (Table S2).

DEMs in response to heat stress and qRT-PCR validation

Expression patterns of the identified miRNAs were analyzed based on the normalized method reported previously [33]. Pearson correlation analysis was performed to evaluate the correlation among the biological replicates for MC and MH samples. As shown in Fig. 3a, the correlation coefficient (R) of the three biological replicates was greater than 0.913 and 0.945 for MC and MH libraries, respectively, suggesting high reproducibility of the biological replicates of the high-throughput sequencing. DEMs from the known and novel miRNAs were screened followed the differential expression analysis criteria, and 35 miRNAs were shown significantly different expression between MC and MH samples, including 26 known miRNAs and 9 novel miRNAs. Compared with the expression in MC library, 21 miRNAs were up-regulated and 14 members were down-regulated (Fig. 3b and 3c). We noted that the 35 DEMs belong to different miRNA families, and some DEMs were shown with large fold change, such as zma-miR160f-5p_R-1, mtr-miR171c and zma-miR168b-3p_R+1 (Table S3). In addition, the same family miRNAs were shown with same expression patterns, for example, each of the five miR164 family members shown up-regulated, while down-regulated expression were detected for the three miR168 family members. To validate the expression patterns of DEMs, 10 DEMs were selected to validate the high-throughput sequencing results by quantitative real-time polymerase chain reaction (qRT-PCR). As shown in Fig. 4, the expression patterns of the 10 selected miRNAs were similar with the sequencing data, suggesting the reliability of sequencing data and can be used for further analysis.

Identification of miRNA targets by degradome sequencing

The target genes of miRNAs were predicted using degradome sequencing technology. We used the CleaveLand and TargetFinder software to detect potentially sliced targets of miRNAs in the MC and MH libraries. Totally, 21136209 and 11630262 raw reads were obtained for the MC and MH libraries, respectively. And a total of 9395171 (44.45%) and 4861092 (41.80%) of the reads were mapped to the maize transcript in MC and MH libraries, respectively (Table S4). The cleavage products were divided into five categories (0-4) according to the relative abundance of the signatures at the target cleavage sites (Fig. 5). In the MC library, a total of 1563 cleavage events were identified, which 683 (43.70%), 72 (4.61%), 437 (27.96%), 29 (1.86%) and 342 (21.88%) were grouped into category 0, 1, 2, 3, and 4, respectively. In the MH library, a total of 1386 cleavage events were identified, with 596 (43.00%), 97 (6.70%), 459 (33.12%), 9 (0.65%) and 225 (16.23%) cleavage events in category 0, 1, 2, 3, and 4, respectively (Fig. 6).

Of these cleavage events, a total of 174 (527 transcripts) target genes were predicted to be cleaved by 115 miRNAs (104 known miRNAs and 11 novel miRNAs). We found that many of the target genes encode diﬀerent families of transcription factors, such as auxin response factor (ARF), homeobox-leucine zipper protein (HD-Zip), squamosa-promoter binding protein (SBP), NAC domain-containing protein and GRAS protein family (Table S5). Many miRNAs have multiple target genes, and the same family miRNA always have the same type of target genes. For example, the targets of zma-miR160 family encode ARF transcription factors; the targets of zma-miR171 family encode GRAS transcription factors, which might suggest the conserved regulation mechanism of miRNAs. Of the 35 DEMs, degradome sequencing showed that only 17 members had target genes, including 16 known miRNAs and 1 novel miRNAs (Table S3).

GO and KEGG enrichment analysis of miRNA targets

To better understand the potential biological functions of the 174 targets identified by degradome sequencing, GO and KEGG enrichment analysis were performed for these target genes. GO enrichment analysis includes three ontologies of molecular function (MF), cellular component (CC), and biological process (BP) on the basis of the ontological descriptions of GO terms. The 30 most significantly GO terms were shown in Fig. 7. In the molecular function category, the most abundant of targets were enriched in GO terms of DNA binding (GO:0003677), sequence-specific DNA binding transcription factor activity (GO:0003700), and protein dimerization activity (GO:0046983). In the cellular component category, the most abundant of targets were enriched in nucleus (GO:0005634) and nucleolus (GO:0005730). In the biological process category, GO terms with most genes were regulation of transcription, DNA-templated (GO:0006355) and auxin-activated signaling pathway (GO:0009734). Particularly, we also found some genes were enriched in GO term response to hormone (GO:0009725), heat acclimation (GO:0010286), response to hydrogen peroxide (GO:0042542) and etc., which might have potential roles in response to heat stress.

KEGG analysis indicated that a total of 42 pathways were enriched for the 174 target genes, and 30 pathways were shown in Fig. 8 based on the significance of qvalue. Among the 30 pathways, 4 of them were significantly enriched pathways (qvalue<0.05), including ubiquinone and other terpenoid-quinone biosynthesis (ko00130), selenocompound metabolism (ko00450), monobactam biosynthesis (ko00261), sulfur metabolism (ko00920). Of these enriched pathways, the metabolic pathways (ko01100) pathway was the most enriched pathway with 18 genes, followed by biosynthesis of secondary metabolites (ko01110) and ribosome (ko03010) pathways, which contain 10 and 7 members, respectively. In addition, 4 and 3 members were enriched in plant hormone signal transduction (ko04075) pathway and spliceosome (ko03040), and one gene was enriched in MAPK signaling pathway-plant (ko04016) pathway.

Regulatory network of heat-responsive miRNAs in maize

As mentioned above, one miRNA can function on several different target genes, and one gene also can be targeted by one or more miRNAs. To explore the potential regulatory mechanisms of miRNAs under heat stress in maize, regulatory network was constructed using Cytoscape software based on the relationships of miRNAs, target genes and the enriched GO terms and KEGG pathways. As shown in Fig. 9, 11 up-regulated DEMs and 5 down-regulated DEMs were shown to involve in the network, respectively. It can be seen that the target genes were mainly enriched in several GO terms, such as regulation of transcription, DNA-templated (GO:0006355), nucleus (GO:0005634), DNA binding (GO:0003677), and the mainly enriched pathways were metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), and ribosome (ko03010). The relationships of miRNA, targets and the enriched GO terms and pathways were clearly presented in the network. Thus, the regulatory network provided an important foundation for further identification of important miRNAs and target genes in responsive to heat stress.

High temperature has become one of the most serious environmental stresses that affecting plant growth and development. For example, high temperature occurred in the flowering period can lead to serious yield reduction in maize production. Thus, improving the adaptability to high temperature has become a new focus in the process of crop molecular breeding. Plant has developed elaborate regulatory mechanisms at morphological, physiological and molecular levels to cope with heat stress [34, 35], and many studies have paid more attention on gene expression and its roles in response to abiotic stresses [36–39]. MiRNAs are a group of non-coding sRNAs that play essential roles in plant growth, development, and stress responses [20–23, 40, 41]. In our previous study, we showed that inbred line CM1 play high tolerance to heat stress, and revealed the transcriptomic divergence of the maize hybrid An’nong 591, its maternal line (CB25), and paternal line (CM1) under heat stress [32], but the heat-responsive miRNAs and their genetic regulation mechanisms remain poorly known.

In this study, a total of 340 miRNAs (215 known and 125 novel) were identified in the MC and MH libraries. The known miRNAs have been further divided 40 different families, and many miRNA families have been previously reported to be involved in the regulation of stress responses, such as miR156, miR167, miR169, and miR396 [11]. Our result indicated that the length distribution of sRNAs and the known miRNAs varied from 18 to 25 nt, and the 21 nt sequences were the most frequent length and 22 nt ranked in the second. For the novel miRNAs, the 24 nt miRNA was the most abundant (64.00%), and the length distribution ranges from 18 to 24 nt. Some studies reported that the 24 nt sRNAs are more than the 21 nt sRNAs in plants. It is well known that distinct sizes of sRNAs are produced by different members of the DCL family, for example, CL1 is mainly responsible for the biogenesis of 21 nt miRNAs, while DCL3 is involve in the production of 24 nt miRNAs. Thus, functional diversity of different DCL members may responsible for this difference [42, 43]. Compared with the expression of the known miRNAs, no members were shown with high expression level among the 125 novel miRNAs (Table S1 and S2). Among the 35 DEMs, we noted that 8 miRNAs, including zma-miR164f-5p_R-1, zma-miR394a-5p, zma-miR168a-5p, zma-miR168b-3p_R + 1, zma-miR160f-5p_R-1, zma-miR164a-5p, tae-MIR9774-p3_1ss4AG, and zma-miR160a-5p, exhibited high expression in the MC or MH library. This finding was consistent with previous studies that novel miRNAs often expressed at lower levels [19, 44]. The expression patterns of the DEMs were also compared with previous studies. In our study, five conserved miRNA families, including miR156, miR159, miR167, miR168, and miR169 among the DEMs, exhibited down-regulated expression, and the same expression patterns of four families (miR156, miR159, miR167, and miR169) were found among rice heat-responsive miRNAs. However, different expression patterns were also found. For example, miR164 and miR399 were up-regulated in our study while down-regulated expression was detected for the two families in rice [45]. In addition, miR160 exhibited up-regulated expression in rice and our study, while down-regulated expression in response to heat stress was reported in other species [45–47]. It can be concluded that different plant species and heat treatment method are possibly responsible for the differences among these studies. On the other hand, these observations also suggested the complex regulatory mechanisms of plant heat response mediated by miRNAs, although many miRNA families are evolutionarily conserved in different species.

MiRNAs act as negative regulators to regulate the expression of mRNA at the post-transcriptional level [11, 12], and the identification of target genes is an important step to explore miRNA function. In the present study, the target genes of miRNAs were detected by degradome sequencing. A total of 174 (527 transcripts) target genes were predicted to be cleaved by 115 miRNAs, and no targets were identified for the other 225 miRNAs by degradome sequencing. We noted that many of the target genes encode transcription factor families. For example, SPL transcription factors were targeted by miR156 family, and GRAS transcription factors were targeted by miR171 family, suggesting the conserved function and regulation mechanism of the miRNAs. Increasing studies indicated that transcription factors have important roles in plant growth and development, and response to abiotic stresses [9, 48]. In Arabidopsis, NAC transcription factor NTL4 has been demonstrated to participate in a reactive oxygen species (ROS)-mediated positive feedback loop that modulates programmed cell death under heat stress conditions [49]. Chao et al. (2017) reported that overexpression of Arabidopsis SPL transcriptional factors SPL1 or SPL12 can enhanced the thermotolerance in transgenic plants at reproductive stage [50]. Recently, Wang et al. (2019) demonstrated the association of four HD-Zip I members with genetic variations in perennial ryegrass responses to heat stress through investigating the expression patterns between heat-sensitive and heat-tolerant lines [51]. We noted that two heat shock proteins were targeted zma-miR396e-5p and PC-3p-184258_13, respectively, and were only detected in the MH library. Heat shock proteins are a group of molecular chaperones and their expression can be rapidly induced by various environmental stresses, which have been demonstrated to play crucial roles in protecting proteins and cells from damage when under stressed conditions[52–54]. KOG annotation indicated that these genes involved in the posttranslational modification, protein turnover, and chaperones, suggesting the potential roles of zma-miR396e-5p and PC-3p-184258_13 in the regulation of heat response.

GO and KEGG enrichment analysis were performed for the 174 target genes to further explore the biological function of miRNAs. GO analysis indicated that the target genes mainly involved in DNA binding, DNA binding transcription factor activity, protein dimerization activity, nucleus, regulation of transcription, DNA-templated, and auxin-activated signaling pathway. In the GO biological process, we should note that some genes were enriched in the GO terms of response to hormone, response to hydrogen peroxide, heat acclimation, and protein stabilization. KEGG pathway analysis indicated that target genes were mainly involved in the metabolic pathways, biosynthesis of secondary metabolites, ribosome, ubiquinone and other terpenoid-quinone biosynthesis, and plant hormone signal transduction pathways based on the enriched gene number. Some target genes were involved in some important pathways. For example, 3 genes were found to be involved in the spliceosome pathway. Consistent with many studies, our previous study indicated that heat stress can result in a significant increase in the number of alternative splicing (AS) events [32], which can be regulated by different plant developmental stages and environmental stress [55, 56]. In addition, one gene was targeted by a DEM (zma-miR169i-3p_R + 1_1) with down-regulated expression, and involved in the mitogen-activated protein kinase (MAPK) signaling pathway that has been demonstrated to play important roles in plant growth and development, and stress response [57]. These findings might suggest the potential functions of these targets in response to heat stress. To further explore the important target genes and the regulatory mechanisms of miRNAs, the regulatory network was constructed based on the directed relationships of miRNA, target genes, the enriched GO terms and KEGG pathways. A total of 16 DEMs, and their enriched terms and pathways were highlighted in the regulatory network. The regulatory network provides an important foundation for identification of important miRNA and target genes to explore the regulation mechanism of heat-responsive miRNAs in maize.

In this study, high-throughput sequencing was performed to identify heat-responsive miRNA in the leaves of maize inbred line CM1. We identified 340 miRNAs in the MC and MH libraries, including 215 known miRNAs and 125 novel miRNAs. A total of 35 members were shown significantly different expression between control conditions and heat treatment. Target genes of the miRNAs were predicted by degradome sequencing, and 174 genes (527 transcripts) were found to be targeted by 115 miRNAs. GO and KEGG pathway enrichment analysis were performed for the predicted targets, and regulatory network was further constructed based on their relationships of miRNAs and targets. Our results provided valuable information to elucidate the regulatory mechanisms of heat-responsive miRNAs in maize.

Plant materials, heat treatment and RNA isolation

The elite maize inbred line CM1 was obtained from Maize Engineering Technology Research Center of Anhui Province (provided by Professor Qing Ma and Beijiu Cheng), which is a paternal line of maize hybrid An’nong 591. Seedlings of the CM1 were grown in a greenhouse at 28 °C/23 °C (day/night) with a 16-h light/8-h dark photoperiod. To identify heat-responsive miRNAs, seedlings were subjected to 42 °C treatment for 24 h at the three leaves stage. After heat treatment, the third leaves were collected, immediately frozen in liquid nitrogen and stored at −80 °C for RNA isolation. Total RNA was isolated from each sample using a total RNA purification Kit (LC Sciences, USA) following the manufacturer’s protocols. The quantity and purity of the total RNAs were checked by an Agilent Bioanalyzer 2100 (Agilent, USA) and RNA 6000 Nano LabChip Kit (Agilent, USA), and sample with a minimum RNA integrity number (RIN) value > 7.0 was used for further research. Three biological replicates were performed for heat treatment group (MH) and control condition group (MC), respectively.

Small RNA library construction, deep sequencing and data analysis

Small RNA libraries were constructed using a TruSeq Small RNA Sample Prep Kit (Illumina, USA) following the manufacturer’s protocol. 1 µg of total RNA sample was used for library construction, and ligated to 3′ adaptors and 5′ adaptors (Illumina, USA), respectively. Subsequently, the first-strand cDNA was synthesized by reverse transcription PCR (RT-PCR), and the resulting cDNAs were used to PCR amplification using the adaptor primers. After the gel puriﬁcation of PCR products, the sRNA libraries were sequenced on an Illumina Hiseq2500 (Illumina, USA) platform to generate 1×50 bp single-end reads. Raw reads of sequencing data was firstly processed with ACGT101-miR (v4.2) software (LC Sciences, USA) by removing low-quality reads, adapter sequences, junk reads, and read length < 18 nt and >25 nt. The remaining reads were aligned to mRNA, RFam (http://rfam.xfam.org) and Repbase (https://www.girinst.org/repbase/) databases, and then common RNA families (rRNA, tRNA, snRNA, and snoRNA), repeats, and sequences matching maize mRNA were ﬁltered. The last valid sequences with length in 18–25 nt were mapped to specific species precursors in miRBase 22.0 (http://www.mirbase.org/) [58] by BLAST search to identify known miRNAs and novel 3p- and 5p- derived miRNA candidates. One nucleotide mismatch at both 3′ and 5′ ends was allowed in the alignment. Finally, the unmapped sequences were aligned to the maize genome (AGPv4) to identify potential novel miRNAs, and secondary structures were predicted using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi).

Identification of DEMs in response to heat stress

To identify the DEMs between control condition and heat treatment, expression levels of all miRNAs in the six libraries were normalized according to previous study [33]. Differentially expressed sequence was analyzed using a t-test. MiRNAs were considered to be significantly differentially expressed members followed the criteria: |log₂(Fold change)| >1 and p-value <= 0.05. Normalized counts of the expressed miRNAs were transformed by log₂ (norm+1), and used to drawn a heat map with R package pheatmap.

qRT-PCR analysis of miRNAs

To validate the expression of miRNAs by high-throughput sequencing, qRT-PCR analysis was performed. Total RNA used for library construction was used for first-strand cDNA synthesis using Poly (A)-tailing reaction with Mir-X™ miRNA First-Strand Synthesis Kit (Clontech, USA). All primers used in this study were listed in Table S6. qRT-PCR reactions were performed according to our previous studies [32]. The 5S rRNA was used as internal control for the expression analysis of miRNAs [32, 59]. The relative expression level of the miRNAs was calculated by the 2^-ΔΔCT method [60]. Three biological replicates were performed for each treatment.

Identification of miRNA targets by degradome sequencing

In this study, two degradome libraries (MC and MH) were constructed from the mixed total RNAs of the three biological samples, which used for sRNA library construction. The method used for degradome library construction was previously described by Ma et al. [61] with some modiﬁcations. Poly-(A) RNA was purified using poly-T oligo-attached magnetic beads. The 5' adapters were ligated to the 5' end of the 3′ cleavage product of the mRNA, and used to synthesize first-strand cDNA and followed by PCR ampliﬁcation. Finally, 1×50 bp single-end sequencing was performed using an Illumina Hiseq2500 (Illumina, USA) platform. After sequencing, raw reads were filtered to remove adaptor sequences and low quality reads. CleaveLand v3.0 [62] and TargetFinder (http://targetfinder.org/) were used for detect potentially sliced targets of miRNAs. And t-plots were generated to analyze the degradation patterns based on the distribution of signatures and abundance along the transcripts, and divided into five categories (0-4).

GO and KEGG enrichment analysis and regulatory network construction

GO and KEGG enrichment analysis were performed to investigate the biological function of miRNA target genes. GOSeq R package (v1.34.1) [63] was used to identify the functional categories based on molecular function, biological processes and cellular component. Pathway enrichment analysis was performed using KEGG pathway database (http://www.genome.ad.jp/kegg/) [64]. To explore the regulatory network of miRNAs in response to heat stress, the relationships of miRNAs and the target genes enriched in the GO terms and KEGG pathways were used to construct the network using Cytoscape (v3.7.2) software with the similar method as described previously [65]. MiRNAs, targets, enriched GO terms and pathways were represented as nodes, and the directed relationships among them were represented as edges.

AGO: Argonaute; ARF: Auxin response factor; AS: Alternative splicing; BP: Biological process; CC: Cellular component; DCL1: Dicer-like 1; DEM: Differently expressed miRNA; GO: Gene ontology; HD-Zip: Homeobox-leucine zipper protein; KEGG: Kyoto encyclopedia of genes and genomes; MAPK: Mitogen-activated protein kinase; MF: Molecular function; miRNA: microRNA; nt: Nucleotide; qRT-PCR: Quantitative real-time polymerase chain reaction; pre-miRNA: precursor miRNAs; pri-miRNA: primary miRNA; RIN: RNA integrity number; RISC: RNA-induced silencing complex; ROS: Reactive oxygen species; SBP: Squamosa-promoter binding protein; sRNA: Small RNA; TF: Transcription factor.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable.

Availability of data and materials

The generated raw reads were submitted to the Sequence Read Archive (SRA) of NCBI under accession number PRJNA643647.

Competing interests

The authors declare that they have no competing interests.

Funding

This research was supported by the National Key Research and Development Program (2017YFD0101205) and the Science and Technology Major Project of Anhui Province (18030701180). The funders played no role in the study design, review and interpretation of data, or manuscript writing.

Author contribution statement

YZ and QM conceived and designed the research. QYW, TCC, LJX, and JL performed the data analysis. QYW and XGZ conducted the gene expression experiments; YZ and GMH performed the construction of regulatory network. YZ drafted the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We would like to thank Shanshan Xie for her help in this study.

Schnable PS, Ware D, Fulton RS, Stein JC, Wei F, Pasternak S, et al. The B73 maize genome: complexity, diversity, and dynamics. Science. 2009;326(5956):1112-1115.
Tuteja N. Abscisic Acid and abiotic stress signaling. Plant Signal Behav. 2007;2(3):135-138.
Yamaguchi-Shinozaki K, Shinozaki K. Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu Rev Plant Biol. 2006;57:781-803.
Shinozaki K, Yamaguchi-Shinozaki K. Gene networks involved in drought stress response and tolerance. J Exp Bot. 2007;58(2):221-227.
Nakashima K, Tran LS, Van Nguyen D, Fujita M, Maruyama K, Todaka D, et al. Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J. 2007;51(4):617-630.
Liu SX, Wang XL, Wang HW, Xin HB, Yang XH, Yan JB, et al. Genome-wide analysis of ZmDREB genes and their association with natural variation in drought tolerance at seedling stage of Zea mays L. Plos Genet. 2013;9(9):e1003790.
Casaretto JA, El-Kereamy A, Zeng B, Stiegelmeyer SM, Chen X, Bi YM, et al. Expression of OsMYB55 in maize activates stress-responsive genes and enhances heat and drought tolerance. BMC Genomics. 2016;17:312.
Lee DK, Jung H, Jang G, Jeong JS, Kim YS, Ha SH, et al. Overexpression of the OsERF71 transcription factor alters rice root structure and drought resistance. Plant Physiol. 2016;172(1):575-588.
Chen WJ, Zhu T. Networks of transcription factors with roles in environmental stress response. Trends Plant Sci. 2004;9(12):591-596.
Hoang XLT, Nhi DNH, Thu NBA, Thao NP, Tran LSP. Transcription factors and their roles in signal transduction in plants under abiotic stresses. Curr Genomics. 2017;18(6):483-497.
Covarrubias AA, Reyes JL. Post-transcriptional gene regulation of salinity and drought responses by plant microRNAs. Plant Cell Environ. 2010;33(4):481-489.
Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281-297.
Krutzfeldt J, Poy MN, Stoffel M. Strategies to determine the biological function of microRNAs. Nature Genetics. 2006;38:S14-S19.
Kurihara Y, Watanabe Y. Arabidopsis micro-RNA biogenesis through Dicer-like 1 protein functions. P Natl Acad Sci USA. 2004;101(34):12753-12758.
Mallory AC, Elmayan T, Vaucheret H. MicroRNA maturation and action--the expanding roles of ARGONAUTEs. Curr Opin Plant Biol. 2008;11(5):560-566.
Baumberger N, Baulcombe DC. Arabidopsis ARGONAUTE1 is an RNA Slicer that selectively recruits microRNAs and short interfering RNAs. Proc Natl Acad Sci U S A. 2005;102(33):11928-11933.
Lee RC, Feinbaum RL, Ambros V. The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell. 1993;75(5):843-854.
Wang R, Yang Z, Fei Y, Feng J, Zhu H, Huang F, et al. Construction and analysis of degradome-dependent microRNA regulatory networks in soybean. BMC Genomics. 2019; 20(1):534.
Li H, Peng T, Wang Q, Wu Y, Chang J, Zhang M, et al. Development of incompletely fused carpels in maize ovary revealed by miRNA, target gene and phytohormone analysis. Front Plant Sci. 2017;8:463.
Sorin C, Declerck M, Christ A, Blein T, Ma L, Lelandais-Briere C, et al. A miR169 isoform regulates specific NF-YA targets and root architecture in Arabidopsis. New Phytol. 2014;202(4):1197-1211.
Liu N, Wu S, Van Houten J, Wang Y, Ding B, Fei Z, et al. Down-regulation of AUXIN RESPONSE FACTORS 6 and 8 by microRNA 167 leads to floral development defects and female sterility in tomato. J Exp Bot. 2014;65(9):2507-2520.
Balyan S, Kumar M, Mutum RD, Raghuvanshi U, Agarwal P, Mathur S, et al. Identification of miRNA-mediated drought responsive multi-tiered regulatory network in drought tolerant rice, Nagina 22. Sci Rep. 2017;7(1):15446.
Cao X, Wu Z, Jiang FL, Zhou R, Yang ZE. Identification of chilling stress-responsive tomato microRNAs and their target genes by high-throughput sequencing and degradome analysis. BMC Genomics. 2014; 15(1):1130.
Xin M, Wang Y, Yao Y, Xie C, Peng H, Ni Z, et al. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 2010;10:123.
Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65(8):2119-2135.
Zhang X, Zou Z, Gong P, Zhang J, Ziaf K, Li H, et al. Over-expression of microRNA169 confers enhanced drought tolerance to tomato. Biotechnol Lett. 2011;33(2):403-409.
Yang CH, Li DY, Mao DH, Liu X, Ji CJ, Li XB, et al. Overexpression of microRNA319 impacts leaf morphogenesis and leads to enhanced cold tolerance in rice (Oryza sativa L.). Plant Cell Environ. 2013;36(12):2207-2218.
Denver B, Ullah H. miR393s regulate salt stress response pathway in Arabidopsis thaliana through scaffold protein RACK1A mediated ABA signaling pathways. Plant Signal Behav. 2019;14:1600394.
Wang Y, Sun F, Cao H, Peng H, Ni Z, Sun Q, et al. TamiR159 directed wheat TaGAMYB cleavage and its involvement in anther development and heat response. PLoS One. 2012;7(11):e48445.
Yue E, Liu Z, Li C, Li Y, Liu Q, Xu JH. Overexpression of miR529a confers enhanced resistance to oxidative stress in rice (Oryza sativa L.). Plant Cell Rep. 2017;36(7):1171-1182.
Dong CH, Pei HX. Over-expression of miR397 improves plant tolerance to cold stress in Arabidopsis thaliana. J Plant Biol. 2014;57(4):209-217.
Zhao Y, Hu F, Zhang X, Wei Q, Dong J, Bo C, et al. Comparative transcriptome analysis reveals important roles of nonadditive genes in maize hybrid An'nong 591 under heat stress. BMC Plant Biol. 2019;19(1):273.
Li X, Shahid MQ, Wu J, Wang L, Liu X, Lu Y. Comparative small RNA analysis of pollen development in autotetraploid and diploid rice. Int J Mol Sci. 2016;17(4):499.
Zandalinas SI, Mittler R, Balfagon D, Arbona V, Gomez-Cadenas A. Plant adaptations to the combination of drought and high temperatures. Physiol Plant. 2018;162(1):2-12.
Ravichandran S, Ragupathy R, Edwards T, Domaratzki M, Cloutier S. MicroRNA-guided regulation of heat stress response in wheat. BMC Genomics. 2019;20(1):488.
González-Schain N, Dreni L, Lawas LM, Galbiati M, Colombo L, Heuer S, et al. Genome-wide transcriptome analysis during anthesis reveals new insights into the molecular basis of heat stress responses in tolerant and sensitive rice varieties. Plant Cell Physiol. 2016;57(1):57-68.
Wang B, Liu C, Zhang D, He C, Zhang J, Li Z. Effects of maize organ-specific drought stress response on yields from transcriptome analysis. BMC Plant Biol. 2019;19(1):335.
Ueda A, Kathiresan A, Bennett J, Takabe T. Comparative transcriptome analyses of barley and rice under salt stress. Theor Appl Genet. 2006;112(7):1286-1294.
Cui JH, Ren GZ, Qiao HY, Xiang XD, Huang LS, Chang JH. Comparative transcriptome analysis of seedling stage of two sorghum cultivars under salt stress. J Plant Growth Regul. 2018;37(3):986-998.
Sunkar R, Chinnusamy V, Zhu JH, Zhu JK. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007;12(7):301-309.
Ding D, Zhang LF, Wang H, Liu ZJ, Zhang ZX, Zheng YL. Differential expression of miRNAs in response to salt stress in maize roots. Ann Bot. 2009;103(1):29-38.
Huang X, Feng JJ, Wang R, Zhang HS, Huang J. Comparative analysis of microRNAs and their targets in the roots of two cultivars with contrasting salt tolerance in rice (Oryza sativa L.). Plant Growth Regul. 2019;87(1):139-148.
Rogers K, Chen X. Biogenesis, turnover, and mode of action of plant microRNAs. Plant Cell. 2013;25(7):2383-2399.
Liu HJ, Qin C, Chen Z, Zuo T, Yang XR, Zhou HK, et al. Identification of miRNAs and their target genes in developing maize ears by combined small RNA and degradome sequencing. BMC Genomics. 2014;15:25.
Li J, Wu LQ, Zheng WY, Wang RF, Yang LX. Genome-wide identification of microRNAs responsive to high temperature in rice (Oryza sativa) by high-throughput deep sequencing. J Agron Crop Sci. 2015;201(5):379-388.
Pan Y, Niu MY, Liang JS, Lin EP, Tong ZK, Zhang JH. Identification of heat-responsive miRNAs to reveal the miRNA-mediated regulatory network of heat stress response in Betula luminifera. Trees. 2017;31(5):1635-1652.
Wang RH, Xu L, Zhu XW, Zhai LL, Wang Y, Yu RG, et al. Transcriptome-wide characterization of novel and heat-stress-responsive microRNAs in radish (Raphanus Sativus L.) using next-generation sequencing. Plant Mol Biol Rep. 2015;33(4):867-880.
Huang GT, Ma SL, Bai LP, Zhang L, Ma H, Jia P, et al. Signal transduction during cold, salt, and drought stresses in plants. Mol Biol Rep. 2012;39(2):969-987.
Lee S, Lee HJ, Huh SU, Paek KH, Ha JH, Park CM. The Arabidopsis NAC transcription factor NTL4 participates in a positive feedback loop that induces programmed cell death under heat stress conditions. Plant Sci. 2014;227:76-83.
Chao LM, Liu YQ, Chen DY, Xue XY, Mao YB, Chen XY. Arabidopsis transcription factors SPL1 and SPL12 confer plant thermotolerance at reproductive stage. Mol Plant. 2017;10(5):735-748.
Wang J, Zhuang LL, Zhang J, Yu JJ, Yang ZM, Huang BR. Identification and characterization of novel homeodomain leucine zipper (HD-Zip) transcription factors associated with heat tolerance in perennial ryegrass. Environ Exp Bot. 2019;160:1-11.
Kadota Y, Shirasu K. The HSP90 complex of plants. Biochim Biophys Acta. 2012;1823(3):689-697.
Wang W, Vinocur B, Shoseyov O, Altman A. Role of plant heat-shock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 2004;9(5):244-252.
Campbell JL, Klueva NY, Zheng HG, Nieto-Sotelo J, Ho TD, Nguyen HT. Cloning of new members of heat shock protein HSP101 gene family in wheat (Triticum aestivum (L.) Moench) inducible by heat, dehydration, and ABA(1). Biochim Biophys Acta. 2001;1517(2):270-277.
Laloum T, Martin G, Duque P. Alternative splicing control of abiotic stress responses. Trends Plant Sci. 2018;23(2):140-150.
Filichkin S, Priest HD, Megraw M, Mockler TC. Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Curr Opin Plant Biol. 2015;24:125-135.
Danquah A, de Zelicourt A, Colcombet J, Hirt H. The role of ABA and MAPK signaling pathways in plant abiotic stress responses. Biotechnol Adv. 2014;32(1):40-52.
Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36(Database issue):D154-158.
Nie Z, Ren Z, Wang L, Su S, Wei X, Zhang X, et al. Genome-wide identification of microRNAs responding to early stages of phosphate deficiency in maize. Physiol Plant. 2016;157(2):161-174.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method. Methods. 2001;25(4):402-408.
Ma ZR, Coruh C, Axtell MJ. Arabidopsis lyrata small RNAs: Transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell. 2010;22(4):1090-1103.
Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009;25(1):130-131.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.
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. 2008;36(Database issue):D480-484.
Xu YJ, Zhu SW, Liu F, Wang W, Wang XW, Han GM, et al. Identification of arbuscular mycorrhiza fungi responsive microRNAs and their regulatory network in maize. Int J Mol Sci. 2018;19(10):3201.

Download PDF

Journal Publication

published 03 Nov, 2021

Read the published version in Plant Growth Regulation →

Version 1

posted

You are reading this latest preprint version

Identification and Characterization of Heat-Responsive miRNAs and their Regulatory Network in Maize

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1