Unraveling the microRNA of Caragana korshinskii along a precipitation gradient on the Loess Plateau, China, using high-throughput sequencing

Drought remains one of the main factors that negatively affect plant growth and development. Caragana korshinskii is widely distributed on the Loess Plateau, China, where it mediates soil and water loss and helps prevent desertification. However, little is known about the stress response mechanisms of C. korshinskii in water-starved environments. MicroRNAs (miRNAs) have been implicated in the regulation of plant responses to several types of biotic and abiotic stress. Here, we describe the miRNAs of wild C. korshinskii from Huangling, Yulin, and Dalad Banner, which occur along a precipitation gradient. Using next-generation sequencing technology, we obtained a total of 13 710 681, 15 048 945, and 15 198 442 reads for each location, respectively; after filtering and BLAST analysis, 490 conserved miRNAs and 96 novel miRNAs were characterized from the sRNAome data, with key functions determined using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses. We also designed stem-loop qRT-PCR to validate the expression patterns of 5 conserved miRNAs (miR390, miR398, miR530, miR2119, and miR5559) that obviously responded to water stress in plants grown both under natural and experimental water stress conditions and found that the expression levels of miR2119 and miR5559 were negatively correlated with their predicted target genes. This study is the first to identify miRNAs from wild C. korshinskii and provides a basis for future studies of miRNA-mediated gene regulation of stress responses in C. korshinskii.


Introduction
MicroRNAs (miRNAs) are small, non-coding RNAs of~22 nucleotides (nt) in length that play pivotal regulatory roles in both animals and plants [1]. In fact, miRNAs are indispensable players in many types of gene regulation and function; they act by binding to complementary messenger RNA (mRNA) from specific target genes, thus causing mRNA degradation and gene silencing [2]. In addition, miRNAs constitute the majority of gene regulatory molecules a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 involved in regulatory pathways and influence the output of a variety of protein-coding genes. A large amount of evidence also indicates that miRNAs are involved in various aspects of plant growth and development, such as leaf morphogenesis [3], phase transition [4], and adaptive responses to abiotic and biotic stress [5,6].
Among the many types of environmental stress, drought remains one of the most important factors limiting plant growth and development. For example, the Loess Plateau, which is located in the semi-arid regions of northwest China, is a fragile ecosystem that is affected by drought, as well as by other natural factors. Global climate change has recently heightened drought in this region, and both grassland degradation and desertification now threaten to destroy the ecosystem [7,8]). Therefore, investigating drought-resistant plants of the Loess Plateau may reveal effective ways to prevent the spread of desertification, which can then be used to facilitate the long-term protection of the ecosystem's balance and stability.
Caragana korshinskii (Fabaceae), widely known as a "life-saving plant," is a perennial leguminous shrub that is distributed in the arid and semi-arid regions of northwest China [9,10]. It is a high-quality forage plant, perfectly adapted to drought stress. It is also a pioneer plant and functions to control soil erosion and desertification of the grasslands of the Loess Plateau. Investigating the molecular mechanisms of the plant's drought resistance is still an attractive and popular target for research, owing to its potential uses in agricultural production, vegetation rehabilitation, and ecological restoration [11]. However, few studies have investigated the contribution of miRNA to the adaptive mechanisms of C. korshinskii.
In recent years, the HiSeq high-throughput sequencing technology has been widely applied to investigating novel species-specific or low-abundance miRNAs [12,13]. The technology has been applied to multiple model plants such as rice [12] and tobacco [14], as well as non-model plants, such as Chinese yew [15], peanut [16], and cowpea [17].
We used high-throughput sequencing to analyze the miRNA expression profile of wild C. korshinskii along a precipitation gradient on the Loess Plateau, China. In addition, we aimed to determine the significance of known and novel miRNAs from their expression profiles and to predict their target genes. The present study presents significant evidence for the role of miRNAs in the drought-tolerance of C. korshinskii. Further characterization of the targets of drought-associated miRNAs will help elucidate the response and tolerance of C. korshinskii to drought.

