Skip to main content
Download PDF

MicroRNA-guided regulation of heat stress response in wheat

BMC Genomics volume 20, Article number: 488 (2019) Cite this article

Abstract

Background

With rising global temperature, understanding plants’ adaptation to heat stress has implications in plant breeding. MicroRNAs (miRNAs) are small, non-coding, regulatory RNAs guiding gene expression at the post-transcriptional level. In this study, small RNAs and the degradome (parallel analysis of RNA ends) of leaf tissues collected from control and heat-stressed wheat plants immediately at the end of the stress period and 1 and 4 days later were analysed.

Results

Sequencing of 24 small RNA libraries produced 55.2 M reads while 404 M reads were obtained from the corresponding 24 PARE libraries. From these, 202 miRNAs were ascertained, of which mature miRNA evidence was obtained for 104 and 36 were found to be differentially expressed after heat stress. The PARE analysis identified 589 transcripts targeted by 84 of the ascertained miRNAs. PARE sequencing validated the targets of the conserved members of miRNA156, miR166 and miR393 families as squamosa promoter-binding-like, homeobox leucine-zipper and transport inhibitor responsive proteins, respectively. Heat stress responsive miRNA targeted superoxide dismutases and an array of homeobox leucine-zipper proteins, F-box proteins and protein kinases. Query of miRNA targets to interactome databases revealed a predominant association of stress responses such as signalling, antioxidant activity and ubiquitination to superoxide dismutases, F-box proteins, pentatricopeptide repeat-containing proteins and mitochondrial transcription termination factor-like proteins.

Conclusion

The interlaced data set generated in this study identified and validated heat stress regulated miRNAs and their target genes associated with thermotolerance. Such accurate identification and validation of miRNAs and their target genes are essential to develop novel regulatory gene-based breeding strategies.

Background

Rising global temperatures with increasing occurrences of heat waves negatively impact crop productivity and global food production [1]. Heat stress disrupts membranes and inactivates proteins by unfolding, misfolding and aggregation [2]. Such metabolic imbalances have an adverse impact on all aspects of plant growth and impair several physiological and developmental processes including photosynthesis, respiration, reproduction, grain filling and yield [reviewed by [3] ]. Plants perceive and respond to heat stress using an array of stress signaling pathways to regulate their metabolism and cell function. Perception of heat stress begins at the plasma membrane, with subsequent alterations in membrane fluidity, lipid saturation and protein stability [4]. Changes in membrane fluidity initiate Ca2+ influx and rapid accumulation of reactive oxygen species (ROS) such as H2O2 [4, 5]. Stress-induced Ca2+ influxes and ROS activate mitogen-activated protein kinase (MAPK) and calcium-dependent protein kinase (CDPK) stress signaling pathways [6, 7] to regulate the expression of heat stress transcription factors [8] and WRKY proteins in order to confer thermotolerance [9].

Small RNAs (sRNAs) are a class of endogenous non-coding RNAs that regulate gene expression at transcriptional and post-transcriptional levels. In plants, sRNAs are classified based on their biogenesis and function. MicroRNAs (miRNAs) are primarily 21 nucleotide (nt)-long molecules derived from single stranded precursors with a hairpin structure. miRNA coding genes are transcribed by RNA polymerase II (Pol II) to generate primary miRNAs (pri-miRNAs) that fold into characteristic self-complementary stem-loop secondary structures [10]. The pri-miRNAs are cleaved at least twice by Dicer-like1 (DCL1) endonucleases to yield the precursor-miRNAs (pre-miRNAs) that are further processed to release miRNA/miRNA* duplexes. Mature miRNAs are loaded onto Argonaute 1 (AGO1) to form miRNA-induced silencing complexes (miRISCs) that guide target-specific mRNA cleavage or translational repression [reviewed by [11] ]. In contrast, short interfering RNAs (siRNAs) are transcribed from transposons and repetitive regions by RNA polymerase IV (Pol IV). The single stranded RNA transcripts are copied into double-stranded RNA (dsRNA) molecules by RNA-dependent RNA polymerase 2 (RDR2) that are cleaved by DCL3 to produce the 24-nt siRNAs. These siRNAs are loaded onto AGO4 to induce sequence-specific transcriptional gene silencing (TGS) by recruiting methyltransferases [12].

miRNA-guided post-transcriptional gene silencing (PTGS) has emerged as a mechanism to reprogram the expression of development and stress associated genes. PTGS regulates nucleic acid binding [13, 14], phytohormone biosynthesis [15, 16], chlorophyll biosynthesis [17], flowering [18], root differentiation [19], nutrient homeostasis [20] and stress responses [21]. Earlier studies on Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa) identified a high degree of conservation among miRNAs and their mRNA targets [22]. While some miRNAs and their targets are evolutionarily conserved, many are lineage or even species-specific [11, 23]. The emergence of high-throughput sequencing technologies and genome-wide analysis has facilitated the identification of new miRNA variants from an ever increasing number of plant species [24, 25]. Similarly, computational approaches remain a valuable tool to predict plant miRNA targets by complementary base pairing in cDNA and whole genome sequences [22, 26]. However, experimental validation of miRNA-guided cleavage site is still required considering the high rate of false positive target prediction [27, 28]. Parallel analysis of RNA ends (PARE), also known as degradome analysis, has enabled the massive cloning and sequencing of 5′ ends of cleaved or uncapped mRNA to validate many miRNA-guided cleavage sites [29].

Advancements in sequencing technologies have revealed rapid divergence and high diversity in miRNAs and their targets [24, 25, 30]. Stress responsive miRNAs and their targets have been identified from several cereal crops including wheat, rice, maize and barley (reviewed in [31]). Despite the identification of miRNA from an increasing number of plant species and stress responses, experimental validation of miRNA targets is lacking in polyploid crops including wheat. However, studies in Arabidopsis and rice have experimentally validated a few miRNA targets regulating thermotolerence. For example, in Arabidopsis miR156-overexpressing plants exhibited enhanced tolerance to heat stress and was shown to be required for plant heat stress memory [32]. In Arbidopsis, plants expressing miR398-resistant version of CSD (copper/zinc superoxide dismutase) were more sensitive to heat stress. miR159 is downregulated in response to heat stress in order to regulate MYB-like family transcription factors. Transgenic rice plants overexpressing miR159 were more sensitive to heat stress [33].

In wheat, heat stress responsive miRNAs have been reported from several susceptible and resistant cultivars, primarily investigating early responses within 24 h of treatment and the role of isomiRs and differentially accumulated mature miRNAs among families was largely unnoticed [34,35,36]. Further, miRNA target identification in previous studies [34,35,36,37] have relied on prediction models based on draft genome assemblies. Although advancements in sequencing technologies have benefited the identification of several stress responsive miRNAs, lack of a high quality reference sequence and high confidence annotations have limited previous studies with the identification of miRNAs, isomiRs and validation of their targets.

Previous studies in cotton [38] and cowpea [39] have also shown genotype-specific miRNA expression profiles in responses to stress. In cereal crops such as rice, domestication has led to specific gain and/or loss of miRNA regulated gene expression [40]. As new miRNA variants are spawned and lost frequently, especially in polyploid crops such as wheat, accurate identification and experimental validation of miRNAs and their targets in response to heat stress is essential to develop novel regulatory-RNA-based plant breeding strategies. The recent release of the Chinese Spring wheat reference sequence [41] with high-confidence annotations has enabled the accurate identification of miRNAs and validation of their targets. The gold standard reference sequence can further be utilized to functionally validate miRNA targets.

