Genome-Wide Analysis of MicroRNA Responses to the Phytohormone Abscisic Acid in Populus euphratica

MicroRNA (miRNA) is a type of non-coding small RNA with a regulatory function at the posttranscriptional level in plant growth development and in response to abiotic stress. Previous studies have not reported on miRNAs responses to the phytohormone abscisic acid (ABA) at a genome-wide level in Populus euphratica, a model tree for studying abiotic stress responses in woody plants. Here we analyzed the miRNA response to ABA at a genome-wide level in P. euphratica utilizing high-throughput sequencing. To systematically perform a genome-wide analysis of ABA-responsive miRNAs in P. euphratica, nine sRNA libraries derived from three groups (control, treated with ABA for 1 day and treated with ABA for 4 days) were constructed. Each group included three libraries from three individual plantlets as biological replicate. In total, 151 unique mature sequences belonging to 75 conserved miRNA families were identified, and 94 unique sequences were determined to be novel miRNAs, including 56 miRNAs with miRNA* sequences. In all, 31 conserved miRNAs and 31 novel miRNAs response to ABA significantly differed among the groups. In addition, 4132 target genes were predicted for the conserved and novel miRNAs. Confirmed by real-time qPCR, expression changes of miRNAs were inversely correlated with the expression profiles of their putative targets. The Populus special or novel miRNA-target interactions were predicted might be involved in some biological process related stress tolerance. Our analysis provides a comprehensive view of how P. euphratica miRNA respond to ABA, and moreover, different temporal dynamics were observed in different ABA-treated libraries.

