Experiments were conducted in 2013 and 2014 using standard-type “Nagafu 2” and spur-type mutant “Yanfu 6” apple (Malus domestica) trees that were grafted on Malling “M.26” rootstocks and planted in 2009. Trees were located in the experimental plots of the Yangling National Apple Improvement Center, Yangling, China (34.31° N, 108.04° E), trained and managed using standard horticultural practices. The meteorological conditions of the experiment site in 2013 and 2014 are described in Figure S1. Three biological replicates consisting of five trees of each genotype were used, for a total of 15 trees.

Apical portions (5–8 mm below and including the shoot meristem) of YF and CF shoots (without any leaves or petiole tissue), were collected at 45 days after bud break (DABB) in 2013 and used for deep miRNA sequencing. Similar apical portions of shoots were also collected at 65, 85, 105, 125, and 145 DABB in 2013 for miRNA and mRNA expression analysis. Approximately 30 shoot samples of each biological replicate of YF and CF were collected. Leaves near the shoot tips were also collected for plant hormone measurements. All samples were immediately frozen in liquid nitrogen and stored at −80°C for subsequent RNA and hormone extraction.

Six individual shoots on each of six trees were marked in each genotype to assess shoot elongation on a weekly basis in 2013 and 2014. The number of internodes present on each shoot were also recorded. Internode lengths were estimated by dividing the overall shoot length by the number of internodes.

Approximately 0.5 g of leaf samples were used to measure IAA, GA, ZR, and ABA. The plant hormones were extracted using the method described by Dobrev and Kamınek (2002) and Dobrev et al. (2005), quantifed by the High Pressure Liquid Chromatography (Waters, USA) and enzyme-linked immunosorbent assay.

YF and CF trees were sprayed with 500 mg·L⁻¹ GA₄₊₇ at the end of April 2014. Shoot tips, as previously described, were subsequently sampled at 0, 2, 4, 6, and 8 days after spraying. Both shoot length and internode length were also recorded.

Total RNA was isolated using the CTAB method (Chang et al., 1993). After ligating adaptors to the 5′ and 3′ RNA, the purified RNAs were reverse transcribed using primers with partial complementarity to the adaptors. The DNA pool amplified from the first-strand cDNA was subsequently sequenced using a Hiseq 2000 (Illumina, San Diego, CA) sequencer located at the Beijing Genomics Institute (BGI), Shenzhen.

The YF sample for miRNA sequencing was also subjected to degradome sequencing in order to better determine the miRNA/mRNA pairing relationship. For the construction of the degradome libraries, a 5′ RNA adapter with a Mme I recognition site was ligated to the 5′ terminus of purified polyA RNA. The 5′ ligated products were reverse transcribed using five PCR cycles, digested with Mme I, and ligated to a 3′ double DNA adapter. Lastly, the ligated products were amplified by PCR and gel-purified for deep sequencing.

The miRNA and degradome sequencing raw data were available at NCBI SRA (BioProject ID: PRJNA361540).

The original sequence data from the sequencing machine were used for basic analysis. Filtering the low quality reads and contaminants sequence (adaptor, polyA, reads shorter than 18 nt, insert tag, reads without 3′ primer) and discarding the remaining reads with lengths smaller than 18 nt, the clean reads were used for the advanced analysis. The length distribution of the clean reads were summarized to analyze the compositions of small RNA sample. Unique tags and total tags existence among the sample were summarized to verify whether the uniformity of the two samples on the whole of sequencing was good.