We previously showed the temporal expression profile of wheat cv Glenlea with a pronounced differential expression of miRNA occurring immediately after heat stress until 3 days post treatment [37]. Here, wheat cv Chinese Spring was used to generate the interlaced data set of small RNA and PARE. The recently released wheat reference sequence [41] of cv Chinese Spring was used to enhance the accuracy of miRNA and isomiRs annotation and to validate heat stress regulated miRNAs target genes associated with thermotolerance. Heat stress regulated miRNAs were highly altered immediately after heat stress to regulate thermotolerance but expression levels largely returned to control levels during the 4 day recovery period following the end of the stress. Differential accumulation of isomiRs indicates highly regulated post transcriptional modifications dictating the accumulation of specific mature miRNAs in response to heat stress. Degradome sequencing revealed chloroplast and mitochondria localized pentatricopeptide repeat-containing protein and mitochondrial transcription termination factor-like protein targets that were previously not known to undergo miRNA-guided cleavage.

Methods

Plant materials and growth conditions

Triticum aestivum cv Chinese Spring seeds were obtained from Dr. Eric Kerber (Agriculture and Agri-Food Canada, Winnipeg, Canada). A single Chinese Spring plant was selfed to produce the seeds for this experiment. A total of 144 plants were grown in a PG-40 growth cabinet (Conviron Technologies, Winnipeg, Canada) under long-day conditions of 16 h of light (300 μmol m− 2 s− 1) at 18 °C and 8 h of darkness at 16 °C. At the boot stage, a set of 72 plants were exposed to a 37 °C heat stress for 5 days while the other 72 plants remained under the control conditions mentioned above. The heat stressed plants were amply watered twice per day during the 5-day stress period to avoid a confounding effect with drought stress. After the 5-day stress exposure, the plants were returned to the control conditions. The experiment was setup as a completely randomized design with four biological replicates, two treatments (control and heat) with 18 plants per treatment and replicate.

Sampling and RNA extraction

Leaf tissue was collected from all heat-stressed and control individual plants immediately at the end of the stress period, i.e., time point zero day after treatment (0 DAT), and at one and four days after treatment (1 and 4 DAT). Leaf tissue from individual plants was immediately frozen in liquid nitrogen. Total RNA was isolated from all 432 leaf tissue samples using TRI-reagent (Ambion, Naugatuck, CT) as previously described [37]. Total RNA quality and quantity were obtained using an Agilent 2100 Bioanalyzer with the RNA 6000 Nano chip (Agilent Technologies, Santa Clara, CA). Equimolar amounts of RNA from the 18 individual leaf samples representing a replicate, treatment and time point were pooled to create the 24 total RNA samples representing the four biological replicates, two treatments and three time points. These 24 total RNA samples were used to construct sRNA and degradome libraries as described below.

Small RNA library preparation, data processing and annotation of miRNA

sRNA libraries were constructed from low molecular weight (LMW) RNA isolated from 72 μg of total RNA as described by Lu and Souret (2010) [42]. Briefly, sRNA purified from LMW RNA was size-selected on a 15% polyacrylamide gel. A gel slice, corresponding to the 16 to 30 nt size range was eluted to obtain the sRNA fractions that were subsequently quantified on an Agilent 2100 Bioanalyzer with the Small RNA kit. The TruSeq® small RNA sequencing adapters (Illumina, San Diego, CA) were ligated to the size-selected sRNA samples and libraries were constructed using Illumina TruSeq® small RNA kit following the manufacturer’s instructions. The libraries were quantified on an Agilent 2100 Bioanalyzer using a High Sensitivity DNA kit. Six pools, each containing equimolar amounts of 12 random libraries and where each library was present in three different pools, were resolved on a 6% polyacrylamide gel from which the 145–160 bp fragments, corresponding to the adapter-ligated sRNA, were size-selected. The ethanol precipitated library pools were resuspended and quantified on an Agilent 2100 Bioanalyzer with the DNA 1000 kit. The library pools were normalized to 2 nM and sequenced on a MiSeq for 50 cycles using a MiSeq Reagent Kit v2 (Illumina). A total of six runs were performed to sequence the six pools where each library was sequenced in three separate runs.

A tool chain for processing small RNA reads [37] was used with slight modifications. miRNA annotation was carried out following the guidelines of [43] and A Kozomara and S Griffiths-Jones [44] for mature miRNA and precursor evidence, respectively. Distinct tags represented by ≥10 reads per million (RPM) in at least one library were mapped to the pre-miRNA loci identified from the hexaploid bread wheat genome (IWGSC RefSeq v1.0., 2018) using Bowtie2 [45]. Distinct tags with perfect matches were extracted for hairpin structure prediction using RNAfold [46]. Star sequences were predicted using a custom Python script.

Only the distinct tags that showed evidence of biosynthesis based on matching precursor that could fold into secondary hairpin structure with <− 0.2 kcal/mol/nt were retained. In addition, presence of the complementary star sequence on the opposite strand with a 2-bp offset on the 3′-end of the hairpin structure was assessed. Distinct tags represented by ≥10 RPM in at least one library were also matched to previously annotated miRNA and star sequences [47]. Distinct tags with strong evidence of biosynthesis, including precursor and star sequence supports, were selected for differential expression analysis and to identify miRNA target using degradome data.

Read counts were grouped to define treatment combination and normalized [48] to account for differences in library sequencing depth and RNA composition using edgeR [49]. MDS plots were constructed to check for outliers. After estimating dispersion a general linear model (GLM) was used with glmLRT function to identify miRNAs differentially expressed across treatments and over time. miRNAs with false discovery rate (FDR) < 0.05 [50] were considered differentially expressed. Hierarchical clustering ananlysis was performed using the hclust function in R. Heat maps were generated with normalized log2 read counts using the maximum algorithm and complete linkage method on ClustVis [51].

Sequence extraction and multiple sequence alignment

The precursor sequences from model plants Arabidopsis thaliana (ath), Brachypodium distachyon (bdi) and cereal crops including Oryza sativa (osa), Sorghum bicolor (sbi), Hordeum vulgare (hvu) and Zea mays (zma) were extracted from miRBase release 22. Multiple sequence alignment was performed using MUSCLE with default parameters.

Degradome

Degradome libraries were constructed from 75 μg of total RNA as previously described [52]. cDNA synthesized from 5′ adapter-ligated poly(A) RNA was digested with MmeI. After 3′ adapter ligation, double-stranded DNA (dsDNA) was purified on a gel and the libraries were PCR-amplified. The 128-bp fragment was size-selected and verified on an Agilent 2100 Bioanalyzer using a High Sensitivity DNA kit. dsDNA concentrations were determined using the Qubit High Sensitivity kit and twelve random libraries were pooled in equimolar amounts. Sequencing of the degradome libraries was performed by the McGill University and Genome Quebec Innovation Centre. The clustering was done on an Illumina cBot and the sequencing was performed on a HiSeq 2000 for 50 cycles. A phiX library was used as control at 10% level. bcl2fastq v1.8.4 (Illumina) was used to demultiplex the sample reads and generate the fastq files.

Raw reads were quality-checked and adapter sequences were removed using a custom Perl script. Reads ranging from 19 to 21 nt were extracted using Geneious (Biomatters Ltd., Auckland, NZ) and were converted to distinct tags using Tally [53]. The distinct tags were mapped to the coding sequences of wheat (IWGSC RefSeq v1.0, 2018) using sPARTA [54]. The previously annotated miRNA sequences were queried for target prediction using the sPARTA’s miRferno target prediction module. miRNA targets and their abundance were extracted and target plots were constructed using the R package ggplot2. miRNA targets with a threshold P value < 0.05 and read abundance > 4 were retained for further analysis.

Prediction of heat stress associated regulatory networks

