High-throughput sequencing reveals microRNAs and their targets in response to drought stress in wheat (Triticum aestivum L.)

Abstract

Drought is a major factor that influences the productivity of crops including wheat. However, drought tolerance is a complex trait involving many genes and cellular pathways. To help understand the complex regulatory mechanisms governing drought tolerance in wheat, we performed small RNA sequencing analysis in a highly drought-tolerant wheat variety (Triticum aestivum L. cv. XF 20). Experimental plants were cultivated in the presence of PEG6000 to induce drought stress, while control plants were grown in the absence of PEG6000. Sequencing results confirmed the expression of 199 previously reported miRNAs. We also identified 32 novel miRNAs in this species. The expression of 37 miRNAs was significantly affected by PEG6000-induced drought stress (27 upregulated, 10 downregulated), including 5 novel miRNAs. Degradome sequencing showed that 222 target gene mRNAs were cleaved by these miRNAs. The target genes identified by psRNATarget and degradome were all enriched in gene ontology categories. There were significant differences in the number of target genes regulated by differentially expressed miRNAs, especially in signal transducer, response to stress and antioxidant pathways. MiRNA targets that are known to be related to drought and/or stress response in plants included 6-phosphogluconate dehydrogenase, epoxide hydrolase, sucrose synthase and phytochrome and flowering time regulatory protein. Our results shed light on the molecular mechanisms of drought tolerance in wheat, and identify drought-responsive miRNAs and their targets in a drought-resistant variety. These results could be exploited to improve wheat and related plants, which is expected to be vital as agriculture addresses the consequences of climate change.

Keywords:

• XF 20
• degradome
• drought tolerance
• miRNAs
• RNA sequencing
• wheat

Introduction

Plants have developed remarkable capabilities to adapt to extreme environments, such as high temperature, drought and salt stress [1]. Drought is a common environmental stress, and dry weather is becoming more frequent with intensification of the greenhouse effect and climate change [2, 3]. This is particularly true in the agricultural areas of northwestern China, which are frequently exposed to drought conditions and dry, hot winds. Drought limits the growth, development and ultimately, the yields of crop plants worldwide [4]. Plants respond to stresses by regulating the expression of specific genes to avoid or minimize cellular damage, thus adapting to stress conditions [5].

Recently, the discovery of microRNAs (miRNAs) revealed a new and additional layer of post-transcriptional gene regulation. MiRNAs are 21 to 24 nucleotides (nt) long, and are noncoding RNAs that can negatively regulate gene expression [6–8]. Plant miRNAs are highly complementary to target messenger RNAs (mRNAs), and miRNA sequences are generally highly conserved between related organisms. Increasing evidence indicates that miRNAs participate in a diverse set of processes such as plant growth, signal transduction, apoptosis, and biotic and abiotic stress responses [9–13]. Drought-stress-related miRNAs have been identified by small RNA sequencing technology in rice [14], wheat [15], sugarcane [16], Medicago truncatula [17], and potato [18]. Upon drought stress, approximately 14 miRNAs were found to be significantly differentially expressed in bread wheat [2], while expression patterns and response varied between root and leaf [19]. One miRNA may regulate many target genes [20, 21], and the biological functions of miRNAs can be deduced by identification and characterization of their target genes. To this end, characterization of the degradome typically includes deep sequencing, bioinformatic analysis, and 5′-rapid amplification of cDNA ends. Degradome analysis has been used for global identification of miRNA–target RNA pairs in numerous plant species, including Arabidopsis [22], rice [23], wheat [24], soybean [25], maize [26] and grapevine [27].

Hexaploid wheat (Triticum aestivum L. AABBDD) is a major food crop, nurturing more than one-third of the world’s population and providing nearly 55% of carbohydrates grown and consumed [28, 29]. Environmental constraints, such as high temperature (or heat stress), drought, and the combination of the two, severely and antagonistically limit wheat yield and quality [5, 7–9]. It is predicted that a global temperature increase of 2 °C above current averages could lead to wheat yield reductions of up to 50% via perturbations in physiological and biochemical processes [30]. Moreover, drought is reported to adversely affect >50% of current global wheat cultivation land, causing considerable yield loss through the inhibition of photosynthesis [6, 31].