Sampling sites, conditions, and plant material
The Loess Plateau is located in the continental monsoon region of China, where the annual precipitation forms a natural gradient from 150 mm in the northwest to 800 mm in the southeast; the majority of the rainfall is concentrated between June and September [18], and the annual evaporation ranges from high in the northwest to low in the southeast [19,20]. Three C. korshinskii experimental sites, Huangling (HuL; 35˚39'N,109˚14'E), Yulin (YuL; 38˚19'N, 1094 3'E) and Dalad Banner (DaB; 40˚14'N, 109˚59'E), were selected along the precipitation gradient across Shaanxi Province and Inner Mongolia, Northwest China (Fig 1). As shown in our previous study of C. korshinskii [21], the annual precipitation, annual evaporation, soil moisture content, and plant leaf water potential of the sites differ significantly, and the positive correlations that leaf water potential shows with annual precipitation (R = 0.968, P < 0.01) and soil moisture content (R = 0.963, P < 0.01) indicate that water stress increases gradually from the south to the north. Moreover, the annual average temperature from the northwest to the southeast, and the soil pH and altitude along the precipitation gradient are known to be similar [22].
Three circular plots (20 m radius) were established at each site, and both plant and soil samples were collected from each site in July 2013. All of the sampled C. korshinskii shrubs were over 20 years old. Leaf samples (~5 g) were collected, and their water potential was determined using a WP4-T dew-point meter (Decagon Devices, Pullman, WA, USA). Upper mature leaflets with healthy extensions were immediately frozen in liquid nitrogen, transported to the laboratory, and stored at -80˚C prior to enzyme activity analysis and total RNA isolation.
We state that no specific permissions were required for these locations/activities. We confirm that the field studies did not involve endangered or protected species.

Small RNA library construction and sequencing
To construct the small RNA (sRNA) library, all three samples were pooled, and total RNA was extracted from the mixture according to the manufacturer's instructions. About 10 mg of sRNA was used for sequencing. The sRNA fractions with lengths of 10-40 nt were isolated using 15% denaturing polyacrylamide gel electrophoresis, ligated to 5 0 and 3 0 adaptors, and reverse-transcribed to cDNA, according to the Illumina protocol. The resulting cDNA libraries were then sequenced by the GENE DENOVO Company (China).

miRNA identification
The raw sequences were first processed using Illumina's Genome Analyzer Pipeline software in order to filter out the adapter sequences, as well as low quality and low-copy sequences. Then, RNA sequences of 18-30 nt in length were mapped to transcriptome data (https://www. ncbi.nlm.nih.gov/sra/?term=SRX1117100), RFam (http://rfam.xfam.org/), and Repbase (http://www.girinst.org/repbase/), in order to identify and filter out mRNA, rRNA, tRNA, snRNA, snoRNA, and repeat sequences. Finally, the remaining unique sequences were compared to the database of published miRNA sequences, miRBase 21 (http://www.mirbase.org/) using BLASTn search to identify conserved miRNAs. A maximum of three mismatches was allowed between the miRNAs of C. korshinskii and those from other plant species, and only perfectly matched sequences were considered conserved [23]. small RNA tags that were not identified from the public databases were aligned with the reference genome. According to their genome positions and hairpin structures predicted by the software MIREAP (version 0.2, http://sourceforge.net/projects/mireap/;) [24], the novel miRNA candidates were identified [25,26].