Degradome validated miRNA targets were extracted and translated to query the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING version 10.5) database [55]. Orthologous proteins matching rice (Oryza sativa ssp. indica L.) with identity > 50% were used to associate coexpression and interactome data. A high confidence of 0.7 was used to construct the regulatory networks.

Data availability and accession numbers

The smallRNA and degradome (PARE) data were submitted to the NCBI Gene Expression Omnibus under accession number GSE113358.

Results

Small RNA profile and annotation of microRNAs

A total of 24 sRNA libraries, obtained from leaf tissues harvested from four biological replicates representing control and heat stressed plants sampled at 0, 1 and 4 DAT were sequenced. After quality check and removal of low quality reads, a total of 55.2 M 18–24 nt-long reads corresponding to 18.7 M distinct tags were retained for miRNA annotation (Additional file 1). The number of high quality raw reads and distinct tags was comparable regardless of the treatment and sampling time (Additional file 2: Figure S1a,b). A size distribution analysis of raw reads showed that the 21- and 24-nt size classes exceeded other size classes (Additional file 2: Figure S1c) and the proportion of distinct tags was skewed towards the 24-nt size class (Additional file 2: Figure S1d).

The adapter trimmed, high quality reads were filtered and annotated using the in-house pipeline described in the methods to identify a total of 202 mature miRNAs belonging to 37 miRNA families (Additional file 1, Additional file 3: Figure S2a). Evidence for complementary miRNA star (miRNA*) sequence on the other arm of the precursor with 2-bp overhang was found for 104 mature miRNAs. The majority (161/202) matched more than one pre-miRNA coding loci across all sub-genomes (Additional file 4). Here, miRNAs derived from a single precursor are termed isomiRs. For example, all five isomiRs of miR528 were processed from the single pre-miRNA miR538–1 transcribed from locus Chr5A:663320619–663,320,737 (Fig. 1a). Similarly, all nine 19–21 nt-long isomiRs of miR9662 were processed from pre-miR9662–1 aligning to Chr6D:95858457–95,858,576 (Fig. 1b). IsomiRs processed from these two pre-miRNAs differed among themselves by minor shifts of one to three nucleotides (Additional file 4; Fig. 1a,b). A search for orthologous pre-miR528 in miRBase (release 22), identified a single miRNA coding locus in rice (Oryza sativa) and sorghum (Sorghum bicolor) while maize (Zea mays) carried two loci on chromosomes 1 and 9, respectively. A similar search for pre-miR9662 failed to identify orthologous precursors.

Fig. 1
figure 1

Foldback structures and sequence conservation among miRNA families. Schematic drawings of the foldback structure of pre-miRNA528–1 (a) and pre-miRNA9662–1 (b) and the multiple alignments of their encoded isomiRs. Alignments of miR156 (c), miR166 (d), miR393 (e) and miR528 (f) precursors from Arabidopsis thaliana (ath), Brachypodium distachyon (bdi), Oryza sativa (osa), Sorghum bicolor (sbi), Hordeum vulgare (hvu), Zea mays (zma) and Triticum aestivum (tae; our study). Identity is color-coded as follows: dark green (100%), light green (80–100%), orange (60–80%) and grey (< 60%)

Full size image

Multiple sequence alignment of orthologous precursors revealed a high degree of sequence conservation primarily in the miRNA/miRNA* region of the foldback structure (Fig. 1c-f). The position of sequence conservation varied among miRNA families in the foldback structures. Orthologous precursors of miR156 were highly conserved at the 5p arm of the foldback structure and all mature miRNAs were derived from the same (Additional file 4; Fig. 1c). Precursors of miRNA166 were highly conserved at the 3p arm but not so at the 5p arm. Of the 23 mature miR166, 20 were derived from the 3p arm. The less diverse and lineage specific families such as miR393 and miR528 were conserved at both 5p and 3p arms of the precursor (Fig. 1e,f).

Among the 37 miRNA families, 24 families comprised more than one mature miRNAs (Additional file 3: Figure S2a). Families miR1135 and miR166 were the most diverse with 27 and 23 mature miRNAs, respectively. miR1137 comprised only 21 and 24-nt size classes, while 58% (7/12) of miR1118 were 24-nt long (Additional file 3: Figure S2b). The conserved 21-nt size class represented 57% of the miRNA distinct tags corresponding to 91% of the raw reads (Additional file 5).

Relative abundance and differential expression of heat stress regulated miRNAs

Of the 37 families, 29 were weakly expressed representing less than 1% of the reads (Additional file 6: Figure S4a). Despite the presence of 27 distinct miRNAs in family miR1135, this family accounted for only 0.5% of the reads. In contrast, family miR166 was the most abundant totalling 79% of the reads, including six members with more than 1000 RPM in all libraries and the remaining 17 ranging from 10 to 1000 RPM (Additional files 7 and 8). miR166–4.1-3p was the most abundant with an average read count of 414,544 RPM. Overall, ten miRNA families had read counts exceeding 1000 RPM which globally represented more than 90% of all reads (Additional file 6: Figure S4b).

Difference in accumulation was also observed among miRNA processed from the same precursor, such as the case with miR528 and miR9662 where isomiRs varied from 33 to 660 RPM and 14 to 13,658 RPM, respectively (Additional file 8; Fig. 2a,b). Analysis of differences in accumulation of mature miRNAs and their complementary miRNAs* showed the abundance of mature miR9662–1.4-3p to be considerably lower (8738 RPM) than its complementary star miR9662–1.3-5p (82,966 RPM) (Additional file 9: Figure S6a). Inversely, although not nearly as pronounced, the abundance of mature miR9662–1.7-3p was lower (1928 RPM) than its complementary star miR9662–1.8-5p (2591 RPM) (Additional file 9: Figure S6b).

Fig. 2
figure 2

Box plots of the relative abundance of isomiRs belonging to miRNA family (a) miR528 and (b) miR9662. The median is shown by a horizontal black line. Boxes represent the 25th and 75th percentiles while error bars indicate the 10th and 90th percentiles

Full size image

An MDS plot constructed using the normalized read counts showed the relative similarities among the biological replicates (Additional file 10: Figure. S7). Of the 202 mature miRNAs, 36 belonging to 18 families were differentially expressed upon heat stress (Additional file 11; Fig. 3a). Cluster dendrogram of differentially expressed miRNA revealed a temporal expression pattern where the heat-stress induced alterations in miRNA levels were more pronounced immediately at the end of the heat stress (0 DAT; Fig. 3b). For the majority of the differentially expressed miRNAs, fold change expression differences were smaller during the recovery phase at 1 and 4 DAT compared to 0 DAT (Fig. 3b). Of the 36 differentially expressed miRNAs, 18 were downregulated immediately after heat stress (0 DAT; Additional file 11). Hierarchical clustering of differentially expressed miRNAs showed that miR528 and miR1135 were downregulated after heat stress (Fig. 3b). In contrast, members of miR169 and miR9772 were upregulated by the end of the 5-day heat stress period (Fig. 3b; Additional file 11). Levels of miR1121–33-5p were reduced in response to heat stress, but unlike others, they did not return to the same level as the control during the recovery period (1–4 DAT; Additional file 11). Members of miR166, miR1118, miR1120, miR1121 and miR1122 were also differentially regulated in response to heat stress (Fig. 3a,b). Comparison between heat 0 and 4 DAT identified nine differentially expressed miRNAs: miR528–1.1-5p, miR528–1.3-3p, miR166–11-5p, miR528–1.2-5p, miR399–18-3p, miR166–13-3p, miR167–1-5p, miR1120–2337-3p, and miR166–5.1-3p (Fig. 3b). Comparison between control 0 and 4 DAT did not reveal any differentially expressed miRNA.