The small RNA tags were mapped to the apple reference genome at GDR (http://www.rosaceae.org) (Jung et al., 2014) using SOAP software to analyze their distribution in the genome and their expression levels. All unique, clean reads, particularly the non-redundant reads, were subjected to further analysis. The clean tags were queried against the known apple miRNAs in the mirBASE 18.0 (http://www.mirbase.org) database in order to identify known microRNAs based on the following criteria: (1) the tags aligned to a miRNA precursor in the mirBASE database with no mismatch, and; (2) based on the first criteria, the tags aligned to mature miRNA in the mirBASE database with at least 16 nt overlap, allowing offsets. The unmatched unique reads in each sample were screened in Genbank in order to annotate non-coding RNA sequences (rRNA, scRNA, snoRNA, snRNA, and tRNA). Differentially expressed miRNAs were identified by comparing the sequencing reads between CF and YF and determining statistically significant differences (P < 0.05 and P < 0.01). miRNAs with more than one normalized read were analyzed by calculating their fold-change and p values. miRNAs with p < 0.05 and fold changes (log₂) less than −1 or greater than 1 were considered to have significantly altered expression.

The psRNATarget server (http://plantgrn.noble.org/psRNATarget/) with default parameters was used to predict putative targets of the identified apple miRNAs. All of the identified apple miRNA were used as a query against “Malus domestica (apple), predicted consensus gene set CDS sequences, version 1” database located on the psRNATarget server.

For the degradome sequencing raw data about 50 nt length obtained from HiSeq sequencing, go through steps likely adaptor, pollution, low quality trimmed to obtain “clean” sequence. After a series of annotation of the degradated sequence, the degradation mRNA obtained go through comparison with miRNA to find mRNA-miRNA pairing.

GO enrichment and KEGG pathway analysis for all of the annotated miRNAs and their targets, using the apple reference genome background, were conducted in order to predict the function of the miRNA and their targets. Blast2GO was employed to store information from the GO and KEGG pathway databases. All of the miRNA targets were queried against the GO protein database (http://www.geneontology.org/) using BLASTX. A combined query was used in order to complete the GO annotation and pathway analysis against the GO and KEGG databases.

Stem-loop RT-qPCR method was used to quantify miRNA expression (Chen et al., 2005). mRNAs of the target genes were also quantified. 5S rRNAs were used as an internal reference of miRNA expression. MdActin and MdEF were used as internal references for mRNA expression. The primers utilized are listed in Tables S12, S13. miRNA and mRNA relative expression were calculated using the comparative ΔΔCT method (Livak and Schmittgen, 2001).

The data were analyzed using Data Processing System (DPS, version 7.05; Zhejiang University, Hangzhou, China) software. Significant differences were determined using a Student's t-test at (p ≤ 0.05).

The shoot length in YF was significantly shorter than in CF (Figures 1A,E). YF also had shorter internodes, that were approximately 66% shorter than CF (Figures 1C,D). Most YF shoots stopped growing between 80 and 105 DABB, which was earlier than CF shoots. After 105 DABB, any continued growth of YF shoots extended their length with a section of shorter internodes.

IAA in YF leaves near the shoot tips was higher than CF leaves at 45 DABB. In contrast, IAA levels in YF were significantly lower than CF at 125 and 145 DABB (Figure 2A). The GA content in YF leaves remained at a consistently low level, relative to CF, throughout the entire growth period (Figure 2B). The ABA content in YF was significantly lower than CF at 65 and 85 DABB (Figure 2C). The ZR content in YF was significantly lower than CF at 45, 65, 125, and 145 DABB (Figure 2D). The changes in the levels of the growth promoting hormones, IAA, GA, and ZR, were in accordance with the growth of shoot, namely, exhibiting a high content in the level during the rapid growth period, and a low level as the growth period waned (Figure 2). The stress hormone ABA exhibited the opposite trend, exhibiting a low level during the active growth period and a high level as growth slowed and terminated.

Two cDNA libraries of small RNAs from CF and YF were sequenced in order to identify miRNAs involved in apple shoot development. The sequencing produced a total of 13,264,874 raw reads from CF and 19,773,007 from YF (Table S1). After removal of the adaptor sequences, filtering out low-quality tags, and omitting contamination resulting from adaptor–adaptor ligations, 13,083,017 (98.95% of the total) CF and 19,580,358 (99.41%) YF clean reads remained, comprising 5,872,194 and 7,680,257 CF and YF unique sequences, respectively (Table S2). The unique sequences were mapped to the apple genome assembly using SOAP, resulting in 3,674,249 (62.57%) and 4,473,261 (58.24%) genome-matched reads for CF and YF, respectively. A total of 143,629 (2.44%) and 278,103 (2.54%) exon RNAs (sense and antisense), 267,025 (5.08%) and 321,431(4.18%) intron RNAs (sense and antisense), 130,796 (2.23%) and 204,833 (2.67%) rRNAs, 1,973,617 (33.61%) and 2,370,059 (30.86%) repeat regions, 4733 (0.08%) and 6054 (0.08%) snRNAs, 1,238 (0.02%) and 1,451 (0.02%) snoRNAs, 11,364 (0.19%), and 15,596 (0.20%) tRNAs, were removed from the CF and YF reads, respectively, This left 3,338,736 (56.86%) and 4,564,590 (59.43%) unannotated CF and YF RNAs, respectively, and a total of 1,056 (0.02%) and 1,149 (0.01%) miRNA reads for CF and YF, respectively. These were utilized as miRNA candidates in subsequent analyses. The read lengths in the two small RNA libraries were assessed (Figure 3A).The lengths ranged of the designated small RNAs ranged from 21 to 24 nt and accounted for more than 90% (approximately 92.9% and 94.6% of the CF and YF sequences, respectively) of all the clean reads in each library. The 24 nt small RNA was the most abundant small RNA length (about 50%). Similar results were reported in apple leaves, flowers, fruit, and roots (Xia et al., 2012), peach (Prunus persica) (Zhu et al., 2012), pear (Pyrus bretschneideri) (Wu et al., 2014), tomato (Lycopersicon esculentum) (Gao et al., 2015), and chickpea (Jain et al., 2014). The next most abundant small RNAs were 21 nt > 22 nt > 23 nt. There was a greater number of 23 nt and 24 nt small RNA in CF than in YF. The number of 21 and 22 nt small RNA were greater in YF than in CF.

The clean reads in each library were queried against the Malus domestica microRNA database in miRBase in order to identify the known miRNAs that were present in CF and YF apple shoot tips. Collectively, one million reads from the two libraries were perfectly matched to the Malus domestica microRNA database. A total of 202 unique sRNAs sequences were assigned to 41 microRNA families (Table S3, Figure 3B).

A comparison of family members revealed that the number of miRNA members varied among the different miRNA families (Figure 3C). The top five miRNA families, mdm-miR156, mdm-miR172, mdm-miR171, mdm-miR167, and mdm-miR399, had more than 10 members. For example, the mdm-miR156 family had 31 members. mdm-miR1511, mdm-miR391, mdm-miR7125, mdm-miR7126, mdm-miR827, and mdm-miR858 only had one member. miRNA156 is conserved in embryophytes (Zhang et al., 2006; Taylor et al., 2014) and represents the largest microRNA family in pear (Wu et al., 2014), which like apple, is a member of the Rosaceae. The differences in the number of members in the different microRNA families may be due to the different evolutionary rates in the microRNA families and to genome duplication events. Two known microRNAs families, mdm-miR7123, and mdm-miR7128, were not found to be present in the libraries (Figure 3C), therefore, these microRNAs may not be expressed in apple shoot tips. Statistical analysis of the read counts for the known miRNA families indicated significant differences in the level of expression among the different miRNAs (Table S3). Read counts of some miRNAs in a family were nearly the same, but some were obviously different. For example, all of the mdm-miR395 family members had the same number of counts. In contrast, mdm-miR172a, mdm-miR172o, and mdm-miR172g-h, members of the mdm-miR172 family (15 members), had a high read count, while mdm-miR172 a-c, mdm-miR172i-k, and mdm-miR172 m-o had a moderate read count, and mdm-miR172l had a low count. Perhaps some of the microRNA clustered in the same area of the genome, and thus had the same cis-regulation elements.

The read counts of about 30% of the known miRNAs were greater than 1,000 (Figure 3D). Another 30% of the known miRNA had read counts between 100 and 1,000. Due to high level of expression of miRNA in the apple shoot tip tissue, these miRNAs could be accurately identified. Approximately 20% of the known miRNA had read counts between 0 and 20.

Among the miRNA in YF and CF shoot tips, mdm-miRNA166 families exhibited the highest level of expression, with approximately 10,000 standardized reads (Figure 3E). This was followed by mdm-miR482a-5p, which had about 4000 standardized reads. The juvenile to adult phase transition related microRNAs, mdm-miR156 and mdm-miR172, and the flowering related, mdm-miR535, mdm-miR168, and mdm-miR167, also showed high expression levels in apple shoot tips.

After the identification of the known miRNAs, the unannotated tags were analyzed using MIREAP software in order to identify novel miRNAs based on the characteristics of the secondary hairpin structure of miRNA precursors. Using the hairpin structures of the precursors, 498 novel miRNAs were identified. The 23 nt length novel miRNAs represented an overwhelming amount 76.8% of the 306 total novel miRNAs identified (Table S4).

The novel miRNA was used a query to blast the miRBASE 18.0. Results indicated that 11 of the 498 novel miRNAs possessed a high significant similarity to known miRNAs from rice, wheat, and other plant species (Figure 4). Some varied only in nucleotide length, which may be due to differential cleavage sites in the pre-miRNA. These miRNAs represent new members of conserved miRNA families, and as a result, were designated as novel miRNA isoforms. The novel miRNAs would be species-specific Malus domestica miRNA.

Lower negative minimal folding free energy (MFE) value for RNA secondary structure, indicates greater stability (Bonnet et al., 2004). The average negative MFE was −51.03 kcal/mol (Table S4). This MFE value is similar to the data for Arabidopsis thaliana miRNA precursors (−57 kcal/mol) (Debat and Ducasse, 2014), chickpea (−57.58 kcal/mol) (Jain et al., 2014), maize (−61.15 kcal/mol) (Li et al., 2013), and pear (−51.48 kcal/mol) (Wu et al., 2014).

The GC content of miRNA has been shown to be related to biological function (Mishra et al., 2009). The GC content of 90% of the novel miRNA ranged between 30 and 70% (Table S4) with an average value of 42%. This average GC is similar to Arabidopsis (44%), pear (45%), and chickpea (44%), but less than grape (50%) and the GC content of the known miRNA of apple (50%).

The U residue is the most prevalent nucleotide at the first 5′ nucleotide site of miRNA, especially 21 nt miRNA (Chen, 2005). About 46% of the 5′ end nucleotides of the 21 nt novel miRNAs were U (Figure S2).

The abundances of miRNA reads in the CF and YF libraries were compared in order to identify miRNAs involved in shoot development and to determine their expression pattern in the two different genotypes with contrasting phenotypes. A total of 42 miRNAs, belonging to 12 known miRNA families, were identified whose expression differed significantly between the CF and YF shoot tips (Figure 5A). Among these, only mdm-miR5225c and mdm-miR7125 were upregulated in YF. In contrast, five members of mdm-miR164, three members of mdm-miR171, six members of mdm-miR172, three members of mdm-miR393, nine members of mdm-miR395, two members of mdm-miR396, six members of mdm-miR399, two members of mdm-miR5225, two members of mdm-miR7124, and mdm-miR858 were all downregulated in YF shoot tips.

Among of the 498 novel miRNAs, the expression of 93 differed significantly between CF and YF. A total of 34 upregulated and 59 were downregulated in YF (Figure 5B). A total of 49 novel miRNAs were expressed only in CF while 28 novel miRNAs were expressed only in YF (Table S4).

Plant miRNAs pair with their complementary target and inhibit translation or degrade the target mRNA. The complementary pairing of the miRNA/mRNA duplex provides a useful approach for target prediction. In order to better understand the role of apple miRNAs in shoot development, psRNATarget software was used to predict potential miRNA target genes. Using previously described criteria, 689 potential targets for 168 of the known miRNAs, and 3,663 targets for the 466 novel miRNAs were predicted (Tables S5, S6). The number of predicted targets for the various miRNAs varied from one to as many as 20, suggesting that these miRNAs may have diverse biological functions. On the other hand, a single gene may also be targeted by several miRNAs. MDP0000237964 encodes a LAC gene that is targeted by mdm-miR397 at two different sites, and two similar CuRO_LCC_plant domains (Figure 6A). miRNA828 and miRNA858 target the MYB family at the R3 domian (Xia et al., 2012; Figure 6B). MDP0000147309 is likely to be targeted by mdm-mi159 and mdm-mi319 (Figure 6C).

Degradome sequencing was conducted in order to better understand the role of miRNA-target regulation in shoot tips. Approximately 23 million total clean reads were generated from the degradome sequencing. About 49.5% of the unique reads mapped to apple mRNA (Table S7). The degradome sequencing results indicated that 246 genes are potentially regulated by 164 miRNAs, including both known and novel miRNAs. The main miRNAs and their target genes verified by degradome sequencing are listed in Tables 1, 2. Two-thirds of the degradome validated known miRNA target genes also identified in the target gene prediction software.

BLASTX querying against the protein database, Gene Ontology (GO), and KEGG pathway analysis were used to annotate the target genes identified in the degradome analysis in order to determine the potential functions of the miRNA target genes. A total of 69 degradome-identified target genes of 21 known miRNA were associated with 238 GO terms (Table S8). The first three enrichment GOs were GO:0005634, GO:0003700, and GO:0006351. GO:0005634, involved with the nucleus in the cellular component had 49 associated target genes. GO:0003700 involved in transcription factor activity in the molecular function component had 45 associated target genes. GO:0006351, involved in transcription, DNA-dependent in the biological process component had 41 associated target genes.

KEGG pathway analysis indicated that the target genes were involved in 21 pathways, in which, Transcription factors, ko03000; Metabolic pathways, ko01100, and; Plant hormone signal transduction, ko04075 were significantly enriched (Table S9).

Twenty-nine of the degradome identified genes targeted by 15 of the novel miRNAs were assigned to 176 GO terms (Table S10). Like the known miRNA, the GO enrichment in novel miRNA were transcription regulation related GO terms: 13 genes in GO:0005634; 9 genes in GO:0006355; 9 genes in GO:0003677, and; 7 genes in GO:0003700. KEGG pathway analysis indicated that the target genes were involved in 35 pathways, in which, Transcription factors, ko03000, and Metabolic pathways: ko01100 were significantly enriched. (Table S11).

YF exhibited a shorter duration of shoot elongation (Figure 1B). The duration of shoot growth is regulated by the shoot apical meristem (SAM). Based on the GO annotation of the target genes and with the KEGG pathway analysis, the potential targets of miR164, miR166, miR171, miR172, and miR482 are involved in meristem development (Table 3). In particular, miR164 participates in the formation of meristem-to-organ boundaries in the SAM (Laufs et al., 2004). The potential targets of miR166, miR159, and miR156 are involved in meristem phase transition. The potential targets of miR167, miR160, and miR159 play a role in auxin response and the auxin signaling pathway may also participate in SAM development.

Internode length is determined by cell division and cell elongation. Based on the GO annotation of the target genes and the KEGG pathway analysis, miRNAs related to cell division and cell differentiation were identified (Table 4). The potential targets of miR159, miR166, miR167, miR171, miR172, miR393, miR858, and miR828 are involved in cell growth. Gibberellin also plays an important role in regulating internode length. The potential targets of miR159 are involved in the gibberellic acid mediated signaling pathway and gibberellin biosynthetic process. While the potential targets of miRNA160 and miRNA167 regulate auxin response to control cell division. Thus, the potential targets of miR167, involved in response to brassinosteroid stimulus and the auxin-mediated signaling pathway, may play a role in internode development. Secondary metabolism affects primary growth by competing for common substrates. The potential targets of miR858, miR828, and miR156 participate in the fatty acid biosynthetic process. Cell walls also affect cell elongation and cell size. The potential targets of miR166 and miR167 are involved in cell wall biogenesis and cell wall organization. Mineral nutrition affects plant physiology and growth and the potential targets of miR169, miR166, miR7125, miR393, and miR395 are involved in mineral element response.

The miRNAs and their targets putatively involved in SAM development that were differentially expressed in YF vs. CF are illustrated in Figure 7. The expression analysis of miRNAs and their associated target genes revealed inversed expression patterns. The expression of miRNA164 decreased from 45 to 145 DABB in both YF and CF while expression of the potential target gene, MdNAC, increased from 45 to 55 DABB, decreased from 85 to 105 DABB, and then increased rapidly after 125 DABB. The expression of miRNA164 in YF was significantly lower than in CF at 45, 65, and 85 DABB. The expression of miRNA166 in CF and YF peaked at 105 DABB. The expression of miRNA166 in YF was significantly higher than in CF at 45, 105, and 125 DABB. The expression of miRNA171 in CF and YF increased from 45 to 65, decreased from 85 to 125, and then increased to a high level at 145 DABB. The potential target of miRNA171, MdHAM, fluctuated in its expression level from 45 to 105 DABB, and was then downregulated after 125 DABB. The expression of miRNA171 in YF was significantly lower than CF from 65 to 145 DABB.

MiRNAs and their targets, which are putatively involved in internode elongation, that were differentially expressed in YF vs. CF are shown in Figure 8. The expression of miRNA159 in YF was significantly lower at 65, 85, and 125 DABB than it was in CF. The potential target, MdMYB65, was downregulated from 45 to 145 DABB. The highest expression level of miRNA167 in both CF and YF was observed during the slow growth period at 80 DABB. The potential target, MdARF6, was downregulated from 45 to 125 DABB, and upregulated at 145 DABB. The expression of miRNA167 in YF was significantly higher than in CF at 65, 85, and 105 DABB. The expression of miRNA396 in YF peaked at 105 DABB. The expression of miRNA396 in YF was significantly higher than in CF at 105, 125, and 145 DABB. The miRNA396 potential target gene, MdGRF, exhibited a high level during the later growth stage. The expression of miRNA159, miRNA167, and miRNA396 in YF and CF shoot tips exhibited their highest level of expression during the period of slow shoot growth. These data suggest that these miRNAs inhibited cell division in developing shoots.

A strong correlation has been observed between final internode length and cell number (Brown and Sommer, 1992). Therefore, the expression of cell cycle and cell elongation related genes were examined (Figure 9). The expression of the cell cycle related genes, MdCYCD, MdCYCB1.2, MdCYCD1.3, and MdRBR1 exhibited high levels of transcript abundance during the period of rapid growth after bud break, and were downregulated during the period of slow growth, followed by a slight upregulation at the later stages of sampling. MdCYCD expression in YF was significantly lower than in CF from 45 to 125 DABB. Expression levels of MdCYCB1;2, MdCYCD1;3, and MdRBR1 genes in YF were significantly lower than in CF during the period of rapid growth just after bud break. The cell elongation genes, MdTCH4 and MdXTH23, exhibited high expression levels during the period of slow growth and in the autumn. MdTCH2 and MdXTH23 expression in YF exhibited a lower expression level than in CF during these periods.

Based on the observed effect of plant hormones on growth, the expression of hormone-related genes was also examined (Figure 10). The expression of the GA2ox gene, which catalyzes the catabolism of bioactive GA or their precursors, was significantly higher in YF than in CF. In addition, the GA signal transduction gene, MdRGA, which negatively regulates GA response, was also more highly expressed in YF than in CF, especially at 80 and 105 DABB. The expression of the cytokinin synthesis gene, MdIPT9, was significantly lower in YF than in CF during the period of rapid growth at 45, 65, 125, and 145 DABB. The MdCKX5 gene, which is involved in the breakdown and catabolism of cytokinin, was upregulated to a higher level in YF than in CF at 45 and 65 DABB. The expression of the ABA synthesis gene, MdNCED3, was significantly lower in YF than in CF at 65 and 85 DABB, but higher in YF than in CF at 105 DABB.

Auxin and GA₄₊₇ was applied to YF and CF trees. Results indicated that the application of GA promoted shoot growth and internode elongation in YF trees (Figure 11A). The terminal shoot started to grow rapidly 1 week after the GA treatment was applied (Figure 11B). In total, internode length in YF trees increased by about 33% after the GA treatment, however, this was still less than the internode length in CF trees (Figure 11A). The auxin application did not stimulate shoot growth in either YF or CF trees (data not shown).

The GA treatment dramatically downregulated miRNA166 expression (Figure 12). The potential target of miRNA166, MdHD-ZIP, was also downregulated during the first 2 days after the GA treatment. After the 6th day, however, MdHD-ZIP was upregulated. During the first 4 days after the GA treatment, the expression of miRNA164 in CF and YF trees was slightly upregulated. Then miRNA164 was downregulated in CF, but upregulated in YF. The expression of miRNA171 in YF trees exhibited little response to the GA treatment. The expression of miRNA171 in CF trees was upregulated on the 8th day after the application of the GA treatment. The expression of miRNA159 in both YF and CF trees showed little response to the GA treatment during the first 6 days, however, on the 8th day, the miRNA159 was upregulated in CF trees but downregulated in YF trees. Relative to the expression levels in CF trees, miRNA396 expression in YF trees slowly went down after the application of the GA treatment. The expression of miRNA167 exhibited little response during the first 2 days after the GA application but was downregulated on the 6th day, and then subsequently upregulated.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.