Target prediction and functional annotation
The target genes of the known miRNAs were predicted as those occurring at the intersection of the MIREAP, miRanda, and TargetScan results. The criteria used for target prediction were based on those suggested by Allen et al. [25] and Schwab et al. [26]: (1) no more than 4 mismatches between the sRNA and the target (G-U bases count as 0.5 mismatches); (2) no more than two adjacent mismatches in the miRNA/target duplex; (3) no adjacent mismatches in positions 2-12 of the miRNA/target duplex (5' of miRNA); (4) no mismatches in positions 10-11 of the miRNA/target duplex; (5) no more than 2.5 mismatches in positions 1-12 of the miRNA/target duplex (5' of miRNA); and (6) minimum free energy of the miRNA/target duplex should be !75% of the minimum free energy of the miRNA bound to its perfect complement.
In order to better understand miRNA target function and classification, as well as the metabolic regulatory networks associated with C. korshinskii miRNAs and their targets, Gene Ontology (GO; http://www.geneontology.org/) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.genome.jp/keg/) were performed for the target genes, as described previously [27].

Verification of miRNAs by quantitative real-time PCR (qRT-PCR)
Stem-loop qRT-PCR was carried out to validate the presence and expression of the identified miRNAs. Total RNA was extracted from leaves using Trizol (TaKaRa, Dalian, China), according to the manufacturer's instructions, and the RNA concentrations and integrity were analyzed using a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis. Finally, the RNA was treated with RNase-free DNase I (TaKaRa) to remove genomic DNA. Sequence-specific forward primers were designed for 9 selected miRNAs (S1 Table), and qRT-PCR was performed using SYBR RealMasterMix (Tiangen Biotech, Beijing, China) and the ABI 7300 real-time PCR detection system (Applied Biosystems, USA), with the following amplification conditions: denaturation at 95˚C for 3 min followed by 40 cycles of denaturation at 95˚C for 15 s and annealing and extension at 61˚C for 30 s. Finally, melting analyses were performed to confirm the absence of false-positive peaks. All reactions were performed in triplicate for each sample, and U6 snRNA was used as an internal reference [28]. The relative expression levels of the miRNAs were calculated using the 2 -ΔΔCt method [29].
Drought treatment of C. korshinskii in laboratory C. korshinskii seeds were obtained from the Loess Plateau and allowed to germinate at 28˚C in the dark. After 2 d, the seeds were individually sown in pots (11 × 35 cm) filled with a 1:1 (v:v) mixture of rich soil and vermiculite and then kept in a growth chamber with a 16-h light/8-h dark photoperiod and 25/18 (±2)˚C day/night temperature cycle [30]. After 3 months, when the seedlings had grown to~25 cm, the plants were assigned to one of three soil moisture levels (75%, 55%, 35%), which corresponded to suitable water availability (CK), mild drought stress (MD), and severe drought stress (SD), respectively [31]. The soil moisture levels were then maintained constant for 20 d, as described previously [32]. After the drought treatment, leaf samples were collected, frozen in liquid nitrogen, and stored at -80˚C. Three replicates were performed for each treatment.

Statistical tests
All data were analyzed using one-way ANOVA in SPSS 17.0 (SPSS Inc., Chicago, IL, USA), with differences of P < 0.05 considered statistically significant. Experimental data are reported as mean values ± SE (n = 3).
The majority of the small RNAs were 20-24 nt in length, with the 24-nt class being the most abundant (64.7%), followed by the 21-nt class (16.35%); many of the small RNAs were annotated as noncoding RNAs (Fig 2). We easily notice all kinds of small RNA, in which the unannotated miRNA accounts for a proportion of the unique RNA, approximately 80%, miRNA accounts for proportion just 5-7% in total sRNAs. These distribution patterns are typical of small RNA distributions in plants such as Medicago truncatula [33] and peanut [34].
The 96 miRNAs ranged from 20 to 23 nt in length, among which 21-nt sequences were the most common (73.96%), followed by 22-nt (19.79%), 20-nt (3.13%), and 23-nt (3.13%) sequences. According to the statistical novel miRNA nucleotide bias at each position, we were able to evaluate the accuracy of the novel miRNA predictions. We found that up to 50% of the novel predicted miRNAs possessed 5 0 U bases, and A bases were observed at a high frequency in the 10 th and 11 th base positions (Fig 4).

Identification of drought stress-related miRNAs
In the present study, 2298 unigene sequences were predicted as the targets of 75 conserved miRNAs, and 1274 unigene sequences were predicted as the targets of 18 novel miRNAs (S4 Table). The predicted target transcripts were known to be involved in a variety of stressrelated metabolic processes, including apoptotic processes, the ER-nucleus signaling pathway, oxidoreductase activity, fatty acid metabolic processes, leaf development and the biosynthesis of jasmonic acid, oxidative stress tolerance, leaf morphogenesis, xylem and phloem pattern formation, and the auxin-activated signaling pathway. However, considering our focus on drought stress, we carefully selected 39 miRNAs for the GO and KEGG analyses, in order to elucidate the key drought-resistance mechanisms of C. korshinskii from the Loess Plateau (Table 2).

Key miRNAs in the C. korshinskii water stress response
Regarding the gradient of moisture stress among C. korshinskii habitats, we compared the results from the laboratory water stress experiment and natural habitat. Combining the GO and KEGG annotation results (Table 2), we further carried out a trend analysis on the 39 miR-NAs (Fig 5), based on a decrease in precipitation from HuL to DaB. Five putative droughtrelated miRNAs (miR390, miR398, miR530, miR2119, and miR5559) were selected for the stem-loop qRT-PCR analysis because their target genes were involved in abiotic stress-related processes, such as oxidative stress, auxin-response, cell death, and senescence. Using qRT-PCR, we measured the expression of these miRNAs in the laboratory-and field-grown samples. We found that the genes targeted by the 5 selected putative drought-related miRNAs (miR390, miR398, miR530, miR2119, and miR5559) exhibited a strong response to changes in soil moisture, and the expression of the five conserved miRNAs was downregulated by water stress under both field and laboratory conditions, which well matched the transcripts per million mapped reads results from the sequencing analysis (Fig 6). Although additional studies are needed to elucidate the function of these miRNAs in the stress-resistance of C. korshinskii, the consistency of our observations in both field-collected and laboratory-grown samples supports the role of the five identified miRNAs in the moisture response of C. korshinskii along the precipitation gradient from HuL to DaB.
Key target genes and miRNA of C. korshinskii collaborate in responding to water stress To further validate the mechanism of the miRNA response to drought, we selected two unauthenticated miRNAs (miR2119 and miR5559) to investigate their target genes (S4 Table). We found that the Bax inhibitor-1 (BI-1) and Clathrin assembly protein (CAP) genes were potential targets of miR2119, whereas the histone regulation (hira) and sieve element occlusion (SEO) genes were potential targets of miR5559. Using qPCR, we found that the expression of the 4 genes in both the field-collected and laboratory-grown samples was significantly upregulated as a result of a water deficit effect, which suggested that there was a negative correlation between miR2119 and miR5559 and their specific target genes (Fig 7). Therefore, the 4 genes (BI-1 and CAP; HIRA and SEO), collaborating with their respective miRNAs (miR2119; miR5559), were involved in the stress response in C. korshinskii, where further study would bring more information on the functions of the 2 miRNAs.

Discussion
Drought stress is one of the major natural constraints on plant growth and development. As an important species for preventing desertification in the Loess Plateau, C. korshinskii is well adapted to water shortage, and even modifies its gene expression and morphology under drought conditions. For instance, our previous work revealed that the size and density of C. korshinskii trichomes are correlated with the precipitation gradient that occurs across the Loess Plateau, China [22]. However, the underlying regulatory mechanism of this species' adaptation to natural drought conditions remains poorly understood. In the present study, we generated a large amount of miRNA expression data from C. korshinskii using the high-throughput Solexa sequencing technology. By deep sequencing C. korshinskii small RNAs and aligning them to known plant miRNAs, we identified 490 known miRNAs, as well as 96 novel miRNAs. MiRNAs have been implicated in the regulation of plant responses to several types of biotic and abiotic stress. Some of the known miRNAs have been detected in other plant groups, whereas others are uni-que to legumes [35]. Thus, our study provides a foundation for further research into the miRNA-mediated gene regulation of stress responses in C. korshinskii.
The main aim of this study was to investigate the miRNAs of C. korshinskii, particularly those related to drought tolerance. By analyzing the functions of putative target genes using the GO analysis, we selected and screened 39 miRNAs with biological functions that were possibly related to drought tolerance ( Table 2). The targets of these miRNAs included genes involved in fatty acid metabolism, apoptotic process, programmed cell death, plant hormone signal transduction, and responses to osmotic and oxidative stress. The variance in miRNA expression among the sampling sites may contribute to the differential regulation of these target genes.
Many miRNAs were down-or upregulated in response to the precipitation gradient (Fig 5), and some of them have been shown to respond to biotic or abiotic stress by previous studies. For example, miR5244 and miR5232 have been reported to accumulate in response to arbuscular mycorrhizal (AM) symbiosis in M. truncatula roots [36]. In addition, in Arabidopsis, miR390, TAS3 tasiRNAs, and their auxin response factor targets constitute an autoregulatory network that quantitatively regulates lateral root growth; lower levels of the miR390 transcript facilitate the expression of auxin response factor genes, thus promoting the development of  DNA binding response to acid; response to stress.

MIR2119
bax inhibitor 1 fatty acid metabolic process; apoptotic process; response to stress; ER nucleus signaling pathway.

MIR2645
heat shock 70 kDa protein 16 oxidoreductase activity, acting on the CHCH group of donors; NAD or NADP as acceptor; response to stimulus.

MIR319
transcription factor GAMYB isoform X1 response to hormone; programmed cell death; regulation of cellular process.

MIR3629
heat shock 70 kDa protein 17 Protein processing in endoplasmic reticulum plastid; oxidoreductase activity,acting on the CHCH group of donors, NAD or NADP as acceptor; response to stimulus.
abscisic acid receptor PYL8 phosphatase inhibitor activity; regulation of phosphoprotein phosphatase activity; abscisic acid activated signaling pathway; Plant hormone signal transduction.
WRKY transcription factor DNA binding; nucleic acid binding; transcription factor activity.
zinc finger CCCH domain containing protein 32 cation binding.

MIR394
F-box only protein 6 leaf morphogenesis; xylem and phloem pattern formation; auxin activated signaling pathway.
mitochondrial inner membrane protease ATP23 endopeptidase activity protein processing.

MIR397
laccases oxidoreductase activity, acting on diphenols and related substances as donors, oxygen as acceptor lignin metabolic process.

MIR398
Oxidative stress Response to oxidative stress.

MIR408
replication protein A Nucleotide excision repair; DNA replication; Mismatch repairorganic cyclic compound binding.
E3 ubiquitin protein ligase Ubiquitin mediated proteolysis; membrane transition; metal ion binding; proteasome mediated ubiquitin dependent protein catabolic process.

MIR472
auxin response factor binding to auxin response elements.

MIR5175
BI1 protein fatty acid metabolic process; apoptotic process; response to stress; ER nucleus signaling pathway.

MIR5232
chlorophyll a/b binding protein Photosynthesis antenna protein splastid envelope; chloroplast thylacation binding; generation of precursor metabolites and energy.

MIR5244
WRKY transcription factor DNA binding; transcription factor activity.  auxin response factor 9 protein binding hormone mediated signaling pathway.

MIR5245
WRKY transcription factor nucleic acid binding; transcription factor activity; DNA binding response to amino acid; response to stress.

MIR5770
peroxiredoxin Q response to oxidative stress.
Calcium dependent protein kinase Plant pathogen interaction protein kinase activity; hormone mediated signaling pathway; phosphate containing compound metabolic process.

MIR6024
abscisic acid 8'hydroxylase oxidoreductase activity; response to fungus; response to light stimulus; sesquiterpenoid metabolic process; response to stress.
NAC domain protein DNA binding.
BI1 protein fatty acid metabolic process; apoptotic process; response to stress; ER nucleus signaling pathway.

MIR6485
ethylene responsive transcription factor ERF054 Plant hormone signal transduction; DNA binding.

MIR7539
MYB transcription factor MYB54 nucleic acid binding; response to UV; response to acid; regulation of cellular metabolic process.
Abscisic acid receptor PYL1 Plant hormone signal transduction; protein dimerization activity; isoprenoid binding regulation of phosphor protein phosphatase activity.
Salt tolerance protein transition metal ion binding; response to abiotic stimulus; response to organonitrogen compound; response to hormone. transcription factor activity; cell death; response to bacterium; response to organonitrogen compound; response to stress; nucleic acid binding.
(Continued) strong root systems and drought tolerance. In addition, miR530 is highly upregulated in response to fungal infection in chickpea [35], and its target genes vary under different conditions and in different plant species and tissues. One of the target genes of miR530 is a tandem zinc knuckle/PLU3 domain encoding gene (At5g43630; TZP), which is involved in the regulation of morning-specific growth in Arabidopsis [37]. miR398 targets two Cu/Zn superoxide dismutases [Cu/Zn Superoxide dismutase 1 (CSD1), and Cu/Zn Superoxide dismutase 2 (CSD2)] in Arabidopsis thaliana [38] and is downregulated in response to both biotic (Pseudomonas syringae) and abiotic stress (ozone and salinity) [39]. Experimental results indicated that CSD1 mRNA levels were found to be negatively correlated with miR398 levels during exposure to ozone, salinity, and biotic stress, which suggests that a link exists between oxidative stress and miR398 [39]. The evidence gathered by another study indicated that miR398 regulates the responses of rice to a wide range of environmental stress and to ethylene, and that miR398 functions by mediating the expression of CSDs (Cu/Zn superoxide dismutases) and cellular levels of reactive oxygen species [40]. MiR319 has been demonstrated to target TCP [for TEOSINTE BRANCHED/CYCLOIDEA/PROLIFERATING CELL FACTORS (PCF)] genes, which are involved in cell differentiation, leaf development, senescence, and the biosynthesis of jasmonic acid [32,41,42]. The miR319-mediated downregulation of target genes (TCPs) in overexpressing omiR319 transgenic plants causes changes in various biological processes, including those associated with water retention capacity, leaf wax synthesis, and salt uptake [32]. The prediction of potential target genes for the conserved and novel miRNAs was based on sequence homology [35], since many miRNAs respond similarly under both biotic and abiotic stress. The existence of crosstalk between the pathways that are involved in regulating stress indicates that miRNAs that respond to the precipitation gradient may be responsible for the adaptive mechanisms of wild C. korshinskii on the Loess Plateau [35]. Moreover, our finding that some key conserved miRNAs respond to water stress in C. korshinskii agrees with previous studies. During drought stress, the levels of miR408 transcripts decrease significantly in rice, whereas they remain elevated in tolerant rice [43]. By predicting the target genes of osa-miR408 in rice, it was determined that osa-miR408 regulates members of the plastocyanin-like protein family. Phytocyanins (PCs) are ancient blue copper proteins that can bind to a single copper atom, function as electron transporters, and are involved in responses to drought stress [43]. MiR408 was also downregulated by drought stress in Populus [44]. MiR394 was identified among the stress-responsive microRNAs of Populus [44,45] and was downregulated under drought conditions [44]. The predicted targets of miR394 are annotated as dehydration-responsive and F-box proteins, which were recently reported to be differentially regulated under stressful conditions and to play significant roles in the abiotic stress response pathway.
The identification of miRNA target genes is critical for understanding their functional implications. Of the miRNAs that were down-or upregulated by the precipitation gradient, miR2119 and miR5559 have only been reported in legumes, and their function is still unclear. By computational analysis of the transcriptome data, BI-1 and CAP were predicted as the targets of the miR2119. BI-1, a widely conserved cytoprotective protein that is localized in the endoplasmic reticulum membrane, is an attenuator of both biotic and abiotic types of cell death [46,47]. Moreover, CAP is a key protein in the formation of clathrin-coated pits, which are an intermediate in the process of clathrin vesicle formation [48]. HIRA and SEO are the predicted target genes of miR5559. HIRA is involved in cell senescence via a pathway that appears to depend on the flux of heterochromatic proteins through promyelocytic leukemia bodies [49], while SEO genes, which encode forisome subunits (structural phloem proteins unique to the sieve elements of legumes), have been identified in M. truncatula and other legumes [48].
To support our analysis of wild C. korshinskii, we designed an indoor experiment to detect the expression of the 5 miRNAs (miR390, miR398, miR530, miR2119, and miR5559) in response to reduced soil moisture. Under water stress, the expression of the 5 miRNAs was consistent with both the sequencing and qPCR analyses of the wild plant samples. Furthermore, the expression of BI-1 and CAP, as well as HIRA and SEO, the predicted target genes of miR2119 and miR5559, respectively, were measured for both the wild and indoor-grown plants. Although we found that the expression levels of the miRNAs and their target genes were negatively correlated, indicating that key miRNAs could be involved in water stress responses, further experiments are needed to elucidate interactions of miRNAs with their predicted targets and their regulatory mechanisms in physiological processes.
The predicted targets for the novel miRNAs were involved in plant growth, development, stress response, hormone signal transduction, and transcription. With a major concern on responding to the precipitation gradient on the Loess Plateau, novel_mir_43 was predicted to target the WRKY transcription factor, which is essential in the pathogen and salicylic acid responses of higher plants, as well as in a variety of other plant-specific reactions [50,51]. One study found that the WRKY transcription factor was differentially regulated in plants that were exposed to abiotic stress factors such as NaCl, polyethylene glycol, cold, and heat stress [52]; another reported that the transcription factor enhanced drought tolerance by retarding leaf wilting and increasing the survival rate of green plant parts [53]. Putative targets were also predicted for novel_mir_49, including the WRKY transcription factor, Na+/H+ antiporter, zinc finger CCCH domain-containing protein, and a series of molecules involved in plant hormone signal transduction, such as the auxin response factor, abscisic acid receptor, and ethylene responsive transcription factor. Because they clearly responded to the precipitation gradient on the Loess Plateau, the validity of these novel miRNAs deserves further investigation, in order to confirm their differential expression, which will provide clues about their biological functions.

Conclusions
Using deep-sequencing technology, we identified 490 known miRNAs and 96 novel miRNAs in C. korshinskii sampled from the Loess Plateau, including 39 miRNAs that may be associated with drought response. After analyzing the response of the miRNAs to a precipitation gradient, we investigated the expression of 5 key miRNAs under both natural and experimental water stress conditions. Our results provide a better understanding of the functional relationships and specific miRNAs involved in the stress-related gene regulatory pathways. Additional work is needed to further elucidate this relationship, utilizing the available unique datasets in C. korshinskii. Further characterization of the targets of drought-associated miRNAs will also help understand the drought response and the tolerance of C. korshinskii, which serves as an important model for stress resistance in plants facing desertification.
Supporting information S1