Fig. 3
figure 3

Differential expression of miRNAs following heat stress. (a) Distribution of the 36 miRNAs differentially expressed in wheat leaf tissue after heat stress into 18 families (FDR < 0.05). (b) Heat map of the 36 differentially expressed miRNAs in response to heat stress. The blue color represents upregulation and red is for downregulation. The dendrogram represents the hierarchical clustering of the control and heat treatments at the three sampling time points of 0, 1 and 4 days after treatment (DAT)

Full size image

Identification of cleaved miRNA targets through sequencing of the degradome

Sequencing of the 24 degradome libraries yielded a total of 404 M reads. After quality check and adapter trimming, 370 M 19–21 nt-long reads were analyzed. A total 84 miRNAs were identified to target 589 transcripts corresponding to 206 protein coding loci (P value < 0.05; PARE abundance > 4; Additional file 12) thus validating several conserved miRNA targets (Table 1). IsomiRs of miR156 and miR159 targeted squamosa promoter-binding-like proteins and a MYB transcription factor, respectively (Fig. 4). Homeobox leucine-zipper proteins were targeted by no less than 20 miRNAs of family miR166. Similarly, stress responsive mitochondrial transcription termination factor-like proteins were targeted by three isomiRs of miR9662.

Table 1 Selected miRNA targets validated through degradome sequencing. Conserved miRNA and their targets are italicized
Full size table
Fig. 4
figure 4

Identification of cleaved miRNA targets. The target plots (t-plots) show read abundance aligned to the indicated transcript. The arrows above the miRNA-transcript alignment indicate the miRNA-directed cleavage site

Full size image

Transcripts associated with superoxide dismutase activity were targeted by both miR398–1.b-3p and miR528–1.b-5p. Kinases and transcripts involved in ubiquitin-mediated proteolysis were regulated by several miRNAs from multiple families including miR1118, miR1120, miR1135 and miR1136. F-box family proteins were targeted by isomiRs of miR1133, miR1135, miR1136 and miR9772. Anthranilate synthase was targeted by isomiRs of family miR1136, miR1137 and miR1439 (Table 1).

Prediction of heat stress associated regulatory networks

Degradome validated miRNA targets queried to STRING DB predominantly associated with redox homeostasis, antioxidant activity and ubiquitination (Fig. 5). Cu-Zn Superoxide dismutases targeted by miR528 and miR398 coexpressed with other Fe/Zn Superoxide dismutases. Cellular and redox homeostasis related network were also observed with pentatricopeptide repeat-containing protein (PPR) involving dihydrofolate reductase/thymidylate synthase. A prominent interaction cluster was observed for F-box associated with ubiquitin family of proteins. Query of mitochondrial transcription termination factor (mTERF) resulted in a similar cluster associating both ubiquitin family and phosphatidylinositol 3- and 4-kinase domain carrying proteins. Predicted interacting proteins and their annotation are listed (Additional file 13).

Fig. 5
figure 5

STRING network view of known and predicted interaction partners. Degradome validated miRNA targets were used as query to search orthologous proteins in rice (Oryza sativa). The confidence view screen shows the known and predicted protein-protein interactions where stronger associations are represented by thicker lines. The bold font indicate the queried protein

Full size image

Discussion

Understanding the ability of plants to sense and respond to abiotic stresses such as heat has implications in plant breeding. In recent years, the prominent role of miRNA-guided post-transcriptional gene regulation has been uncovered in plant development and stress responses [reviewed by [31, 56, 57] ]. In this study, we report on the regulatory role of miRNAs in response to heat stress by analyzing an interlaced data set of sRNA and degradome sequences. With stringent annotation guidelines [43], we identified 202 mature miRNAs, of which, 36 were differentially altered with temporal expression profile upon heat stress. Through degradome analysis, we identified 84 miRNA targeting 206 protein coding loci. As reported for other plant species, a high degree of conservation was observed in miR156 regulating squamosa promoter-binding-like protein, miR159 regulating MYB transcription factor, miR166 regulating homeobox leucine-zipper protein and miR398 regulating superoxide dismutase [22, 31]. In addition to the typical heat stress response target transcripts of superoxide dismutases, peroxidases and an array of transcription factors, we also identified miRNA targets for organelle specific transcripts such as pentatricopeptide repeat-containing proteins and mitochondrial transcription termination factor-like proteins through miRNA-guided cleavage; a phenomenon that has not previously been reported in monocots.

Polyploid crops such as wheat have evolutionary advantages in development and adaptation [reviewed by [58] ]. Polyploid genomes better tolerate gene loss or silencing compared to their diploid counterparts because of their innate functional redundancy [59]. Hybridization and whole genome duplication have been demonstrated to result in expansion of miRNA families and enhance stress tolerance [60, 61]. Here, isomiRs of heat stress responsive miR528 and miR9662 were processed from single pre-miRNAs to regulate ROS and mitochondrial transcription termination factor-like proteins, respectively. The absence of the homoeologous coding loci may be from miRNA gene loss on the homoeologues or, alternatively, from selective gain in allohexaploid wheat. Such cases were however the exception rather than the norm considering that 161 miRNAs aligned to more than one pre-miRNA across all sub-genomes with many being highly redundant. For instance, miR1136–235-3p, targeting the 30S ribosomal protein S3, aligned to 1725 pre-miRNA coding loci across all sub-genomes. Previous studies on allohexaploid wheat and their progenitors have shown non-additive expression of miRNAs to regulate growth and adaptability [62]. Allohexaploid wheat and their progenitors are known to carry sub-genome donor- specific miRNA profiles [63].

Plant pre-miRNAs are variable in size and foldback structure compared to their animal counterparts [23]. The presence of one to three nucleotide bulges, mismatches and large loops in the foldback structure alter DCL recognition and processing resulting in new miRNA variants or isomiRs [24, 64]. In our study, all five isomiRs of miR528 and nine isomiRs of miR9662 were processed from pre-miRNA miR528–1 and miR9662–1, respectively. Despite being processed from the same pre-miRNAs, these isomiRs were not produced in equimolar amounts and only a subset of them was differentially expressed after heat stress. Mature miRNA accumulation is dictated by the rate of transcription, processing and their 3′ end modifications, the latter being typically methylated to be protected [65] from 3′ terminal uridination and subsequent degradation [66]. Previous studies have also shown isomiRs accumulation to vary across tissue types and in response to external stimuli such as stresses [67, 68]. Hence, it is imperative to quantify isomiRs because they differ in function and their expression from a single locus is controlled at the post-transcriptional level.

miRNA genes are transcribed independently based on their own promoter activity, a likely explanation for the differential accumulation of miRNAs of the same family. Although, members of miR166 were transcribed by only 16 miRNA coding loci, they accounted for more than 70% of the reads. IsomiRs of miR166 were also predominant in leaf tissue of durum wheat [69] and soybean [70]. Differences were also observed in accumulation of mature miRNA or guide strand and the complementary miRNA originating from the opposite arm, often called miRNA-star(*) (Additional file 9). The most abundant mature miRNAs were considered the dominant and functional products that regulate PTGS, whereas the miRNA* were considered minor products of the duplex processing [71]. However, miRNA* should not be discounted because it also has the ability to regulate gene expression [72]. The abundance of miRNA* miR9662–1.3-5p and miR9662–1.8-5p was higher than their corresponding miRNAs and both were differentially altered in response to heat stress, indicating a possible functional role in wheat (Additional file 11). However, despite their differential expression, we did not observe accumulation of their cleaved mRNA products in the degradome dataset. Nevertheless, three isomiRs processed from the 3p arm of pre-miR9662–1 were identified to target mitochondrial transcription termination factor-like protein encoding genes, lending credence to their functional role. Studies on Arabidopsis have shown that accurate processing and miRNA guide strand selection are regulated by HYL and CPL1 phosphatases [73]. The loading of miRNA to specific AGO is directed by the 5′ terminal nucleotides and mismatches in the duplex structure [74]. Consequently, it is possible that miRNAs processed from 3p and 5p arms can be preferentially selected to be loaded into different AGO complexes to regulate different transcripts in a tissue and stress specific manner.