In this work, we investigated the differential expression of miRNAs in the wheat cultivar XF 20 in response to drought stress, and identified targets of these miRNAs. We also identified 32 novel miRNAs in wheat. As this variety of wheat is drought-resistant, knowledge of its regulatory and biochemical mechanisms in response to drought may be used in programmes to improve the resistance of other varieties of wheat.

Materials and methods

Plant growth and treatment

The XF 20 wheat cultivar has been used in the identification of factors influencing drought resistance in a national regional trial and was first released by Qingyang Academy of Agricultural Sciences in 1984 [32]. This variety has grade 2 drought resistance, and it is mainly resistant to heat and drought, hot wind (National Infrastructure of Plant Germplasm Resources, Beijing). Supplementary information Table S1 shows the indices of different wheat varieties under drought stress—the drought-resistant phenotype of Xifeng 20 is obvious.

For our experiments, seeds of wheat cultivar XF 20 were surface-sterilized in 1% sodium hypochlorite for 20 min, rinsed in distilled water six times, and soaked in the dark at room temperature overnight. All seedlings were incubated in a growth chamber [24 °C/18 °C (day/night), 14 h/10 h (light/dark), 60% humidity]. The seedlings were irrigated with half-strength Hoagland solution [33] once every other day, and after 21 days they were randomly divided into two groups. For drought stress (treatment group), seedlings were treated with 20% (w/v) PEG-6000 for 6 h (roots were totally covered). The untreated group (control) was incubated as above without exposure to PEG-6000. All experiments were performed in parallel. Roots were collected 6 h after the start of treatment. For the two groups (control and treatment), the roots of six randomly chosen seedlings were pooled to form a biological replicate, snap frozen in liquid nitrogen and stored at −80 °C until further use.

Small RNA library construction and sequencing

Small RNAs were extracted from plant samples from the two groups (control and drought) using TRIzol reagent (Invitrogen, Carlsbad, CA), according to the manufacturer’s instructions. Small RNAs were ligated sequentially to 59- and 39-RNA/DNA chimeric oligonucleotide adaptors. The resulting ligation products were gel-purified by 15% denaturing polyacrylamide gel electrophoresis, and reverse-transcribed to produce cDNAs. cDNAs were sequenced using a Genome Analyzer IIx System (Illumina, USA), according to the manufacturer’s instructions.

Identification of conserved and novel miRNAs

Raw data were first processed by filtering out low-quality reads, trimming the adaptors and removing other noisy reads, to obtain clean reads. Clean reads were aligned to the Rfam database to remove other noncoding RNAs, including rRNAs, tRNAs and snRNAs. The remaining reads were mapped to assembled transcriptome sequences from Brachypodium distachyon [34]. Mapped reads were retrieved and aligned to plant miRNA sequences from miRBase V21 [5] using Bowtie [35], which identified and annotated conserved miRNA genes. To identify novel miRNA genes, first, miRDeep-P retrieved the flanking assembled transcriptome sequences of mapped reads, which were identified as candidate precursors of miRNAs [36]. Second, based on the precursors, novel miRNAs were identified using MIREAP, and their secondary structures were verified by the software RNAfold, as previously reported [37]. In addition, plant miRNA criteria for Arabidopsis [38] and rice [39] were also considered for novel wheat miRNA identification.

All miRNA abundances were evaluated and normalized using the tags per million reads (TPM) method, based on BLAST mapping results. The TPM values were calculated as: TPM = (number of mapped miRNA reads × 10⁶)/number of clean sample reads. The normalised expression was adjusted to 0.01 when miRNA expression (TPM) was zero to avoid negating subsequent calculation of fold-change. The miRNA expression fold-changes between drought-stress and control group plants were computed, and chi-square tests were performed to determine the significance of miRNA expression for each comparison using R software. The miRNAs with fold change (TPM ratios) ≥2 or ≤0.5, and p value ≤0.05, were deemed to be differentially expressed in response to drought stress, as described by Xie et al. [40].

