Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Genome-wide identification and comparative analysis of drought-related microRNAs in two maize inbred lines with contrasting drought tolerance by deep sequencing

  • Xuyang Liu,

    Roles Investigation, Visualization, Writing – original draft

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Xiaojing Zhang,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Baocheng Sun,

    Roles Investigation

    Affiliation Institute of Grain Crops, Xinjiang Academy of Agricultural Sciences, Urumqi, China

  • Luyang Hao,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Cheng Liu,

    Roles Investigation

    Affiliation Institute of Grain Crops, Xinjiang Academy of Agricultural Sciences, Urumqi, China

  • Dengfeng Zhang,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Huaijun Tang,

    Roles Investigation

    Affiliation Institute of Grain Crops, Xinjiang Academy of Agricultural Sciences, Urumqi, China

  • Chunhui Li,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Yongxiang Li,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Yunsu Shi,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Xiaoqing Xie,

    Roles Investigation

    Affiliation Institute of Grain Crops, Xinjiang Academy of Agricultural Sciences, Urumqi, China

  • Yanchun Song,

    Roles Investigation

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Tianyu Wang ,

    Roles Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing

    wangtianyu@caas.cn (TW); liyu03@caas.cn (YL)

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

  • Yu Li

    Roles Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing

    wangtianyu@caas.cn (TW); liyu03@caas.cn (YL)

    Affiliation Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing, China

Abstract

Drought has become one of the most serious abiotic stresses influencing crop production worldwide. Understanding the molecular regulatory networks underlying drought adaption and tolerance in crops is of great importance for future breeding. microRNAs (miRNAs), as important components of post-transcriptional regulation, play crucial roles in drought response and adaptation in plants. Here, we report a miRNome analysis of two maize inbred lines with contrasting levels of drought tolerance under soil drought in the field. Differential expression analysis showed 11 and 34 miRNAs were uniquely responded to drought in H082183 (drought tolerant) and Lv28 (drought sensitive), respectively, in leaves. In roots, 19 and 23 miRNAs uniquely responded to drought in H082183 and Lv28, respectively. Expression analysis of these drought-responsive miRNA-mRNA modules revealed miR164-MYB, miR164-NAC, miR159-MYB, miR156-SPL and miR160-ARF showed a negative regulatory relationship. Further analysis showed that the miR164-MYB and miR164-NAC modules in the tolerant line modulated the stress response in an ABA (abscisic acid)-dependent manner, while the miR156-SPL and miR160-ARF modules in the sensitive line participated in the inhibition of metabolism in drought-exposed leaves. Together, our results provide new insight into not only drought-tolerance-related miRNA regulation networks in maize but also key miRNAs for further characterization and improvement of maize drought tolerance.

1 Introduction

To meet the rising demand resulting from population growth and economic development, crop production should maintain a sustainable increase [1]. However, the increasing environmental stresses resulting from global climate change, such as drought, have become major threats to crop production [2]. Maize (Zea mays L.) is one of the most important crops for food and feed worldwide [3]. Although maize yield has grown over the past decades, maize has also demonstrated increasing sensitivity to drought [4]. Thus, breeding of maize varieties with drought tolerance is an important alternative approach to relieving the threat; however, it requires a better understanding of the molecular regulatory networks of drought tolerance.

microRNAs (miRNAs), which are endogenous, short, single-stranded non-coding RNAs, play key roles in the stress response and various developmental processes in plants as post-transcriptional regulators [5, 6]. miRNAs suppress the expression of target genes by binding to the partially complementary sequences in messenger RNAs (mRNAs) and then cleaving and degrading mRNA molecules or inhibiting translation with the cooperation of RNA-induced silencing complexes (RISCs) [7]. The target genes of miRNAs can be predicted based on bioinformatic analysis of miRNA-mRNA complementary sequences [8]. And the degradome sequencing also provides supporting evidence of interactions between miRNAs and their target genes [9, 10]. Most target genes of miRNAs in plants are transcription factors (TFs, 65.9%), followed by NB-LRR proteins (nucleotide-binding domain and leucine-rich repeat domain, 6.6%), pathogen proteins (2.2%), lncRNAs (long non-coding RNAs, 1.1%) and other proteins (24.2%) [5].

The genome-wide identification of miRNAs that respond to drought has been reported in plants by using microarrays [1115] or by RNA sequencing (RNA-Seq) [1620]. Numerous miRNAs have been reported to respond to drought in different plant species, such as Arabidopsis [21], rice (Oryza sativa L.) [19], wheat (Triticum aestivum L.) [22] and soybean (Glycine max L.) [16]. For example, in Arabidopsis, miR393, miR397b and miR402 have been shown to be up-regulated while miR319c and miR389a were down-regulated under dehydration stress [23]. In another study, 17 miRNAs, including miR164c, miR319b, and miR1861d, were down-regulated under drought in rice at the seedling stage, while 16 miRNAs (miR166h, miR172d, miR408, etc.) were up-regulated [24].

In the maize genome, 321 mature miRNAs derived from 172 precursors, which could be divided into 29 families, were identified [25, 26]. Many miRNAs from different families have been reported to respond to drought in maize [27]. For example, a maize inbred line with drought tolerance was treated by simulated drought (16% polyethylene glycol, PEG) at the seedling stage and miRNA expression analysis showed that 68 miRNAs responded to drought, with 25 up-regulated and 43 down-regulated miRNAs [15]. Two drought-tolerant and two drought-sensitive inbred lines were treated by simulated drought (16% PEG for 16 and 24 hours), and the results of the miRNA microarray analysis identified 303 differentially expressed miRNAs under drought, with 177 up-regulating and 126 down-regulating [28]. It was also found that eight miR169 members could target seven NF-YA (NUCLEAR FACTOR-Y subunit A) genes in maize, and they responded to PEG treatment in reverse correlation patterns in leaves and roots [29, 30].

Previous studies showed that the responses of miRNAs to drought stress are highly complex because the expression patterns of the same miRNAs under drought might vary among plant species, tissues, treatment methods or even different genotypes of the same species [27]. For example, miR172 was up-regulated in Arabidopsis while it was down-regulated in rice in response to drought [14, 31]. Four miRNAs (miR397a/b, miR398b, miR408-5p and miR528-5p) were down-regulated in drought-tolerant rice varieties but up-regulated in a drought-sensitive variety [19]. With respect to drought treatments, rapid soil drying in limited pots and under simulated drought with PEG were shown to vary significantly from real drought conditions in the field [32]. Hence, different drought treatments may lead to contradictory results pertaining to miRNA regulation patterns. For instance, in maize, miR398 was up-regulated under PEG treatment but down-regulated under soil drought [13, 15]. Therefore, it is of great importance to detect the actual regulatory mechanisms of the miRNome under soil drought conditions in the field from the viewpoint of plant breeders.

To gain a better understanding of the regulatory mechanisms performed by miRNAs to improve drought tolerance, we report a miRNome study using two maize inbred lines with contrasting levels of drought tolerance during soil drought in the field. The aims of this study were to identify the responding miRNAs under soil drought, compare the different regulatory patterns between the drought-tolerant and drought-sensitive lines and investigate the candidate miRNAs with potential for genetic improvement of drought tolerance in maize.

2 Materials and methods

2.1 Plant materials, field experimental design and leaf relative water content measurement

Two maize inbred lines, i.e., H082183 and Lv28, were used in this study. H082183 was an inbred line with strong drought tolerance that was obtained by double haploid (DH) technology and the seeds for which were provided by Prof. Jiuran Zhao of the Maize Research Center of Beijing Academy of Agricultural and Forestry Sciences. Lv28 was a drought-sensitive inbred line and a representative of the Luda Red Cob heterotic group in China [33]. The field experimental design has been previously described [34]. Briefly, these two materials were planted in the field at Urumqi (Xinjiang, China; 43.98°N, 87.51°E), which features limited rainfall and an enormous evaporation capacity. For each inbred line, two water treatment and six biological replicates were established with two row plots containing 11 plants for each row. The rows were 3 meters long and spaced 0.6 meters apart. All plants were watered by drip irrigation at the germination stage at sowing. Plastic films for uniform germination were removed at stage V4. Plants in the well-watering treatment were then irrigated every week, while watering was withheld from plants in the drought treatment. To evaluate the physiological and molecular responses of the plants to drought, the last fully expanded leaves of both genotypes and treatments were sampled to measure the leaf relative water content (RWC); moreover, the second-to-last fully expanded leaves (10 cm leaf tips) and roots (10 cm root tips) were sampled and frozen in liquid nitrogen immediately and stored at -80°C for small RNA sequencing. Sampling was performed at 11:00 (GMT +8:00; 9:00 local time), and the RWC measurement was performed according to the Barrs method [35].

2.2 RNA extraction, library construction and Illumina sequencing

