miRNA and Degradome Sequencing Reveal miRNA and Their Target Genes That May Mediate Shoot Growth in Spur Type Mutant “Yanfu 6”

The spur-type growth habit in apple trees is characterized by short internodes, increased number of fruiting spurs, and compact growth that promotes flowering and facilitates management practices, such as pruning. The molecular mechanisms responsible for regulating spur-type growth have not been elucidated. In the present study, miRNAs and the expression of their potential target genes were evaluated in shoot tips of “Nagafu 2” (CF) and spur-type bud mutation “Yanfu 6” (YF). A total of 700 mature miRNAs were identified, including 202 known apple miRNAs and 498 potential novel miRNA candidates. A comparison of miRNA expression in CF and YF revealed 135 differentially expressed genes, most of which were downregulated in YF. YF also had lower levels of GA, ZR, IAA, and ABA hormones, relative to CF. Exogenous applications of GA promoted YF shoot growth. Based on the obtained results, a regulatory network involving plant hormones, miRNA, and their potential target genes is proposed for the molecular mechanism regulating the growth of YF. miRNA164, miRNA166, miRNA171, and their potential targets, and associated plant hormones, appear to regulate shoot apical meristem (SAM) growth. miRNA159, miRNA167, miRNA396, and their potential targets, and associated plant hormones appear to regulate cell division and internode length. This study provides a foundation for further studies designed to elucidate the mechanism underlying spur-type apple architecture.


INTRODUCTION
Apple is an important fruit crop in temperate regions of the world. The use of dwarfing rootstocks and high density planting of apple trees is common in commercial orchards as it increases yields per unit area and reduces the cost of management practices, such as pruning. The use of spur-type scion varieties is one of effective methods used to facilitate high-density planting, along with the use of dwarfing rootstocks. "Fuji" is the major cultivar planted in China, accounting for more than 65% of apple plantings. Several spontaneous "Fuji" spur-type mutants have been identified and utilized in the breeding and planting of "Fuji" apples. The spur-type growth habit in apple is characterized by a reduced number of vegetative shoots and a corresponding increase in the number of fruiting spurs, short distances between nodes, compact and dense growth, and dark green and relatively thick leaves (Lapins and Lapins, 1969). Spur-type apple trees exhibit a high rate of bud break but have weak branching, characterized by rosette foliage and overall reduced levels of vegetative growth. The anatomical structure of spur-type shoots is characterized by vessel elements that are shorter and narrower than in standard-type shoots. The spur-type vessel elements have longer ends (tails), and thicker secondary cell walls (Yang et al., 2000).
Several earlier studies have reported that spur-type apple trees have lower levels of GA than standard-type trees, are less sensitive to applications of GA, and have reduced rates of conversion of GA19 to GA20 under high temperature conditions in summer (Steffens and Hedden, 1992;Song et al., 2012). The majority of shoot elongation is the result of cell expansion, with gibberellins playing a dominant role in this process. GAs promote internode elongation and overall shoot elongation by regulating cell proliferation, differentiation, and expansion. The "green revolution, " semi-dwarf varieties of cereal crops, are defective in GA biosynthesis and/or signaling. Auxin is recognized as the major hormone involved in shoot development (Gallavotti, 2013). Auxin also influences stem elongation and regulates the formation, activity, and fate of meristems. Cytokinins (CK) inhibit the differentiation of cells in the SAM and instead promotes cell division. CK promotes SAM enlargement partly through WUSCHEL, and defects in CK synthesis or perception reduces the size of the SAM size (Leibfried et al., 2006;Kyozuka, 2007). Except for GA, however, the levels of plant hormones in spur-type apples have not been reported or are unclear.
MicroRNAs (miRNAs) are a class of endogenous 20-24 nt non-coding RNAs. In plants, miRNAs are transcribed by RNA polymerase II, 5 ′ capped, spliced, and polyadenylated at the 3 ′ end to create pri-miRNAs. DICER-LIKE (DCL) protein processes the stem-loop structure in pri-miRNAs, creating pre-miRNA, and subsequently miRNA/miRNA * duplexes. A single mature miRNAs is generated by the removal of the miRNA * . The mature miRNA is recruited into a RISC complex (RISC) which degrades mRNA targets or suppresses their translation. Plant miRNAs play a critical role in growth, organ development, phase changes, senescence, fruit ripening, stress response, disease resistance, stress tolerance, and mineral assimilation and translocation.
Some miRNAs in plants have been shown to regulate plant architecture. A point mutation in OsSPL14 perturbs the cleavage of OsSPL14 transcripts by OsmiR156, which affects tiller number and lodging resistance, as and enhances grain yield (Jiao et al., 2010). OsmiR397 down regulates a laccase-like protein OsLAC involved in brassinosteroid sensitivity, and thereby promotes panicle branching . OsMIR444 participates in a miRNA/MADS/TCP/D14 (miMTD) regulation module, in which D14 functions in strigolactone (SL) signaling to control tillering, by suppressing OsMADS57 expression (Guo et al., 2013). In pear trees, miR6390 is involved in the transition from endormancy to ecodormancy by targeting PpDAM transcripts, degrading them and stimulating the release of PpFT2 (Niu et al., 2015). Overexpression of the MbDRB1, a gene involved in miRNA biogenesis in Malus baccata, increased the level of adventitious rooting, curly leaf, and columnar-like architecture (You et al., 2014).
The molecular mechanisms responsible for spur-type architecture is unclear. Few studies have been conducted on the role of miRNAs in the architecture of fruit trees. In the present study, miRNA and degradome sequencing was used to identify miRNA involved in shoot development in Malus domestica "Nagafu 2" and spur bud mutation "Yanfu 6." Morphological indicies, plant hormone levels, miRNA, and target gene expression, and the application of various hormone treatments were also measured and/or used to compare and contrast shoot development in "Yanfu 6" vs. "Nagafu 2." This comprehensive analysis was conducted in order to identify and determine the potential role of miRNA-mediated regulatory network in shoot growth and involvement of miRNA in the architecture of fruit trees and other woody, perennial plants.

