Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.).

Flower induction in apple (Malus domestica Borkh.) trees plays an important life cycle role, but young trees produce fewer and inferior quality flower buds. Therefore, shoot bending has become an important cultural practice, significantly promoting the capacity to develop more flower buds during the growing seasons. Additionally, microRNAs (miRNAs) play essential roles in plant growth, flower induction and stress responses. In this study, we identified miRNAs potentially involved in the regulation of bud growth, and flower induction and development, as well as in the response to shoot bending. Of the 195 miRNAs identified, 137 were novel miRNAs. The miRNA expression profiles revealed that the expression levels of 68 and 27 known miRNAs were down-regulated and up-regulated, respectively, in response to shoot bending, and that the 31 differentially expressed novel miRNAs between them formed five major clusters. Additionally, a complex regulatory network associated with auxin, cytokinin, abscisic acid (ABA) and gibberellic acid (GA) plays important roles in cell division, bud growth and flower induction, in which related miRNAs and targets mediated regulation. Among them, miR396, 160, 393, and their targets associated with AUX, miR159, 319, 164, and their targets associated with ABA and GA, and flowering-related miRNAs and genes, regulate bud growth and flower bud formation in response to shoot bending. Meanwhile, the flowering genes had significantly higher expression levels during shoot bending, suggesting that they are involved in this regulatory process. This study provides a framework for the future analysis of miRNAs associated with multiple hormones and their roles in the regulation of bud growth, and flower induction and formation in response to shoot bending in apple trees.

Introduction

Apple (Malus domestica Borkh.) is an important economic fruit tree worldwide, in which flower induction and flower buds formation occur in two growing seasons. China is the leading country of apple production, with a planting area of 3.1 million hectares and a production of 33 million tons annually. ‘Fuji’ apple cultivar represents ~65% of the apple trees in China (http://faostat.fao.org/). However, the induction and formation of flower buds can be difficult, and this, accompanied by the alternate bearing in young ‘Fuji’ apple trees, represents a serious problem, restricting the development of the apple industry in northern China. Previous studies showed that shoot bending was an effective measure widely used for promoting flower bud formation in many fruit trees, such as pear (Ito et al., 2004) and apple (Lauri and Lespinasse, 2001).

Additionally, the variability of physiological responses to bending produced in apple trees, such as stimulating bud growth, increasing the sink capacity of the bud relative to other tissues and affecting the tree structure (Han et al., 2007), as well as hormone level changes in response to bending shoots, may play an important role in plant growth and flower bud formation (Liu and Chang, 2011). Additionally, tree shoot bending generates a new long‐distance signal that changes growth and development during the trees’ life cycles (Lopez et al., 2014). However, the molecular regulatory mechanisms of bud growth and flower induction in response to shoot bending are not very clear, especially the regulatory process involving microRNAs (miRNAs), and whether these miRNAs actually exist in buds during the flower induction process.

Previous studies indicate that plant hormones play an important role in the regulation of bud growth and flower induction, as well as in the growth and development of woody plants (Curaba et al., 2014; Mutasa‐Gottgens and Hedden, 2009), and take part in growth regulation in response to stress (Tsai and Gazzarrini, 2014). Phytohormones also act as signal molecules in plants to control growth and development through the regulation of gene expression, involving different pathways in response to environmental changes. miRNAs act as important regulators involved in multiple links with hormone crosstalk‐mediating developmental processes, such as flower development, phase transition and stress responses, in plants (Curaba et al., 2014). In addition, phytohormone treatments can regulate the expression of miRNAs, affecting growth and development during their life cycles (Chen et al., 2012). Recently, a large number of miRNAs involved in multiple aspects of biological processes, such as phase transition (Xing et al., 2014), stress response (Chen et al., 2012; Ding et al., 2013), hormone signalling (Curaba et al., 2014) and flower induction and development (Kim and Ahn, 2014; Spanudakis and Jackson, 2014), were identified in plants. In apples, some conserved miRNAs, such as miR156, 172, 160, 319 and 396, along with their expression profiles in different tissues, were studied and indicated that they may play roles in multiple regulatory functions in plant growth and development (Xia et al., 2012b). In addition, potential targets of identified miRNAs in apple (MdSPLs, MdARFs, MdTOE1 and MdF‐box) were responsible for plant growth, development and flowering and stress responses (Varkonyi‐Gasic et al., 2010), and some genes, such as MYB transcription factors (TFs), are highly conserved in plant species and are also involved in improving plant tolerance to multiple abiotic stresses (Cao et al., 2013). This suggests that these miRNAs and their targets may play important roles in the regulation of plant growth in response to environmental stress.

Additionally, molecular regulatory mechanisms of flower induction and the flowering process have been studied in annual (Kurokura et al., 2013) and woody plants (Wang et al., 2011). Some flowering control genes, such as FT, SOC1, LEAF, SPLs and AP1, play important roles in plant flowering through multiple flowering pathways (Lee and Lee, 2010; Mutasa‐Gottgens and Hedden, 2009; Wahl et al., 2013). In addition, miR156‐SPLs and miR172‐AP2 modules were mainly involved in regulating the ageing flowering pathway (Spanudakis and Jackson, 2014), suggesting that flowering genes and relevant miRNAs may act together in plant flower induction and the flowering process.

Therefore, understanding the molecular regulatory mechanisms of apple bud growth and flower induction associated with hormones and miRNAs in response to shoot bending will greatly assist us in solving the problem of the long juvenile phase and poor quality flower buds in apple trees. In this study, we exploited high‐throughput sequencing and degradome analyses to comprehensively identify miRNAs and their targets in response to shoot bending during the flower induction process in apple.

Experimental procedures

Plant material and sample collection

Six‐year‐old ‘Fuji’ apple (Malus domestica Borkh.) trees growing on M.26 rootstocks were planted in the Apple Demonstration Nursery of Yangling Modern Agriculture Technology Park (Northwest Agriculture & Forestry University), Shaanxi Province of China (34° 52′ N, 108° 7′ E), and buds sampled during bud growth and the flower induction process were collected directly into liquid nitrogen, including those of control and shoot‐bending treatments. Buds on the tops of spurs in both control and shoot‐bending treatments (Control: vertical shoots; Shoot bending: shoots bending at 110°) were separately collected nine times as mixed samples from 60 6‐year‐old trees. The size and weight, as well as the hormone levels, indole‐3‐acetic acid (IAA) and abscisic acid (ABA), were measured three times, the early, middle and later stages of flower induction [ES (5th May, 2013), MS (1th June, 2013) and LS (25th June, 2013), respectively]. All samples (six buds samples each from control and shoot‐bending treatments) were immediately frozen in liquid nitrogen and stored at −80 °C for RNA extraction and small RNA libraries construction. In addition, a mixed sample from the six buds was used for degradome sequencing. Total RNA was isolated from each sample using a modified method (Xing et al., 2014). Each treatment, including control and shoot bending, was subdivided into three subgroups consisting of 20 trees each to allow for replications.

Determination of bud changes during flower induction

Approximately 30 buds from control and shoot‐bending treatments at each time were used to determine bud sizes and weights. The bud changes were measured nine times from ES to LS during flower induction in six ‘Fuji’ trees receiving either the control or shoot‐bending treatment.

Statistical analysis of the flowering rate

Six 6‐year‐old ‘Fuji’ apple (Malus domestica Borkh.) trees were used to calculate the flowering rate in control and shoot‐bending treatments. The methods for calculating the flowering rates in both control and shoot‐bending treatments were determined as described in detail by Ito et al. (1999).

Hormonal content analysis

Approximately 0.5 mg fresh weight of buds in each sample from the three developmental stages, ES, MS and LS, of control and shoot‐bending treatments during the flower induction process was used for phytohormone extractions (Dobrev and Vankova, 2012). The detection and analysis were performed using a high‐performance liquid chromatography (Waters 2498/UV; Visible Detector, Shaanxi, China) (Djilianov et al., 2013).

RNA deep sequencing and library construction

Using an Illumina Genome Analyzer, six small RNA libraries and one degradome library were constructed and sequenced by the Biomarker Biotechnology Corporation (Beijing, China). The methods of small RNAs extraction and enrichment can be seen in Xing et al. (2014). First, total RNA extractions from six buds libraries were performed using the RNeasyPlant Mini Kit (Qiagen, Hilden, Germany), and each was collected into RNA pools for small RNA construction. Second, ~16–30 nt gel fragments were selected and ligated to a pair of adapters at the 5′‐ and 3′‐ends using T4 RNA ligase. These small RNAs with adapters were transcribed into cDNA using Super‐Script II Reverse Transcriptase (Invitrogen, Shanghai, China). Third, the cDNA products were used for amplification by polymerase chain reaction (PCR), and then, the purified PCR products were sequenced by the Biomarker Biotechnology Corporation (Beijing, China). In addition, total RNA from a mixed sample [six buds each at the ES, MS, and LS under both control conditions (CES, CMS, and CLS, respectively) and under shoot‐bending conditions (BES, BMS, and BLS, respectively)] was used for degradome library sequencing and construction using a previously described method (Xing et al., 2014).

Prediction and identification of targets of known and novel miRNAs

We predicted the potential targets of known and novel miRNAs using methods previously described in detail by Xing et al. (2014), and we also identified some targets of known and novel miRNAs using degradome sequencing that matched the degraded fragments to the apple genome (Malus domestica Borkh.) after removing ncRNAs, as well as polyN fragments, in the mixed bud sample to reduce interference (Xing et al., 2014). In addition, to predict the potential functions of these targets, we used gene ontology (GO) categories (http://www.geneontology.org/) mainly associated with various biological processes, molecular functions and cellular components to classify the identified target genes.

Venn diagrams of known miRNAs and novel miRNAs and their expression profile analyses

Differentially expressed known and novel miRNAs were analysed using Venn diagrams (Martin et al., 2012) with the VENNTURE software (http://www.irp.nia.nih.gov/branches/lci/nia_bioinformatics_software.html). In addition, we analysed the known and novel miRNAs’ expression profiles over time, as well as between control and shoot‐bending treatments, according to the methods that are described in detail by Botton et al. (2011). Additionally, hierarchical clustering heat maps were generated by the HemI Manuel software 1.0 (http://hemi.biocuckoo.org/down.php) using the log‐FPKM values of each miRNA for control and shoot‐bending treatments.

qRT‐PCR validation of miRNAs, their targets and flowering genes

We used stem‐loop RT‐PCR to determine the miRNA expression levels in buds of control and shoot‐bending treatments during flower induction (Varkonyi‐Gasic et al., 2010). In addition, the methods of determining the expression levels of miRNA targets and flowering genes can be seen in Xing et al. (2014). Briefly, cDNAs of miRNA targets and flowering genes were generated from 2 μg each of the total RNAs of six bud samples in both control and shoot‐bending treatments during the flower induction process using a PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Ohtsu, Japan), and qRT‐PCR was performed using SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) with 10 μL of 2X SYBR® Premix Ex Taq II and 0.8 μL of forward and reverse primers in a 20‐μL system (Takara, Ohtsu, Japan) to determine the expression of flowering genes and miRNA targets (Varkonyi‐Gasic et al., 2010). All primers used in the qRT‐PCR experiments, including those for miRNAs, their targets and flowering genes, are listed in Table S1.

Statistical analysis

For bud sizes and weights, flowering rates, hormone levels and qPCR results in both control and shoot‐bending treatments, a one‐way analysis of variance with Tukey–Kramer multiple comparison tests was performed using DPS software, version 7.0 (DPS version 7.0; Zhejiang University, Hangzhou, China).

Results

Bud growth during flower induction and flower promotion in response to shoot bending in ‘Fuji’ apple trees

Bud growth during the flower induction period was studied in ‘Fuji’ apple trees in response to shoot bending (Figure 1). The length, width and dry weight under control and shoot‐bending treatments increased from ES to LS, but shoot bending significantly increased bud size during this developmental period (Figure 1a–c). In addition, bud growth rate changes showed that the bud sizes changed mainly during the early period of flower bud induction, from ES to MS, and the increase was more significant under the shoot‐bending treatment (Figure 1d–f).

Figure 1

Bud growth and growth rates during flower induction in ‘Fuji’ apple. (a) Length, (b) Width, (c) Fresh weight, (d) Growth rate of bud length, (e) Width and (f) Fresh weight.

Our results showed that shoot bending significantly increased, by ~58%, the flowering rate in ‘Fuji’ apple trees compared with the control (Figure 2).

Figure 2

Flowering rate of ‘Fuji’ apple trees in both control and shoot‐bending treatments.

Hormone content changes in buds during flower induction

The hormone levels were analysed in ‘Fuji’ buds in shoot‐bending and control plants at three time points (ES, MS, and LS) during flower induction (Figure 3). The auxin (AUX) content in both shoot‐bending and control plants was high in the ES and MS, but showed relatively lower levels in the LS stage (Figure 3a). In addition, in the earlier stages of buds growth (ES and MS), the AUX content was significantly higher in the shoot‐bending treatment than in control buds (Figure 3a). The changes in the ABA contents in buds during flower induction between shoot‐bending and control treatments increased from ES to LS (Figure 3d). In the later stages (MS and LS), the ABA content was also significantly higher in the shoot‐bending treatment than in control buds (Figure 3d). The changes in the CK content in buds decreased from ES to LS in both control and shoot‐bending treatments, but displayed a significantly higher level in shoot‐bending treatments than in the control (Figure 3b). Meanwhile, the changes in gibberellic acid (GA) content had a result similar to that of the CK content in buds during the flower induction process, but the GA level in the shoot‐bending treatment was significantly lower than that in the control (Figure 3c).

Figure 3

Hormone content of buds during the flower bud physiological differentiation stage in ‘Fuji’ apple. (a) AUX: Auxin; (b) CK: cytokinin; (c) GA: gibberellin; (d) ABA: abscisic acid.

Construction and sequencing of small RNA and degradome libraries

To determine responsive sRNAs in the shoot‐bending and control treatments during bud growth (ES, MS, and LS), six miRNA libraries were constructed and sequenced (Table S2). A total of 14 095 388, 14 106 904, 15 242 433, 13 259 311, 15 969 504 and 13 523 103 raw reads were generated by the high‐throughput Illumina HiSeq 2000 Sequencing System in CES, CMS, CLS, BES, BMS and BLS, respectively (Table S3), with 6 382 280 (45.3%), 6 627 824 (47.0%), 6 464 149 (42.4%), 5 349 068 (40.3%), 6 531 048 (41.0%) and 5 713 280 (42.2%) unique raw reads from the CES, CMS, CLS, BES, BMS and BLS libraries, respectively (Table S3). In addition, we also constructed a degradome library using total RNA from six mixed samples (CES, CMS, CLS, BES, BMS and BLS) (Table S4). A total of 30 762 927 (93.08%) clean and 373 406 (1.13%) unique reads were generated (Table S4). For the six miRNA libraries, the reads’ size distributions were quite similar (Figure S4), but the read size distribution in the degradome library was significantly different from those of the miRNA libraries (Figure S2). The length of sRNA varied, with 21‐, 22‐, 23‐ and 24‐nt small RNAs forming the major population, and 24 nt was the most dominant length in the six libraries (Figure S1), which was similar to the results obtained from other plants, such as Malus hupehensis and Arabidopsis thaliana (Confraria et al., 2013; Xing et al., 2014). The total raw reads were categorized into 12 different classes, exon _antisense (0.91%), exon_sense (1.42%), intron _antisense (2.31%), intron_sense (2.72%), miRNA (2.50%), rRNA (19.16%), repeat (3.04%), scRNA (0), snRNA (0.02%), snoRNA (0.02%), tRNA (2.06%) and others (65.83%), for the CES library (Table S3), by matching them to the domesticated apple's genome in the Rfam (http://www.sanger.ac.uk/resources/databases/rfam.html) and GenBank (http://www.ncbi.nlm.nih.gov/genbank) databases using the tag2annotation software. This and other library information can been seen in Table S3.

The sequencing data of the six libraries and one degradome library have been deposited in the NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi).