Degradome library construction and target identification

To investigate the potential targets of miRNAs, a degradome library was constructed from drought stress samples, as previously described [41]. In brief, poly(A)-enriched RNAs were isolated, and ligated to an RNA oligonucleotide adaptor containing a 39-MmeI recognition site; the ligated products were used to synthesize first-strand cDNA. A short PCR program (five cycles) was used to amplify the cDNA, and the product was ligated to a double-stranded DNA adaptor before being subjected to gel purification again for PCR amplification. The final cDNA library was purified and sequenced on an Illumina Genome Analyzer IIx by BGI-Shenzhen Co. Ltd.

Adaptor sequences and low-quality sequencing reads were removed from the raw reads, and the clean reads were used to identify potentially cleaved targets based on B. distachyon assembled transcriptome sequences (described below) by the CleaveLand4 pipeline [42]. The psRNATarget tool was also used to predict miRNA targets using a set of default parameters [43]. MiRNA target genes were also predicted from B. distachyon assembled transcriptome sequences. Prediction targets were used to cross-check the degradome sequencing results.

Data availability

The Triticum aestivum (cv. XF 20) small RNA sequences, degradome sequences, and transcriptome sequences have all been deposited in the NCBI SRA database with accession number SRP 166912.

Results and discussion

High-throughput sequencing of small RNA libraries

To identify novel miRNAs expressed in wheat (cv. XF 20), we successfully produced a total of 11,723,790 and 12,034,102 clean reads from control and drought-treated plant cDNA libraries, respectively (summarized in Table 1). The length of the small RNA reads from the two libraries ranged from 18 to 30 nt. The clean reads were mapped to the Triticum aestivum genome assembly; about 79.9% of control and 84.4% of drought-treated plant reads were mapped. To identify conserved miRNAs, the mappable reads were aligned to known plant miRNAs in the miRbase database, using Bowtie with no more than one mismatch. In total, 814,254 and 1,287,494 reads from the control and drought libraries, respectively, were identified as homologous to known miRNAs. The flanking sequences of the remaining mappable reads were retrieved and analyzed by MIREAP using plant default parameters. In total, 32 novel miRNAs were identified in XF 20, most being 21 and 24 nt long (Supplementary information Table S2).

Table 1. Overview of small RNA sequences for two wheat (Triticum aestivum L. cv. XF 20) libraries.

Data type	Control		Drought
	Total reads	Unique reads	Total reads	Unique reads
Clean reads	11723790	1508601	12034102	1357647
Genome	9365098	904379	10151515	918905
Rfam	3935209	166060	3590058	114691
Known miRNA	814254	18198	1287494	27292
Novel miRNA	12972	1025	12512	1124

Note: This table summarizes the number of total and unique reads for both the control and drought-stressed plant cDNA libraries. Clean reads, reads mapped to the genome of Brachypodium distachyon, reads mapped to Rfam, known miRNA and novel miRNA read numbers are listed.

To identify miRNAs involved in the response of XF 20 plants to drought stress, we evaluated and analyzed the miRNA expression in the roots of both the control and drought plant groups. Compared with the control group, miRNAs with a log2 fold-change >1, or < −1, combined with a p value <0.05, were identified as differentially expressed miRNAs. We identified 37 significantly differentially expressed miRNAs in response to drought stress (Table 2). Most (27/37) were downregulated by drought stress; ten ones (10/37) were upregulated. Five of the 37 differentially expressed miRNAs were novel miRNAs, two (upregulated) and three (down-regulated) predicted novel miRNAs were possibly related to drought resistance in wheat.

Table 2. Differentially expressed miRNA genes in wheat (Triticum aestivum L. cv. XF 20) in response to drought stress.