Plant Materials
Experiments were conducted in 2013 and 2014 using standardtype "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.

Measurement of Tree Growth
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.

Hormone Measurements
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.

GA Treatment
YF and CF trees were sprayed with 500 mg·L −1 GA 4+7 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.

Small RNA and Degradome Sequencing
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).

Bioinformatic Analysis of miRNA
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 2 ) less than −1 or greater than 1 were considered to have significantly altered expression.

Prediction of miRNA Targets and Functional Annotation
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
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 . 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).

Statistical Analysis
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).

Growth Characteristics of YF and CF
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.

Plant Hormone Analyses
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.

Small RNA Sequencing Profiles in CF and YF
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 , peach (Prunus persica) , 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-miR172ik, 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 speciesspecific Malus domestica miRNA.
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).

Differentially Expressed miRNAs
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).

Target Prediction of miRNAs
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 ; 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.

GO and KEGG Analysis of the Degradome-Predicted Target Genes in Shoot Tips
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-toorgan 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.

RT-qPCR of Differentially Expressed miRNAs and Their Targets Involved in SAM Development
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.

RT-qPCR of Differentially Expressed miRNAs and Their Targets Involved in Internode Elongation
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.

Expression of Cell Cycle-and Cell Elongation-Related Genes
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.

Expression of Hormone-Related Genes
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.

Effect of GA on Shoot Growth and miRNA Expression
Auxin and GA 4+7 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.