Gene expression can be regulated at transcriptional, posttranscriptional and translational levels, with miRNAs negatively regulating it at the posttranscription level (Obernosterer et al., 2006). With assistance from one of the argonaute proteins (AGO1), miRNAs regulate complementary target mRNA via specific cleavage or through translational inhibition (Phillips et al., 2007;Voinnet, 2009;Rogers and Chen, 2013). In plant, miRNAs are transcribed to form long primary transcripts (pri-miRNA). Then the pri-miRNAs are trimmed, producing miRNA precursors (pre-miRNA) with stem-loop secondary hairpin structure. Then pre-miRNA is cleaved into a double stranded RNA consisting of a miRNA and its complementary sequence miRNA * (Park et al., 2002;Han et al., 2004;Kurihara et al., 2006). The double-stranded RNA is conveyed to the cytoplasm. One of the TWO strands acts as mature miRNA, loaded on the RNA induced silencing complex (RISC) to target mRNA (Bartel, 2004;Jones-Rhoades et al., 2006;Brodersen et al., 2008), meanwhile the other strand, miRNA * , is typically degraded (Jones-Rhoades et al., 2006).
The phytohormone abscisic acid (ABA) plays vital physiological roles in response to abiotic stress. The ABA level is controlled by complex regulatory mechanisms including biosynthesis, catabolism, transport and signal transduction pathways. The 9-cis-epoxycarotenoid dioxygenase 3 (NCED3), a key enzyme for ABA biosynthesis, is induced by drought stress, and it upregulates endogenous ABA levels in overexpressed transgenic plants Thompson et al., 2000;Iuchi et al., 2001;Schwartz et al., 2003). CYP707A is the key enzyme for ABA catabolism. Cyp707a3-1 mutant induces endogenous ABA level and reduces transpiration rate, thereby resulting in high tolerance to drought stress in Arabidopsis (Kushiro et al., 2004). Pyrabactin resistance1/PYR1-like/regulatory components of ABA receptors (PYR/PYL/RCAR protein family), the type 2C protein phosphatases (PP2Cs), and subfamily 2 of the SNF1-related kinases (SnRK2s) are a major breakthrough in the field of ABA signaling (Ma et al., 2009;Melcher et al., 2009;Nishimura et al., 2009;Park et al., 2009;Santiago et al., 2009;Yin et al., 2009;Gonzalez-Guzman et al., 2012). PP2C has a negative role in ABA signaling (Miyazono et al., 2009). Under stress conditions, the abundance of PP2C transcripts increase (Rubio et al., 2009;Szostkiewicz et al., 2010). With the help of ABA, SnRK2s are released from the PP2C inhibition, and activate their downstream targets (Klingler et al., 2010). Srk2d srk2e srk2i (srk2d/e/i) triple mutants show reduced tolerance to drought stress and highly enhanced insensitivity to ABA (Fujita et al., 2009). PYR/PYL/RCAR protein family, ABA receptor, play a major role in regulation of seed germination and establishment, basal ABA signaling required for vegetative and reproductive growth, stomatal aperture, and transcriptional response to the hormone (Gonzalez-Guzman et al., 2012).
Recent evidence indicates that miRNA and ABA affect each other, and that the expression levels of some miRNAs are regulated by ABA. For example, miR159, miR169, miR172 are regulated by ABA in the embryogenic callus of Japanese larch (Larix leptolepis) . In maize roots, the entire miR169 family is downregulated by ABA (Luan et al., 2014).
Under ABA treatment, the expression level of miR169a decreases, and the target of miR169, NFYA5 is upregulated (Li et al., 2008;Ni et al., 2013). The cistronic miRNA pair, miR842, and miR846, is a product of alternative splicing regulated by ABA in the roots of Arabidopsis. ABA regulating the alternative splicing, leading to reduce the expression of miR846, along with the accumulation of its target jacalin At5g28250 (Jia and Rock, 2013a,b).
Conversely, miRNAs change the sensitivity of plants to ABA. In Arabidopsis, miR160 overexpression reduces ABA sensitiveness during germination (Liu et al., 2007), and also causes abnormal root morphology, that leads to the lack of gravitropic responses in root tips and the promotion of more adventitious roots (Wang et al., 2005). In Arabidopsis, miR172b overexpression increased sensitivity to ABA and osmotic stress during a specific postgerminative stage (Zou et al., 2013).
In addition, exogenous ABA also influences some miRNA expression, with miRNA regulating the downstream genes of ABA. Overexpressing miR168a lead to ABA hypersensitivity and drought tolerance, while the loss-of-function mutant miR168a-2 displays ABA hyposensitivity and drought hypersensitivity. Both the precursor and mature of miR168 are induced by ABA . ABA also positively regulates both mature miR394 and precursor miR394a/b in Arabidopsis. Although ABA negatively regulates LCR, which is the target of miRNA394, the overexpression of miR394a/b leads to ABA hypersensitivity and ABA-associated phenotypes, whereas LCR overexpressing plants show ABA resistant phenotypes. Moreover, the overexpression of miR394a/b plants accumulate higher levels of ABA-induced hydrogen peroxide and superoxide anion radicals compared to wild-type and LCR-overexpressing plants (Song et al., 2013). MiR159 is induced by ABA, and it targeted MYB33 and MYB101 germinating seeds of Arabidopsis; the two MYB transcription factors are positively regulated by ABA. Overexpression of miR159a reduces sensitivity to ABA (Reyes and Chua, 2007). However, previous research has generally focus on the relationship between a specific single miRNA and ABA, with most regulation studied in the seed or root.
Populus euphratica (P. euphratica), which exhibits remarkable tolerance to environmental stresses, is among a few tall tree that can survives in saline and alkaline area. It is a model woody plant for studying the molecular mechanisms of abiotic stress responses (Ye et al., 2009;Lv et al., 2014). Previous studies show miRNA participate in the drought and salt stress responses of P. euphratica, and the high-throughput sequencing has widely been used in miRNA research of P. euphratica (Li et al., 2011Si et al., 2014), while other species poplus also involved in stress response (Ren et al., 2012(Ren et al., , 2013(Ren et al., , 2015Chen et al., 2012aChen et al., ,b, 2015Shuai et al., 2013). Hence, screening for ABA-responsive miRNA in P. euphratica may help elucidate the responses of woody plants to ABA and thus the mechanisms underlying such responses to abiotic stress; however, the relationship between the mechanism of ABA and miRNA regulation at a genomic level has not been reported. In this study, we attempt to provide new insight into this issue by identifying more relationships between ABA and miRNA regulation mechanisms in P. euphratica leaves.

Plant Material and Growth Conditions
Uniformly grown 1-year-old P. euphratica acquired from the Xinjiang Uygur Autonomous Region of China, were planted in individual 5 L pots containing a loam soil and placed in a greenhouse at Beijing Forestry University. Each container contained three individual plants of similar height. They were well watered and grown under natural conditions for 4 months (from April to July). The relative soil moisture content (RSMC) was measured using a FieldScout TM TDR 300 Soil Moisture Meter (Spectrum Technologies, Aurora, IL, USA). For all seedlings, the RSMC was controlled within 70-75%. To ensure that all leaves received a similar amount of ABA, an aqueous solution of ABA was used to water the P. euphratica. For the treated groups, 1 L 300 µM ABA solution was applied to each pot; and for the control group, pure water instead of ABA solution was used. Every pot was placed on a tray to prevent the solution from flowing away. Each group included three libraries from three individual plantlets as biological replicate. Leaves collected 1 day after adding the ABA solution were considered the shortterm ABA treatment and named SL1, SL2, and SL3. Those collected 4 days after adding the ABA solution were considered the long-term ABA treatment and named LL1, LL2, and LL3. The control groups were named CL1, CL2, and CL3. Leaf tissues were collected at about 10:00 h, and all 8-12th (count from the apex) mature leaves, from independent plants in each group were collected. Prior to collecting, the photosynthetic rate, intercellular CO 2 concentration, stomatal conductance and transpiration rate were measured using a Li-6400XT portable photosynthesis system (Li-Cor Inc., Lincoln, NE, USA). Two leaves in one plantlet were random selected, and measured with three technical repeat. And all three biological replicate for each group were measured. All of the collections were immediately frozen in liquid nitrogen, and stored at −80 • C until RNA extraction.

RNA Extraction, Libraries Construction, and Sequencing
Using the modified CTAB method, total RNA enriched with small RNA (sRNA) was isolated from the P. euphratica leaves in the nine collections (Jaakola et al., 2001). RNA quality and integrity were checked using a 2100 Bioanalyzer with the RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries were generated using NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB, USA.) following manufacturer's recommendations. All the nine libraries from the three groups were for sRNA sequencing. High-throughput sequencing was performed using Illumina HiSeq technology was according to the manufacturer's protocol (Illumina, San Diego, CA, USA).

Analysis of Small RNA Sequencing Data
From the raw sequence reads obtained from the sRNA sequencing, we first removed low quality reads, including those shorter than 18 nt or longer than 30 nt, those with more than 10 nt single nucleotide repeats, or more than 10% N, and those with 5 ′ adapter contaminants, or without a 3 ′ adapter or the insert tag. Then the 3 ′ and 5 ′ adapters were removed to obtain clean reads without any mismatches, which were mapped to P. euphratica genome (Ma et al., 2013) using bowtie software (http://sourceforge.net/projects/bowtiebio/files/) without any mismatch.

sRNA Reads Annotation and miRNA Identification
All mapped reads were annotated as follows. First, mapped reads were annotated as conserved miRNAs, which were previously discovered, which were registered in miRBase (Release 21) for Populus trichocarpa by BlastN algorithm with both mature and hairpin were without any mismatches. The remaining reads were annotated as non-coding RNA (i.e., tRNAs, rRNAs, scRNAs, snRNAs, and snoRNAs). The sequences were collected from the GenBank (http://www.ncbi.nlm.nih.gov/genbank/) and Rfam (11.0 release, http://rfam.xfam.org/) database. The similarity was investigated using the BlastN algorithm. The RepeatMasker was used to remove the repeat-associated RNAs (Repbase v.18.07, http://www.girinst.org/). Then nat-siRNAs were removed (P. trichocarpa in PlantNATsDB, http://bis.zju.edu.cn/pnatdb/). The remaining sRNA exactly matched the mRNA exons and introns in the P. euphratica genome (Ma et al., 2013). The remaining reads were used to predict novel miRNA utilizing miREvo and miRdeep2, based on secondary structure, the Dicer cleavage site, and the minimal folding free energies (Friedlander et al., 2012;Wen et al., 2012).

Differential Expression Analysis of miRNA Response to ABA
The expression levels of miRNA between the two groups were compared to determine which miRNAs were differentially expressed. The frequency of miRNA read counts was normalized as transcripts per million (TPM): normalized expression = (number of miRNA reads/total number of clean reads) * 1,000,000 . Raw data were used with the "DESeq2" library in the R statistical software package for this analysis (Love et al., 2014). The P-values were adjusted, and P-adjusted < 0.05 was considered to indicate significantly different expression (Benjamini and Hochberg, 1995).

miRNA Targets Prediction and Function Analysis
Conserved and novel miRNA target were predicted by the psRNA Target Server (http://plantgrn.noble.org/psRNATarget/) with default parameters (Dai and Zhao, 2011). All discovered miRNA targets and differentially expressed miRNA targets were classified based on gene ontology (GO) performed using the online agriGO program (Du et al., 2010). miRNA targets prediction and GO classification analysis were based on the P. euphratica genome (Ma et al., 2013).

Real-Time Quantitative PCR Analysis of miRNAs and Predicted Targets
Ten miRNAs were randomly selected for real-time quantitative polymerase chain reaction (qPCR) for each comparison. The RNAs were extracted from each sample using CTAB method (Jaakola et al., 2001). The miRNA First-Strand cDNA Synthesis Kit (Aidlab Biotechnologies, Beijing, China), which is based on a poly-adenylation protocol was used for mature miRNA reverse transcription, and miRNA Real-Time PCR assay kit (Aidlab Biotechnologies, Beijing, China) was used for real-time qPCR. Real-time qPCR was performed using a total reaction volume of 20 µL, which contained 0.5 µL of diluted cDNA, 0.8 µM primer mix, 10.0 µL of 2× miRNA qPCR Mix, and 8.7 µL ddH 2 O, which were performed using an ABI StepOnePlus TM instrument (Applied Biosystems, Foster City, CA, USA). Amplification reactions were performed as follows: 95 • C for 10 s, 60 • C for 20 s and 72 • C for 30 s All reactions were performed in triplicate. P. euphratica 5.8 s rRNA was used as internal control for miRNA (Lu et al., 2008), and the 2 − CT method was applied to calculate the relative changes in gene expression from real-time qPCR experiments (Livak and Schmittgen, 2001). All primers used for real-time qPCR are listed in Supplementary Data 7.
The expression analyses of several target genes were also examined using qRT-PCR. The FastQuant RT Kit (with gDNase) (Tiangen Biotech CO., LTD, Beijing, China), was used predicted targets mRNA reverse transcription, and SuperReal PreMix Plus (SYBR Green) (Tiangen Biotech CO., LTD, Beijing, China) was used for real-time qPCR according to the manufacturer's instructions. Real-time qPCR was performed using a total reaction volume of 20 µL, which contained 0.8 µL of diluted cDNA, 1.2 µM primer mix, 10.0 µL of 2× SuperReal PreMix Plus, and 8.0 µL ddH 2 O, which were performed using an ABI StepOnePlus TM instrument (Applied Biosystems, Foster City, CA, USA). Amplification reactions were performed as follows: 95 • C for 10 s, 60 • C for 20 s and 72 • C for 30 s. All reactions were performed in triplicate. P. euphratica UBQ was used as internal control for mRNA . And the 2 − CT method was applied to calculate the relative changes in gene expression from real-time qPCR experiments (Livak and Schmittgen, 2001). All primers used for real-time qPCR are listed in Supplementary Data 7.

Accession Number
Sequencing data obtained in this work have been submitted to the Sequence Read Archive the accession number SRP077948.

Physiological Characterization of P. euphratica in Response to ABA
One-year-old P. euphratica seedlings were exposed to soil with ABA. The photosynthetic rate, transpiration rate and stomatal conductance were greatly affected after 1 day of treatment, but had recovered by day 4. Compared to control, the photosynthetic rate decreased 55% on day 1 but was 1.1 times the control level on day 4 (Figure 1). Stomatal conductance decreased 80% after 1 day of treatment and recovered to 82% of the control level after 4 days of treatment. Similarly, the transpiration rate decreased 73% after 1 day of treatment, but recovered to 68% of the control level after 4 day of treatment. The intercellular CO 2 concentration generally decreased over the 4 days compared to controls, and was 9 and 14% lower, respectively, after 1 and 4 days of treatment. In this study, the stomatal conductance decreased as the photosynthetic rate and transpiration rate fell after 1 day of treatment and then returned to normal, indicating that stomatal conductance was a key limitation to the photosynthetic rate in an early ABA response, although other factors, such as the activity of Rubisco, photosystem I (PSI), and photosystem II (PSII) were all essential factors for photosynthesis, may also inhibited photosynthesis under ABA conditions. The impact of the ABA treatment on photosynthesis in the early stage appeared to be recoverable.

Overview of sRNA Libraries Data for High-Throughput Sequencing
To identify and characterize conserved and novel miRNAs in P. euphratica, nine sRNA libraries from CLs, SLs, and LLs were constructed. Through high-throughput sequencing, more than 10 million raw reads were obtained. After filtering the lowquality reads, and adaptor and contaminant sequences, ∼90% of the clean reads remained ( Table 1). Those unique sequences were then perfectly mapped to the P. euphratica genome (no mismatches allowed) (Ma et al., 2013), and the results showed that in most libraries over 50% of the total sRNA matched the P. euphratica genome perfectly. Among these libraries, nearly 60% of the total sequences were matched in all of the SL libraries, which was the highest proportion ( Table 1).
The size distribution of the sequencing of unique reads from the nine libraries was similar, and the majority of sRNAs were 21-24 nt in length (Figure 2). Specifically, 24 nt length reads were the most abundant in most libraries (except CL1), and over 40% were 24 nt sRNAs. The 21 nt class was the second most abundant. The sRNAs with different sizes perform different functions: 21 nt sRNAs usually mediate posttranscriptional gene silencing, while 24 nt sRNAs typically perform gene silencing mediated by RNAdependent DNA methylation and heterochromatin maintenance (Zhang and Zhu, 2011;He et al., 2014;Lewsey et al., 2016).

Conserved miRNAs Discovered in P. euphratica
Conserved miRNA, non-coding RNA, repeat-associated RNA, nat-siRNA, exons, introns, and unknown sequences were identified successfully (Supplementary Data 1).The known conserved miRNAs were first identified, with no mismatches. Presently 253 unique mature sequences (include both 5p and3p), belonging to 136 families registered in miRBase (Release 21) for P. trichocarpa, and here 151 unique mature sequences, belonging to 75 miRNA families, were expressed in at least one of the nine libraries. Detailed information about the known miRNAs is shown in Supplementary Datas 2, 3. Seventy-two conserved miRNA families were identified, including 91 unique miRNAs belonging to 27 miRNA families and highly conserved (identified in more than 10 species of angiosperm, Figure 3), and 60 miRNAs belonging to 45 families were specific to Populus. The average TPM of the highly conserved and Populusspecific miRNA were 53,197.01 and 8,540.66, respectively. The abundance of highly conserved miRNA families was significantly higher than the Populus-specific miRNA.

Analysis of Novel miRNAs from P. euphratica
All novel miRNA candidates were obtained from miREvo and miRdeep2, and confidently met the novel miRNA requirements: (i) sequencings represented both the miRNA and miRNA * , and (ii) in miRNA * -deficient cases, the novel miRNA came from multiple and independent libraries (Xu et al., 2013). Ninety-four unique sequences were identified as novel miRNAs, while for 56 of them we found the corresponding miRNA * sequences ( Table 2 and Supplementary Data 4). The length of mature novel miRNA sequences were 18-24 nt. Mfold was utilized to predict the stem-loop structure of precursors of all non-conserved miRNAs (Zuker, 2003), and the information was provided in Figure 4 and Supplementary Figure  1. The pre-miRNA length ranged from 56 to 292 nt. The negative minimal folding free energies (MFEs) varied from −127.0 to −23.4 kcal/mol; the average value was -49.47, which was much less than that of the tRNA (−27.5 kcal/mol) and rRNA (−33 kcal/mol) (Bonnet et al., 2004). The minimal folding free energy index (MFEI) values ranging from 0.66 to 2.93. Seventy-one miRNAs (∼74%) were above 0.85, which was a key characteristic in distinguishing pre-miRNAs from other sRNAs (Zhang et al., 2006). The MFEI of all pre-miRNAs was significantly higher than that of tRNAs (0.64), rRNAs (0.59), and mRNAs (0.62-0.66). All data indicated that pre-miRNAs possessed highly stable hairpin structures. A nucleotide bias tendency indicated that in 46 of 94 cases, the first nucleotide was U in 48% of the novel non-conserved miRNAs, in agreement with that of conserved miRNAs in other plants (Voinnet, 2009). These results agree with that miRNAs were preferentially loaded onto AGO1-containing RISC, and preferentially contain a U at the 5 ′ -end ( Table 2; Mi et al., 2008).
The expression levels of these novel miRNA showed a large TPM range (Supplementary Data 3). Seven of 94 novel miRNA were above 10,000 TPM for the sum of all nine libraries, with the most abundant being miR-n1 at an amazing 623,863.3 TPM. More than half of these miRNAs (∼78%) were discovered to be less than 1000 TPM, which is consistent with the finding that novel miRNAs are often expressed at a lower level compared to conserved miRNAs. The average TPM of all Populus-specific microRNAs, including known Populus-specific microRNAs and novel miRNA was 8,611.81. The abundances of those consistent with speciesspecific conserved miRNAs were expressed lower than those of highly conserved miRNAs (Zhao et al., 2010).

Differentially Expressed miRNAs Response to ABA
To compare the expression of miRNA in different samples, the "DEseq2" library in the R statistical software package was performed on the raw counts for all of the conserved and novel miRNAs. The analysis of the differential expression of miRNAs showed that 31 known miRNAs and 31 novel miRNAs under ABA were significantly differentially expressed in P. euphratica ( Table 3 and Supplementary Data 5).
Results showed that of miRNA expressed after 1 day of treatment, eight conserved and six novel miRNAs were upregulated, whereas five conserved and eight non-conserved miRNA were downregulated, compared to that control. Among them miR-n91 and peu-miR408-3p were upregulated, and peu-miR6421 and miR-n87 was downregulated. Similarly, after 4 days of treatment, eight conserved and six novel miRNAs were upregulated, whereas seven conserved and nine novel miRNAs were downregulated, compared to   the control. For example, miR-n87 and peu-miR6462-3p were downregulated, and peu-miR6421 and miR-n9 were upregulated. After 4 days of treatment, 11 conserved and 10 novel miRNAs were upregulated, whereas 12 conserved and 11 novel miRNAs were downregulated, compared to those after 1 day of treatment. For example, miR-n19 was downregulated, and peu-miR6462-3p and peu-miR6421 were down-regulated. In comparing the three groups with each other, only five miRNAs (miRNA-n2, miRNA-n91, peu-miR164a, peu-miR395a, and peu-miR6421) were all significantly differentially expressed. Compared to the control libraries, nine miRNAs were significantly differentially expressed in both the 1-and 4-day ABA treatment libraries. Compared to the control, the expression levels of three miRNAs were all significantly increased in the 1 day of treatment libraries, and the other six miRNAs decreased for 4 days. Between SL vs. CL and LL vs. SL, nine miRNAs were significantly differentially expressed. Three of them first ascended, then descended, while the others were opposite. Compared to the 4 days of treatment libraries, both in the control and 1 day of treatment libraries, 11 miRNAs were significantly differentially expressed. Five of them were upregulated, and other six miRNAs were all downregulated.

Confirmation of miRNAs by Real-Time-qPCR
Real-time qPCR analysis was utilized to confirm the expression patterns for the significantly differentially expressed miRNA. Ten miRNAs were selected for each comparison (SL/CL, LL/CL, and LL/SL) with three experimental and three biological replicates to validate and measure the sequencing results. The comparison expression levels of these miRNAs between real-time qPCR and sequencing analyses were consistent (Figure 5).

Prediction and Annotation of Target Genes of miRNAs
To understand the functions of miRNA responses to ABA, the prediction and identification the function of their targets is a crucial step of in analysis. The psRNA Target Server was used to predicted target of miRNAs. In total of 4132 targets genes were predicted for known and novel miRNAs. 1584 were for novel miRNAs, and 2548 were for conserved miRNAs (Supplementary Data 6).
Generally, several targets were regulated by a single miRNA. MiRNA-n61 targeted 155 transcripts, the most targets of all discovered miRNA. Considering conserved miRNAs only, miRNA482a.1 target 129 transcripts. Whereas a single gene could be targeted by multiple miRNA. CCG015515.1 was targeted by 12 members of the miRNA169 family and two novel miRNAs, miR-n62 and miR-n91. Furthermore, one gene being targeted by several miRNAs from at least two miRNA families was not unusual. In particular, CCG030854.1 was targeted by six miRNAs, peu-miR478b, peu-miR478d, peu-miR481a, peu-miR481b, peu-miR6421, and peu-miR7812, which were belonged to four different miRNA famines.
To better understand the functions of these genes, GO analysis was employed to classify target genes based on their involvement in biological processes, cellular components and molecular functions. Findings showed that 2179 of all predicted target genes could be categorized into totally 706 GO terms, including 260 biological processes, 365 molecular functions and 81 cellular components. For all differential expressed miRNA targets 280 GO terms, included 139 biological processes, 116 molecular functions and 25 cellular components. Compared to all discovered miRNA targets, some GO terms had different proportions. The secondary level GO terms for all of differentially expressed miRNA targets, all miRNA targets and the reference genome was determined (Figure 6). In this study, the cellular components, the proportions of differentially expressed miRNA targets of "organelle" (GO: 0043226), "cell part" (GO: 0044464), and "cell" (GO: 0005623) were higher compared to those of all discovered miRNAs. In molecular functions, that of  "binding" (GO: 0005488), "transporter activity" (GO: 0005215) and "electron carrier activity" (GO: 0009055) were also higher than that for all discovered miRNA here. Similarly, in biological processes, those of "death" (GO: 0016265), "response to stimulus" (GO: 0050896), and "immune system process" (GO: 0002376) were higher than that for all of the discovered miRNAs.

High-Throughput Sequencing of P. euphratica
Compared with other deep sequencing studies for P. euphratica before, the P. euphratica genome was utilized as the reference genome instead of P. trichocarpa genome in this study. Finding additional new P. euphratica-specific miRNAs was desirable here. Consistently, 94.7% (89/94) of novel miRNAs were not identified in other P. euphratica studies before, where P. trichocarpa genome was used as reference genome. Only 5 out of 94 novel miRNAs (miRNA-n5, miRNA-n7, miRNA-n33, miRNA-n34, and miRNA-n82) had been found Si et al., 2014; Table 4). And more novel P. euphratica-specific miRNAs, which do not exist in the P. trichocarpa genome, were successfully discovered. All of the novel miRNAs discovered had not been reported in other previously studies in P. trichocarpa or registered in miRBase (Puzey et al., 2012;Shuai et al., 2013). Our results indicate that it is a more useful approach to discover P. euphratica-specific miRNAs based on the P. euphratica genome used here. Several conserved ABA-responsive miRNAs in this study, also significant differentially expressed under drought or salt in previous studies in P. euphratica Si et al., 2014). For example, the expression level of two miRNAs, peu-miR1446 and peu-miR319a, significant changed under salt treatment; four miRNAs, peu-miR399c, peu-miR403-3p, peu-miR472a, and peu-miR6421 were all expression significant changed to respond to drought, other nine miRNAs, including peu-miR156a, peu-miR160a-3p, peu-miR164a, peu-miR168a-3p, peu-miR171c-3p/d-3p, peu-miR395, peu-miR408, peu-miR472b, FIGURE 5 | Verification of selected miRNAs by real-time quantitative PCR. Differentially expressed miRNAs identified by sequencing were confirmed by real-time qPCR, and their expression levels were compared among the three groups. The expression level of miRNA in deep sequences was performed with R statistical software package, which was"DESeq2" library used with raw date. *P < 0.05, **P < 0.01. and peu-miR475, were all differentially expressed under both salt and drought. These observations revealed that some conserved miRNAs respond to multiple forms of stress, and that complex miRNA effects are involved in resistance to abiotic stress in P. euphratica.
miRNAs have been shown to play an important role in abiotic stress responses in plants (Ferdous et al., 2015), and ABA is a core signal for abiotic stress responses (Lee and Luan, 2012). The analysis regarding the mechanism between ABA and miRNA at a genomics-level had not been reported previously. In total, 79 conserved miRNAs and 31 novel miRNAs were differentially expressed in response to ABA. Many conserved plant miRNAs regulatory mechanisms have been previously reported, which were also discovered here. For example, a negative regulator in responses abiotic stress, no apical meristem (NAM), was targeted by miR164, which was upregulated by ABA (Souer et al., 1996;Wang et al., 2009;Puranik et al., 2012;Fang et al., 2014). In addition, miR395 targeted two different groups of genes, ATP sulfurylases and SULTR2;1, which are both involved in sulfate translocation and assimilation (Jones-Rhoades and Bartel, 2004;Allen et al., 2005). Targets of all discovered miRNAs and differentially expressed miRNAs were both predicted and subjected to a functional analysis. Predicted target genes of differentially expressed miRNA were more focus on response to ABA addition, with the help of further study that may provide us insight into the molecular mechanism underlying ABA in P. euphratica.

Photosynthesis and Stomata Movement
Photosynthesis was one of the most sensitive physiological processes responsive to ABA. ABA significantly influenced the transcriptional abundance of genes involved in photosynthesis (Yamburenko et al., 2013;Mou et al., 2015). Protection of the photosynthesis apparatus is very important for stress tolerance (Xiao et al., 2015). Plants can survive severe drought by maintaining the photosynthesis ability to avoid carbon-starvation (McDowell et al., 2008). Similar findings were recorded indicating that P. euphratica maintained high photosynthetic rates under moderate drought stress levels (Chen et al., 2006;Tang et al., 2013). In the present study, we found that the photosynthetic rate of P. euphratica first decreased and then increased. At 4 days of treatment, the photosynthesis rate of the treated seedlings was higher than that of the control. The internal CO 2 concentration, however, always decreased suggesting that after 1 day of treatment, stomatal limiting was the most important factor for photosynthesis and transpiration rate, corroborating previous studies (Chaves et al., 2009). After 4 days of treatment, the increase of photosynthesis consumed more CO 2 , and then loaded internal CO 2 concentration continued to decline. Even when the stomata were still mildly closed, with the transpiration rate lower than of control, the photosynthetic rates still remained at high levels. ABA is a phytohormone, and exogenously supplied ABA can be considered simulating the stress that caused by ABA increasing, suggesting that P. euphratica could maintain a high photosynthetic ability to reserve the energy and materials necessary for stress adaptation.
Although variance is correctly used it might be confusing as it is a term used in statistics, variation or changes would be more adequate. ROP2 inhibited ABA-induced stomatal closure (Lemichez et al., 2001;Hwang et al., 2011;Miyawaki and FIGURE 6 | GO analysis of miRNA putative target genes. GO annotation categorized all of the predicted miRNA target genes and differentially expressed miRNA target genes into biological processes, molecular functions, and cellular components. Yang, 2014). Two ROPs were the predicted target of peu-miR6462a-3p. The expression of peu-miR6462a-3p decreased in response to ABA, while ROPs were up-regulated (Figure 7J), which would lead to ABA-induced stomatal closure. That agreed with the phenomenon in this study that stomatal conductance recovered after 4 days of treatment. CAXs involved in calcium (Ca 2+ ) transport and homeostasis (Conn et al., 2011;Manohar et al., 2011;Punshon et al., 2012). CAX was a predicted target of miR-n10. After ABA treatment 4 days, miR-n10 was downregulated by ABA, while CAX was upregulated ( Figure 7H). In addition, due to CAX genes playing a critical role in sequestering Ca 2+ into vacuole (Barkla et al., 2008), the increasing of CAX resulted in stomatal opening, which is agreed with the stomatal conductance in this study. And it can been speculated that the regulation mechanism between miRNA and stomata movement may exist, which need further exploration.

miRNA Take Part in ABA Synthesis and Metabolism
ABA is one of the final products of carotenoids Wang et al., 2013). Carotenoid oxygenase was one of the peu-miR408-5p predicted targets. As peu-miR408-5p was upregulated by ABA after being treated for 1 day, while carotenoid oxygenase was downregulated (Figure 7E), which would influenced the synthesis of ABA. The group A PP2C interacted with SnRK2. And without ABA, SnRK2 is inhibited  A: Si et al. (2014); B: Li et al. (2013). by PP2C; with ABA, PP2C binds to the receptor to release the inhibition SnRK2 (Umezawa et al., 2009;Ng et al., 2011;Soon et al., 2012). MiRNA-n87 was downregulated by ABA. PP2C was a predicted gene of miRNA-n87, while PP2C was upregulated ( Figure 7D). PP2C interacted with SnRK2, and ABA inhibited the reaction-the feedback regulation for ABA and PP2C. It can be speculated that miRNAs may be involved in balance the ABA level and metabolism, and more studies will be needed in the future to address this issue.
miRNA Involved in the Crosstalk Between ABA and Other Phytohormone Overexpression of Cb5s confers lower ethylene sensitivity (Chang et al., 2014). ABA negatively regulated ethylene production (Dong et al., 2015). The two stress-induced hormones ABA and ethylene interacted each other in stomata movement (Wilkinson and Davies, 2010). Cb5 was predicted target of peu-miR408-5p, which increased significantly after ABA treatment 1 day, and Cb5 was increased in expression ( Figure 7E). As Cb5 may affected ethylene signaling, suggesting that miRNAs may facilitate another way for ABA to inhibit ethylene signaling. Auxin response factor (ARF), the repressor of indole-3-acetic acid (IAA), decreased in expression in response to ABA. Exogenous IAA increased sensitivity to ABA in Arabidopsis. Furthermore, overexpression of miR160 reduced sensitivity to ABA (Liu et al., 2007). MiR160 appears to promote auxin activity by suppressing the levels of the ARF (Nizampatnam et al., 2015). ARF was the target of miR160 in plant (Liu et al., 2007;Turner et al., 2013;Damodharan et al., 2016;Tian et al., 2016). Therefore, miR160 targeted ARF to decrease ARF-mediated IAA-induced ABA hypersensitivity. In this study, six targets of peu-160b-5p involving ARF were downregulated by ABA at 4th day, corroborating the findings of previous studies. In general, miRNA may involves a very complex set of the crosstalk between ABA and other phytohormone.

miRNA Regulated SNARE Interactions in Vesicular Transport
In this study, peu-miR6421 expression was first increase and then decreased to response to ABA. Sec22, a synaptobrevin, was one of the predicted target genes that participated in soluble N-ethylmaleimide sensitive factor attachment protein receptor (SNARE) interactions in the vesicular transport pathway (Supplementary Figure 2). In addition, Gos1 and VAMP7 were the predicted targeted of miR482a.1 and miR475a-5p, respectively. They were both upregulated by ABA ( Figure 7G) and also involved in the vesicular transport pathway. Sec22 played an essential role in early secretory traffic between the ER and the Golgi (El-Kasmi et al., 2011). In Arabidopsis, VAMP721 and VAMP722 protein levels were downregulated by ABA, leading to slow down plant growth; (Yi et al., 2013). A reasonable speculation was that this inhibited expression by miRNA, so further study was needed.

miRNA Involved in Stress-Related Genes Regulation
Di19s were upregulated by the supply of ABA (Li et al., 2010b), and it induced sensitivity to ABA and tolerance to stress in Arabidopsis and rice (Li et al., 2010a;Qin et al., 2014;Wang et al., 2014). Di19 was one predicted targets of miR-n87, which was reduced by ABA. While Di19 was induced ( Figure 7D) to improve stress resistance, an assumption was consistent with that of pervious reports. Spermine was a part of polyamines. Free spermine accumulation was showed as a particular metabolic feature of being under long-term salt stress (Maiale et al., 2004). Meanwhile polyamines had been shown to be an important part of plant responses to improve stress resistance (Shi and Chan, 2014;Tiburcio et al., 2014). Spermine synthase was a predicted target of miRNA408-5p, which was significant upregulated by ABA at 1 day, while spermine synthase was inhibited ( Figure 7E), which resulted in free spermine accumulation to affect stress resistance. But more evidence was needed to support the interaction between miRNAs and the predicted targets.

CONCLUSIONS
We constructed nine sRNA libraries based on P. euphratica leaves for high-throughput sequencing. In total, 151 unique mature sequences belonging to 75 conserved miRNA families were identified. Meanwhile, 56 novel miRNAs of 94 sequences were discovered. Among them, the expression levels of 31 conserved miRNAs belonging to 22 families were significantly different. Confirmed by real-time qPCR, the expression profiles of miRNAs and their predicted target genes were complementary. Based on function analysis, we suggest may play critical roles in maintaining a high photosynthetic ability to facilitate adaptation to stress. And involved several pathways and cellular processes that help this plant to cope with stresses. How individual genes participate in stomatal closure, photosynthesis and other processes involves a very complex set of mechanisms. Our results provide a foundation for further analyses of plant miRNA responses to ABA, and provide new insight into the mechanism underlying the role of ABA in the abiotic stress response and other biological processes, in P. euphratica.

AUTHOR CONTRIBUTIONS
HD, XX, and WY designed experiments; HD, XL, and CL carried out experiments; HD and YA analyzed experimental results. HD and XX wrote the manuscript.

ACKNOWLEDGMENTS
This research was supported by grants from the Ministry of Science and Technology of China (2015BAD07B01, 2013AA102701), the National Natural Science Foundation of China (31270656, 31570308), Program for Changjiang Scholars and Innovative Research Team in University (IRT13047), the 111 Project (B13007). We would like to thank Prof. Jianquan Liu of Lanzhou University for providing P. euphratica genomic annotation information.