The temporal expression profile of miRNAs in wheat cv Chinese Spring after heat stress is similar to our previous report on wheat cv Glenlea [37]. However, the relative accumulation of miRNAs varied between the two cultivars. For example, miR398 and miR169 were abundant in cv Glenlea whereas miR1135 and miR166 were the most abundant in cv Chinese Spring. These variations can occur as a consequence of shear genetic and regulatory differences among cultivars as reported in cowpea [39], cotton [38] and rice [75]. sRNA and degradome sequencing confirmed that few miRNA families and their targets were generally conserved across many plant species. miRNA families 156, 159, 166, 393 and 398 targeting squamosa promoter-binding-like proteins, MYB transcription factors, homeobox leucine-zipper proteins, transport inhibitor response proteins and superoxide dismutases, respectively, are highly conserved across dicots and monocots [31, 76]. In addition to conserved miRNA families and their targets, monocot specific miR528 was differentially altered after heat stress to regulate antioxidant activity. Heat stress induced H2O2 has been shown to function as a key signalling molecule to regulate the expression of heat shock proteins [7, 77]. The transient downregulation of miR528 and validation of the guided cleavage of the superoxide dismutase target through degradome sequencing exemplify a miRNA-directed stress response. A recent study has also shown that suppression of miR528 enhances antiviral defense in rice [78]. Triticum-specific miR9772 was validated to target F-box family proteins, implicating the regulating ubiquitin-mediated proteasome degradation. In Arabidopsis, F-box family proteins are regulated by miR393 and miR394 to promote ubiquitination and proteasome degradation [22]. Overexpression of wheat F-box protein TaFBA1 confers oxidative and drought stress tolerance with enhanced antioxidant activity [79]. Anthranilate synthase, mediating the biosynthesis of chorismate to tryptophan, was targeted by miR1137–229-3p. In plants, tryptophan serves as a precursor for the synthesis of phytohormone auxins, phytoallexins, glucosinolates and alkaloids to control an array of developmental processes and to confer biotic and abiotic stress tolerances [80]. Heat stress associated regulator network revealed distinct clusters linking biosynthesis of tryptophan, phytohormones, antioxidant activity and ubiquitination. Although, six transport inhibitor response proteins were identified as targets of miR393, only one was differentially expressed with the expected inverse correlation. Sub-genome specific transcriptional level gene regulation and relative rates of biogenesis and decay of miRNA and their targets are possible [81].

miRNA-guided gene expression was also observed among organelle specific regulatory proteins. For instance, miR2020 targeted chloroplast and mitochondria localized a pentatricopeptide repeat-containing protein (PRP) that facilitates processing, splicing, editing, stability and translation of RNAs [82]. In rice, loss of pentatricopeptide repeat-containing proteins (WSL1) led to enhanced sensitivity to salinity, sugar and abscisic acid [83]. Similarly, mitochondrial transcription termination factor-like protein (mTERF) was targeted by miR9662. PRP and mTEFR carry tandem degenerate helical repeats and direct RNA splicing [84]. Studies on Arabidopsis and maize have shown that mTERF regulates responses to high light [85], heat stress [86] and salinity [87].

Conclusion

The interlaced data set generated in this study identified and validated heat stress regulated miRNAs and their target genes associated with thermotolerance. We showed that miRNAs were highly altered immediately after heat stress to regulate thermotolerance but expression levels largely returned to control levels during the 1–4 day period following the end of the stress. Degradome sequencing revealed that the target genes of miRNA156, miR159, miR166 and miR398 are conserved among wheat, other cereal crops and dicot plants. Chloroplast and mitochondria localized pentatricopeptide repeat-containing and mitochondrial transcription termination factor-like proteins were shown to be regulated through miRNA-guided cleavage, a phenomenon that has not been reported to date in monocots. A thorough understanding of miRNA-mediated epigenetic changes constitutes an additional source of variation that can be capitalized upon for developing new breeding strategies targeting expression rather than structural changes.

Availability of data and materials

The smallRNA and degradome (PARE) data were submitted to the NCBI Gene Expression Omnibus under accession number GSE113358.

Abbreviations

DAT:

days after treatment

miRNA:

microRNA

nt:

nucleotide

PARE:

Parallel Analysis of RNA Ends

RPM:

reads per million