miRNA name	Control	Drought	Fold	p value
miR165a-5p	3551.50	1688.12	−1.07	**
miR419	167.01	1.00	−7.38	**
miR482a-5p	1192.36	1.00	−10.22	**
miR530b	50.15	1.00	−5.65	**
miR533e	158.40	1.00	−7.31	**
miR858a	1.00	21.52	4.43	**
miR1023a-3p	47.08	20.28	−1.22	**
miR1120c-5p	1.00	11.72	3.55	**
miR1127b-3p	1.00	29.25	4.87	**
miR1140	649.02	1.00	−9.34	**
miR1918	23.54	1.00	−4.56	**
miR2111c	192.69	1.00	−7.59	**
miR4391	10.41	1.00	−3.38	**
miR5037c	766.99	102.29	−2.91	**
miR5054	90.33	27.17	−1.73	**
miR5060	52.80	1.00	−5.72	**
miR5374-5p	237.12	1.00	−7.89	**
miR5377	127.26	1.00	−6.99	**
miR5384-3p	24.99	7.23	−1.79	**
miR5662	935.11	105.87	−3.14	**
miR5800	35.40	13.13	−1.43	**
miR6117	1.00	125.81	6.98	**
miR6135f	1003.69	1.00	−9.97	**
miR7513	64.31	27.34	−1.23	**
miR7696a-3p	1.00	7367.56	12.85	**
miR7731-3p	37.53	1.00	−5.23	**
miR7810	541.72	181.57	−1.58	**
miR8651	1.00	305.63	8.26	**
miR8714	27.21	6.15	−2.15	**
miR9674a-3p	14.59	1.00	−3.87	**
miR9674b-3p	1.00	9.64	3.27	**
miR9763	1.00	1733.16	10.76	**
NmiRNA07	33.52	1.00	−5.07	**
NmiRNA08	20.73	1.00	−4.37	**
NmiRNA15	19.45	1.00	−4.28	**
NmiRNA23	1.00	71.38	6.16	**
NmiRNA25	1.00	11.30	3.50	**

Note The list includes all the miRNAs that are significantly differentially expressed between control and drought-stressed plants. Mean expression levels are listed in transcripts per kilobase million for each group, and the fold-change in expression in drought-stressed plants is shown relative to control plants. The table is ordered by miRNA name. The chi-square test was used to calculate p values.

** p < 0.01 indicates significant differences.

Identification of drought-responsive miRNA targets

To determine the function of drought-responsive miRNAs in XF 20 wheat, degradome sequencing was used to identify the miRNA targets. A total of 20,665,443 clean reads were obtained from the degradome cDNA libraries; most were 20- or 21-nt long. The reads were aligned to wheat transcriptome sequences: 17,591,286 reads mapped to the assembled transcriptome sequences. The CleaveLand pipeline was used for identification of target mRNAs; 222 potential target mRNAs that are predicted to be cleaved by known and novel miRNAs were identified from the CleaveLand degradome library analysis (Supplementary information Figure S1). To validate target genes identified through degradome sequencing, miRNA and transcript sequences were submitted to psRNA. There were 90 (90/222, 44.6%) targets from the degradome that were present in the psRNA target results list. These results indicated that degradome sequencing is an effective method to identify cleaved targets of miRNAs.

Analysis of target genes

To assess the functions of differentially expressed miRNAs, the target gene list was BLAST searched against a combined database of Arabidopsis, rice and B. distachyon. Considering all expressed miRNAs, a total of 2,905 target mRNAs were identified; 512 of these target mRNAs were targets of differentially expressed miRNAs. These target mRNAs were analyzed by gene ontology (GO) (Figure 1). Relative to the targets of all expressed miRNAs, GO categories ‘signal transducer’ and ‘response to stress’ showed a higher number of differentially expressed miRNAs than would be expected; the category ‘antioxidant’ contained fewer differentially expressed miRNAs than would be expected. Thus, in drought conditions, differentially expressed miRNAs reduced the inhibition of antioxidant genes, increasing the expression of antioxidant genes to resist drought, while the opposite occurred for signal transduction genes.

Figure 1. Gene ontology analysis of miRNA target genes identified in drought stress of wheat (Triticum aestivum L. cv. XF 20). Note: The digits on the left y-axis and right y-axis show the percentage and number of miRNA targets, respectively. Targets of all expressed miRNAs are shown in cyan; targets of all differentially expressed miRNAs are shown in pink.