DISCUSSION
YF is a natural spur mutant. YF trees have short internodes and reduced shoot length. Whether or not spur-type apple trees are also a fruiting characteristics, has been a source of debate, especially since the spur type phenotype of siblings derived from the spur-type apple exhibiting great variability (Fideghelli et al., 2003). Regardless, spur-type apple cultivars undoubtedly moderately reduced tree size as compared to standard-type trees (Costes et al., 2010). Unlike dwarfing rootstocks that reduce the proportion of lateral buds that develop into long shoots during the early development of trees (Lauri and Trottier, 2006;Seleznyova et al., 2008), YF produces a dwarf tree architecture by reducing the duration of shoot vegetative growth and shortening internode length.
In the present study, YF trees were demonstrated to contain lower levels of GA, ZR, IAA, and ABA. An application of GA promoted cell division in the shoot apical meristem and an increase in internode elongation in YF trees (Figures 11A,B). In contrast, YF trees exhibited no response to IAA. Previous studies have reported that the application of 6-benzylaminopurine and GA to apple nursery trees increased shoot elongation to a much greater extent than applying 6-benzylaminopurine alone (Zhang et al., 2011;Doric et al., 2015). Hence, among the hormones, GA, ZR, IAA, and ABA, it appears that GA plays a major role in YF trees.
The miRNA sequencing conducted in the present study confirmed the presence of 202 known miRNAs expression in apple shoot tips and discovered an additional 498 novel miRNAs. A total of 42 known miRNAs and 93 novel miRNAs were differentially expressed in the shoot tips of CF and YF trees. The differential expression of the miRNAs appeared to be related to the differences in shoot development observed in the two genotypes. This supports the premise that miRNAs are involved in the regulation of shoot architecture in YF trees.
Internode length and the number of nodes are two essential components that determine plant height. Compared to internode length, however, the impact of shoot apical meristem activity FIGURE 10 | Expression of hormone related genes. *Indicates P < 0.05, Student's t-test; Error bars indicate SE. has a greater effect on shoot development. For example, the shoot apical meristem is responsible for terminating cell division when trees enter dormancy or when the vegetative meristem is converted to a floral meristem. The number of nodes is also determined by continued cell divisions in the shoot apical meristem. Therefore, factors affecting SAM development would also affect the number of nodes produced on a shoot during an annual growth cycle. Our results indicate that miRNA164, miRNA166, and miRNA171 are potentially involved in SAM development and differentially expressed in YF vs. CF apple shoot tips. A degradation of the NAC-domain transcription factors, CUC1 and CUC2 (Mallory et al., 2004;Raman et al., 2008) by miRNA164 has been reported to constrain the expansion of the boundary domain by degrading the NAC-domain transcription factors, CUC1 and CUC2 (Mallory et al., 2004;Raman et al., 2008). In Arabidopsis, CUC1, CUC2, and CUC3 function redundantly in initiating SAM and establishing organ boundaries (Aida et al., 1997). Abolishing miR164 regulation of CUC2 resulted in progressive enlargement of the boundary domain (Raman et al., 2008). The shoot tips in YF trees exhibited abnormal shoot tips, with the dome of the apex appearing a little larger than the apical dome in CF trees. This phenotype is also present in apples with columnar architecture and in some dwarfing rootstocks (Petersen and Krost, 2013). The mir164abc mutant frequently displayed extremely short internodes, small floral organs, and reduced fertility. Thus, miRNA164 may be involved in abnormal shoot and internode development in YF trees.
Among all of the miRNAs, miRNA166 had the highest level of expression in the shoot tips. miRNA166 potentially inhibited SAM development by targeting a set of HD-ZIP transcription factors that are important for maintaining pluripotency in shoot apical meristems (Williams et al., 2005). AGO10 binds miRNA165/166 to prevent its association with the broadly expressed AGO1 and thus prevents the degradation of HD-ZIP transcription factors in the SAM (Zhou et al., 2015). YF trees have a high proportion of spur shoots. Many of the lateral shoots in YF trees stopped elongating soon after bud break. Thus, miRNA166 may be involved in apple regulating the spur shoot habit in apple trees. miRNA171 potentially targets the HAM subfamily genes, SCL6-II, SCL6-III, and SCL6-IV, within the GRAS gene family. miRNA171 affects the maintenance of SAM indeterminacy by regulating the shoot apical meristem WUS-CLV feedback loop (Schulze et al., 2010). The expression level of miRNA171 was significantly lower in YF trees than in CF trees. The height of Arabidopsis plants overexpressing miR171c was significantly higher than wild-type plants (Long et al., 2010). Upregulation of OsMIR171c in a delayed heading mutant rice genotype exhibited pleiotropic phenotypic defects, including a prolonged vegetative phase and a delayed heading date (Fan et al., 2015). These observations correspond to the reduced growth in YF shoots, FIGURE 13 | Tentative model for the association of miRNAs with the phytohormone crosstalk regulation the spur type apple tree architecture. Diagram represents a schematic drawing of an apple shoot apex. Arrows and inhibition lines linking the gene names represent positive and negative interactions, respectively. Up and down arrows behind the gene names represent up and down regulation, respectively. relative to shoot growth in CF trees, and the earlier date in which growth ceased in YF trees, relative to CF trees.
As a result, it appears that the final internode phenotype results from an interplay between cell division and elongation that progresses in an inverse directions. Cell number has a more predominant role in regulating the growth and development of internodes because of cell length increases only two to three times, but cell numbers increase by 10 to 30 times (Pallardy, 2008). Analysis of the expression pattern of cell cycle genes further supports the premise that the SAM in YF trees has a lower rate of cell division than in CF trees (Figure 12). The expression of MdCYCD, MdCYCB1;2, MdCYCD1;3, and MdRBR1 was significantly lower in YF shoot tips than in CF shoot tips. The period of high expression of the cell cycle genes corresponded to the period of rapid shoot growth; suggesting that cell division may play an instrumental role in shoot growth, and more so than cell elongation. Ripetti et al. (2008) also found that genetic variation in internode length could be primarily attributed to cell numbers, while cell length played more of a secondary role. The miRNA sequencing and GO annotation in the present study identified several miRNAs related to cell division. miRNA159 potentially deregulates its target genes, MYB33 and MYB65, in vegetative tissues, and inhibit growth by reducing cell proliferation (Alonso-Peral et al., 2010). The expression of miRNA159 was significantly higher in CF shoot tips than in YF shoot tips. The miRNA159 loss of function mutant, mir159ab, has an enlarged shoot apical meristem and smaller lateral organs (Alonso-Peral et al., 2010). miRNA159 and its targets have also been reported to be involved in GA-mediated flowering (Gocal et al., 2001;Achard et al., 2004). The expression of miRNA159 in YF shoot tips exhibited little response to the application of GA during the first 6 days after spraying GA. Interestingly, miR159a and miR159b also remained unchanged after the application of GA in wild-type Arabidopsis (Alonso-Peral et al., 2010). The potential target gene of miR159a, MYB65, however, was downregulated after the application of GA in Malus domestica. Thus, GA maybe regulate the expression of the target of miRNA159 to produce the miRNA159 phenotype. miRNA396 negatively regulates its potential target gene, GROWTH RESPONDING FACTOR (GRF), to affect cell division in shoot meristems, leaves, and roots. miRNA396 appears to act as a strict repressor of the GA and CK pathways, and negatively regulates cell proliferation by decreasing cell division activity and the expression of cell cycle-related genes . The potential downregulation of the GA pathway genes by miR396 corresponds to the dwarf phenotype observed in miR396 overexpressing lines (Liu et al., 2009). The high level of expression of miRNA396 in YF shoot tips may inhibit cell division and internode elongation. miRNA396 was downregulated after GA application. Thus, GA may affect downstream GA response by downregulating miRNA396. miRNA167 regulates vegetative and reproductive growth by controlling the expression of auxin response factors 6 and 8 (ARF6/8) in plants. Overexpression of miRNA167a results in twisted leaves, short inflorescences, and arrested flower development; thereby fully producing the mature phenotype of arf6arf8 plants (Wu et al., 2006). The expression of miRNA167 was significantly higher in YF shoot tips than in CF shoot tips. GA treatment downregulated miRNA167 at first, and then subsequently induced it; resulting in upregulating. But GA also increased the expression of the target gene, ARF6. DELLA has been reported to interact with ARF6 and inhibit ARF6 DNA-binding to modulate ARF6 gene expression (Oh et al., 2013). Hence, GA may promote cell elongation mechanisms in YF trees that involve miRNA167-ARF-mediated responses.
Based upon the findings of the current study, we present a model to explain how plant hormones and miRNA are involved in regulating the shoot development in spur-type YF trees (Figure 13). Low levels of IAA, GA, and CK are associated with reduced shoot growth in YF trees. GA appears to play a major role in regulating shoot development. Low levels of GA in YF shoots may reduce the number of cell divisions in the SAM, relative to the SAM in CF shoots, thus causing shoot growth to cease much earlier in the YF trees as compared to in the CF trees. In the model, low expression of miRNA164, miRNA171, and high expression of miRNA166 and their potential target genes affected SAM development and reduced the duration of shoot growth in YF trees. Of the three miRNAs, the regulation of miRNA166 by GA appeared to play a major role in shoot tips elongation in YF trees. The high levels of expression of miRNA167 and miRNA396, and the low expression of miRNA159, in YF trees appears to have inhibited cell division and reduced internode length. The results of the present study, and the hypothetical model presented, contributes to a better understanding of shoot development in perennial plants; especially those that exhibit a spur-type growth architecture.

AUTHOR CONTRIBUTIONS
CS, DZ, JM, and MH participate in the experimental design. CS, JZ, LZ, and WL participate in material sampling, field measurements and the laboratory data measurement. CS, BZ, YL, and GL participate in the laboratory data measurement. CS, DZ, JM, and MH participate in the paper writing and discussion. All authors reviewed the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 00441/full#supplementary-material Figure S1 | The meteorological condition of the experiment site. Figure S2 | Size distribution of novel miRNAs and the identity of the first nucleotide.
Table S1 | Quality evaluation of CF and YF sequencing data.