Identification of F1 plants with different drought stress tolerances
We identified six F1 progenies (0–293, 0–106, 0–334, 7–275, 8–134, and 0–330) as drought-tolerant plants and four (0–297, 0–131, 1–338, and 8–11) as drought-sensitive plants (Fig. 1 and Table 1) based on their recovery rate after rehydration. Plants with recovery rates >50% were defined as drought-tolerant plants and those with recovery rates <50% were defined as drought-sensitive plants.
The phenotype of F1 progenies of R3 (M. × domestica) × M. sieversii is depicted after drought treatment and rehydration. F1 progenies with drought-sensitive or drought-tolerant plants were propagated by tissue culture and transplanted to soil for an additional 10 months. The plants were then treated with drought stress for 12 days and rehydrated for 30 days
Overview of sRNA sequencing
Agarose electrophoresis analysis suggested a decent quality of the extracted RNA for each biological replicate of F1 progeny under control or drought treatment (Supplementary Fig. S1a). Analysis by agarose electrophoresis and Agilent 2100 also showed a decent quality of the pooled RNA (Supplementary Fig. S1b and c). In addition, the quality of all libraries was good enough for downstream analysis (Supplementary Table S1). We sequenced a total of 368.88 million (M) raw and 356.58 M clean reads, with a mean of 29.71 M (97%) clean reads per library (Table 2). For further analysis, we isolated 264.21 M clean reads with lengths ranging from 18 to 30 nt (74%) and 67.21 M unique reads (25%) (Supplementary Table S1). All libraries had a similar distribution of lengths (~72% for redundant reads and 84% for unique reads, Fig. 2). The 24 nt reads were the most abundant in redundant reads (34%) and unique reads (52%), followed by the 20 nt redundant reads (14%) and 23 nt unique reads (10%) (Fig. 2). Among the 18–30 nt sRNA, 83% redundant reads and 71% unique reads were successfully mapped to the apple reference genome (Supplementary Table S1).
Identification of known and novel miRNAs in apple
To identify novel miRNAs in apple plants, a total of 1.3 M redundant reads (0.57%) were mapped onto 274 hairpin structures, identifying 253 candidates, among which 123 had miRNA* sequences (Supplementary Tables S2 and S4). As miRNA* is a common tool to verify new miRNAs46, we designated the 123 candidates with miRNA* sequences as novel miRNAs and the remaining 130 without miRNA* sequences as novel miRNA candidates. The lengths of the 253 novel apple miRNAs ranged from 18 to 24 nt, with a dominant distribution of 24 nt (211 miRNAs), followed by 21 nt (18 miRNAs), 22 nt (12 miRNAs), 23 nt (4 miRNAs), 18 nt (3 miRNAs), 19 nt (3 miRNAs), and 20 nt (2 miRNAs) (Supplementary Table S4).
The identified known and novel miRNAs were classified into 40 families. Among these families, the largest was miR156, with 31 members, followed by miR172 (15 members), miR171_1 (14 members), and miR167_1 (10 members). Only two novel miRNAs were classified into miRNA families: miRn234 into the miR159 family and miRn85 into the miR169_2 family (Fig. 3 and Supplementary Table S5).
In addition to the known and novel miRNA tags, small RNA reads were also mapped to rRNA, tRNA, snRNA, snoRNA, repeat sequences, NAT-siRNA, and genes. The most abundant redundant reads were mapped onto genes (17%), followed by repeat sequences (16%), rRNAs (8.1%), and NAT-siRNAs (7.2%). Unique reads were also most abundant in repeat sequences (26%), followed by genes (7.8%), rRNAs (1.3%), and NAT-siRNAs (1.2%). Among unique and redundant reads, the minimum RNA classes were snoRNA (0.11% and 0.22%, respectively), snRNA (0.080% and 0.10%, respectively), and tRNA (nearly 0% for both) (Supplementary Table S2).
Expression analyses of the miRNAs in apple
The miRNAs identified in drought-sensitive or drought-tolerant plants in each treatment were classified into six categories of expression based on normalized read counts: no (0 reads), extra-low (0–10 reads), low (10–100 reads), moderate (100–1000 reads), high (1000–10,000 reads), and extra-high (>10,000 reads) expression (Fig. 4). Highly similar percentages of identified miRNAs fell into the same category in each treatment (Fig. 4a). The largest percentage of known miRNAs fell into the low category, followed by the moderate and extra-low categories, whereas the lowest percentages fell into the no expression category (Fig. 4b and Supplementary Table S6). MiR1511 was the most abundant miRNA, with over 60,000 normalized reads in each treatment. The majority of novel miRNAs had low abundances in each treatment (Fig. 4c and Supplementary Table S6). For example, miRn-54 had 0.75 normalized reads in a library of untreated control drought-tolerant plants. Moreover, different members of the same miRNA families displayed different expression levels. For example, miR156h had more than 300 normalized reads, but miR156t had <10 (Supplementary Table S6). Notably, 20 miRNAs were absent in at least one treatment, including 12 known and 8 novel miRNAs (Supplementary Table S6). In drought-tolerant plants, miR3627 a/b/c, miRn-50, miR-61, and miRn-158 were not present under drought treatment, whereas miR399 e/f/g/h, miRn-67, and miRn-203 were not detected under control conditions. In drought-sensitive plants, miR828a and miR828b were absent under both conditions, whereas miRn-11 and miRn-78 were not present under drought. In addition, miR169c and miR169d were only detected in drought-sensitive plants under drought, whereas miR7128 and miRn-54 were only detected in drought-tolerant plants under control conditions (Supplementary Table S6).
The expression levels of miRNAs in drought-tolerant or drought-sensitive plants under control or drought treatment are provided in Figs. 5 and 6, and Supplementary Table S7. Thirty-nine percent of miRNAs displayed significant expression levels (p-value < 0.05) and 93 differentially expressed miRNAs were identified in the drought response. After drought treatment, 21 miRNAs were downregulated and 14 miRNAs were upregulated in drought-sensitive plants, whereas 32 miRNAs were downregulated and 29 miRNAs were upregulated in drought-tolerant plants (Fig. 5a, b). Among these 93 differentially expressed miRNAs, only 3 novel miRNAs (miRn15, miRn85, and miRn246) were responsive to drought in both plant groups (Fig. 6 and Supplementary Table S7). Specifically, miRn15 was upregulated and miRn85 was downregulated in both plants in response to drought, whereas miRn246 was downregulated in drought-tolerant plants but upregulated in drought-sensitive plants in response to drought. Moreover, drought treatment induced most novel miRNAs (76%) and repressed most known miRNAs (88%) in drought-tolerant plants, but upregulated (71%) and downregulated (71%) most known miRNAs in drought-sensitive plants (Fig. 5d and Supplementary Table S7).
a Distribution of miRNA differences among the four comparison groups. b Differentially expressed miRNAs induced or repressed by drought. c Differentially expressed miRNAs in drought-tolerant vs. drought-sensitive plants. d Drought-induced or -repressed novel and known miRNAs in both phenotypes. e High abundant or low abundant miRNAs in drought-tolerant plants under drought treatment or control conditions compared with drought-sensitive plants. S-C: drought-sensitive plants under control conditions; S-D: drought-sensitive plants under drought treatment; T-C: drought-tolerant plants under control; T-D: drought-tolerant plants under drought treatment. Upward arrows indicate upregulation or higher abundance and downward arrows indicate downregulation or lower abundance
When comparing drought-sensitive and drought-tolerant plants, 151 miRNAs were differentially expressed under treatment, including 119 miRNAs under control and 67 miRNAs under drought (Fig. 5a, c). These 67 differentially expressed miRNAs represent a unique response of miRNAs in drought-tolerant plants and may be important miRNAs for improving the drought tolerance of apple rootstocks. Among these 67 miRNAs, 35 were differentially expressed under both normal and drought conditions, including 23 lower abundance and 12 higher abundance miRNAs in drought-tolerant plants (Fig. 5c). Among these miRNAs, miRn-157, miRn-158, and miRn-101 were classified as lower abundance in drought-sensitive plants, whereas miR395 a/b/c/d/e/f/g/h/i were classified as higher abundance (Supplementary Table S7).
In addition, the expression levels of 22 miRNAs that were similar in both plant groups under normal conditions were either reduced or increased after drought treatment, including miR156 aa/p/q/r/s/x/y/z, miR390 a/b/c/d/e/f, miRn-134, miRn-249, and miR5225 a/b (Supplementary Table S7). Of these 22 miRNAs, nine were responsive to drought in drought-sensitive plants and 13 were responsive to drought in drought-tolerant plants (Figs. 5a and 6).
Under control conditions, 69 miRNAs displayed higher abundance in drought-tolerant plants, whereas 50 miRNAs displayed lower abundance in drought-tolerant plants (Fig. 5c). In drought-tolerant plants relative to drought-sensitive plants, 35 miRNAs were higher in abundance, whereas 32 miRNAs were lower in abundance after drought treatment (Fig. 5c). Under both normal and drought conditions, most known miRNAs were more abundant (91% in control and 74% in drought treatment), whereas most novel miRNAs were less abundant (90% in control and 72% in drought treatment) in drought-tolerant plants than in drought-sensitive plants (Fig. 5e and Supplementary Table S7).
Target prediction and GO analysis
We predicted 2754 miRNA target pairs, including 1090 unique protein-coding mRNAs. Among the 1090 protein-coding mRNAs, the majority of genes encoded transcription factors (SPL transcription factors, C3HC4-type RING finger), resistance-associated proteins (TIR-NBS-LRR class disease resistance protein), transport proteins (sulfate transporter protein), and enzymes, such as kinases (receptor kinase), transferases (S-adenosyl-l-methionine-dependent methyltransferases), and phosphatases (RNA polymerase II C-terminal domain phosphatase) (Supplementary Table S8).
GO analysis of 1090 protein-coding mRNAs suggested that enriched biological processes included cellular process (GO:0009987), metabolic process (GO:0008152), and response to stimulus (GO:0050896), whereas binding (GO:0005488) and catalytic activity (GO:0003824) were the most abundant classes in the molecular function category (Fig. 7). The most significant cellular components were cell (GO:0005623), cell part (GO:0044464), and organelle (GO:0043226).
The x-axis is the GO term targeting gene enrichment, whereas the y-axis represents the percentage of target genes mapped by the term. The query list (yellow bar) represents the percentage of target genes mapped to the GO term in all genes of the input list. The background/reference list (dark bar) represents the percentage of all GO term genes from all reference database genes
We also evaluated the potential functions of 67 miRNA-targeted genes that were differentially expressed between both plant groups under drought. We selected a total of 1213 miRNA target pairs for GO enrichment analysis, including 390 protein-coding mRNAs. Hierarchical graphs displayed significantly enriched GO terms in biological process and molecular function categories, rather than the cellular component. Within the biological process category, programmed cell death (GO:0012501), defense response (GO:0006952), postembryonic development (GO:0009791), regulation of cellular process (GO:0050794), and biological regulation (GO:0065007) were significantly enriched (Fig. 8). Within the molecular function terms, receptor activity (GO: 0004872), ubiquitin-protein ligase activity (GO: 0004842), transferase activity (GO:0016740), ion binding (GO:0043167), and protein binding (GO: 0005515) were significantly enriched (Fig. 8).
Expression verification of miRNAs and their targets by qRT-PCR
We used stem-loop qRT-PCR to validate sequencing results by determining miRNA expression levels in both plant groups under control and drought treatments. From the high-throughput sequencing results, we selected eight miRNAs (four conserved and four novel) for stem-loop qRT-PCR, including miR156, miR395, miR408a, miR5225, miRn-101, miRn-157, miRn158, and miRn-249 (Fig. 9). Primer analysis suggested that all stem-loop qRT-PCR primers had high efficacy and were thus suitable for miRNA expression detection (Supplementary Fig. S2a and b). The stem-loop qRT-PCR results confirmed the expression of seven of the eight miRNAs (with the exception of miR156), indicating the reliability of our small RNA-seq data. The expression of miR156, miR395, and miR408a was much higher in drought-tolerant plants than in drought-sensitive plants under control or drought conditions. MiR408a expression was induced by drought in both plant groups, whereas miR156 and miR395 were not. Under control conditions, the expression of miR5225 and miRn-249 was similar between drought-tolerant and drought-sensitive plants; however, under drought conditions, the increased expression of both miRNAs in drought-tolerant plants was observed. In addition, the expression of both miR5225 and miRn-249 was only induced by drought in drought-tolerant plants. MiRn-101, miRn-157, and miRn-158 had similar expression patterns with lower expression in drought-tolerant plants under both control and drought conditions (Fig. 9). The expression pattern of miR156 by stem-loop qRT-PCR was not consistent with the small RNA-seq results, possibly because stem-loop qRT-PCR detected all miR156 family members (Figs 6 and 9).
The expression of four known miRNAs (miR156, miR395, miR408a, and miR5225) and four novel miRNAs (miRn-101, miRn-157, miRn-158, and miRn-249) was verified by stem-loop qRT-PCR, with MDH (malate dehydrogenase) as an internal reference gene. MDH reverse primer and miRNA-specific reverse transcription primer were mixed, and the mixed primer and total RNA (same as sRNA-seq) were used for cDNA synthesis. Error bars indicate ±SD (n = 3, from three technical replicates). Three biological replicates were performed. Student’s t-test was performed for data analysis and statistically significant differences were indicated by *p-value < 0.05, **p-value < 0.01, and ***p-value < 0.001
We next analyzed the expression patterns of miRNA-targeted genes by using high-efficiency primers (Supplementary Fig. S3a and b), including seven miR156 targets, three miR395 targets, and one novel miRNA (miRn-249) targeted gene. The qRT-PCR results demonstrated that three predicted target genes of miR395 (RNA polymerase II C-terminal domain phosphatase, WRKY33, and receptor kinase) had expression patterns opposite that of miR395 in both plant groups under control and drought conditions. This finding indicated a negative correlation with miR395, supporting the potential of these genes as miR395 targets (Fig. 10). Zinc finger (C3HC4-type RING finger) was predicted as the target of miRn-249 under drought, with an expression pattern that negatively correlated with that of miRn-249 in both plant groups under drought. The expression patterns of SPL genes were more complicated. SPL_4 and SPL_6a displayed expression patterns opposite that of miR156 under control or drought, whereas the expression of SPL_5 and SPL_13a was only negatively correlated with that of miR156 transcripts under drought. Thus, both SPL genes may only be targeted by miR156 under drought conditions. SPL_2, SPL_9, and SPL_12 had the same expression pattern as miR156 under control and drought, indicating that these three genes may not be the targets of miR156 in apple rootstocks (Fig. 10).
Eight miR156 target genes (SPL_2, SPL_4, SPL_5, SPL_6a, SPL_9, SPL_12, and SPL_13a), three miR395 target genes (RNA polymerase II C-terminal domain phosphatase, WRKY33, and receptor kinase), and one novel miRNA (miRn-249) target gene (Zinc finger (C3HC4-type RING finger)) were analyzed by qRT-PCR. Oligo-dT primers and total RNA (same as sRNA-seq) were used for cDNA synthesis, with MDH (malate dehydrogenase) as an internal reference gene. Error bars indicate ±SD (n = 3, from three technical replicates). Three biological replicates were performed. Student’s t-test was performed for data analysis and statistically significant differences were indicated by *p-value < 0.05, **p-value < 0.01, and ***p-value < 0.001
Function verification of miRNAs in response to osmotic stress
To verify the function of these miRNAs, we selected two miRNAs, miR156p and miRn-249, and generated transgenic calli carrying 35S:MIR156p or 35S:MIRn-249. As drought stress could cause osmotic stress in plants, we treated the WT and transgenic calli with 150 mM mannitol. Stem-loop qRT-PCR analysis demonstrated that miR156p was upregulated ~400-fold in the transgenic calli compared with that in WT calli and 10-fold for miRn-249 in miRn-249 transgenic calli (Fig. 11a). When treated with 150 mM mannitol for 3 weeks, the relative growth rate of miRn-249 and miR156p transgenic calli was higher than that of WT (Fig. 11b), suggesting positive roles for these two miRNAs in osmotic stress tolerance. We also detected the expression levels of the target genes of these two miRNAs. Our results showed that four target genes of miR156 (SPL_6a, SPL_9, SPL_12, and SPL_13a) were repressed by elevated miR156p expression (Fig. 11c). However, the expression of SPL_2 and Zinc finger (C3HC4-type RING finger), which were predicted as the targets of miR156 and miRn-249, respectively, was not reduced in miR156p or miRn-249 transgenic calli (Supplementary Fig. S4).
a Expression levels of miR156p or miRn-249 in overexpressing calli under control and 150 mM mannitol treatment for 6 h. b Morphology of the calli after culture for 3 weeks on MS medium (upper) or MS medium supplemented with 150 mM mannitol (middle), and relative growth rate between them (under, error bars indicate ±SD, n ≥ 10). c Expression levels of the miR156p target genes under control and 150 mM mannitol treatment for 6 h. Stem-loop qRT-PCR and qRT-PCR used MDH (malate dehydrogenase) as an internal reference gene, error bars indicate ±SD (n = 3, from three technical replicates) and three biological replicates were performed. Student’s t-test was performed for data, and statistically significant differences were indicated by *p-value < 0.05, **p-value < 0.01, and ***p-value < 0.001