According to the leaf RWC results, the samples at 27 and 46 days after drought (DAD), which corresponded to moderate drought and severe drought, respectively, were chosen for sequencing. The samples at 27 and 46 days after drought under the well-watered condition were also selected as the moderate drought control and severe drought control. For each genotype and treatment, the leaf and root samples from two individual plants with the closest RWC in six replicates were selected as two biological replicates for sequencing. In total, 32 samples of the two genotypes (H082183 and Lv28), two tissues (leaf and root), two treatments (drought and well-watered), two drought levels (moderate and severe drought) and two biological replicates were used for miRNA and mRNA sequencing. The mRNA sequencing results have been published in our previous studies [34] (SRA accession number: SRP102142, SRP119805).

For miRNA sequencing, the total RNA of 32 samples was extracted using TRIzol (Invitrogen, USA) according to the manufacturer’s protocol. The purity of the total RNA was examined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA). The concentration and integrity of the RNA were measured with a Qubit 2.0 fluorometer (Thermo Fisher Scientific, USA) and Bioanalyzer 2100 (Agilent Technologies, USA), respectively. The RNA samples with OD260/280 > 1.8 and OD260/230 > 1.0 and RNA concentration > 250 ng/μL and RIN > 8.0 were selected for the following small RNA library construction. Small RNA sequencing libraries were constructed using the Next Ultra Small RNA Sample Library Prep Kit for Illumina (NEB, USA). Briefly, small RNA molecules were ligated with 3’ and 5’ adapters by T4 RNA Ligase. Subsequently, RNA samples were reverse-transcribed and amplified by PCR. Then, the cDNA was separated by 6% polyacrylamide gel electrophoresis (PAGE), and the short molecules were recovered and purified. The Illumina HiSeq 2500 (Illumina, USA) platform was employed to sequence the 32 small RNA samples at Biomarker Technologies (China).

2.3 Alignment and mapping of the small RNA sequence data