References

  1. Lesk C, Rowhani P, Ramankutty N. Influence of extreme weather disasters on global crop production. Nature. 2016;529:84–7.

    Article  CAS  PubMed  Google Scholar 

  2. Sharma SK, De Los Rios P, Christen P, Lustig A, Goloubinoff P. The kinetic parameters and energy cost of the Hsp70 chaperone as a polypeptide unfoldase. Nat Chem Biol. 2010;6:914–20.

    Article  CAS  PubMed  Google Scholar 

  3. Kaushal N, Bhandari K, Siddique KHM, Nayyar H. Food crops face rising temperatures: an overview of responses, adaptive mechanisms, and approaches to improve heat tolerance. Cogent Food & Agriculture. 2016;2:1134380.

    Article  Google Scholar 

  4. Saidi Y, Finka A, Muriset M, Bromberg Z, Weiss YG, Maathuis FJ, Goloubinoff P. The heat shock response in moss plants is regulated by specific calcium-permeable channels in the plasma membrane. Plant Cell. 2009;21:2829–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Konigshofer H, Tromballa HW, Loppert HG. Early events in signalling high-temperature stress in tobacco BY2 cells involve alterations in membrane fluidity and enhanced hydrogen peroxide production. Plant Cell Environ. 2008;31:1771–80.

    Article  CAS  PubMed  Google Scholar 

  6. Sangwan V, Orvar BL, Beyerly J, Hirt H, Dhindsa RS. Opposite changes in membrane fluidity mimic cold and heat stress activation of distinct plant MAP kinase pathways. Plant J. 2002;31:629–38.

    Article  CAS  PubMed  Google Scholar 

  7. Volkov RA, Panchuk II, Mullineaux PM, Schöffl F. Heat stress-induced H2O2 is required for effective expression of heat shock genes in Arabidopsis. Plant Mol Biol. 2006;61:733–46.

    Article  CAS  PubMed  Google Scholar 

  8. Evrard A, Kumar M, Lecourieux D, Lucks J, von Koskull-Döring P, Hirt H. Regulation of the heat stress response in Arabidopsis by MPK6-targeted phosphorylation of the heat stress factor HsfA2. PeerJ. 2013;1:e59.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. Li S, Zhou X, Chen L, Huang W, Yu D. Functional characterization of Arabidopsis thaliana WRKY39 in heat stress. Mol Cell. 2010;29:475–83.

    Article  CAS  Google Scholar 

  10. Xie Z, Allen E, Fahlgren N, Calamar A, Givan SA, Carrington JC. Expression of Arabidopsis miRNA genes. Plant Physiol. 2005;138:2145–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Cui J, You C, Chen X. The evolution of microRNAs in plants. Curr Opin Plant Biol. 2017;35:61–7.

    Article  CAS  PubMed  Google Scholar 

  12. Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genetics. 2014;15:394–408.

    Article  CAS  PubMed  Google Scholar 

  13. Palatnik JF, Wollmann H, Schommer C, Schwab R, Boisbouvier J, Rodriguez R, Warthmann N, Allen E, Dezulian T, Huson D, et al. Sequence and expression differences underlie functional specialization of Arabidopsis microRNAs miR159 and miR319. Dev Cell. 2007;13:115–25.

    Article  CAS  PubMed  Google Scholar 

  14. Kruszka K, Pacak A, Swida-Barteczka A, Nuc P, Alaba S, Wroblewska Z, Karlowski W, Jarmolowski A, Szweykowska-Kulinska Z. Transcriptionally and post-transcriptionally regulated microRNAs in heat stress response in barley. J Exp Bot. 2014;65:6123–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Mallory AC, Bartel DP, Bartel B. MicroRNA-directed regulation of Arabidopsis AUXIN RESPONSE FACTOR17 is essential for proper development and modulates expression of early auxin response genes. Plant Cell. 2005;17:1360–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Wu MF, Tian Q, Reed JW. Arabidopsis microRNA167 controls patterns of ARF6 and ARF8 expression, and regulates both female and male reproduction. Development. 2006;133:4211–8.

    Article  CAS  PubMed  Google Scholar 

  17. Ma Z, Hu X, Cai W, Huang W, Zhou X, Luo Q, Yang H, Wang J, Huang J. Arabidopsis miR171-targeted scarecrow-like proteins bind to GT cis-elements and mediate gibberellin-regulated chlorophyll biosynthesis under light conditions. PLoS Genet. 2014;10:e1004519.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. 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 

  19. Sorin C, Declerck M, Christ A, Blein T, Ma L, Lelandais-Brière C, Njo MF, Beeckman T, Crespi M, Hartmann C. A miR169 isoform regulates specific NF-YA targets and root architecture in Arabidopsis. New Phytol. 2014;202:1197–211.

    Article  CAS  PubMed  Google Scholar 

  20. Zhao M, Ding H, Zhu JK, Zhang F, Li WX. Involvement of miR169 in the nitrogen-starvation responses in Arabidopsis. New Phytol. 2011;190:906–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Guan Q, Lu X, Zeng H, Zhang Y, Zhu J. Heat stress induction of miR398 triggers a regulatory loop that is critical for thermotolerance in Arabidopsis. Plant J. 2013;74:840–51.

    Article  CAS  PubMed  Google Scholar 

  22. Jones-Rhoades MW, Bartel DP. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell. 2004;14:787–99.

    Article  CAS  PubMed  Google Scholar 

  23. Cuperus JT, Fahlgren N, Carrington JC. Evolution and functional diversification of MIRNA genes. Plant Cell. 2011;23:431–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Chavez Montes RA. de Fatima Rosas-Cardenas F, De Paoli E, Accerbi M, Rymarquis LA, Mahalingam G, Marsch-Martinez N, Meyers BC, Green PJ, de Folter S. sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs. Nat Commun. 2014;5:3722.

    Article  PubMed  CAS  Google Scholar 

  25. You C, Cui J, Wang H, Qi X, Kuo LY, Ma H, Gao L, Mo B, Chen X. Conservation and divergence of small RNA pathways and microRNAs in land plants. Genome Biol. 2017;18:158.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Dai X, Zhao PX. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011;39:W155–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Pinzon N, Li B, Martinez L, Sergeeva A, Presumey J, Apparailly F, Seitz H. MicroRNA target prediction programs predict many false positives. Genome Res. 2017;27:234–245.

    Article  PubMed  CAS  Google Scholar 

  28. Seitz H. Issues in current microRNA target identification methods. RNA Biol. 2017;14:831–834.  

    Article  PubMed  PubMed Central  Google Scholar 

  29. German MA, Luo S, Schroth G, Meyers BC, Green PJ. Construction of parallel analysis of RNA ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nature Prot. 2009;4:356–62.

    Article  CAS  Google Scholar 

  30. Smith LM, Burbano HA, Wang X, Fitz J, Wang G, Ural-Blimke Y, Weigel D. Rapid divergence and high diversity of miRNAs and miRNA targets in the Camelineae. Plant J. 2015;81:597–610.

    Article  CAS  PubMed  Google Scholar 

  31. Liu H, Able AJ, Able JA. SMARTER de-stressed cereal breeding. Trends Plant Sci. 2016;21:909–25.

    Article  CAS  PubMed  Google Scholar 

  32. Stief A, Altmann S, Hoffmann K, Pant BD, Scheible WR, Baurle I. Arabidopsis miR156 regulates tolerance to recurring environmental stress through SPL transcription factors. Plant Cell. 2014;26:1792–807.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wang Y, Sun F, Cao H, Peng H, Ni Z, Sun Q, Yao Y. TamiR159 directed wheat TaGAMYB cleavage and its involvement in anther development and heat response. PLoS One. 2012;7:e48445.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Kumar RR, Pathak H, Sharma SK, Kala YK, Nirjal MK, Singh GP, Goswami S, Rai RD. Novel and conserved heat-responsive microRNAs in wheat (Triticum aestivum L.). Funct Integr Genomics. 2015;15:323–48.

    Article  CAS  PubMed  Google Scholar 

  35. Xin M, Wang Y, Yao Y, Xie C, Peng H, Ni Z, Sun Q. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 2010;10:123.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Qin D, Wu H, Peng H, Yao Y, Ni Z, Li Z, Zhou C, Sun Q. Heat stress-responsive transcriptome analysis in heat susceptible and tolerant wheat (Triticum aestivum L.) by using wheat genome Array. BMC Genomics. 2008;9:432.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Ragupathy R, Ravichandran S, Mahdi MSR, Huang D, Reimer E, Domaratzki M, Cloutier S. Deep sequencing of wheat sRNA transcriptome reveals distinct temporal expression pattern of miRNAs in response to heat, light and UV. Sci Rep. 2016;6:39373.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Yin Z, Li Y, Yu J, Liu Y, Li C, Han X, Shen F. Difference in miRNA expression profiles between two cotton cultivars with distinct salt sensitivity. Mol Biol Rep. 2012;39:4961–70.

    Article  CAS  PubMed  Google Scholar 

  39. Barrera-Figueroa BE, Gao L, Diop NN, Wu Z, Ehlers JD, Roberts PA, Close TJ, Zhu J-K, Liu R. Identification and comparative analysis of drought-associated microRNAs in two cowpea genotypes. BMC Plant Biol. 2011;11:127.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Swetha C, Basu D, Pachamuthu K, Tirumalai V, Nair A, Prasad M, Shivaprasad PV. Major domestication-related phenotypes in Indica Rice are due to loss of miRNA-mediated laccase silencing. Plant Cell. 2018;30:2649–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. International Wheat Genome Sequencing Consortium (IWGSC 2018)Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science. 2018;361:eaar7191.

  42. Lu C, Souret F: High-throughput approaches for miRNA expression analysis. In: Methods Mol Biol. Vol. 592, 2009/10/06 edn; 2010: 107–125.

  43. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, Cao X, Carrington JC, Chen X, Green PJ, et al. Criteria for annotation of plant microRNAs. Plant Cell. 2008;20:3186–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42:D68–73.

    Article  CAS  PubMed  Google Scholar 

  45. 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  CAS  Google Scholar 

  46. Lorenz R, Bernhart SH, Honer Zu Siederdissen C, Tafer H, Flamm C, Stadler PF, Hofacker IL. ViennaRNA Package 2.0. Algorithms for Molecular Biology. 2011;6:26.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Sun F, Guo G, Du J, Guo W, Peng H, Ni Z, Sun Q, Yao Y. Whole-genome discovery of miRNAs and their targets in wheat (Triticum aestivum L.). BMC Plant Biol. 2014;14:142.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  50. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57:289–300.

    Google Scholar 

  51. Metsalu T, Vilo J. ClustVis: a web tool for visualizing clustering of multivariate data using principal component analysis and heatmap. Nucleic Acids Res. 2015;43:W566–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zhai J, Arikit S, Simon SA, Kingham BF, Meyers BC. Rapid construction of parallel analysis of RNA end (PARE) libraries for Illumina sequencing. Methods. 2014;67:84–90.

    Article  CAS  PubMed  Google Scholar 

  53. Davis MP, van Dongen S, Abreu-Goodger C, Bartonicek N, Enright AJ. Kraken: a set of tools for quality control and analysis of high-throughput sequence data. Methods. 2013;63:41–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Kakrana A, Hammond R, Patel P, Nakano M, Meyers BC. sPARTA: a parallelized pipeline for integrated analysis of plant miRNA and cleaved mRNA data sets, including new miRNA target-identification software. Nucleic Acids Res 2014;42:e139-e139.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45:D362–d368.

    Article  CAS  PubMed  Google Scholar 

  56. Li C, Zhang B. MicroRNAs in control of plant development. J Cell Physiol. 2016;231:303–13.

    Article  CAS  PubMed  Google Scholar 

  57. 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 

  58. Wendel JF, Jackson SA, Meyers BC, Wing RA. Evolution of plant genome architecture. Genome Biol. 2016;17:37.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Song Q, Zhang T, Stelly DM, Chen ZJ. Epigenomic and functional analyses reveal roles of epialleles in the loss of photoperiod sensitivity during domestication of allotetraploid cottons. Genome Biol. 2017;18:99.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  60. Liu B, Sun G. microRNAs contribute to enhanced salt adaptation of the autopolyploid Hordeum bulbosum compared with its diploid ancestor. Plant J. 2017;91:57–69.

    Article  CAS  PubMed  Google Scholar 

  61. Kenan-Eichler M, Leshkowitz D, Tal L, Noor E, Melamed-Bessudo C, Feldman M, Levy AA. Wheat hybridization and polyploidization results in deregulation of small RNAs. Genetics. 2011;188:263–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Li A, Liu D, Wu J, Zhao X, Hao M, Geng S, Yan J, Jiang X, Zhang L, Wu J, et al. mRNA and small RNA transcriptomes reveal insights into dynamic homoeolog regulation of allopolyploid heterosis in nascent hexaploid wheat. Plant Cell. 2014;26:1878–900.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Alptekin B, Budak H. Wheat miRNA ancestors: evident by transcriptome analysis of a, B, and D genome donors. Funct Integr Genomics. 2017;17:171–87.

    Article  CAS  PubMed  Google Scholar 

  64. Bologna NG, Schapire AL, Zhai J, Chorostecki U, Boisbouvier J, Meyers BC, Palatnik JF. Multiple RNA recognition patterns during microRNA biogenesis in plants. Genome Res. 2013;23:1675–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Yu B, Yang Z, Li J, Minakhina S, Yang M, Padgett RW, Steward R, Chen X. Methylation as a crucial step in plant microRNA biogenesis. Science. 2005;307:932–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Ren G, Chen X, Yu B. Uridylation of miRNAs by HEN1 SUPPRESSOR1 in Arabidopsis. Curr Biol. 2012;22:695–700.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Newman MA, Mani V, Hammond SM. Deep sequencing of microRNA precursors reveals extensive 3′ end modification. RNA. 2011;17:1795–803.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Wyman SK, Knouf EC, Parkin RK, Fritz BR, Lin DW, Dennis LM, Krouse MA, Webster PJ, Tewari M. Post-transcriptional generation of miRNA variants by multiple nucleotidyl transferases contributes to miRNA transcriptome complexity. Genome Res. 2011;21:1450–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Giusti L, Mica E, Bertolini E, De Leonardis AM, Faccioli P, Cattivelli L, Crosatti C. microRNAs differentially modulated in response to heat and drought stress in durum wheat cultivars with contrasting water use efficiency. Funct Integr Genomics. 2017;17:293–309.

    Article  CAS  PubMed  Google Scholar 

  70. Arikit S, Xia R, Kakrana A, Huang K, Zhai J, Yan Z, Valdes-Lopez O, Prince S, Musket TA, Nguyen HT, et al. An atlas of soybean small RNAs identifies phased siRNAs from hundreds of coding genes. Plant Cell. 2014;26:4584–601.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Hutvagner G. Small RNA asymmetry in RNAi: function in RISC assembly and gene regulation. FEBS Lett. 2005;579:5850–7.

    Article  CAS  PubMed  Google Scholar 

  72. Yang JS, Phillips MD, Betel D, Mu P, Ventura A, Siepel AC, Chen KC, Lai EC. Widespread regulatory activity of vertebrate microRNA* species. RNA. 2011;17:312–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Manavella Pablo A, Hagmann J, Ott F, Laubinger S, Franz M, Macek B, Weigel D. Fast-forward genetics identifies plant CPL phosphatases as regulators of miRNA processing factor HYL1. Cell. 2012;151:859–70.

    Article  CAS  PubMed  Google Scholar 

  74. Mi S, Cai T, Hu Y, Chen Y, Hodges E, Ni F, Wu L, Li S, Zhou H, Long C, et al. Sorting of small RNAs into Arabidopsis argonaute complexes is directed by the 5′ terminal nucleotide. Cell. 2008;133:116–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Mondal TK, Ganie SA. Identification and characterization of salt responsive miRNA-SSR markers in rice (Oryza sativa). Gene. 2014;535:204–9.

    Article  CAS  PubMed  Google Scholar 

  76. Zhang B, Wang Q. MicroRNA-based biotechnology for plant improvement. J Cell Physiol. 2015;230:1–15.

    Article  PubMed  CAS  Google Scholar 

  77. Choudhury FK, Rivero RM, Blumwald E, Mittler R. Reactive oxygen species, abiotic stress and stress combination. Plant J. 2017;90:856–67.

    Article  CAS  PubMed  Google Scholar 

  78. Wu J, Yang R, Yang Z, Yao S, Zhao S, Wang Y, Li P, Song X, Jin L, Zhou T, et al. ROS accumulation and antiviral defence control by microRNA528 in rice. Nat Plants. 2017;3:16203.

    Article  CAS  PubMed  Google Scholar 

  79. Zhou S-M, Kong X-Z, Kang H-H, Sun X-D, Wang W. The involvement of wheat F-box protein gene TaFBA1 in the oxidative stress tolerance of plants. PLoS One. 2015;10:e0122117.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  80. Radwanski ER, Last RL. Tryptophan biosynthesis and metabolism: biochemical and molecular genetics. Plant Cell. 1995;7:921–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Mukherji S, Ebert MS, Zheng GX, Tsang JS, Sharp PA, van Oudenaarden A. MicroRNAs can generate thresholds in target gene expression. Nat Genet. 2011;43:854–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Manna S. An overview of pentatricopeptide repeat proteins and their applications. Biochimie. 2015;113:93–9.

    Article  CAS  PubMed  Google Scholar 

  83. Tan J, Tan Z, Wu F, Sheng P, Heng Y, Wang X, Ren Y, Wang J, Guo X, Zhang X, et al. A novel chloroplast-localized pentatricopeptide repeat protein involved in splicing affects chloroplast development and abiotic stress response in rice. Mol Plant. 2014;7:1329–49.

    Article  CAS  PubMed  Google Scholar 

  84. Hammani K, Barkan A. An mTERF domain protein functions in group II intron splicing in maize chloroplasts. Nucleic Acids Res. 2014;42:5033–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Meskauskiene R, Wursch M, Laloi C, Vidi PA, Coll NS, Kessler F, Baruah A, Kim C, Apel K. A mutation in the Arabidopsis mTERF-related plastid protein SOLDAT10 activates retrograde signaling and suppresses 1O2-induced cell death. Plant J. 2009;60:399–410.

    Article  CAS  PubMed  Google Scholar 

  86. Kim M, Lee U, Small I, des Francs-Small CC, Vierling E. Mutations in an Arabidopsis mitochondrial transcription termination factor-related protein enhance thermotolerance in the absence of the major molecular chaperone HSP101. Plant Cell. 2012;24:3349–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Robles P, Micol JL, Quesada V. Mutations in the plant-conserved MTERF9 alter chloroplast gene expression, development and tolerance to abiotic stress in Arabidopsis thaliana. Physiol Plant. 2015;154:297–313.

    Article  CAS  PubMed  Google Scholar 

  88. Cui L-G, Shan J-X, Shi M, Gao J-P, Lin H-X. The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants. Plant J. 2014;80:1108–17.

    Article  CAS  PubMed  Google Scholar 

  89. Swida-Barteczka A, Pacak A, Jarmolowski A, Kruszka K, Nuc P, Szweykowska-Kulinska Z, Alaba S, Karlowski W, Wroblewska Z. Transcriptionally and post-transcriptionally regulated microRNAs in heat stress response in barley. J Exp Bot. 2014;65:6123–35.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  90. Kinoshita N, Wang H, Kasahara H, Liu J, Macpherson C, Machida Y, Kamiya Y, Hannah MA, Chua NH. IAA-ala Resistant3, an evolutionarily conserved target of miR167, mediates Arabidopsis root architecture changes during high osmotic stress. Plant Cell. 2012;24:3590–602.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Liu Q, Zhang Y-C, Wang C-Y, Luo Y-C, Huang Q-J, Chen S-Y, Zhou H, Qu L-H, Chen Y-Q. Expression analysis of phytohormone-regulated microRNAs in rice, implying their regulation roles in plant hormone signaling. FEBS Lett. 2009;583:723–8.

    Article  CAS  PubMed  Google Scholar 

  92. Zhou L, Liu Y, Liu Z, Kong D, Duan M, Luo L. Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa. J Exp Bot. 2010;61:4157–68.

    Article  CAS  PubMed  Google Scholar 

  93. Ganie SA, Dey N, Mondal TK. Promoter methylation regulates the abundance of Osa-miR393a in contrasting rice genotypes under salinity stress. Funct Integr Genomics. 2016;16:1–11.

    Article  CAS  PubMed  Google Scholar 

  94. Jagadeeswaran G, Saini A, Sunkar R. Biotic and abiotic stress down-regulate miR398 expression in Arabidopsis. Planta. 2009;229:1009–14.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We are grateful to Drs. Pamela Green and Vinay Nagarajan of University of Delaware for their valuable suggestions on optimizing small RNA library construction. We thank Dr. Sateesh Kagale of National Research Council Canada in Saskatoon for his insights into data analysis. The contribution from many members of the Cloutier’s laboratory during tissue sampling is acknowledged.