Among the identified target genes, many have reported roles in plant responses to drought and other biotic and abiotic stresses. To give examples, we considered genes whose mRNAs were found to be targets of the differentially expressed miRNAs in our CleaveLand degradome library analysis. Novel miRNAs 08 and 15 were identified in this study (Supplementary information Table S2) and were downregulated in wheat cv. XF 20 in response to induced drought stress (Table 2). These miRNAs were predicted to target mRNA of genes for disease resistance proteins PRM1-like and RPP13-like. In Arabidopsis, these proteins provide resistance to the bacterial and pathogens Pseudomonas syringae and Hyaloperonospora parasitica (formerly Peronospora parasitica), respectively [44, 45]. Notably, NmiRNA15 also targeted mRNA encoding sucrose synthase. Sucrose is the main carbohydrate transported to developing plant organs. In soybean nodules, sucrose synthase was found to decline dramatically in activity and content in response to drought, but the nodule content of sucrose increased [46]. In tall fescue, the sucrose content of leaves increased in response to drought stress [47]. In wheat, it was recently shown that sucrose synthase gene expression was significantly higher in a drought-tolerant cultivar than in a drought-sensitive one. More sucrose accumulated in the tolerant cultivar under drought stress [48]. Similarly, the drought-tolerant cultivar JD17 accumulated sucrose on prolonged drought-stress [49]. Thus, sucrose levels appear to be crucial to the response of wheat to drought, and our results indicate that one response mechanism may involve regulation of sucrose synthase via NmiRNA15.

NmiRNA15 also targeted phytochrome and flowering time regulatory protein (PFT1, also called Mediator 25), a multifunctional nuclear protein involved in response to light and hormone signalling, defence against fungus, and abiotic stress response (including to drought) [50, 51]. This protein positively and negatively regulates drought response, for example, through interactions with the transcription factors dehydration responsive element binding protein 2 and ethylene response factor 1 [51].

Among the previously reported miRNAs that were found in this study to be differentially regulated in response to drought-stress (Table 2), miR482a-5p has been implicated in the response to pathogens in mulberry trees [52]. MiR482 was downregulated in response to drought in T. aestivum cv. Sivas 111/33 [2]. In our analysis, miR482a-5p, as well as miR9674a-3p, was predicted to target mRNA encoding a protein annotated as a stripe rust resistance protein. Stripe rust is a serious disease of wheat caused by the fungus Puccinia striiformis var. tritici. MiR2111 has been reported to play a role in phosphate homeostasis and plant survival [53]. Here, notably, miR2111c targeted 6-phosphogluconate dehydrogenase (6-PGDH), one of the key enzymes of the oxidative pentose phosphate pathway. In rice, transcript levels and enzyme activity of 6-PGDH were upregulated in response to drought, cold and high salinity [54]. Enhanced 6-GDPH activity has been shown to be related to aluminium tolerance in wheat [55]. MiR5054, which was also downregulated in response to drought stress in our study (Table 2), was predicted to target mRNA encoding epoxide hydrolase (EH). EH, which is involved in lipid metabolism, was previously found to be upregulated in a drought-resistant variety of wild emmer wheat (Triticum turgidum ssp. dicoccoides), the progenitor of cultivated wheat [56]. In Arabidopsis, an EH is induced in stems and leaves in response to drought [57].

Conclusions

We report the identification of differentially expressed miRNAs in drought-stressed wheat plants (Triticum aestivum L. cv. XF 20). This variety of wheat is particularly resistant to drought conditions, especially heat and dry wind. Our results shed light on the molecular mechanisms that influence drought tolerance in wheat. Targets of the drought-responsive miRNAs include many factors with well-characterized roles in stress responses, including drought. These miRNAs and their target genes may be exploited via breeding programmes or genetic engineering to improve drought tolerance in other varieties of wheat and related plants.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

This work was supported by the University Research Project of Gansu Province under grant number 2015A-145 and Longdong University Doctoral Scientific Research Fund under grant number XYBY1406. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.