The sequence read length was set to 50 nucleotides (nt). Raw reads with low sequencing quality, and with lengths shorter than 18 nt or longer than 30 nt, were eliminated to produce clean data. The clean read data were deposited in the NCBI Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra, accession number: SRP119798). To identify known non-coding RNAs and repetitive sequences, clean reads were aligned to ribosome RNA (rRNA) sequences in Silva (http://www.arb-silva.de) [36], transfer RNA (tRNA) sequences in GtRNAdb (http://lowelab.ucsc.edu/GtRNAdb) [37], small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA) sequences in Rfam (http://rfam.xfam.org) [38], and repetitive sequences in Repbase (http://www.girinst.org/repbase) [39], using Bowtie software [40]. The clean reads were also aligned to the maize reference genome (Zea mays AGPv3.24, ftp://ftp.ensemblgenomes.org/pub/plants/release-24/fasta/zea_mays) [3] using miRDeep2 [41].

2.4 Identification of known and novel miRNAs

miRDeep2 [41] was used to align the clean reads to miRBase release 22 (http://www.mirbase.org) [26] to identify known miRNAs. The identification of novel miRNA in maize was performed according to the 2018 criteria for plant miRNA annotations [42]. First, the reads obtained by small RNA sequencing were mapped to the reference genome by miRDeep2, and the secondary structure of miRNA precursors was predicted by the RNAfold program [43]. The resulting secondary structure was analyzed to verify characteristic miRNAs based on the following criteria: (1) miRNA-miRNA* duplex formed with two-nucleotide 3’ overhangs and a foldback size less than 300 nucleotides; (2) no more than five mismatched positions, only three of which were nucleotides in asymmetric bulges; (3) an miRNA length ≥ 20 and ≤ 24; and (4) novel miRNAs met all criteria in at least three small RNA sequencing libraries.

For the candidate novel miRNAs, a stricter criterion was used to obtain credible prediction results, which stated that the novel miRNAs should contain the target sequence, as confirmed by degradome sequencing. The novel miRNAs that fit all the criteria were also confirmed by stem-loop reverse transcription PCR. Small RNAs of the maize inbred lines B73, H082183 and Lv28 were extracted using the miRcute miRNA isolation kit (Tiangen Biotech, China) and reverse-transcribed to cDNA with stem-loop reverse transcription primers using the miRcute plus miRNA first-strand cDNA synthesis kit (Tiangen Biotech, China). PCR amplification was performed using 2×EasyTaq PCR SuperMix (Transgen Biotech, China). The sequences of primers were listed in (Table A in S1 File).

The pre-miRNA and mature miRNA sequences of maize (Zea mays), Arabidopsis thaliana, rice and sorghum (Sorghum bicolor L.) were downloaded from miRbase release 22 (http://www.mirbase.org) [26] and used in the phylogenetic analysis performed with MEGA5 [44].

2.5 microRNA expression analysis

For each miRNA, gene expression was calculated and normalized using the TPM algorithm [45]. The TPM value was used for principal component analysis (PCA) and cluster analysis of the miRNome profiling of each sample. DESeq [46] was used to analyze the differentially expressed genes. To detect the drought-responsive miRNAs of each genotype, differential expression analysis was performed for the same genotype subjected to different watering treatments (drought or well-watered conditions), with four comparison groups, i.e., H082183 under moderate drought vs. well-watered control, Lv28 under moderate drought vs. well-watered control, H082183 under severe drought vs. well-watered control, and Lv28 under severe drought vs. well-watered control. The criteria for differential expression were established as |log2(fold change)| > 1 and FDR < 0.05.

2.6 Prediction and annotation of microRNA target genes

The target genes of miRNA were predicted by TargetFinder software [8] with default parameters. The function and annotation of target genes was analyzed using the Blast to Swiss-Prot (http://www.uniprot.org/) [47], NCBI non-redundant protein sequence (Nr) (ftp://ftp.ncbi.nih.gov/blast/db/) [48] and Pfam databases (http://pfam.xfam.org/) [49]. The drought-responsive genes, obtained by mRNA sequencing in our previous studies [34], were calculated using the same methods with miRNAs, with criteria of |log2(fold change)| > 1 and FDR < 0.05 (some genes with |log2(fold change)| > 0.6 and FDR < 0.05 were also considered differentially expressed genes).

2.7 Identification of co-expressed genes and binding sequence of transcription factor genes

Clustering analysis of mRNA genes was performed to identify the genes co-expressed with miRNA target genes. The k-means cluster analysis was performed using R, with k = 25 and n of the cluster = 20. The co-expressed genes with miRNA target transcription factor genes were used for the binding sequence analysis. The motif sequence scan of cis-regulatory elements of co-expressed genes was performed by FIMO [50] using the 2-k sequences in the promoter regions, with a criterion of P < 0.0001. The over-representation of cis-regulatory elements in promoters was identified by Fisher’s exact test, P = F(x, y, n, N), where x was the number of a considered cis-elements in an individual promoter, y was the total number of cis-elements in an individual promoter, n was the number of a considered cis-element in promoters of all genes in maize, and N was the total number of cis-elements in promoters in all genes in the maize genome. The significance level of over-representation was defied as P < 0.001. Gene Ontology (GO) annotation analysis was performed with the Singular Enrichment Analysis (SEA) tool of agriGO (http://bioinfo.cau.edu.cn/agriGO/) [51], using Zea mays AGPv3.30 as the reference. The significance of enriched GO terms was calculated by the Fisher statistical test, and the FDR value were obtained by Bonferroni multi-test adjustment. The threshold of significant enriched GO terms was set to FDR < 0.05.

2.8 Degradome sequencing

Two bulked (mixed) samples of H082183 and Lv28 were used for degradome sequencing. The H082183 and Lv28 were planted in the greenhouse and watered every three days. After three weeks of germination, drought treatment was performed by withholding water for 20 days, while the plants of the well-watered treatment were watered normally every three days. The leaves and roots under drought and well-watered treatments of H082183 and Lv28 were sampled and mixed for degradome sequencing. The degradome libraries were constructed following the method of Ma et al. [52]. Briefly, the cleaved mRNAs with 5’ monophosphates were ligated to adapters containing a MmeI recognition site. The first strand cDNA was produced using oligo(dT) and PCR. The PCR product was digested with MmeI. After that, a double DNA adapter was ligated to digested PCR product with T4 DNA ligase and amplified by PCR. Then, single-end sequencing of degradome libraries was performed on an Illumina HiSeq 2500 (Illumina, USA). The CleaveLand pipeline [9] was used to identify potentially cleaved targets.

2.9 Gene expression under ABA treatment

The H082183 were planted in the greenhouse. After three weeks of germination, the plants were irrigated with 100 μM ABA (abscisic acid). The leaves and roots of H082183 were sampled just before ABA treatment (0 h) and at 6 h, 12 h and 24 h after ABA treatment. Three replicates, which contained three plants in each, were used for subsequent RNA extraction and gene expression analysis.

For the qRT-PCR of mRNA, reverse transcription of mRNA to cDNA was performed using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech, China). qRT-PCR was performed using the PowerUp SYBR Green Master Mix (Applied Biosystems, USA) on an ABI QuantStudio3 (Applied Biosystems, USA). Each 20 μL PCR contained 1 μL cDNA, 10 μL PowerUp SYBR Green Master Mix (2X), 0.5 μL forward primer, 0.5 μL universal reverse primer and 8 μL ddH2O. The thermal cycling parameters were as follows: 50°C for 2 minutes, 95°C for 2 minutes, 40 cycles of denaturing at 95°C for 15 seconds and annealing extension at 60°C for 1 minute, followed by melting curve analysis. The maize GAPDH gene was used as an internal reference. The expression level was calculated by 2-ΔΔCt method [53]. Each sample was also performed with three replicates. The sequences of the primers used in the qRT-PCR are listed in (Table A in S1 File).

2.10 Quantitative real-time RT-PCR of miRNAs

To confirm the miRNA expression obtained by sequencing, quantitative real-time RT-PCR (qRT-PCR) of five miRNAs (miR166e, miR172b-3p, miR319b-3p, miR408b-3p and miR528a-5p) in leaves and five miRNAs (miR1432-5p, miR159b-3p, miR172b-3p, miR398a-5p and miR399j-5p) in roots of H082183 and Lv28 under drought was performed. The small RNAs of 16 samples under drought were extracted using the miRcute miRNA Isolation Kit (Tiangen Biotech, China). The first-strand cDNA fragments were synthesized using the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (Tiangen Biotech, China). Briefly, a poly(A) tail was added to miRNAs by E. coli Poly(A) Polymerase. Then, Oligo(dT)-Universal Tag primers were used to generate reverse-transcript cDNA of the miRNAs. qRT-PCR was performed using the miRcute miRNA qPCR Detection Kit (SYBR Green) (Tiangen Biotech, China) on an ABI7500 (Applied Biosystems, USA). Each 20 μL PCR contained 1μL cDNA, 10 μL 2×miRcute miRNA Premix (with SYBR&ROX), 0.4 μL forward primer, 0.4 μL universal reverse primer and 8.2 μL ddH2O. The thermal cycling parameters were as follows: 95°C for 2 minutes, 40 cycles of denaturing at 95°C for 20 seconds and annealing extension at 60°C for 34 seconds, followed by melting curve analysis. The maize U6 gene was used as an internal reference. The expression level was calculated by 2-ΔΔCt method [53]. Each sample was performed with three replicates. The expression levels of miRNAs in H082183 and Lv28 under moderate drought or severe drought were transformed to log2(H082183/Lv28) for comparison with the sequencing results.

3 Results

3.1 Physiological response of two maize inbred lines to drought stress in the field

Two maize inbred lines, Lv28 and H082183, were grown, and leaf relative water contents (RWCs) were tested under well-watered and soil drought conditions in the field. The results showed that the leaf RWCs of both genotypes under drought were significantly (P < 0.001) lower than those of their well-watered controls at 27 and 46 days after drought (Fig 1A). At 27 days after drought, the leaf RWCs of H082183 and Lv28 was 89.30% and 84.31% under drought, respectively, and showed significant difference (P < 0.001). However, the leaf RWCs of both H082183 and Lv28 decreased to approximately 82%, and showed no significant difference at 46 days after drought.

thumbnail
Fig 1. The change in leaf relative water contents under drought and the sample relationship based on the miRNome results.

(A) The leaf relative water contents (RWCs) of H082183 and Lv28 under drought and well-watered treatments. **: P < 0.01; *: P < 0.05; NS: P > 0.05. (B) Cluster analysis of miRNome of leaf samples. (C) Cluster analysis of miRNome of root samples. (D) Principal component analysis (PCA) of leaf samples. (E) PCA of root samples. HMD: H082183 under moderate drought, HSD: H082183 under severe drought, LMD: Lv28 under moderate drought, LSD: Lv28 under severe drought. HMC: H082183 under moderate drought control, HSC: H082183 under severe drought control, LMC: Lv28 under moderate drought control, LSC: Lv28 under severe drought control.

https://doi.org/10.1371/journal.pone.0219176.g001

3.2 Sequencing and alignment

To understand drought-responsive and tolerance-related miRNAs in maize, we selected the leaf and root samples of both genotypes at 27 and 46 DAD, representative of moderate and severe drought, for sequencing. In addition to the well-watered controls of moderate drought and severe drought, small RNA libraries of 32 samples were sequenced. In total, 743.8 million raw short reads were obtained, with 18.6 million raw reads per sample on average (Table B in S1 File). Reads of low quality and substandard sequence length (fewer than 18 nucleotides or greater than 30 nucleotides) were removed. Therefore, a total of 533.6 million clean reads were acquired, with a mean of 13.3 million clean reads per sample. The Q30 scores of all samples ranged from 94.54% to 99.13%, with an average of 97.32%, indicating that the sequence data were of reliable quality. The clean reads were aligned to the maize reference genome, and the percentage of mapped reads in each library’s clean reads ranged from 13.26% to 29.60%, with a mean value of 18.37% (Table C in S1 File). The clean reads were mapped to different small RNA databases using Bowtie to differentiate the types of small RNAs. The most abundant small RNA was rRNA, with an average percentage of 43.73%, followed by tRNA (7.77%), repeated sequences (1.60%) and snoRNA (0.17%) (Table D in S1 File). The unannotated reads were aligned to miRBase (release 22) to identify the known miRNAs in maize. It was found that 278 and 275 known miRNA genes were expressed in leaf and root samples, respectively (Table E in S1 File). Those known miRNAs could be classified into 28 miRNA families (Fig A in S1 File). In addition, the length distribution of the expressed known miRNAs was highly enriched in 21, 20 and 22 nucleotides (nt) (Fig B in S1 File).

Two bulked (mixed) samples of H082183 and Lv28 were used for degradome sequencing to confirm the predicted target genes of the miRNAs. A total of 30.6 and 25.4 million reads were obtained by degradome sequencing in H082183 and Lv28, respectively (Table F in S1 File). In addition, 99.59% and 99.55% of the degradome sequencing reads could be mapped to the maize reference genome (Zea mays AGPv3.24), and covered 36417 and 37862 transcripts in H082183 and Lv28, respectively.

3.3 Novel miRNAs in maize

A total of 34 novel miRNAs were identified based on the 2018 criteria for plant miRNA annotations [42] (Table G in S1 File). The lengths of novel miRNAs were enriched in 21 and 24 nt. Six of the novel miRNAs were conserved with other plant species (Fig C in S1 File). Three novel miRNAs (named Novel_2_38, Novel_10_88 and Novel_10_147 according to their positions on chromosomes) were confirmed by degradome sequencing (Fig 2). The Novel_2_38, Novel_10_88 and Novel_10_147 miRNAs could cleave the 33873, 4231, and 1783 site of GRMZM2G308707, AC190628.4_FG007, and GRMZM2G064212, respectively. The phylogenetic results showed that the novel miRNAs Novel_2_38, Novel_10_88 and Novel_10_147 were close to osa-miR1905, osa-miR5534 and sbi-miR6229, respectively (Fig D in S1 File). However, the mature sequences of these three novel miRNAs did not have conservative seed regions with their relative miRNAs in other species, indicating these novel miRNAs were not annotated by homologies in plants. In addition, the Novel_2_38, Novel_10_88 and Novel_10_147 were confirmed by stem-loop reverse transcription and PCR amplification in different maize lines.

thumbnail
Fig 2. The second structure and target sequence of predicted novel miRNAs.

The stem-loop structure and miRNA-miRNA* duplex of Novel_2_38 (A), Novel_10_88 (B) and Novel_10_147 (C). The target sequence and cleavage read count by degradome sequencing of Novel_2_38-GRMZM2G308707 (D), Novel_10_88-AC190628.4_FG007 (E) and Novel_10_147-GRMZM2G064212 (F).

https://doi.org/10.1371/journal.pone.0219176.g002

3.4 Expression profiles of miRNAs

The miRNA expression profiles of each of the two biological replicates samples were all significantly correlated (r > 0.95, P < 0.01) (Fig E in S1 File). Principal component analysis (PCA) showed the samples were closest to their corresponding biological samples in both leaves and roots (Fig 1D and 1E). The clustering of different leaf samples showed that the samples were first divided into two groups by genotypes and then separated by treatment (Fig 1B), while the root samples were initially clustered by drought treatments and subsequently by different genotypes (Fig 1C).

3.5 Drought-responsive miRNAs

The drought-responsive miRNAs in HMD, LMD, HSD and LSD were acquired by comparing the expression profiles of HMC-HMD, LMC-LMD, HSC-HSD and LSC-LSD. In leaves, there were three (miR397, miR398, miR1432) and five miRNAs (miR399s and miR827s) that were commonly down- and up-regulated under drought in both H082183 and Lv28, respectively (Fig 3A), while seven (miR164, miR169, miR398, miR528s and miR408s) and four miRNAs (miR395 and miR399s) were exclusively down- and up-regulated in response to drought, respectively, in leaves of H082183. Additionally, 18 (miR156s, miR160s, miR162, miR164s, miR171s and miR399s) and 16 (miR156, miR159, miR166s, miR168s, miR171, miR172s and miR444s) miRNAs were down- and up-regulated, respectively, under drought in leaves of Lv28. In roots, four (miR156s and miR168) and five miRNAs (miR166s and miR399s) were commonly down- and up-regulated under drought in both H082183 and Lv28, respectively (Fig 3B). In addition, four (miR398s and miR171) and fifteen (miR159s, miR167s, miR390s and miR399s) miRNAs were exclusively down- and up-regulated in roots of H082183, and twenty-one (miR156s, miR167 and miR172s) and two miRNAs (miR399 and miR827) were uniquely down- and up-regulated under drought in roots of Lv28, respectively. Interestingly, one member of the miR156 family (miR156e-3p) was up-regulated in roots of H082183 but down-regulated in roots of Lv28. However, no predicted novel miRNAs showed a response to drought in leaves or roots.

thumbnail
Fig 3. Drought-responsive miRNAs in leaves and roots.

(A) The drought-responsive miRNAs in leaves. (B) The drought-responsive miRNAs in roots. The arrows in the Venn plot indicate down- or up-regulation of miRNAs. MD: moderate drought, SD: severe drought.

https://doi.org/10.1371/journal.pone.0219176.g003

3.6 Target genes of drought-responsive miRNAs

Among all the miRNAs expressed in this study, 211 were predicted to target 262 mRNAs in maize leaves or roots (Table H in S1 File). The expression analysis revealed several target miRNA-mRNA modules that had negative expression patterns. In the leaves of H082183, miR164 was down-regulated under moderate drought, while its target genes GRMZM2G096358 (MYB) and GRMZM2G114850 (NAC) were up-regulated under MD (Fig 4). In the roots of H082183, three miRNA-mRNA pairs were negatively regulated under drought, including miR159 and two MYB genes (GRMZM2G004090 and GRMZM2G028054), miR390 and one leucine-rich repeat protein (LRR) gene (GRMZM2G304745), miR398 and one selenium-binding protein (SBP) gene (GRMZM2G103812). In the drought-sensitive line Lv28, miR160 was up-regulated under drought, while its target auxin response factor (ARF) genes AC207656.3_FGT002, GRMZM2G390641 and GRMZM2G159399 were down-regulated under drought. The miR156 displayed the opposite regulation patterns in leaves and roots of Lv28. The up-regulation of miR156 negatively modulated the squamosa promoter binding like protein (SPL) genes GRMZM2G065451, GRMZM2G113779, GRMZM2G414805, GRMZM2G126018, GRMZM2G371033 and GRMZM2G067624. However, miR156 was down-regulated in the roots of Lv28, and the target SPL genes (GRMZM2G113779, GRMZM2G160917, GRMZM5G806833, GRMZM2G307588, GRMZM2G460544 and GRMZM2G061734) were up-regulated under drought. These miRNA-mRNA modules were confirmed by degradome sequencing (Fig F in S1 File).

thumbnail
Fig 4. Regulation networks of drought-responsive miRNAs and their target mRNAs.

The left and right panel include the miRNA-mRNA expression in Lv28 and H082183 respectively. (A) GRMZM2G065451, (B) GRMZM2G113779, (C) GRMZM2G414805, (D) GRMZM2G126018, (E) GRMZM2G371033, (F) GRMZM2G067624, (G) AC207656.3_FGT002, (H) GRMZM2G390641, (I) GRMZM2G159399, (J) GRMZM2G096358, (K) GRMZM2G114850, (L) GRMZM2G113779, (M) GRMZM2G160917, (N) GRMZM5G806833, (O) GRMZM2G307588, (P) GRMZM2G460544, (Q) GRMZM2G061734, (R) GRMZM2G004090, (S) GRMZM2G028054, (T) GRMZM2G304745, (U) GRMZM2G103812. The y-axis shows the log2(fold change) of miRNAs or mRNAs under drought. The x-axis presents the water treatment conditions. ND: well-watered treatment, MD: moderate drought, SD: severe drought.

https://doi.org/10.1371/journal.pone.0219176.g004

3.7 Regulatory networks of drought-responsive miRNAs

Transcription factor (TF) genes regulate gene expression by binding to the cis-element in promoter regions. To explore the regulatory networks of the miRNA-TF pairs, a clustering analysis was performed to detect the genes co-expressed with the target mRNA genes (Fig G in S1 File). For example, 341 genes were co-expressed with the miR164-targeted MYB (GRMZM2G096358) and NAC (GRMZM2G114850) genes in leaves of H082183. Additionally, 364 genes were co-expressed with the miR156-targeted MYB (GRMZM2G028054) gene in roots of H082183. Furthermore, the co-expressed genes with MYB, NAC, ARF or SPL binding motifs in their promoter regions were selected for GO enrichment analysis. The results revealed that the putative modulated genes of the miR164-MYB and miR164-NAC modules mostly involved the GO terms “response to stimulus”, “response to stress” or “response to water deprivation” in the leaves of H082183 (Table 1). In contrast, the miR156-SPL and miR160-ARF modules in the leaves of Lv28 were both related to “photosynthesis”. In roots, the miR159-MYB modules of H082183 were involved in “membrane”, “plasma membrane” and “response to stimulus”. The miR156-SPL modules were involved in “response to stress”, “response to stimulus” and “response to abiotic stimulus”. These results suggested that the miR164-MYB and miR164-NAC modules in the tolerant line H082183 were involved in the regulation of stress-responsive genes, while the miR156-SPL and miR160-ARF modules in the sensitive line Lv28 participated in inhibition of metabolism and development in leaves. Further analysis revealed 12 and 8 genes that were co-expressed with miR164-MYB and miR164-NAC and that had significantly over-represented MYB and NAC binding motif sequences in the 2-k promoter regions, respectively (Fig H in S1 File).

thumbnail
Table 1. GO enrichment of genes co-expressed with miRNA target transcription factors.

https://doi.org/10.1371/journal.pone.0219176.t001

The MYB and NAC transcription factors are known to regulate the drought stress response via abscisic acid (ABA) signaling pathways [54, 55]. Thus, the expression levels of the miR164-MYB, miR164-NAC and miR159-MYB modules were analyzed under ABA treatment in H082183. In leaves, miR164 was down-regulated, but the target genes GRMZM2G096358 (MYB) and GRMZM2G114850 (NAC) were up-regulated under ABA treatment (Fig 5). These results were consistent with the expression patterns of miR164-MYB and miR164-NAC modules under drought (Fig 4). However, the negative regulatory relationship of miR159-MYB was lost in roots of H082183 under ABA treatment, which were all down-regulated (Fig I in S1 File).

thumbnail
Fig 5. The expression of miR164 and target genes in leaves under ABA treatment.

(A) miR164. (B) GRMZM2G114850 (NAC). (C) GRMZM2G096358 (MYB).

https://doi.org/10.1371/journal.pone.0219176.g005

3.8 Confirmation of miRNA expression by qRT-PCR

A total of five miRNAs in leaves and five miRNAs in roots were selected for qRT-PCR to confirm the expression of miRNAs obtained by sequencing. The relative expression of miRNAs under drought were calculated as the log2(fold change) for comparison with the sequencing results. Correlation analysis showed that the gene expression changes between the two genotypes generated by qRT-PCR and sequencing were significantly correlated (r = 0.935, P = 3.397×10−10) (Fig J in S1 File).

4 Discussion

4.1 Drought-responsive miRNAs in plants

miRNAs, as post-transcriptional regulators, play significant roles in the stress responses of plants. Therefore, a few studies involving the genome-wide identification of drought-responsive miRNAs in various plants have been reported, including Arabidopsis [21], rice (Oryza sativa L.) [14, 19], and wheat (Triticum aestivum L.) [22]. In this study, many miRNAs were detected to respond to drought in leaves and roots. The drought-responsive miRNAs in different plant species were thus compared (Table I in S1 File).

In the leaves of maize, miR156, miR160, miR162, miR164, miR171, miR395, miR399 and miR827 were up-regulated under drought, while miR159, miR164, miR166, miR168, miR169, miR171, miR172, miR397, miR398, miR408, miR444, miR528 and miR1432 were down-regulated. Similarly, miR160, miR162, miR395, miR827 were also up-regulated, while miR166, miR172, miR397, miR827 and miR1432 were down-regulated in Arabidopsis, rice or wheat [14, 19, 21, 22]. However, certain miRNA families showed different regulatory patterns in plants, i.e., miR159, miR168, miR169, miR398, miR399, miR408, and miR528. For instance, miR159 was down-regulated in maize but up-regulated in rice and wheat under drought [14, 22]. miR168 was downregulated in rice and maize but up-regulated in Arabidopsis under drought [14, 21]. miR399 was up-regulated in maize and wheat but down-regulated in rice under drought [19, 22]. In addition, the members of three miRNA families (miR156, miR164 and miR171) showed different regulatory patterns under drought in maize, while miR156 and miR171 also displayed opposite regulations in two individual studies of rice under drought [14, 19].

In roots of maize, miR156, miR159, miR166, miR167, miR390, miR399 and miR827 were up-regulated, while miR156, miR159, miR167, miR168, miR171, miR172 and miR398 were down-regulated, under drought. However, when compared with other plants, only miR398 was conservatively down-regulated in the roots of maize and wheat under drought [22]. miR172 was down-regulated in the roots of maize and leaves of rice and wheat, while it was up-regulated in the roots of wheat under drought [14, 22]. miR399 was up-regulated in the roots of maize and leaves of wheat, but down-regulated in roots of wheat [19, 22]. Furthermore, the members of three miRNA families (miR156, miR159 and miR167) showed altered regulation in maize roots, while miR156 displayed opposite regulations in rice, and miR159 showed a contrasting response between leaves and roots of wheat under drought [14, 19, 22].

4.2 Comparison of miRNome generated under field drought conditions and other drought treatment methods

In maize, the genome-wide identification of drought-responsive miRNAs has been reported under simulated drought treatment (PEG) [15, 28] or by withholding water in limited pots in a greenhouse [13, 55]. In this study, drought-induced miRNome using different drought treatment methods were compared, and the regulatory patterns of miRNAs in response to the different drought treatment methods showed great diversity (Table J in S1 File). For example, miR159 and miR166 were down-regulated in maize leaves under soil drought in the field or in a greenhouse pot but up-regulated under PEG treatment [13, 28, 55]. miR160 was up-regulated in maize leaves under drought in the field but down-regulated under rapid drought in a limited pot [13]. miR156 was down-regulated under moderate drought but up-regulated under severe drought in the field; however, miR156 was up-regulated and down-regulated at 16 and 24 h after PEG treatment, respectively [15, 28]. miR399 was up-regulated in both leaves and roots under drought in the field but down-regulated under PEG treatment [15]. This finding indicated an inconformity between the effects of soil drought and drought simulated by PEG (dehydration stress) on the expression of mRNAs. Therefore, it is necessary to identify the drought-induced miRNome or transcriptome in the field to determine the actual expression regulatory networks.

4.3 miR164 contributed to drought tolerance in leaves of the tolerant line H082183

Transcript regulation at the gene expression level is a major part of the stress responses and has been reported to play a role in the stress tolerance of plants [56]. It is notable that, in plants, the target genes of miRNAs mostly encode transcription factors that play important roles in modulating gene expression at the transcriptomic level [5]. The interaction of miRNAs and transcription factors could regulate the drought stress response and tolerance via different signaling and physiological pathways, such as the ABA (abscisic acid) response, auxin signaling, and osmotic production and antioxidant production [57].

In the present study, the miRNome results suggested the putative role of the miR164-MYB and miR164-NAC modules in stress response regulatory pathways in leaves, and their expression was regulated by ABA. The miR164-NAC modules have been reported related to drought tolerance in rice. For example, Fang et al. found that the expression of six NAC genes, which were conservatively targeted by miR164 (OsmiRNA targeted NAC1-6, OMTN1-6), were significantly reduced in leaves under drought. The over-expression transgenic plants of OMTN genes showed increased sensitivity to drought stress, which displayed earlier leaf rolling and wilting [58]. In addition, miR164 could target ZmNAC1 and negatively regulate its expression, and over-expression of ZmNAC1 contributed to the increased lateral root growth in maize [59]. In wheat, the interaction between miR164 and TaNAC21/22 was confirmed, and this module could regulate the resistance to stripe rust [60].

Furthermore, the NAC and MYB transcription factor genes played important roles in signal transduction of the drought stress response in plants. Over-expression transgenic Arabidopsis of a sorghum NAC gene SbSNAC1 resulted in improved drought tolerance with a significantly higher survival rate than wild type plants under drought [61]. The expression of ZmSNAC1 was induced by dehydration treatment, and over-expression of ZmSNAC1 in Arabidopsis enhances tolerance to dehydration stress [62]. The expression of ZmNAC111 was significantly associated with drought tolerance in maize at the seedling stage [63]. In rice, the expression of OsMYB3R-2 was induced by drought, salt and cold stresses, whose over-expression of transgenic Arabidopsis led to increased tolerance to drought and salt stresses [64]. The wheat MYB gene TaMYBsm1 showed functions in improving drought tolerance, the over-expression of which in transgenic Arabidopsis led to higher germination rates under drought stress [65]. Additionally, both MYB and NAC TFs were participated in the ABA signal pathway [66]. For instance, the overexpression of TaNAC29 enhanced the expression of some core genes of the ABA signal pathway, such as ABI5, SAG13 and SAG113 [67]. And the overexpression of TaMYB3R1 increased the expression of RD29A, RD29B and ABF3, which were also involved in the ABA signal pathway [68].

4.4 miR156 inhibited metabolism and development in leaf but promoted it in root under drought in the sensitive line Lv28

The modules of miR156 and SBL TFs have been shown to regulate many important traits in plants. For example, the OsmiR156 and OsSPL14 module controlled the ideal architecture in rice, and mutation of the OsmiR156 target site in OsSPL14 reduced the tiller number and enhanced the grain yield [69]. The SPL TFs have also been shown to regulate grain size, grain quality, panicle branch, and plant height in rice [70, 71]. In maize, miR156 could regulate two SPL TFs, tga1 and tsh4. tga1 (teosinte glume architecture1) has been determined to be a key gene in the domestication of maize from teosinte [72], and the module of miR156 and tga1 has been shown to regulate the phase transition, tillering and inflorescence architecture [73]. The SBL TF gene tsh4 (tasselsheath4) has been shown to play an essential role in meristem boundary establishment [74]. In Arabidopsis, miR156 was up-regulated under drought, and its overexpression resulted in a significantly higher survival rate of transgenic plants under drought stress [75]. Moreover, the miR160-ARF module also plays important role in development in plants. For example, in soybean, overexpression of miR160 in roots results in the hypersensitivity to auxin and cytokinin, resulting in an inhibition of development of the symbiotic nodule [76]. It is interesting that, in our study, the expression of miR156 was induced in leaf bur decreased in root under drought in the drought-sensitive line Lv28, which caused an own-regulation of targeted SPL genes in leaf and up-regulation of the genes in root under drought (Fig 4). Additionally, the putative down-stream genes regulated by miR156-SPL modules of Lv28, which showed co-expression with miR156-SPL genes and had squamosal binding sequences in their promoter region, were enriched in photosynthetic carbon metabolism in leaves, while that genes were enriched in stress response function in roots. These results suggested that inhibition of development in shoots and enhancement of development or activation of stress response signal transduction in roots modulated by the miR156-SPL module are strategies for drought stress adaption in maize.

4.5 miRNAs with potential in drought tolerance improvement

Some drought responsive miRNAs have shown potentials in drought tolerance improvement. For instance, the deletion of miR169a by CRISPR/Cas9 could enhance drought tolerance in Arabidopsis, whose transgenic plants showed a high survival rate under drought [77]. The overexpression of OsmiR393 down-regulated OsTIR1 and OsAFB2 and led to higher sensitivity to drought and salt than the wild type in rice [78]. miR398 has been shown responsive to drought in Arabidopsis, rice, pea (Pisum sativum L.), and Medicago truncatula [7981]. miR398 was also considered a key miRNA involved in drought tolerance in rice, as it has showed opposite responsive patterns between the tolerant and sensitive cultivars under drought [19]. miR398 conservatively targeted CSD1 and CSD2 encoding copper/zinc superoxide dismutases (Cu/Zn-SODs), COX5b encoding a subunit of mitochondrial cytochrome c oxidase and CCS1 encoding the copper chaperone for SOD [82]. The interaction of miR398 and its target genes was also involved in responses to other abiotic stresses, such as heat, salt and oxidative stress [8385]. In the present study, miR398b was specifically down-regulated in the leaves and roots of H082183 under drought, implying that the miRNA might also be used in stress tolerance improvement. In addition, the miR160-ARFs module could regulate leaf development via the auxin signaling pathways, and miR165/166-HD-ZIP IIIs module could confer drought tolerance through the ABA signaling pathways in Arabidopsis [86]. The double mutant plants of miR160 and miR165/166 showed defects in drought tolerance and leaf development [86]. Our results revealed that the miR160-ARF module was responsive to drought in the leaves of sensitive line Lv28. In Arabidopsis, the double mutant of miR159a/b had larger and longer roots than the wild type, and the expressions of target genes MYB33, MYB65 and MYB101 were increased [87]. In the present study, the miR159-MYB modules showed responsive to drought in roots of H082183. Thus, detailed studies on miR159 are also needed to explore the potential in the improvement of drought tolerance.

5 Conclusion

In conclusion, the drought-associated miRNome of maize were analyzed by examining the drought-tolerant and drought-sensitive inbred lines grown in the field by deep sequencing. A total of 53 and 52 drought responsive miRNAs were identified in leaves and roots, respectively. Particularly, several miRNA-mRNA target modules showed uniquely responsive to drought in the two lines with contrasting drought tolerance. The miR164-MYB and miR164-NAC modules, regulated by ABA, were drought responsive specifically in H082183. The miR156-SPL module was uniquely responsive to drought in Lv28, with opposite regulation patterns in leaves and roots. In addition, our study detected three novel miRNAs with high confidence, which merit further validation and biological function studies. This study expands our knowledge of drought-induced molecular regulatory mechanisms modulated by miRNAs in maize and extends the list of candidate genes for molecular breeding for drought tolerance.

Supporting information

S1 File.

Additional File (.zip) included Figs A-J and Tables A-J. Fig A. Percentage of expressed miRNAs in each family. (A) The expressed miRNAs of different families in leaves. (B) The expressed miRNAs of different families in roots. Fig B. Distribution of length of known miRNAs that expressed in the present study. Fig C. The conservation analysis of predicted novel miRNAs. The mature sequence of predicted novel miRNAs in the present study and known miRNAs in plants were aligned. (A) Novel_1_44595. (B) Novel_1_47179. (C) Novel_3_27520. (D) Novel_3_29234. (E) Novel_9_5030. (F) Novel_1-_3438. zma: Zea mays. osa: Oryza sativa. bdi: Brachypodium distachyon. tae: Triticum aestivum. mtr: Medicago truncatula. ath: Arabidopsis thaliana. ptc: Populus trichocarpa. vvi: Vitis vinifera. rco: Ricinus communis. Fig D. The phylogenetic tree and sequence of three novel miRNAs. The phylogenetic trees of Novel-2-38 (A), Novel-10-88 (B) and Novel-10-147 (C) were shown. The sequence of relative miRNAs in different species of Novel-2-38 (D), Novel-10-88 (E) and Novel-10-147 (F). ath: Arabidopsis thaliana, osa: Oryza sativa, sbi: Sorghum bicolor, zma: Zea mays. (G) The novel mature sequence of miRNAs were cloned using stem-loop methods in B73 (B), H082183 (H) and Lv28 (L). The maker (M) was 100-bp. Fig E. Correlation result of miRNome in each replicate. The TPM value of all known miRNAs in each replicate was calculated as log10(TPM+1) for the correlation analysis. The correlation coefficient r and the P value are labelled in each replicate sample. H and L indicate the two maize inbred lines H082183 and Lv28, respectively. MD and SD indicate moderate and severe drought, respectively. MC and SC indicate well-watered controls of moderate and severe drought, respectively. 1 and 2 indicate the two replicates. Fig F. The T-plot of drought responsive miRNAs and their target genes. (A) GRMZM2G065451, (B) GRMZM2G113779, (C) GRMZM2G414805, (D) GRMZM2G371033, (E) AC207656.3_FGT002, (F) GRMZM2G390641, (G) GRMZM2G159399, (H) GRMZM2G096358, (I) GRMZM2G114850, (J) GRMZM2G113779, (K) GRMZM5G806833, (L) GRMZM2G061734, (M) GRMZM2G004090, (N) GRMZM2G028054, (O) GRMZM2G304745, (P) GRMZM2G103812. Fig G. The cluster analysis of genes. (A) GRMZM2G096358 (MYB) and GRMZM2G114850 (NAC) in leaves of H082183. (B) GRMZM2G390641 (ARF) in leaves of Lv28. (C) GRMZM2G159399 (ARF) in leaves of Lv28. (D) AC207656.3_FG002 (ARF), GRMZM2G065451 (SPL) and GRMZM2G113779 (SPL) in leaves of Lv28. (E) GRMZM2G067624 (SPL) and GRMZM2G371033 (SPL) in leaves of Lv28. (F) GRMZM2G004090 (MYB) in roots of H082183. (G) GRMZM2G028054 (MYB) in roots of H082183. (H) GRMZM2G061734 (SPL) and GRMZM2G160917 (SPL) in roots of Lv28. (I) GRMZM2G460544 (SPL) and GRMZM5G806833 (SPL) in roots of Lv28. Fig H. Co-expressed genes of miR156 targeted MYB and NAC genes that had over-presented transcription factor binding sequence in promoter regions. The co-expressed genes of NAC and MYB genes that had significant over-presented MYB (A) or NAC (B) binding motif sequence in promoter regions were listed. The MYB (A) and NAC (B) binding sequences in promoter regions were plotted. Fig I. The expression of miR159 and target genes in roots under ABA treatment. (A) miR159. (B) GRMZM2G004090 (MYB). (C) GRMZM2G028054 (MYB). Fig J. Correlation result of miRNA expression by qRT-PCR and sequencing. Five miRNAs in leaves and five miRNAs in roots were selected to perform qRT-PCR validation under moderate drought and severe drought in H082183 and Lv28. Table A. Primers used in the present study. Table B. Sequencing data of all 32 small RNA libraries. a H: H082183; L: Lv28; MD: moderate drought; SD: severe drought; MC: well-watered control of moderate drought; SC: well-watered control of severe drought; 1 and 2: two replicates. Table C. Result of reference genome alignment of small RNA reads. a H: H082183; L: Lv28; MD: moderate drought; SD: severe drought; MC: well-watered control of moderate drought; SC: well-watered control of severe drought; 1 and 2: two replicates. Table D. Alignment result of different small RNA databases. a H: H082183; L: Lv28; MD: moderate drought; SD: severe drought; MC: well-watered control of moderate drought; SC: well-watered control of severe drought; 1 and 2: two replicates. Table E. Numbers of expressed miRNAs in each sample. a H: H082183; L: Lv28; MD: moderate drought; SD: severe drought; MC: well-watered control of moderate drought; SC: well-watered control of severe drought; 1 and 2: two replicates. Table F. The data output of degradome sequencing. Table G. Predicted novel miRNAs. a The number of libraries that had both miRNA and miRNA* reads. Table H. Expression of all miRNA-mRNA pairs in leaves and roots. Table I. Drought responsive miRNAs in different species. a miR156e-3p was down-regulated under moderate drought, while miR156b-5p and miR156d-5p were up-regulated under severe drought in Lv28. b miR164h-5p was down-regulated under moderate drought in H082183, but miR164c-3p and miR164h-3p were up-regulated under both moderate and severe drought in H082183. c miR171g-3p was down-regulated under moderate drought, but miR171h-3p and miR171k-3p were up-regulated under severe drought in Lv28. d miR156d-3p, miR156f-3p and miR156g-3p were down-regulated, while miR156e-3p was up-regulated under severe drought in H082183. e miR159a-3p, miR159b-3p, miR159f-3p, miR159j-3p and miR159k-3p were up-regulated under moderate drought in H082183, but miR159c-3p, miR159c-5p, miR159d-3p and miR159d-5p were down-regulated under moderate drought in Lv28. f miR167a-5p, miR167b-3p, miR167b-5p, miR167c-5p and miR167d-5p were up-regulated under severe drought in H082183, while miR167j-3p was down-regulated under moderate drought in Lv28. Table J. Drought responsive miRNAs by using different drought treatment methods in maize. a This line present the maize inbred lines used in experiment, in which 81565, 87–1, HKI-1532 and H082183 were drought tolerant lines, and 21ES, Dan340, V-372 and Lv28 were drought sensitive lines. b This line present the sampling time point after drought. The h and d present hours and days after drought treatment. c miR171h-3p and miR171k-3p were up-regulated, while miR171g-3p was down-regulated under drought. d miR156e-3p was up-regulated, while miR156d-3p, miR156f-3p and miR156g-3p were down-regulated under drought.

https://doi.org/10.1371/journal.pone.0219176.s001

(ZIP)

Acknowledgments

Thanks go to Prof. Jiuran Zhao from Maize Research Center of Beijing Academy of Agricultural and Forestry Sciences for providing seeds of H082183.

References

  1. 1. Brown ME, Funk CC. Climate. Food security under climate change. Science. 2008;319(5863):580–581. pmid:18239116
  2. 2. Takeda S, Matsuoka M. Genetic approaches to crop improvement: responding to environmental and population changes. Nat Rev Genet. 2008;9(6):444–457. pmid:18475268
  3. 3. 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. pmid:19965430
  4. 4. Lobell DB, Roberts MJ, Schlenker W, Braun N, Little BB, Rejesus RM, et al. Greater sensitivity to drought accompanies maize yield increase in the U.S. Midwest. Science. 2014;344(6183):516–519. pmid:24786079
  5. 5. Tang J, Chu C. MicroRNAs in crop improvement: fine-tuners for complex traits. Nat Plants. 2017;3:17077. pmid:28665396
  6. 6. Aydinoglu F, Lucas SJ. Identification and expression profiles of putative leaf growth related microRNAs in maize (Zea mays L.) hybrid ADA313. Gene. 2019;690:57–67. pmid:30597233
  7. 7. Bonnet E, Van de Peer Y, Rouze P. The small RNA world of plants. New Phytol. 2006;171(3):451–468. pmid:16866953
  8. 8. Bo X, Wang S. TargetFinder: a software for antisense oligonucleotide target site selection based on MAST and secondary structures of target mRNA. Bioinformatics. 2005;21(8):1401–1402. pmid:15598838
  9. 9. 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. pmid:19017659
  10. 10. Gao S, Yang L, Zeng HQ, Zhou ZS, Yang ZM, Li H, et al. A cotton miRNA is involved in regulation of plant response to salt stress. Sci Rep. 2016;6:19736. pmid:26813144
  11. 11. Zhao B, Liang R, Ge L, Li W, Xiao H, Lin H, et al. Identification of drought-induced microRNAs in rice. Biochem Biophys Res Commun. 2007;354(2):585–590. pmid:17254555
  12. 12. Ding D, Zhang L, Wang H, Liu Z, Zhang Z, Zheng Y. Differential expression of miRNAs in response to salt stress in maize roots. Ann Bot. 2009;103(1):29–38. pmid:18952624
  13. 13. Wei L, Zhang D, Xiang F, Zhang Z. Differentially expressed miRNAs potentially involved in the regulation of defense mechanism to drought stress in maize seedlings. Int J Plant Sci. 2009;170(8):979–989.
  14. 14. 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(15):4157–4168. pmid:20729483
  15. 15. Li J, Fu F, An M, Zhou S, She Y, Li W. Differential expression of microRNAs in response to drought stress in maize. J Integr Agri. 2013;12(8):1414–1422.
  16. 16. Li H, Dong Y, Yin H, Wang N, Yang J, Liu X, et al. Characterization of the stress associated microRNAs in Glycine max by deep sequencing. BMC Plant Biol. 2011;11:170. pmid:22112171
  17. 17. Barrera-Figueroa BE, Gao L, Wu Z, Zhou X, Zhu J, Jin H, et al. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice. BMC Plant Biol. 2012;12:132. pmid:22862743
  18. 18. Yin F, Gao J, Liu M, Qin C, Zhang W, Yang A, et al. Genome-wide analysis of water-stress-responsive microRNA expression profile in tobacco roots. Funct Integr Genomics. 2014;14(2):319–332. pmid:24664280
  19. 19. Cheah BH, Nadarajah K, Divate MD, Wickneswari R. Identification of four functionally important microRNA families with contrasting differential expression profiles between drought-tolerant and susceptible rice leaf at vegetative stage. BMC Genomics. 2015;16(1):692.
  20. 20. Hackenberg M, Gustafson P, Langridge P, Shi BJ. Differential expression of microRNAs and other small RNAs in barley between water and drought conditions. Plant Biotechnol J. 2015;13(1):2–13. pmid:24975557
  21. 21. Liu HH, Tian X, Li YJ, Wu CA, Zheng CC. Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA. 2008;14(5):836–843. pmid:18356539
  22. 22. Akdogan G, Tufekci ED, Uranbey S, Unver T. miRNA-based drought regulation in wheat. Funct Integr Genomics. 2016;16(3):221–233. pmid:26141043
  23. 23. Sunkar R, Zhu JK. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell. 2004;16(8):2001–2019. pmid:15258262
  24. 24. Zhang F, Luo X, Zhou Y, Xie J. Genome-wide identification of conserved microRNA and their response to drought stress in Dongxiang wild rice (Oryza rufipogon Griff.). Biotechnol Lett. 2016;38(4):711–721. pmid:26667133
  25. 25. Zhang L, Chia JM, Kumari S, Stein JC, Liu Z, Narechania A, et al. A genome-wide characterization of microRNA genes in maize. PLoS Genet. 2009;5(11):e1000716. pmid:19936050
  26. 26. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42(Database issue):D68–73. pmid:24275495
  27. 27. Ferdous J, Hussain SS, Shi BJ. Role of microRNAs in plant drought tolerance. Plant Biotechnol J. 2015;13(3):293–305. pmid:25583362
  28. 28. Wang YG, An M, Zhou SF, She YH, Li WC, Fu FL. Expression profile of maize microRNAs corresponding to their target genes under drought stress. Biochem Genet. 2014;52(11–12):474–493. pmid:25027834
  29. 29. Luan M, Xu M, Lu Y, Zhang L, Fan Y, Wang L. Expression of zma-miR169 miRNAs and their target ZmNF-YA genes in response to abiotic stress in maize leaves. Gene. 2015;555(2):178–185. pmid:25445264
  30. 30. Luan M, Xu M, Lu Y, Zhang Q, Zhang L, Zhang C, et al. Family-wide survey of miR169s and NF-YAs and their expression profiles response to abiotic stress in maize roots. PLoS ONE. 2014;9(3):e91369. pmid:24633051
  31. 31. Jones-Rhoades MW, Bartel DP, Bartel B. MicroRNAs and their regulatory roles in plants. Annu Rev Plant Biol. 2006;57:19–53. pmid:16669754
  32. 32. Blum A. Genomics for drought resistance—getting down to earth. Funct Plant Biol. 2014;41(11):1191–1198.
  33. 33. Li C, Li Y, Sun B, Peng B, Liu C, Liu Z, et al. Quantitative trait loci mapping for yield components and kernel-related traits in multiple connected RIL populations in maize. Euphytica. 2013(193):303–316.
  34. 34. Zhang X, Liu X, Zhang D, Tang H, Sun B, Li C, et al. Genome-wide identification of gene expression in contrasting maize inbred lines under field drought conditions reveals the significance of transcription factors in drought tolerance. PLoS ONE. 2017;12(7):e0179477. pmid:28700592
  35. 35. Barrs HD, Weatherley PE. A re-examination of the relative turgidity technique for estimating water deficit in leaves. Aust J Biol Sci. 1968(15):413–428.
  36. 36. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41(Database issue):D590–596. pmid:23193283
  37. 37. Chan PP, Lowe TM. GtRNAdb: a database of transfer RNA genes detected in genomic sequence. Nucleic Acids Res. 2009;37(Database issue):D93–97. pmid:18984615
  38. 38. Nawrocki EP, Burge SW, Bateman A, Daub J, Eberhardt RY, Eddy SR, et al. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 2015;43(Database issue):D130–137. pmid:25392425
  39. 39. Bao W, Kojima KK, Kohany O. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 2015;6:11. pmid:26045719
  40. 40. Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics. 2010; Chapter 11:Unit-11.7.
  41. 41. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52. pmid:21911355
  42. 42. Axtell MJ, Meyers BC. Revisiting criteria for plant microRNA annotation in the era of big data. Plant Cell. 2018;30(2):272–284. pmid:29343505
  43. 43. Lorenz R, Bernhart S, Siederdissen C, Tafer H, Flamm C, Stadler P, et al. ViennaRNA Package 2.0. Algorithm Mol Biol. 2011;6:26.
  44. 44. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetic analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28(10):2731–2739. pmid:21546353
  45. 45. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of miRNA genes. PloS ONE. 2007;2(2):e219. pmid:17299599
  46. 46. Anders S, McCarthy DJ, Chen YS, Okoniewski M, Smyth GK, Huber W, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013;8(9):1765–1786. pmid:23975260
  47. 47. Bateman A, Martin MJ, O'Donovan C, Magrane M, Alpi E, Antunes R, et al. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2017;45(Database issue):D158–169. pmid:27899622
  48. 48. Pruitt KD, Tatusova T, Maglott DR. NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 2005;33(Database issue):D501–504. pmid:15608248
  49. 49. Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(Database issue):D279–285. pmid:26673716
  50. 50. Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27(7):1017–1018. pmid:21330290
  51. 51. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(Web Server issue):W64–70. pmid:20435677
  52. 52. 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. pmid:20407023
  53. 53. Shinozaki K, Yamaguchi-Shinozaki K, Seki M. Regulatory network of gene expression in the drought and cold stress responses. Curr Opin Plant Biol. 2003;6(5):410–417. pmid:12972040
  54. 54. Wang H, Wang H, Shao H, Tang X. Recent advances in utilizing transcription factors to improve plant abiotic stress tolerance by transgenic technology. Front Plant Sci. 2016;7:67. pmid:26904044
  55. 55. Aravind J, Rinku S, Pooja B, Shikha M, Kaliyugam S, Mallikarjuna MG, et al. Identification, characterization, and functional validation of drought-responsive microRNAs in subtropical maize inbreds. Front Plant Sci. 2017;8:941. pmid:28626466
  56. 56. Shinozaki K, Yamaguchi-Shinozaki K. Gene networks involved in drought stress response and tolerance. J Exp Bot. 2007;58(2):221–227. pmid:17075077
  57. 57. Ding Y, Tao Y, Zhu C. Emerging roles of microRNAs in the mediation of drought stress response in plants. J Exp Bot. 2013;64(11):3077–3086. pmid:23814278
  58. 58. 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. pmid:24604734
  59. 59. Li J, Guo G, Guo W, Guo G, Tong D, Ni Z, et al. miRNA164-directed cleavage of ZmNAC1 confers lateral root development in maize (Zea mays L.). BMC Plant Biol. 2012;12:220. pmid:23171309
  60. 60. Feng H, Duan X, Zhang Q, Li X, Wang B, Huang L, et al. The target gene of tae-miR164, a novel NAC transcription factor from the NAM subfamily, negatively regulates resistance of wheat to stripe rust. Mol Plant Pathol. 2014;15(3):284–296. pmid:24128392
  61. 61. Lu M, Zhang DF, Shi YS, Song YC, Wang TY, Li Y. Expression of SbSNAC1, a NAC transcription factor from sorghum, confers drought tolerance to transgenic Arabidopsis. Plant Cell Tiss Org. 2013;115(3):443–455.
  62. 62. Lu M, Ying S, Zhang D, Shi Y, Song Y, Wang T, et al. A maize stress-responsive NAC transcription factor, ZmSNAC1, confers enhanced tolerance to dehydration in transgenic Arabidopsis. Plant Cell Rep. 2012;31:1701–1711. pmid:22610487
  63. 63. Mao H, Wang H, Liu S, Li Z, Yang X, Yan J, et al. A transposable element in a NAC gene is associated with drought tolerance in maize seedlings. Nat Commun. 2015;6:8326. pmid:26387805
  64. 64. Dai X, Xu Y, Ma Q, Xu W, Wang T, Xue Y, et al. Overexpression of an R1R2R3 MYB gene, OsMYB3R-2, increases tolerance to freezing, drought, and salt stress in transgenic Arabidopsis. Plant Physiol. 2007;143(4):1739–1751. pmid:17293435
  65. 65. Li MJ, Qiao Y, Li YQ, Shi ZL, Zhang N, Bi CL, et al. A R2R3-MYB transcription factor gene in common wheat (namely TaMYBsm1) involved in enhancement of drought tolerance in transgenic Arabidopsis. J Plant Res. 2016;129(6):1097–1107. pmid:27542160
  66. 66. Zhao Y, Gao J, Im Kim J, Chen K, Bressan RA, Zhu JK. Control of plant water use by ABA induction of senescence and dormancy: an overlooked lesson from evolution. Plant Cell Physiol. 2017;58(8):1319–1327. pmid:28961993
  67. 67. Huang Q, Wang Y, Li B, Chang J, Chen M, Li K, et al. TaNAC29, a NAC transcription factor from wheat, enhances salt and drought tolerance in transgenic Arabidopsis. BMC Plant Biol. 2015;15:268. pmid:26536863
  68. 68. Cai H, Tian S, Dong H, Guo C. Pleiotropic effects of TaMYB3R1 on plant development and response to osmotic stress in transgenic Arabidopsis. Gene. 2015;558(2):227–234. pmid:25560188
  69. 69. Jiao Y, Wang Y, Xue D, Wang J, Yan M, Liu G, et al. Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat Genet. 2010;42(6):541–544. pmid:20495565
  70. 70. Miura K, Ikeda M, Matsubara A, Song XJ, Ito M, Asano K, et al. OsSPL14 promotes panicle branching and higher grain productivity in rice. Nat Genet. 2010;42(6):545–549. pmid:20495564
  71. 71. Si L, Chen J, Huang X, Gong H, Luo J, Hou Q, et al. OsSPL13 controls grain size in cultivated rice. Nat Genet. 2016;48(4):447–456. pmid:26950093
  72. 72. Wang H, Nussbaum-Wagler T, Li B, Zhao Q, Vigouroux Y, Faller M, et al. The origin of the naked grains of maize. Nature. 2005;436(7051):714–719. pmid:16079849
  73. 73. Chuck G, Meeley R, Irish E, Sakai H, Hake S. The maize tasselseed4 microRNA controls sex determination and meristem cell fate by targeting Tasselseed6/indeterminate spikelet1. Nat Genet. 2007;39(12):1517–1521. pmid:18026103
  74. 74. Chuck G, Whipple C, Jackson D, Hake S. The maize SBP-box transcription factor encoded by tasselsheath4 regulates bract development and the establishment of meristem boundaries. Development. 2010;137(8):1243–1250. pmid:20223762
  75. 75. Cui LG, Shan JX, Shi M, Gao JP, Lin HX. The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants. Plant J. 2014;80(6):1108–1117. pmid:25345491
  76. 76. Turner M, Nizampatnam NR, Baron M, Coppin S, Damodaran S, Adhikari S, et al. Ectopic expression of miR160 results in auxin hypersensitivity, cytokinin hyposensitivity, and inhibition of symbiotic nodule development in soybean. Plant Physiol. 2013;162(4):2042–2055. pmid:23796794
  77. 77. Zhao Y, Zhang C, Liu W, Gao W, Liu C, Song G, et al. An alternative strategy for targeted gene replacement in plants using a dual-sgRNA/Cas9 design. Sci Rep. 2016;6:23890. pmid:27033976
  78. 78. Xia K, Wang R, Ou X, Fang Z, Tian C, Duan J, et al. OsTIR1 and OsAFB2 downregulation via OsmiR393 overexpression leads to more tillers, early flowering and less tolerance to salt and drought in rice. PLoS ONE. 2012;7(1):e30039. pmid:22253868
  79. 79. Jagadeeswaran G, Saini A, Sunkar R. Biotic and abiotic stress down-regulate miR398 expression in Arabidopsis. Planta. 2009;229(4):1009–1014. pmid:19148671
  80. 80. Trindade I, Capitao C, Dalmay T, Fevereiro MP, Santos DM. miR398 and miR408 are up-regulated in response to water deficit in Medicago truncatula. Planta. 2010;231(3):705–716. pmid:20012085
  81. 81. Jovanovic Z, Stanisavljevic N, Mikic A, Radovic S, Maksimovic V. Water deficit down-regulates miR398 and miR408 in pea (Pisum sativum L.). Plant Physiol Biochem. 2014;83:26–31. pmid:25064597
  82. 82. Zhu C, Ding Y, Liu H. MiR398 and plant stress responses. Physiol Plant. 2011;143(1):1–9. pmid:21496029
  83. 83. Jia X, Wang WX, Ren L, Chen QJ, Mendu V, Willcut B, et al. Differential and dynamic regulation of miR398 in response to ABA and salt stress in Populus tremula and Arabidopsis thaliana. Plant Mol Biol. 2009;71(1–2):51–59. pmid:19533381
  84. 84. Lu X, Guan Q, Zhu J. Downregulation of CSD2 by a heat-inducible miR398 is required for thermotolerance in Arabidopsis. Plant Signal Behav. 2013;8(8):e24952. pmid:23733060
  85. 85. Leng X, Wang P, Zhu X, Li X, Zheng T, Shangguan L, et al. Ectopic expression of CSD1 and CSD2 targeting genes of miR398 in grapevine is associated with oxidative stress tolerance. Funct Integr Genomics. 2017;17(6):697–710. pmid:28674744
  86. 86. Yang TX, Wang YY, Teotia S, Wang ZH, Shi CN, Sun H, et al. The interaction between miR160 and miR165/166 in the control of leaf development and drought tolerance in Arabidopsis. Sci Rep. 2019;9(1):2832. pmid:30808969
  87. 87. Xue T, Liu Z, Dai X, Xiang F. Primary root growth in Arabidopsis thaliana is inhibited by the miR159 mediated repression of MYB33, MYB65 and MYB101. Plant Sci. 2017;262:182–189. pmid:28716415