Identification of known miRNAs and their expression profiles

To identify the known miRNAs in apple (Malus domestica Borkh.), the sRNAs in the six bud libraries were queried using BLASTN to known mature plant miRNAs of M. domestica in the miRBase 18.0 (http://www.mirbase.org) and plant miRNA (http://bioinformatics.cau.edu.cn/PMRD) databases (Xing et al., 2014). A total of 195 known miRNAs that belonged to 41 miRNA families were identified (Figure 4 and Data S1). Among them, the number of miRNA members within different families is very different (Figure 4). A majority of the 39 known miRNA families had several members, such as mdm‐miR156 (27 members), 164 (6 members), 171 (14 members), 172 (15 members) and 396 (7 members), and six of the known miRNA families, mdm‐miR1511, 391, 7125, 7126, 827 and 858, had only one member (Figure 4).

Figure 4

Numbers of identified miRNAs in known miRNA families in ‘Fuji’ apple ( Borkh.).

In the six buds libraries, known miRNAs expression levels were determined by their read content's frequencies (Data S1). The information on the mature sequences of identified known miRNAs can be seen in Data S2. The expression levels of known miRNAs in each stage of control and shoot‐bending treatments were classified into eight categories based on their reads counts: no expression (0 reads), lowest (1–9 reads), lower (10–49 reads), low (50–99 reads), moderate (100–499 reads), high (500–999 reads), higher (1000–9999 reads) and highest (>10 000 reads) (Figure 5). Our results showed that the percentages of known miRNAs that fell into the eight categories in each library were highly similar (Figure 5). The largest percentages of known miRNAs fell into the lowest (1–9 reads) category at 28.35%, 33.51%, 32.47%, 31.96%, 34.54% and 36.08% in the CES, CMS, CLS, BES, BMS and BLS libraries, respectively, and the lowest percentages of reads fell into the highest (>10 000 reads) category in each library at 3.09%, 3.61%, 2.58%, 3.09%, 3.09% and 2.58% in CES, CMS, CLS, BES, BMS and BLS, respectively (Figure 5).

Figure 5

Expression levels of identified known miRNAs with their read content's frequencies in each library.

A Venn analysis showed that the number of differentially expressed known miRNAs that changed over time was 129 miRNAs (CMS vs CES, CLS vs CMS, CLS vs CES, BMS vs BES, BLS vs BMS and BLS vs BES) (Figure 6), and among them, the expression levels of 42 miRNAs decreased, including mdm‐miR1511 and 171f/n/i/k, and those of 20 miRNAs increased, including mdm‐miR172a, 164a and 7121 g (Figure 6a–c). In addition, there were 124 differentially expressed known miRNAs between control and shoot bending (BES vs CES, BMS vs CMS and BLS vs CLS) (Figure 6), and among them, the expression levels of 68 miRNAs decreased, including mdm‐miR171b, 167e and 477b, and 27 increased, including mdm‐miR5225a, 166 h and 396c, in shoot‐bending compared with the control treatments (Figure 6d–f).

Figure 6

Venn diagrams of differentially expressed known miRNAs in ‘Fuji’ apple buds. (a) Decreased expression of miRNAs over time; (b) increased expression of miRNAs over time; (c) decreased expression of miRNAs between shoot‐bending and control treatments; (d) increased expression of miRNAs between shoot‐bending and control treatments; (e) Venn analysis of decreased and increased expression levels of miRNAs over time; (f) Venn analysis of decreased and increased expression levels of miRNAs between shoot‐bending and control treatments; (g) Venn analysis of decreased expression levels of miRNAs over time and differentially expressed known miRNAs between shoot‐bending and control treatments; (h) Venn analysis of both decreased and increased expression levels of miRNAs over time and differentially expressed known miRNAs between shoot‐bending and control treatments; (i) Venn analysis of increased expression levels of miRNAs over time and differentially expressed known miRNAs between shoot‐bending and control treatments.

Additionally, the Venn analysis showed that the expression levels of 33 miRNAs, including mdm‐miR535b, 160b, 391, 7124a, 398a and 393c, decreased over time and in shoot‐bending compared with the control treatments (Figures 6g and 7a). However, the expression levels of just five miRNAs, mdm‐miR162a, 164b, 171a, 167f and 172b, decreased over time but increased in the shoot‐bending treatment compared with the control (Figures 6g and 7e). A total of 11 miRNAs, including mdm‐miR156x/z, 164a, 395a and 7121c, had increased expression levels over time, but decreased expression levels in the shoot‐bending treatment compared with the control (Figures 6i and 7b). Only seven miRNAs, including mdm‐miR5225a, 399c and 160a, had increased expression levels in the shoot‐bending treatment compared with the control and had increased expression levels over time during flower induction (Figures 6i and 7e). In addition, 22 miRNAs, including mdm‐miR394b, 396a/b, 398c and 162b, had decreased expression levels in the shoot‐bending compared with the control treatments (Figures 6g and 7c). By contrast, 14 miRNAs, including mdm‐miR7120b, 7121a and 398b, showed increased expression levels in the shoot‐bending compared with the control treatments (Figures 6i and 7f). Another 22 miRNAs, including mdm‐miR2118a, 159b, 167b/d and 399b, were differentially expressed both over time and in the treatments (shoot‐bending and control) (Figures 6h and 7g).

Figure 7

Hierarchical clustering of differentially expressed known miRNAs over time, and between shoot‐bending and control treatments in ‘Fuji’ apple buds.

Identification of putative novel miRNAs and their expression profiles

The domesticated apple (Malus domestica Borkh.) was used as a reference genome to predict potential novel miRNAs. The cut‐off to eliminate miRNAs was <10 reads per million in at least one library. We identified 137 putative unique novel miRNAs in each library (Data S3). The information on the mature sequences of identified novel miRNAs can be seen in Data S4. For novel miRNAs in each library, the expression level, based on their read content's frequencies, were classified into six categories: no expression (0 reads), lower (1–9 reads), low (10–49 reads), moderate (50–99 reads), high (100–1000 reads) and higher (>1000 reads) (Figure 8). In our results, the percentage of novel miRNAs that fell into the six categories were very similar in each library (Figure 8). The highest percentage of novel miRNAs fell into the lower (1–9 reads) category at 47.45%, 40.15%, 51.82%, 57.66%, 35.77% and 53.28%, followed by those that fell into the low (10–49 reads) category at 38.69%, 45.99%, 32.85%, 27.74%, 47.45% and 30.66%, for CES, CMS, CLS, BES, BMS and BLS, respectively (Figure 8). However, the lowest expression level in each library was the higher (>1000 reads) category at 0.73% for CES, CMS, CLS and BLS; 1.46% for BMS; and 0% for BES (Figure 8).

Figure 8

Expression levels of identified novel miRNA with their read content's frequencies in each library.

The results of a Venn analysis showed that 38 differentially expressed novel miRNAs had decreased expression levels over time (CMS vs CES, CLS vs CMS, CLS vs CES, BMS vs BES, BLS vs BMS and BLS vs BES), including novel‐m1736‐3p, m0829‐3p, m0697‐3p and m2042‐5p (Figure S3a). However, 34 differentially expressed novel miRNAs, including novel‐m1424‐5p, m0546‐5p, m0559‐3p, m1618‐3p and m2000‐3p, had increased expression levels over time (Figure S3b). Additionally, 24 differentially expressed novel miRNAs, including novel‐m1913‐5p, m1900‐5p, m1336‐3p, m0017‐3p, m1618‐3p, m1977‐5p and m0163‐5p, were expressed at significantly higher levels in the shoot‐bending treatment than in the control (Figure S3c). However, 23 other differentially expressed novel miRNAs, including novel‐m1401‐5p, m0893‐5p, m0088‐3p, m0490‐3p, m1161‐3p, m0173‐3p, m0899‐5p, m0829‐3p and m0968‐5p, were expressed at significantly higher levels in the control than in the shoot‐bending treatments (Figure S3d).

The hierarchical clustering of 31 differentially expressed novel miRNAs in the control and shoot‐bending treatments of ‘Fuji’ buds during flower induction resulted in five major clusters (Figure 9). The six novel miRNAs, novel‐m0088‐5p, m0893‐5p, m0530‐3p, m0173‐3p, m2000‐3p and m1953‐3p, forming cluster 1 displayed higher expression levels in the shoot‐bending treatment than in the control (Figure 9), but the novel miRNAs within cluster 5, novel‐m1913, m1064‐5p, m1317‐3p, m1736‐3p, m1018‐3p and m1727‐3p, displayed the opposite result, having higher expression levels in the control than in the shoot‐bending treatment during flower induction (Figure 9). In addition, the expression levels of the other 19 novel miRNAs, which belonged to the three remaining clusters (2, 3 and 4), displayed relatively low expression levels in each library under both control and shoot‐bending treatments during bud growth (Figure 9).

Figure 9

Hierarchical clustering of differentially expressed novel miRNAs in ‘Fuji’ apple buds.

Targets of known and novel miRNAs

To explore the identified known and novel miRNAs in diverse biological processes and obtain insights into the regulatory mechanisms behind shoot bending's promotion of flower induction and development in ‘Fuji’ apple trees, we identified targets of known and novel miRNAs by degradome sequencing (Table 1 and Data S5; Table 2 and Data S6). A total of 195 differentially expressed known miRNAs belonging to 41 miRNA families, including miRNA1511, 156, 159, 160, 393 and 396, were detected in control and shoot‐bending treatments during flower induction (Table 1). The detailed distribution information on the known miRNAs and their target cleavage sites are shown in Table 1 and Data S7. In addition, 40 targets of 31 differentially expressed novel miRNAs were detected in each library (Table 2). Meanwhile, the detailed distribution information on the novel miRNAs and their target cleavage sites are shown in Table 2 and Data S8. In most cases, multiple targets can be regulated by a single miRNA. For example, mdm‐miR160 regulated ARF16 and ARF17; mdm‐miR393 regulated AFB2, AFB3 and TIR1 genes; and mdm‐miR396 regulated GRF TFs, including GRF1, 2, 4, 5, 7, 8 and 9 (Table 1). An example of the cleavage information on mdm‐miR396 and its GRFs targets can be seen in Figure 10. The known miRNAs regulated some TF targets, such as SPL2 and SPL9 (mdm‐miR156), MYB5 (mdm‐miR858), GRFs (mdm‐miR396), TCP2 and TCP4 (mdm‐miR319), NF‐YA1 and NF‐YA7 (mdm‐miR169), and MYB65 and MYB101 (mdm‐miR159), while others were regulatory proteins, such as WLIM1 protein (mdm‐miR398), NAC036 protein (mdm‐miR7121) and PDL2 protein (mdm‐miR164) (Table 1).

Table 1

Potential targets of the identified known miRNAs from ‘Fuji’ apple ( Borkh.) buds and their GO biological process

Known miRNAs	Targets
	Gene ID	Description	GO biological process	GO ID	Alignment Range	Cleavage Site
1511	MDP0000320775	Protein of unknown function (DUF502) AT2G20120.1 (COV1)	Stem vascular tissue pattern formation	GO:0010222	370–390	381
156a‐z	MDP0000297978	Squamosa promoter binding protein‐like 9(SPL9)	Vegetative to reproductive phase transition of meristem	GO:0010228	794–814	805
	MDP0000155354	Squamosa promoter binding protein‐like 2(SPL2)	Regulation of timing of transition from vegetative to reproductive phase	GO:0048510	1136–1155	1145
	MDP0000249364	Acyl‐CoA synthetase 5 (AT1G62940.1)	Anther development	GO:0048653	2714–2733	2724
160a/b/c/d/e	MDP0000232116	Auxin response factor 17(ARF17)	Auxin‐mediated signalling pathway	GO:0009734	1302–1322	1313
	MDP0000221322	Auxin response factor 16(ARF16)	Response to auxin stimulus	GO:0009733	1597–1616	1607
166a/b/c/d/e/f/g/h/i	MDP0000943529	Homeobox‐leucine zipper family protein (ATHB14)	Meristem initiation	GO:0010014	587–607	598
	MDP0000240313	Cystathionine beta‐synthase (CBS) protein	Biological process unknown	GO:0008150	1992–2012	2003
168a/b	MDP0000772373	General regulatory factor 2(GRF2)	BR‐mediated signalling pathway	GO:0009742	424–444	435
169a/b/c/d	MDP0000296077	Nuclear factor Y, subunit A1(NF‐YA1)	Regulation of timing of transition from vegetative to reproductive phase	GO:0048510	1631–1651	1642
	MDP0000164531	Nuclear factor Y, subunit A7(NF‐YA7)	Regulation of timing of transition from vegetative to reproductive phase	GO:0048510	914–934	925
171a/b/f/g/k/l/m/n/j	MDP0000151144	GRAS family transcription factor (ATHAM3)	Cell differentiation	GO:0030154	1377–1397	1388
	MDP0000820534	RNA‐binding KH domain‐containing protein(AT5G53060.1)	Nucleic acid binding	GO:0003676	1377–1397	1388
	MDP0000784909	GRAS family transcription factor (ATHAM1)	Cell division	GO:0051301	1227–1247	1235
393a/b/c	MDP0000173674	Cyclin family protein(AT5G45190.1)	Flower development	GO:0009908	1704–1725	1716
	MDP0000249678	RNA polymerase III subunit RPC82 family protein(AT3G49000.1)	Endomembrane system	GO:0012505	39–60	50
	MDP0000268652	Auxin signalling F‐box 3(AFB3)	Auxin binding	GO:0010011	1503–1522	1513
	MDP0000498419	F‐box/RNI‐like superfamily protein (TIR1)	Response to auxin stimulus	GO:0009733	1524–1543	1534
	MDP0000125975	F‐box/RNI‐like superfamily protein (TIR1)	Response to auxin stimulus	GO:0009733	1506–1525	1516
	MDP0000203334	Auxin signalling F‐box 2(AFB2)	Auxin binding	GO:0010011	1860–1879	1870
394a/b	MDP0000303022	DNA‐binding storekeeper protein‐related transcriptional regulator	Transcription regulator activity	GO:0030528	636–654	646
	MDP0000310920	Unknown protein(AT3G54750.3)	Biological process unknown	GO:0008150	712–730	721
397a/b	MDP0000306215	GCR2‐like 1(GCL1)	Catalytic activity	GO:0003824	820–839	830
399a/b/c/e/f/g/h/i/j	MDP0000253476	Unknown protein	Unknown protein	Unknown	233–253	244
482b/c	MDP0000190531	LRR and NB‐ARC domains‐containing disease resistance protein	Defence response	GO:0006952	583–604	595
	MDP0000296741	Plant invertase/pectin methylesterase inhibitor superfamily	Cell wall modification	GO:0042545	59–78	69
5225a/b	MDP0000253714	Plasmodesmata‐located protein 2(PDLP2)	Plasmodesma	GO:0009506	583–604	595
	MDP0000270176	Thioredoxin superfamily protein(AT5G42850.1)	Cell redox homoeostasis	GO:0045454	944–966	955
535a/b	MDP0000207199	Transmembrane proteins 14C(AT3G43520.1)	Chloroplast envelope	GO:0009941	624–643	635
7121a/b/c	MDP0000480581	NAC domain‐containing protein 36(NAC036)	Negative regulation of cell size	GO:0045792	379–399	390
159a/b/c	MDP0000147309	MYB domain protein 65(MYB65)	Response to salicylic acid stimulus	GO:0009751	954–973	964
	MDP0000179306	myb domain protein 101(MYB101) ABA responsive gene	Pollen tube growth	GO:0009860	1077–1096	1087
	MDP0000319634	myb‐like HTH transcriptional regulator family protein(DUO1)	Pollen sperm cell differentiation	GO:0048235	587–606	597
	MDP0000308617	myb domain protein 101(MYB101) ABA responsive gene	Pollen tube growth	GO:0009860	1089–1108	1099
164a/b/c/d/e/f	MDP0000269941	PLC‐like phosphodiesterase family protein (PDL2)	Cell tip growth	GO:0009932	24–44	34
	MDP0000121265	NAC domain‐containing protein 100(TNAC5 OR NAC100)	Multicellular organismal development	GO:0007275	645–664	655
	MDP0000911724	NAC domain‐containing protein 100(TNAC5 OR NAC100)	Multicellular organismal development	GO:0007275	645–664	655
	MDP0000298182	NAC domain‐containing protein 1(ac021,ANAC022,NAC1)	Auxin‐mediated signalling pathway	GO:0009734	618–637	628
	MDP0000528658	NAC domain‐containing protein 1(ac021,ANAC022,NAC1)	Shoot apical meristem specification	GO:0010072	528–547	538
319a/b/c	MDP0000243495	TCP family transcription factor 4(EE35,TCP4)	Cell differentiation	GO:0030154	1540–1559	1550
	MDP0000442611	TCP family transcription factor 4 (EE35,TCP4)	Cell differentiation	GO:0030154	841–860	852
	MDP0000328318	TCP family transcription factor 4 (EE35,TCP4)	Cell differentiation	GO:0030154	988–1007	999
	MDP0000184743	TCP family transcription factor 4 (EE35,TCP4)	Cell differentiation	GO:0030154	1489–1508	1490
	MDP0000916623	TCP family transcription factor 4 (EE35,TCP4)	Positive regulation of development	GO:0045962	861–881	872
	MDP0000287069	TEOSINTE BRANCHED 1, cycloidea and PCF transcription factor 2(TCP2)	Cell differentiation	GO:0030154	1180–1199	1190
	MDP0000920127	TEOSINTE BRANCHED 1, cycloidea and PCF transcription factor 2(TCP2)	Cell differentiation	GO:0030154	1186–1205	1197
396 a/b/c/d/e/f/g	MDP0000125282	Growth‐regulating factor 1 (GRF1)	Leaf development	GO:0048366	776–795	786
	MDP0000215583	Growth‐regulating factor 1 (GRF1)	Transcription activator activity	GO:0016563	755–774	756
	MDP0000808163	Growth‐regulating factor 4 (GRF4)	Transcription activator activity	GO:0016563	546–567	556
	MDP0000274400	Growth‐regulating factor 5 (GRF5)	Transcription activator activity	GO:0016563	351–372	361
	MDP0000243533	Growth‐regulating factor 5 (GRF5)	Transcription activator activity	GO:0016563	348–369	358
	MDP0000814056	Growth‐regulating factor 7 (GRF7)	Transcription activator activity	GO:0016563	630–651	641
	MDP0000261112	Growth‐regulating factor 8 (GRF8)	Transcription activator activity	GO:0016563	312–333	322
	MDP0000322576	Growth‐regulating factor 8 (GRF8)	Transcription activator activity	GO:0016563	600–621	610
398a/b/c	MDP0000178745	GATA type zinc finger transcription factor family protein(WLIM1)	Actin filament binding	GO:0051015	935–954	946
	MDP0000530255	Ctr copper transporter family (AT2G26975.1)	Copper ion transport	GO:0006825	33–54	44
	MDP0000217046	Ctr copper transporter family (AT2G26975.1)	Copper ion transport	GO:0006825	24–45	35
172a‐o	MDP0000181606	Related to AP2.7 (RAP2.7,TOE1)	Vegetative to reproductive phase transition of meristem	GO:0010228	1250–1270	1261
	MDP0000200319	Related to AP2.7 (RAP2.7,TOE1)	Organ morphogenesis	GO:0009887	1262–1282	1273
	MDP0000296716	Related to AP2.7 (RAP2.7,TOE1)	Organ morphogenesis	GO:0009887	1656–1576	1267
	MDP0000163645	Related to AP2.7 (RAP2.7,TOE1)	Organ morphogenesis	GO:0009887	1262–1282	1273
	MDP0000137561	Integrase‐type DNA‐binding superfamily protein (AP2,FL1,FLO2)	Flower development	GO:0009908	1637–1657	1648
2111a/b	MDP0000399538	Galactose oxidase/kelch repeat superfamily protein (AT3G27150.1)	Biological process unknown	GO:0008150	228–248	238
	MDP0000892979	Galactose oxidase/kelch repeat superfamily protein (AT3G27150.1)	Biological process unknown	GO:0008150	198–218	208
	MDP0000895750	Galactose oxidase/kelch repeat superfamily protein (AT3G27150.1)	Biological process unknown	GO:0008150	219–239	229
2118a/b/c	MDP0000119199	Heat shock transcription factor A6B (AT‐HSFA6B,HSFA6B)	Regulation of transcription	GO:0006355	757–776	768
408a/b/c/d	MDP0000150953	Plantacyanin (ARPN)	Unknown	Unknown	6–25	16
	MDP0000124552	Plantacyanin (ARPN)	Unknown	Unknown	297–316	307
7124a/b	MDP0000319328	P‐loop containing nucleoside triphosphate hydrolases superfamily protein (AT5G60930.1)	Microtubule‐based movement	GO:0007018	728–748	738
7125c	MDP0000273257	Zinc transporter 1 precursor (ZIP1) AT3G12750.1	Response to zinc ion	GO:0010043	133–153	143
	MDP0000243546	Iron‐sulphur cluster biosynthesis family protein (AT5G03900.2)	Biological process unknown	GO:0008150	1179–1198	1189
828a/b	MDP0000578193	myb domain protein 66 (ATMYB66,MYB66,WER,WER1)	Transcription regulator activity	GO:0030528	312–331	322
	MDP0000253904	myb domain protein 5(ATMYB5,MYB5)	Proanthocyanidin biosynthetic process	GO:0010023	325–346	335
	MDP0000124555	myb domain protein 66 (ATMYB66,MYB66,WER,WER1)	Transcription regulator activity	GO:0030528	369–388	379
	MDP0000650225	myb domain protein 66 (ATMYB66,MYB66,WER,WER1)	Transcription regulator activity	GO:0030528	307–328	318
	MDP0000931057	high response to osmotic stress 10(HOS10,MYB8)	Response to salicylic acid stimulus	GO:0009751	367–388	378
	MDP0000167314	myb domain protein 82 (AtMYB82,MYB82)	Regulation of transcription	GO:0006355	343–364	354
858	MDP0000253904	myb domain protein 5(ATMYB5,MYB5)	Proanthocyanidin biosynthetic process	GO:0010023	293–312	303
	MDP0000184538	myb domain protein 3(ATMYB3,MYB3)	Response to abscisic acid stimulus	GO:0009737	299–318	309
	MDP0000318013	Duplicated homeodomain‐like superfamily protein (ATMYB123)	Regulation of transcription	GO:0006355	233–252	243
	MDP0000437717	Duplicated homeodomain‐like superfamily protein (ATMYB123)	Proanthocyanidin biosynthetic process	GO:0010023	293–312	303
	MDP0000226215	myb domain protein 5(ATMYB5,MYB5)	Regulation of transcription	GO:0006355	293–312	303
	MDP0000887107	myb domain protein 12 (ATMYB12,MYB12,PFG)	Flavonoid biosynthetic process	GO:0009813	290–209	300

Table 2

Potential targets of the identified novel miRNAs from ‘Fuji’ apple ( Borkh.) buds and their reads counts in each library

Novel miRNA	Mature sequence	Targets		Control			Bending
		Gene ID	Description	CES	CMS	CLS	BES	BMS	BLS
DEGs
novel‐m0994‐5p	UUGGAAGCAGAUUUGCAAAUC	MDP0000800338	Translation initiation factor 2, small GTP‐binding protein(FUG1)	16	0	0	0	0	0
		MDP0000245786	NAD(P)‐binding Rossmann‐fold superfamily protein(AT4G24050.1)
novel‐m0088‐3p	UUGUUCAGCUUGAAGAUUCCU	MDP0000162729	Pyridoxal phosphate (PLP)‐dependent transferases superfamily protein	0	40	2	24	33	7
novel‐m0918‐5p	CAUCUCUGUUUUUCUCUGGCAG	MDP0000255472	Protein kinase superfamily protein(AT1G53050.1)	26	25	5	46	69	1
		MDP0000239152	SBP (S‐ribonuclease binding protein) family protein(AT4G35070.1)
novel‐m0893‐5p	AACUGGCUGAAAGAGGACUGC	MDP0000281445	Chaperone DnaJ‐domain superfamily protein(AT5G49580.1)	1	34	17	51	29	12
novel‐m0851‐5p	AUUUUCAACCGUCGGAUCGGAGAA	MDP0000250737	Response regulator 3(ARR3)	181	146	88	142	216	64
novel‐m1401‐5p	CGGUGCGGCUGCUUCGUUCG	MDP0000202468	Protein kinase superfamily protein(AT5G18910.1)	4	1	0	72	1	1
		MDP0000230379	Pentatricopeptide repeat (PPR) superfamily protein(AT4G13650.1)
		MDP0000252034	Protein of unknown function (DUF740)
novel‐m2047‐3p	UUUCUGGAUUUCUGCUGCUCC	MDP0000319833	Major facilitator superfamily protein(AT2G48020.1)	21	27	21	21	52	17
novel‐m0624‐5p	AAAGGGAAUUGUUAUUAGCGCUCC	MDP0000245102	DNA‐directed DNA polymerases(REV1)	6	18	14	9	4	7
novel‐m1913‐5p	UUGGUUUUAGCUUUGGAAUCC	MDP0000263746	Protein of unknown function (DUF3741) (AT5G02390.1)	22	18	7	3	5	4
novel‐m1736‐3p	CUUGGAUCGAAUGGAGCUCC	MDP0000552725	myb‐like HTH transcriptional regulator family protein(DUO1)	28	11	16	11	5	4
novel‐m2000‐3p	UGCCGAUGUGUGGAUCUCCCACA	MDP0000753788	Ribosomal protein S13A(PFL2)	0	1	5	5	18	5
novel‐m1953‐3p	CGACGACUCGGACUUAAGCCG	MDP0000870778	Malectin/receptor‐like protein kinase family protein(FER)	1	4	25	8	39	33
novel‐m0017‐3p	CGUCUUGUGAUUUCAUUGCAAUA	MDP0000173207	K+ efflux antiporter 3(KEA3)	12	3	4	5	19	2
novel‐m0163‐5p	ACGAGUCGUGUUGUGAAUGGCUUC	MDP0000195055	Transducin/WD40 repeat‐like superfamily protein(AT2G20330.1)	182	113	81	226	169	100
		MDP0000388416	annexin 4(ANNAT4)
		MDP0000497339	Transducin/WD40 repeat‐like superfamily protein(AT2G20330.1)
novel‐m0770‐3p	UCUUUGACGUAGGGUUUGUGG	MDP0000286528	P‐glycoprotein 2(PGP2)	10	10	9	2	26	15
		MDP0000284229	TLD domain‐containing nucleolar protein(AT2G05590.2)
		MDP0000715898	UDP‐Glycosyltransferase superfamily protein(AT3G21790.1)
novel‐m1317‐3p	UCUGGAGAGAGGGUUUCAUCC	MDP0000202292	WRKY DNA‐binding protein 35(WRKY35)	8	13	22	0	0	0
novel‐m0014‐5p	UUUAGGGAGGUUCUGCUGAUG	MDP0000149064	ABC‐2 type transporter family protein	1	0	0	2	10	1
novel‐m0056‐3p	UUUGAUUUCGAUUUCUUACGUGC	MDP0000159011	myb domain protein 6(MYB6)	0	0	2	0	1	4
novel‐m0096‐3p	AACCCUUUAGCCGAAACUUGGACC	MDP0000158257	Transducin/WD40 repeat‐like superfamily protein(TTG)	0	2	1	0	0	0
novel‐m0118‐5p	UUGUGCCAGCCUCUGUAUUUC	MDP0000772420	Expansin A8(EXPA8)	0	1	0	0	0	0
novel‐m0173‐3p	CAACAGCUGGGAUUGAAACUU	MDP0000237114	Associated molecule with the SH3 domain of STAM 1(AMSH1)	0	1	0	8	7	22
novel‐m0267‐3p	AGAGGAGGACGUCGCCGGAGA	MDP0000146449	SOS3‐interacting protein 4(SNRK3.22)	1	0	0	1	1	1
novel‐m0530‐3p	GCCUAUGAAGGACGUGUUGCU	MDP0000152242	Triosephosphate isomerase(TP1)	0	1	1	2	1	1
novel‐m1018‐3p	UCUCUGUUGCUGUGAACUAC	MDP0000013380	Isopentenyltransferase 3(IPT3)	3	0	1	0	0	0
novel‐m1026‐3p	AGUAGCUCUGGAUCGGCGACGUGU	MDP0000120398	Auxin‐responsive family protein	0	0	1	1	1	1
novel‐m1064‐5p	UCAUGUCGGGCAAAACAGAGC	MDP0000286933	Cytochrome P450 superfamily protein(CYP75B1)	1	3	0	0	0	0
novel‐m1285‐3p	UCAGGGCGAGCUUAAAUUUGC	MDP0000198561	Carotenoid cleavage dioxygenase 1(NCED1)	1	0	0	1	1	0
novel‐m1316‐5p	AUUACCUUUCUAAUAUAAGUG	MDP0000201058	Beta‐galactosidase 3(BGAL3)	0	0	5	0	0	0
novel‐m1563‐3p	CUAGAUCUGGGGCCUGUUUUCUCC	MDP0000201058	Beta‐galactosidase 3(BGAL3)	0	0	0	0	4	0
novel‐m1727‐5p	GCAUUUGACGGAGAUGACGAU	MDP0000256514	WRKY DNA‐binding protein 48(WRKY48)	1	1	0	0	0	1
novel‐m1791‐5p	AAUCGACUUGGUCAAUUGCCU	MDP0000308379	Cytokinin response factor 4(CRF4)	0	0	0	0	5	0
novel‐m1973‐5p	UGCCUAAAGGGUUUGAAGGAUU	MDP0000327940	Light harvesting complex photosystem II subunit 6(CP24)	1	3	0	0	0	0

Figure 10

The cleavage information of mdm‐miR396 and its GRF targets. (a) Multiple alignment of mdm‐miR396 target site in GRF genes; (b) multiple alignment of mdm‐miR396 family members; (c) an example of miR393 targeting GRF genes.

The targets of novel miRNAs included TFs and regulatory proteins (Tables 2 and S5). For example, novel‐m1736‐3p targeted genes encoding a myb‐like HTH transcriptional regulator family protein (DUO1), novel‐m0014‐5p regulated an ABC‐2 type transporter family protein, novel‐m1064‐5p regulated a cytochrome P450 superfamily protein (CYP75B1) and novel‐m0056‐3p targeted a TF of the MYB6 gene (Table 2). In addition, some novel miRNAs regulated targets associated with plant hormones. Novel‐m1285‐3p targeted a gene encoding a carotenoid cleavage dioxygenase 1 (NCED1) involved in ABA synthesis, and novel‐m1791‐5p regulated a cytokinin (CK) response factor 4 (CRF4) gene involved in CK signalling. Additionally, novel‐m1018‐3p and novel‐m0851‐5p regulated the isopentenyltransferase 3 (IPT3) involved in CK synthesis and the response regulator 3 (ARR3) associated with CK signalling, respectively (Table 2). Two novel miRNAs, novel‐m1317‐3p and m1727‐5p, regulated WRKY35 and WARK48 genes involved in stress responses in plants (Table 2). Meanwhile, several novel miRNAs targeted genes involved in sugar signalling and metabolism, such as BGAL3 (novel‐miR1316‐5P), SOS3‐interacting protein 4 (SNRK3.22) (novel‐m0267‐3p) and triosephosphate isomerase (TP1) (novel‐m0267‐3p) (Table 2).

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with AUX in buds in response to shoot bending

To examine the expression levels of miRNAs and their targets in buds during different developmental stages (ES, MS and LS), both in control and shoot‐bending treatments, as well as to confirm the miRNA sequencing results, we examined the expression levels of three miRNAs and six targets associated with IAA in response to shoot bending and time changes by qRT‐PCR (Figure 11). In early bud growth (ES and MS), the expression level of mdm‐miR396 is significantly higher in the control than in the shoot‐bending treatment. It increased from ES to MS in both treatments, but showed a lower expression level in the LS of flower induction, and a relatively higher expression level in the shoot‐bending treatment (Figure 11). While two genes (mdm‐miR396's targets), GRF7 and GRF8, displayed higher expression levels in the shoot‐bending than in the control treatment in ES and MS, they showed an opposite result in the LS of flower induction (Figure 11). Mdm‐miR160 showed a higher expression level in the LS of flower induction and displayed a significantly higher expression level in the shoot‐bending than in the control treatment during MS and LS (Figure 11). The expression level of the mdm‐miR160 target ARF16 decreased from ES to LS and showed a higher expression level in the control than in the shoot‐bending treatment during the flower induction process (Figure 11). Another target of mdm‐miR160, ARF17, only had a similar result in the LS of bud growth (Figure 11). The expression of mdm‐miR393 was significantly increased in the shoot‐bending treatment compared with the control during bud growth (Figure 11), and its targets (TIR1 and AFB3) showed opposite results over the same time period (Figure 11). The miRNAs and targets that are involved in IAA signalling, and their levels in buds, may participate in the regulation of apple flower bud formation in response to shoot bending.

Figure 11

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with auxin in buds in response to shoot bending.

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with ABA in buds in response to shoot bending

Additionally, we examined, using qRT‐PCR, the expression levels of three miRNAs and four targets associated with ABA in response to shoot bending and time changes (Figure 12). The expression of mdm‐miR319 increased from ES to LS and displayed a higher expression level in the shoot‐bending treatment than in the control in the later stages of bud growth (MS and LS) (Figure 12). However, their target, TCP, showed a decrease in expression from ES to LS in both control and shoot‐bending treatments, and in the shoot‐bending treatment during flower induction (Figure 12). Mdm‐miR164b displayed a decrease in expression from ES to LS, but the target, NAC1, increased at the same time in both control and shoot‐bending treatments. The expression level of mdm‐miR159b decreased from ES to LS and showed a higher expression level in the control than in the shoot‐bending treatment (Figure 12). Its targets, MYB65 and MYB101, increased in expression under the shoot‐bending treatment compared with under control conditions during the flower induction process (from ES to LS) and also increased from ES to LS in both control and shoot‐bending treatments (except MYB101 during LS) (Figure 12). The miRNAs and targets that are involved in ABA signalling, as well as their levels in buds, may participate in the regulation of apple flower bud formation in response to shoot bending.

Figure 12

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with ABA in buds in response to shoot bending.

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with flowering control in response to shoot bending

In our results, the expression of mdm‐miR156ab is significantly decreased in the shoot‐bending treatment compared with the control during the flower induction process (Figure 13). However, two mdm‐miR156ab target genes, SPL2 and SPL9, displayed relatively higher expression levels in the shoot‐bending treatment compared with the control in later stages of bud growth (MS and LS) (Figure 13). In addition, the expression of mdm‐miR172b decreased from ES to LS, but was increased in shoot‐bending compared with the control treatments during flower induction (Figure 13). However, the targets, AP2 and AP2‐like genes, had significantly lower expression levels in the shoot‐bending treatment than in the control, and the AP2 gene's expression increased from ES to LS in both shoot‐bending and control treatments (Figure 13). These two miRNAs and their targets may be involved in the ‘Fuji’ apple trees flower induction in response to shoot‐bending treatment.

Figure 13

Identification by qRT‐PCR of differentially expressed miRNAs and their targets associated with flowering control in response to shoot bending.

Expression of flowering control genes in buds in response to shoot bending

The expression levels of flowering control genes were analysed in ‘Fuji’ apple buds between control and shoot‐bending treatments during flower induction (Figure 14). The expression of the flowering‐time transcription factor FD gene in the control treatment increased from ES to LS, and it showed a higher expression level in the shoot‐bending treatment in MS and LS. The expression level of the FD gene in the shoot‐bending treatment was significantly higher than in the control treatment in the ES, MS and LS during bud growth (Figure 14). The expression level of the FLOWERING LOCUS T (FT) gene increased from ES to MS in both control and shoot‐bending treatments and showed higher expression levels in the shoot‐bending than in the control treatment during the two stages (Figure 14). The expression level of the MADS‐box SOC1 gene decreased from ES to LS, but had a significantly higher expression level in the shoot‐bending treatment than in the control treatment in the MS and LS of bud growth (Figure 14). By contrast, the expression level of the AFL1 (apple FLORICAULA/LFY) gene increased from ES to LS and had a significantly higher expression level in the shoot‐bending treatment than in the control treatment from ES to LS during flower induction (Figure 14). Additionally, the expression level of the APETALA1 (AP1) gene also increased from ES to LS and showed a higher expression level in the shoot‐bending treatment than in the control treatment in earlier stage of ES and MS, but an opposite result was seen in LS (Figure 14).

Figure 14

Identification by qRT‐PCR of flowering control gene expression levels in ‘Fuji’ apple buds in response to shoot bending.

Discussion

Bud growth, flower induction and flowering rate are affected by shoot bending

Bud growth and flower induction in apple (Malus domestica Borkh.) trees, as a woody perennial plant, play important roles in the life cycles (Kurokura et al., 2013). Our results showed that shoot bending significantly promoted bud growth (the length, width and fresh weight per bud) during the flower induction process (ES to LS) compared with the control treatment (Figure 1a–c). In addition, the growth, as measured by bud length, width and fresh weight, was mainly concentrated in the early stages, from ES to MS, of the flower induction process in both shoot‐bending and control treatments (Figure 1d–f). Additionally, the treatment of shoot bending in apple trees accelerated bud growth also in this growth stage (ES to LS) (Figure 1d–f). Previous studies have shown that bending shoots significantly increased bud growth in plants such as apple (Lauri et al., 2008), pear (Ito et al., 2004) and rose (Kim and Lieth, 2004). Our results also showed that the flowering rate was significantly increased in the shoot‐bending treatment compared with the control, and the growth rate increased by ~58% (Figure 2). Similar results can be seen in other studies (Grochowska et al., 2004; Lauri and Lespinasse, 2001). Thus, the results indicated that shoot‐bending increased the flowering rate by stimulating the bud growth much earlier compared with the control treatment during flower induction in ‘Fuji’ trees.

Hormone contents changed in buds during flower induction in response to shoot bending

Our results showed that IAA, CK, GA and ABA levels in buds during flower induction displayed significant changes (Figure 3) and that the IAA, GA and CK contents decreased in buds from ES to LS, but the ABA content had an opposite trend during this process (Figure 3). Meanwhile, studies have shown that the ABA content increases during flower bud differentiation, especially during flower organogenesis, over its level in the previous periods (Guo et al., 2012). Other studies also showed that CK was mainly involved in the initiation of flower bud induction in plants (Lee and Lee, 2010), and that ABA was as an important hormone involved in the regulation of seed development, and floral and phase transitions in responses to environmental stresses (Tsai and Gazzarrini, 2014). In addition, in our study, the CK content in buds was significantly higher in the shoot‐bending treatment than in the control treatment in the ES and MS, and the ABA content had a similar result in the MS and LS (Figure 3b,d). Studies have shown that the CK content in lateral buds was higher in bent shoots than in vertical shoots and formed a greater numb er of flower buds in bending shoots. The acceleration of the flowering rate in bending shoots is accompanied by notable increases in the cytokine level in buds compared with the vertical shoots (Ito et al., 1999). According to our results and previous studies, shoot bending promotes flower bud formation in ‘Fuji’ apple trees by increasing CK and ABA levels during the flower induction process, and both hormones play an active role in flower induction and flower bud formation (Han et al., 2007; Susawaengsup et al., 2011) (Figure 3). However, our results indicated that the GA level in buds was significantly lower in the shoot‐bending treatment than in the control during flower induction (Figure 3c) and that GA may play a negative role in flower bud formation. A similar result can be seen in olive (Ulger et al., 2004), apple (Wilkie et al., 2008) and ‘Satsuma’ mandarin (Koshita and Takahara, 2004; Koshita et al., 1999). Studies showed that IAA played an important role in the regulation of bud growth and flower induction (Balzan et al., 2014). Our results showed that the IAA content in buds was significantly higher in the shoot‐bending treatment than in the control in the ES and MS of flower induction, but that both shoot‐bending and control treatments had a relatively lower level in the LS of bud growth (Figure 3a). Moreover, studies showed that the stimulated bud growth in bent shoots was related to the increased sink capacity of the bud relative to the adjacent shoot tissues and the relatively higher IAA level (Han et al., 2007). In addition, during the vegetative to inflorescence bud stages, IAA increased by 581% in the shoot tips, the most significant change during flowering (Zhang et al., 2014).

Identification of known and novel miRNAs, as well as their expression profiles

In this study, six sRNA libraries of ‘Fuji’ apple buds under shoot‐bending and control treatments during the flower induction process (ES, MS and LS) were constructed using high‐throughput sequencing. In total, 195 known miRNAs that belonged to 41 miRNA families were identified (Figure 4 and Data S1). The majority were conserved in plants such as A. thaliana (Meyers et al., 2008) and apple (Xia et al., 2012b). Some known miRNAs family had more than one member, such as mdm‐miR156 (27 members), 171 (14 members) and 172 (15 members), with the mdm‐miR156 family having the most members (Figure 4). Similar results can be seen in Xing et al. (2014). In addition, the expression levels of these identified known miRNAs in the six libraries between shoot‐bending and control treatments were divided into eight different categories based on their read content's frequencies, with almost no differences among the libraries (Figure 5). This result implied that these differentially expressed miRNAs were involved in the regulation of bud growth and flower induction in response to environmental signals, which was similar to the results of other studies (Ding et al., 2013; Zhu et al., 2011).

In our study, some known miRNAs, including mdm‐miR535b, 160b, 156, 172, 398a and 393c, in buds had significantly different expression levels over times in both shoot‐bending and control treatments during the flower induction process (Figure 6a–c), suggesting that bud growth regulation and flower induction in ‘Fuji’ apple trees involved in these known miRNAs (Table 1). Additionally, these differentially expressed miRNAs, including mdm‐miR156, 159, 160, 393, 164 and 172, may be involved in biological functions that regulate plant development and growth, hormone synthesis and metabolism, phase transition, flower induction and reproductive growth (Huang et al., 2012; Liu et al., 2010; Wollmann et al., 2010; Xing et al., 2014). Meanwhile, our results indicated that some miRNAs, such as mdm‐miR394b, 396a/b, 398c and 162b, had decreased expression levels in the shoot‐bending compared with the control treatment (Figures 6g and 7c), while others, such as mdm‐miR5225a, 399c and 160a, had increased expression levels in the shoot‐bending treatment (Figures 6i and 7e). These miRNAs, which are differentially expressed between shoot‐bending and control treatments, may also be involved in hormone regulation, flower induction and stress responses (Ding et al., 2013; Liu et al., 2010) (Table 1).

In addition, ~137 putative unique novel miRNAs were identified in each library (Data S2), which were divided into six categories based on their expression patterns (Figure 8). The difference expression patterns among these putative novel miRNAs suggested that they may play different roles in regulating plant growth and development, which involve various biological processes (Cuperus et al., 2011). And 38, including novel m1736‐3p, m0829‐3p and m0697‐3p, and 34, including novel m1424‐5p, m0546‐5p and m0559‐3p, differentially expressed novel miRNAs had decreasing and increasing expression levels, respectively, over time (Figure S3a and S3b), suggesting that these novel miRNAs may take part in regulating apple bud growth and flower induction during the physiological differentiation phase. However, 24 differentially expressed novel miRNAs, such as novel‐m1913‐5p, m1900‐5p and m1336‐3p, had increased expression levels and 23, such as novel‐m1401‐5p, m0893‐5p, m0088‐3p and m0490‐3p, had decreased expression levels in response to the shoot‐bending treatment (Figure S3c and S3d), suggesting that these novel miRNAs may be involved in the regulation of ‘Fuji’ bud growth and flower induction in response to shoot bending.

The GO analysis revealed that the targets of known and novel miRNAs were mainly associated with growth, hormone‐mediated signalling pathway, meristem structural organization, phase transition and flower development in plants (Tables 1 and 2). Among them, SPLs as mdm‐miR156 targets and AP2 as a mdm‐miR172 target played important roles in the regulation of the vegetative to reproductive growth transition and the transition to flowering, respectively (Wu et al., 2009; Yu et al., 2012). Several targets of mdm‐miR160 and miR393, such as ARF16, ARF17, TIR1, AFB2 and AFB3, are associated with the AUX‐mediated signalling pathway and the response to AUX stimulus, which is mainly involved in bud growth, flower induction and regulation, and flower development (Kim and Ahn, 2014; Liu et al., 2010) (Table 2). NAC1, as a mdm‐miR164 target (Table 1), may be involved in the regulation of multiple biological processes, such as organ development, meristem structural organization and hormone signalling responses (Guo et al., 2005; Liang et al., 2014). The targets of mdm‐miR159, MYB65 and MYB101 (Table 1), are involved in regulating flower induction and floral development through the GA pathway (Mutasa‐Gottgens and Hedden, 2009) and also take part in bud growth and phase transition in response to ABA‐signalling pathways (Ortiz‐Morea et al., 2013). In addition, the target TCP4 is involved in plant growth, flowering‐time regulation and flower induction (Schommer et al., 2012), as well as playing an important role in controlling cell proliferation in plants (Schommer et al., 2014). GRFs, as important TFs in plants (targets of mdm‐miR396) (Table 1), are involved in multiple aspects of cell division and proliferation (Casadevall et al., 2013) and flower development (Baucher et al., 2013; Liang et al., 2014).

Our results also showed that the targets of novel miRNAs included TFs and regulatory proteins in plants (Tables 2 and S5). Among them, NCED1 is a target of novel‐m1285‐3p that is involved in ABA synthesis and CRF4 (target of novel‐m1791‐5p), IPT3 (target of novel‐m1018‐3p) and ARR3 (target of novel‐m0851‐5p) are associated with the CK signalling and synthesis processes (Table 2), suggesting that these novel miRNAs are associated with the plant hormones that play important roles in bud growth and the flower induction process (Spanudakis and Jackson, 2014). In addition, the targets of novel‐m1317‐3p and m1727‐5p, the WRKY35 and WARK48 genes, respectively, are involved in stress responses in plants (Table 2), suggesting that shoot bending uses miRNAs in its regulation of flower induction and formation, and in response to environmental stresses (Zhu et al., 2011).

miRNAs associated with IAA and their targets regulate bud growth and flower bud induction in response to shoot bending

Our results showed that the IAA content in buds increased in response to shoot‐bending compared with the control in the ES and MS of bud growth (Figure 3) and that GRF7 and GRF8 (targets of mdm‐miR396b) (Figure 2), as growth‐regulatory factors, were expressed at higher levels in the shoot‐bending treatment at the same time (Figure 11). However, the changes in mdm‐miR396 levels had opposite results in shoot‐bending and control treatments (Figure 11), suggesting that IAA and GRFs (GRF7 and GRF8) had positive roles in bud growth regulation and flower induction during the early stages, while mdm‐miR396 had a negative regulatory role. In addition, they also took part in the cell division and proliferation associated with bud growth. A similar result can be seen in Casadevall et al. (2013). Indeed, studies showed that the IAA level in the shoot apical meristem was essential for floral induction and differentiation (Hou and Huang, 2005) and that GRFs play important roles in flower development (Liang et al., 2014). However, our results showed that mdm‐miR160a expression associated with the AUX‐mediated signalling pathway was increased in response to the shoot‐bending treatment (Figure 11), but the expression level of the mdm‐miR160 target ARF16 was decreased in response to the shoot‐bending treatment. Additionally, it was decreased in both the shoot‐bending and control treatments from the ES to LS (Figure 11), suggesting that the novel miRNA and its targets were involved in regulating ‘Fuji’ bud growth and flower induction in response to shoot bending. In addition, it has been reported that ARFs are involved in regulating bud growth (Nagpal et al., 2005) and also influence the ABA sensitivity of plants in response to environmental stress (Nonogaki, 2008). Meanwhile, mdm‐miR393a, which is involved in the AUX response and signalling, had an increased expression level in shoot‐bending compared with the control treatments, but its two targets, TIR1 and AFB3, had a decreased expression level in the shoot‐bending treatment (Figure 11), suggesting that the IAA signalling involved in bud growth and flower induction is regulated by miRNAs. In addition, studies showed that miR393 and its targets, TIR1 and AFB2, are associated with the regulation of flower development in plants in response to stress (Xia et al., 2012a), and they play important roles in multiple biological processes by regulating the ABA and IAA signalling pathways involved in biotic and abiotic stress responses (Bian et al., 2012; Navarro et al., 2006).

miRNAs associated with ABA and their targets regulate bud growth and flower bud induction in response to shoot bending

In this study, the ABA level increased in both shoot‐bending and control treatments from ES to LS and displayed a higher level in the shoot‐bending treatment than in the control (Figure 3). Additionally, more flower buds formed under the shoot‐bending treatment (Figure 2), suggesting that the increase in the flowering rate in ‘Fuji’ apple trees in response to shoot bending may involve ABA signalling and responses during bud growth and flower induction. In addition, it was reported that miRNAs play crucial roles in a broad range of developmental processes in plants, such as germination, bud growth and flower induction, in response to the stress associated with ABA and AUX (Martin et al., 2010). Our results showed that the expression levels of these miRNAs, including mdm‐miR164, 159, 319, and their targets, NAC1, MYB65, MYB101 and TCP4, display great differences between the shoot‐bending and control treatments during bud growth and flower induction (Table 1 and Figure 12), suggesting that the flower induction in ‘Fuji’ apple trees in response to shoot bending may involve these regulatory miRNAs and their associated hormones. Studies have shown that NAC TFs play significant roles in plant development and stress signalling (Jensen and Skriver, 2014) and also play an active role in bud growth and initiating the shoot apical meristem, which causes the axillary meristems to develop into flowers (Raman et al., 2008). A similar result was seen in our study, in which a higher flowering rate and NAC1 expression level existed in the shoot‐bending treatment (Figures 3 and 12). In addition, down‐regulated AUX signals were caused by up‐regulated NAC1 mRNA expression levels regulated by miR164 (Guo et al., 2005), and the decreased expression levels of targets associated with IAA signalling and response were found in the shoot‐bending treatment (Figure 12), indicating that this miRNA and target are also involved in regulating flower induction in ‘Fuji’ apple trees by down‐regulating IAA signalling and response.

It has been reported that the miR159‐ and miR319‐regulated MYB and TCP TFs pattern involved in flower induction and floral formation is associated with hormones such as ABA and GA (Rubio‐Somoza and Weigel, 2013). Indeed, mdm‐miR319a displayed a higher expression level in shoot‐bending than in control, but TCP4 (target of mdm‐miR319) had an opposite result between shoot‐bending and control treatments (Figure 12). Additionally, the over‐expression of TCP4 resulted in an increase in the miR396 expression level and a corresponding decrease in GRF expression levels (Rodriguez et al., 2010), which was similar to the results of our study (Figures 11 and 12), indicating that mdm‐miR319a and its TCP4 target play important roles in the regulation of bud growth, cell proliferation and flower induction, which is associated with miR396 and its targets. In addition, miR159 and its target GAMYB play roles in regulating ABA‐signalling pathways during bud outgrowth (Ortiz‐Morea et al., 2013), and higher expression levels of MYB65 and MYB 101 were found in the shoot‐bending compared with the control treatment (Figure 12). Additionally, the expression of mdm‐miR159b produced opposite results in both treatments (Figure 12), suggesting that these miRNAs and their targets are involved in regulating bud growth and the flower induction process in ‘Fuji’ apple trees. It was reported that the deregulation of MYB33 and MYB65 in vegetative tissues inhibits growth by reducing the cell proliferation that is involved in GA signalling (Alonso‐Peral et al., 2010) and that miR159 and its targets are also involved in the regulation of plant floral induction in response to GA signalling (Mutasa‐Gottgens and Hedden, 2009). In addition, the increasing levels of miR159 caused a reduction in LEAFY transcript levels associated with up‐regulated GAMYB activity (Achard et al., 2004), and a similar result can be seen in our study (Figure 12). These results involve multiple regulatory processes in response to shoot bending, giving us a better understanding of the role played by miRNAs in the hormone crosstalk involved in bud growth and flower induction in ‘Fuji’ apple trees.

Shoot bending alters the expression of flowering control genes and miRNAs in buds during flower induction

Previous studies showed that the regulation of the floral transition involved two evolutionarily conserved miRNAs, miR156 and miR172, in both annual and perennial plants (Wang et al., 2011; Yu et al., 2012). Indeed, our results showed that the expression of miR172b was significantly higher in the shoot‐bending treatment than in the control treatment, which had almost no expression (Figure 13), but the AP2 target displayed an opposite result in both the shoot‐bending and control treatments (Figure 13). Additionally, relatively higher expression levels of SPL2 and SPL9 (targets of mdm‐miR156z) were found in the shoot‐bending treatment (Figure 13), suggesting that the partners of mdm‐miR156 and 172, as well as their targets, are involved in regulating growth and flower induction in response to shoot bending.

It was reported that FD and FT, as interdependent partners and long‐distance signals involved in plants flowering (Wahl et al., 2013), acted at the shoot apex to promote the floral transition and to initiate floral development by up‐regulating AP1, a floral meristem identity gene (Abe et al., 2005). Additionally, the overexpression of MdFT1 in apple caused earlier flowering than in wild‐type plants (Kotoda et al., 2010). Indeed, our results showed the formation of more flower buds and higher expression levels of FT, FD and AP1 genes in the shoot‐bending treatment than in the control treatment (Figures 2 and 15). Studies showed that SOC1 played a key role in integrating multiple flowering signals, such as hormone and age‐related signals (Lee and Lee, 2010). In addition, SOC1 was up‐regulated in the age‐dependent flowering pathway by SPL9 and miR156 (Porri et al., 2012), and SOC1 was induced at the shoot apex along with AGL24 to activate LEAFY gene expression and then together they induced AP1 gene expression (Liu et al., 2008). Our results also showed that the expression levels of SOC1 and AP1, as well as other flowering‐related genes (FT, SPL2, and SPL9), were relatively higher in shoot‐bending compared with the control treatments (Figures 13 and 14), suggesting that shoot bending promoted flower induction and flowering in ‘Fuji’ apple trees by up‐regulating these flowering genes, which are involved in multiple pathways.

Figure 15

Hypothetical model for the regulation of flower induction in apple by miRNAs association with phytohormone crosstalk.

In conclusion, a comprehensive study on apple (Malus domestica Borkh.) miRNAs related to bud growth, and flower induction and formation in response to shoot bending were performed, and a hypothetical model for the regulation of flower induction by miRNAs associated with phytohormone crosstalk can be seen in Figure 15.

Figure S1 The length distribution of small RNAs in each apple ( Borkh.) bud library.

Figure S2 The length distribution of small RNAs in the ‘Fuji’ apple (a Borkh.) degradome library.

Figure S3 Venn diagrams of differentially expressed novel miRNAs in ‘Fuji’ apple buds. (a) miRNA with decreased expression levels over time; (b) miRNAs with increased expression levels over time; (c) miRNAs with decreased expression levels between shoot‐bending and control treatments; (d) miRNAs with increased expression levels between shoot‐bending and control treatments.

Table S1 List of primers used in this study.

Table S2 Summary of the sequencing data in each library.

Table S3 Detailed information on the sRNA in each library.

Table S4 Summary of the degradome library sequencing data.

Table S5 GO analysis showing that targets of differentially expressed novel miRNAs from ‘Fuji’ apple ( Borkh.) buds were involved in biological processes.

Click here for additional data file.

Data S1 Identified known miRNAs in ‘Fuji’ apple (a Borkh.) and their read counts.

Data S2 Mature sequences of identified known miRNAs in ‘Fuji’ apple (a Borkh.).

Data S3 Identified novel miRNAs in ‘Fuji’ apple (a Borkh.) and their read counts.

Data S4 Mature sequences of identified novel miRNAs in ‘Fuji’ apple (a Borkh.).

Data S5 Targets of identified known miRNAs in ‘Fuji’ apple ( Borkh.).

Data S6 Targets of identified novel miRNAs in ‘Fuji’ apple ( Borkh.).

Data S7 Detailed information on distinct transcripts targeted by identified known miRNAs detected in ‘Fuji’ apple ( Borkh.).

Data S8 Detailed information on distinct transcripts targeted by identified novel miRNAs detected in ‘Fuji’ apple ( Borkh.).

Click here for additional data file.