Funding

Funding from the Canadian Crop Genomics Initiative (project J-000072) and the Canadian Wheat Alliance (project J-001584) was provided to Cloutier through competitive grants. The funding bodies were not involved in the design of the study and collection, analysis, and interpretation of data, nor in writing the manuscript.

Author information

Authors and Affiliations

  1. Agriculture and Agri-Food Canada, Ottawa Research and Development Centre, 960 Carling Avenue, Ottawa, Ontario, K1A 0C6, Canada

    Sridhar Ravichandran, Tara Edwards & Sylvie Cloutier

  2. Plant Science Department, University of Manitoba, Winnipeg, Manitoba, Canada

    Raja Ragupathy

  3. Department of Computer Science, University of Manitoba, Winnipeg, Manitoba, Canada

    Michael Domaratzki

  4. Present address: Agriculture and Agri-Food Canada, Lethbridge Research and Development Centre, Lethbridge, Alberta, Canada

    Raja Ragupathy

Authors
  1. Sridhar Ravichandran
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Raja Ragupathy
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Tara Edwards
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Michael Domaratzki
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Sylvie Cloutier
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

SC and RR conceived and performed the experiment and devised the strategy for data analysis. RR, MD and SR developed various modules for the bioinformatics tool chain and its implementation on different platforms. TE performed miRNA library construction and sequencing. RR and SR performed the data analyses. SC, RR and SR co-wrote the manuscript. All authors reviewed and contributed to the manuscript.

Corresponding author

Correspondence to Sylvie Cloutier.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interest.

Additional information

Publisher’s Note

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

Additional files

Additional file 1:

Table S1 Sequencing and processing summary of the 24 small RNA libraries. (XLSX 19 kb)

Additional file 2:

Figure S1. Overview of small RNA profile. a Distribution of raw reads by treatment and sampling time point. b Distribution of distinct tags by treatment and sampling time point. The median is shown by a horizontal black line. Boxes represent the 25th and 75th percentiles while error bars indicate the 10th and 90th percentiles. c Size distribution of raw reads after size selection. d Size distribution of distinct tags prior to removal of other non-coding RNA and chloroplast sequences. (TIF 895 kb)

Additional file 3:

Figure S2. miRNAs families. (a) Number of isomiRs in the 37 miRNA families identified in heat stress and control wheat leaves. (b) Relative abundance of size classes within the same conserved 37 miRNA families. (TIF 473 kb)

Additional file 4:

Table S2. Conserved miRNA families, precursor sequence, structure, minimum free energy and star sequence information. (XLSX 54 kb)

Additional file :5

Figure S3. Size distribution of miRNA families represented as percentage of raw reads and distinct tags. (TIF 151 kb)

Additional file 6:

Figure S4. Relative abundance of miRNAs families expressed in heat stress and control wheat leaves. (a) Normalized raw read percentages of miRNA families. (b) Relative abundance of the normalized read counts based on the levels of expression in reads per million (RPM). (TIF 490 kb)

Additional file 7:

Figure S5. Box plots of the relative abundance of mature miRNAs belonging to family miR166. The median is shown by a horizontal black line. Boxes represent the 25th and 75th percentiles while error bars indicate the 10th and 90th percentiles. (TIF 271 kb)

Additional file 8:

Table S3. Normalized read counts of annotated miRNAs for each of 24 small RNA libraries corresponding to two treatments (control and heat), four replicates and three sampling time points (0, 1 and 4 days after treatment). (XLSX 71 kb)

Additional file 9:

Figure S6. Abundance of miRNAs cleaved from the 5p and 3p ends of pre-miRNA9662. (a) miR9662–1.3-5p/miR9662–1.4-3p and (b) miR9662–1.8-5p/miR9662–1.7-3p. (TIF 218 kb)

Additional file 10:

Figure S7. MDS plot displaying the relative similarities among the biological replicates. (TIF 148 kb)

Additional file 11:

Table S4. List of miRNA differentially expressed after heat stress. (XLSX 17 kb)

Additional file 12:

Table S5. List of miRNA targets identified by degradome sequencing. (XLSX 598 kb)

Additional file 13:

Table S6. List of predicted interaction partners for degradome validated miRNA targets identified in rice (Oryza sativa) using STRING DB. (XLSX 15 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ravichandran, S., Ragupathy, R., Edwards, T. et al. MicroRNA-guided regulation of heat stress response in wheat. BMC Genomics 20, 488 (2019). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-019-5799-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-019-5799-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us