Transcriptomic Analysis Reveals Hormonal Control of Shoot Branching in Salix matsudana

: Shoot branching is regulated by axillary bud activities, which subsequently grow into branches. Phytohormones play a central role in shoot branching control, particularly with regard to auxin, cytokinins (CKs), strigolactones (SLs), and gibberellins (GAs). To further study the molecular basis for the shoot branching in Salix matsudana , how shoot branching responds to hormones and regulatory pathways was investigated, and potential genes involved in the regulation of shoot branching were identified. However, how these positive and inhibitory processes work on the molecular level remains unknown. RNA-Seq transcriptome expression analysis was used to elucidate the mechanisms underlying shoot branching. In total, 102 genes related to auxin, CKs, SLs, and GAs were differentially expressed in willow development. A majority of the potential genes associated with branching were differentially expressed at the time of shoot branching in S. matsudana , which have more number of branching. These findings are consistent with the growth and physiological results. A regulatory network model was proposed to explain the interaction between the four hormones that control shoot branching. Collectively, the results presented here contribute to a more comprehensive understanding of the hormonal effects on shoot branching in S. matsudana . In the future, these findings will help uncover the interactions among auxin, SLs, CKs, and GAs that control shoot branching in willow, which could help improve plant structures through the implementation of molecular techniques in targeted breeding.


Introduction
Shoot branching is a highly plastic developmental process that greatly affects plant architecture. During this process, axillary buds in the axil of the leaves are formed and subsequently outgrow and develop into branches. Plant hormones are major determinants that control outgrowth of axillary buds [1,2]. Three classes of plant hormones, auxin, cytokinins (CKs), and strigolactones (SLs), interact and regulate shoot branching [3][4][5][6]. Another hormone, gibberellins (GAs), also plays a role in shoot branching control [7,8].
Shoot branching is regulated by a network of hormones, in which auxin plays a pivotal role. Two hypotheses, the second messenger hypothesis (i.e., CKs and SLs) and auxin transport canalization-based model, suggest that auxin indirectly inhibits shoot branching [9][10][11]. Auxin moves down the main stem to prevent new branches from forming, which is dependent on the polar auxin transport stream (PATS). PIN efflux carrier family-mediated auxin transport members are key components in the control of shoot branching, especially PIN-FORMED1 (PIN1). The PIN1 expression pattern follows the flow of auxin in the stem. Previously, an AUXIN RESISTANT1 (AXR1)-mediated role of auxin in shoot branching regulation was observed after axillary meristem formation in Arabidopsis [12]. Furthermore, indole-3-acetic acid (IAA) has been implicated in the primary shoot apex and inhibits axillary bud growth [13]. Moreover, the quadruple knock-out mutant, yuc1,2,4,6 (auxin biosynthesis genes YUCCA), exhibited increased branching [14].
Unlike auxin, CKs are a class of hormones that directly promote branching and are transported from roots to shoots [5,15]. The long-distance supply of CK enhances the development of an axillary bud into a branch. The isopentyl transferase (IPT) gene plays an important role in the control of CK levels in plants. Transgenic plants expressing the IPT gene, which is involved in CK biosynthesis, exhibit elevated CK levels that are often accompanied by increased branching phenotype [16]. Tanaka et al. [17] discovered that the expression level of IPT increased at the node after decapitation in pea (Pisum sativum), and therefore, proposed that lateral bud outgrowth was due to locally increased CK accumulation. Arabidopsis response regulators (ARRs), which typically act as CK signaling regulators, also play essential roles in shoot branching [18,19]. However, IPT-mediated CK synthesis in bud initiation apparently does not require the ARRs involved in bud activation [18]. Additionally, CKs act as a second messengers that relay the auxin signal to lateral buds [20]. Auxin inhibits CK synthesis in the main stem, which is dependent on the AXR1 gene in Arabidopsis [21], as well as regulates CK levels by the regulation of IPT and cytokinin oxidase (CKX) expression to control shoot branching [17,22]. SLs regulate plant architecture by directly repressing shoot branching in pea, rice (Oryza sativa), petunia (Petunia hybrida), and Arabidopsis thaliana [23,24]. Lateral bud growth is inhibited by SLs that are transported upward from the roots to shoots. Shoot branching affects SLs supported by the SL-dependent interactions between MAX2/D3/RMS4 and/or DWARF14 (D14)/DAD2 [25]. An F-box protein encoded by MAX2 and an α/β hydrolase protein encoded by D14 are locally required in aboveground plant to suppress shoot branching [26][27][28]. Studies on branching mutants have demonstrated that MAX2 and D14 are both related to the downstream SL pathway [23,29,30]. MAX2 and D14 mutants exhibit a highly branched phenotype [25,27,29]. It has also been proposed that SLs do not operate alone, but collaboratively with other hormones [31].
SLs have been proposed to function as secondary messengers for auxin that inhibit shoot branching by dampening auxin transport canalization. In Arabidopsis and chrysanthemum (Dendranthema grandiflorum), the ability of SLs to inhibit bud activation depends on the presence of competing auxin sources [32]. Furthermore, PIN1 is a plausible downstream target of the SL pathway in shoot branching regulation, as SLs repress PIN1 expression and accumulation [33,34]. The effects of SL on PIN1 depletion are dependent on MAX2 [25]. In addition to the control of auxin transport by SLs, auxins are major regulators of SL biosynthesis [4]. The relationship between auxin and SLs is the converse of auxin and CK [4,35]. CKs act as SL antagonists in lateral bud outgrowth regulation [36]. The combined effects of increased SL and reduced CK levels are proposed to lead to bud growth inhibition. SL also affects CK levels. Notably, in a previous study, the CK content of xylem sap in SL-defective mutants was greatly reduced [37]. Moreover, IPT1 expression increased in SL mutants of pea [36]. Young et al. [38] thus suggested that auxin and CKs are dominant regulators of decapitation-induced branching, while SLs are more important in intact plants.
GAs also regulate shoot branching in many plants, including Jatropha curcas, Populus, and sweet cherry (Prunus avium) [39,40]. Semi-dwarf rice mutants, which are disturbed in GA biosynthesis, exhibit a more branched shoot phenotype similar to SL mutants [41]. GA-deficient mutants displayed even greater shoot branching in rice, pea and Arabidopsis [42,43]. Conversely, GA is a positive regulator of axillary bud outgrowth in J. curcas and papaya (Carica papaya) [7]. Moreover, GA biosynthesis is crucial for axillary bud formation and activation in hybrid aspen [8,39]. In Populus, the overexpression of GA2ox produced an increased number of branches [39]. These findings indicate that the mechanism controlling shoot branching may depend on the balance of GA biosynthesis, catabolism, and receptor abundance. Auxin regulates GA biosynthesis by upregulating the expression of GA20ox or GA3ox, or by downregulating the expression of GA2ox [39,44]. GAs and CKs are positive regulators of shoot branching and collectively regulate axillary bud outgrowth in J. curcas, similar to the relationship between SLs and auxin [40]. GAs also interact with or without SLs to regulate shoot branching [7,24,41]. GAs regulate shoot branching together with SLs in J. curcas [7]. However, GA-and SL-deficient double mutants displayed stronger branching than single mutants in pea, suggesting that GAs act independently of SLs to repress branching [45].
Willows (Salix spp.) are one of the most important bio-energetic and ornamental plants worldwide [2,46]. Two willow species, S. matsudana (SM) and S. matsudana var. pseudo-matsudana (SMP), which exhibit diverse shoot branching numbers, were investigated in this study. The underlying molecular mechanisms of how auxin, through its interactions with CKs, SLs, and GAs, regulates shoot branching in willow remains elusive. Therefore, in this study, the effects of hormones on shoot branching at the molecular level were investigated by RNA-Seq coupled with digital gene expression (DGE) at different developmental stages in the two willow species. The aim of this study was to depict hormonal (i.e., auxin, CKs, SLs, and GAs) control of shoot branching in willow and identify potential genes and pathways, focusing specifically on these four hormones and their interactions.

Plant Materials
Current-year branch samples of S. matsudana and S. matsudana var. pseudo-matsudana were harvested during the growing season of 2012 in Beijing, China (39°59′18″N, 116°14′16″E). Plants at two developmental stages, when buds developed into branches on 16 April 2012 (S1) and at the time shoot branching and axillary buds developed into branches on 17 July 2012 (S2), were selected to study shoot branching in SM and SMP. Each branch was divided into two portions, the stem apex and basal branches, which were collected separately. The branches at S1 were too short, therefore, only the stem apex were harvested. Samples from three trees at each stage were harvested, frozen in liquid nitrogen, and stored at −80 °C for RNA extraction.
The separation was achieved on an ACQUITY UPLC HSS T3 column (100 × 2.1 mm, 1.7 µm; Waters). The elution gradient was carried out using a binary solvent system (solvent A, 0.1% formic acid in water; solvent B, 0.1% formic acid in acetonitrile) at a 0.4 mL/min constant flow rate. Auxin (IPA and tryptophan), CKs (iPR), and GAs (GA1, GA3, and GA4), and their internal standards were scanned in negative mode, while IAA, CKs (tZ, tZR, and iP), and SLs were analyzed in positive mode. Hormones were separated by UPLC and analyzed by ESI-MS/MS in the MRM mode. UPLC-ESI-MS/MS assays were repeated three biologically, with each repetition having three replicates. The Student's t-test were used to determine significant differences in hormone content. Statistical significance was set at p = 0.05.

Sample Preparation for Transcriptome and DGE Profiling Analysis
Two branch samples at two developmental stages of the two willow species were selected for transcriptome analysis. For each sample, an equal number of branches from three individual trees were pooled for RNA extraction using TRIZOL reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. RNA-Seq libraries were prepared using an Illumina TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA). The library was sequenced on an Illumina HiSeq TM 2000 platform (Illumina, San Diego, CA, USA). Transcript quantification from RNA-Seq data was performed in a previous study [46]. In each DGE library, tag sequences were mapped to the reference sequences generated in the RNA-Seq-based transcriptome. Gene expression levels were estimated by RSEM for each sample [47].

Screening of Differentially Expressed Genes (DEGs)
The expression of each gene was normalized to the expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM). Quantitative real-time PCR (qRT-PCR) analysis was performed to validate the expression levels of DEGs from RNA-Seq and DGE data in a previous study [46]. Based on the quantification of transcript results, nine pairwise comparisons in six DGE libraries of the two willow species at two developmental stages using two different branch portions were analyzed. Differential expression analysis of the two samples was performed using the DEGseq R package. DEGs were obtained with the thresholds: p-value ≤ 0.05 and an absolute value of the log2 (fold change) ≥ 1. The context likelihood of relatedness algorithm method was used to construct a gene regulatory network using the hormone-related DEGs [48].

Growth and Endogenous Hormone Analysis
SM and SMP buds developed into branches from S1. In SM, part of the branches stopped growing, while others (56.7%) germinated several axillary buds and developed into new branches from S2. In SMP, branches continued to grow and part of the branches (36.7%) germinated small axillary buds and developed into new branches. Additionally, the average number of shoot branching in SM (3.4) was greater than SMP (1.9).
To determine the effects of hormones on shoot branching, the four endogenous hormone contents of auxin, CKs, SLs, and GAs in the willow branches at the two development stages were compared. Hormone concentration measurements revealed that IAA levels in SM significantly increased in the stem apex compared to SMP ( Figure 1A), but decreased in the basal branches. Moreover, at S2, the IAA levels in the stem apex were significantly higher than the basal branches of SMP. Measurements of the IPA contents in the stem apex increased by 176.86% in SMP compared to the basal branches, but decreased by 42.4% in SM. Furthermore, the IPA levels at S1 were significantly higher than S2 in willows. In SM, the tryptophan levels were significantly higher than SMP. The tryptophan contents significantly increased in the stem apex at S2 compared to the basal branches in willows. The total contents of active forms of CK, including tZ and iP, were significantly lower in SM than SMP at S2 in basal branches ( Figure 1B). The content of tZ and iP in the stem apex at S2 were significantly higher than in the basal branches of SM. However, the tZR content was significantly higher in SM than SMP. The contents of iPR in SM were significantly higher at S1 and lower at S2 compared to SMP. The SLs levels in SMP decreased by 11.71% in the stem apex at S1, but increased by 41.43% in the stem apex and by 118.28% in the basal branches at S2 compared to SM ( Figure 1C). The content of SLs at S1 were higher than at S2 in SM. Furthermore, the GA1 contents in branches of SM were lower than SMP ( Figure 1D). The contents of GA3 in the stem apex of SM were significantly higher than SMP. In particular, the GA4 content markedly increased compared to GA1 and GA3. The contents of GA4 in the stem apex of SM at S2 were higher than SMP and significantly higher at S2 than at S1.

DGE Analysis and Mapping Reads to the Transcriptome
Six DGE libraries were constructed from two willow species at two different developmental stages and separated into two branch portions ( Table 1). The total number of clean reads per library ranged from 5.6-7.2 million. In order to control the quality of the libraries, the Q30 and GC content of the sequences were inspected. In each DGE library, the tag sequences were mapped to the reference sequences generated in the RNA-Seq based transcriptome [46]. Results revealed a substantial proportion of clean reads that matched (78.12%-81.47%). In total, 70,504 mapped genes were detected in the six libraries. For functional annotation, 52,578 unigenes (74.57%) were successfully annotated in at least one database (Nr, Nt, Pfam, COG, Swiss-Prot, KEGG, and GO). Q30 percentage is proportion of nucleotides with quality value larger than 30 in reads. GC percentage is proportionof guanidine and cytosine nucleotides among total nucleotides.

Screening of DEGs and Functional Annotation Analysis
To better understand the hormonal control of shoot branching in willow, plant hormones related to auxin, CKs, SLs, and GAs were used for further analysis. To identify DEGs, nine pairwise comparisons were analyzed and 5,816 DEGs were identified. Among them, 102 hormone-related DEGs were found (Table S1). The largest number of DEGs was identified in SM2B vs. SM1T (62) and SM2T vs. SMP2T (61), while the fewest were found in SM2B vs. SM2T (3).
To identify the biological functions involved in the hormonal control of shoot branching, the 102 hormone-related DEGs were mapped in the GO database. Only the biological processes were used for GO analysis. Of the 102 genes, 89 (87.25%) were assigned to 23 GO terms associated with hormone regulation (Figure 2). Many biological processes associated with responses to hormonal signals, biosynthesis, transport, and metabolism were identified. In the biological process category, auxin-activated signaling pathway (57.30%), response to auxin (31.46%), and auxin polar transport (17.98%) comprised the majority.

Regulatory Patterns of Hormone-Related DEGs
Based on expression differences, the 102 hormone-related DEGs were used for the expression pattern analysis. The number of DEGs in the stem apex between SM and SMP at S2 and between S1 and S2 in SM was greater than other pairwise comparisons (Figure 3). These results were in accordance with the growth observations, in which, the timing of shoot branching occurred at S2 and the number of branching of SM was higher than SMP.
The proportion of genes related to auxin was remarkably high. Based on the GO terms in the biological processes category associated with auxin biosynthesis, signals, transport, and metabolism, 65 auxin-related genes were confirmed ( Figure 3A). Most genes were down-regulated in SM (SM vs. SMP). The number of DEGs in the stem apex between SM and SMP at S2 was the greatest. Notably, four PIN1 gene members (c60720_g1, c61704_g5, c62767_g1, and c65489_g7) were identified. The four genes were down-regulated in SM at S2.
For SLs, only three genes (MAX2A, c59815_g2; KAI2, c56849_g1; DAD2, c59671_g1) were differentially expressed ( Figure 3C). DAD2 is an ortholog of D14 and KAI2 is a paralogue of D14 in Arabidopsis. D14 and MAX2 are involved in SL signal transduction. DAD2 was up-regulated in SM (SM vs. SMP) at S2 and down-regulated in the stem apex of SMP at S2 (S2 vs. S1). However, MAX2A and KAI2 were not differentially expressed between the two different development stages in SM and SMP or at the same stage in the stem.  Table S1.
We further analyzed the expression patterns of whole genes involved in GA biosynthesis (GA20ox) and degradation (GA2ox). Three GA20ox (c55294_g1, c55695_g1, and c57581_g1) and two GA2ox (c59623_g1 and c60122_g1) gene members were identified. Three GA20ox were up-regulated in the stem apex of SM and up-regulated in the basal branches of SMP at S2 (Figure 3D), except in the basal branches of SM, while GA2ox was down-regulated in SM and in the basal branches of SMP at S1. Additionally, the DELLA protein GAI was down-regulated in SM at S2.

Transcriptional Networks Associated with Hormone Genes
To explore the relationship between the hormone-related genes and shoot branching, we performed a gene regulation network analysis (Figure 4, Table S2). We set score above 4.5 to identify the tightly regulated relationship between all pairs of genes [48]. Hormone-related transcription factors (TFs) were analyzed as source genes. Our results revealed that most TFs were auxin-(52.78%) and GAs-related (30.56%) genes. TFs with the highest connectivity (score >5.0) included CKs (RPN12A and AHK4), an auxin-binding protein (ABP20), and GAs (SCL3 and GID1B). The CK-related TF, RPN12A, may help control the degradation of one or more factors that repress CK signaling and plays an important role in balancing cell expansion during shoot development. It is also involved in the response to auxin and was down-regulated in SM and at S2. Therefore, the increased branching in SM may be a result of the crosstalk between CK and auxin. Further detailed analyses on the function of these genes in the regulation network could provide novel insights into the hormonal control of shoot branching. Non-hormonal genes involved in the control of shoot branching were beyond the scope of this study.

Discussion
Many hormones are essential for shoot branching, in particular, auxin, CKs, SLs and GAs. In this study, we conducted a comprehensive analysis of the hormone-related genes involved in the control of shoot branching in two willow species. Auxin was the first hormone identified to play a key role in the control shoot branching. Our results showed that the amounts of IAA in the stem apex was enhanced in SM compared to SMP (Figure 1A), and the difference in this effect was observed in the basal branches. High IAA contents inhibit axillary buds outgrowth and may partly explain the dormancy of axillary buds [49,50]. Additionally, PIN genes were down-regulated in the stem apex and were always accompanied by enhanced auxin levels. Compared to SMP, the four PIN1 genes were down-regulated in SM at S2. Similarly, other PIN genes (PIN3 and PIN5) were down-regulated at S2. Previously, AXR1 mutations were found to result in increased shoot branching in Arabidopsis [12,13], causing changes in SCF TIR1/AFB dependent auxin-responsive gene expression [4]. High auxin level might also induced AXR1 gene expression [50]. However, AXR1 was not differentially expressed in two willows or at the two developmental stages. Moreover, TIR1 was down-regulated in the basal branches of SM at S2.
CKs directly promote shoot branching by activating axillary buds. Compared to SMP, the amounts of tZ and iPR decreased in SM, while iP increased (Figure 1). Additionally, tZR, tZ, and iPR decreased at S2 in SM, while iP increased. Previously, branching mutants were found to exhibit greatly reduced amounts of the ZR, Z, and iPR in the xylem sap compared to wild-type plants [37]. These results indicated that the precursors (tZR and iPR) were not effectively converted into active CKs forms (iP) in SMP due to the lower expression levels of LOG1, a CK biosynthesis related genes [51]. In this study, compared to SMP, the expression levels of IPT genes (IPT3 and IPT5) increased in the stem apex of SM (Figure 3). Compared to S1, these genes were up-regulated in SM at S2. These results suggested that CK may regulate shoot branching by locally promoting CK synthesis in the stem apex. In a previous study, IPT3 and IPT5 loss-of-function mutants exhibited reduced CK levels and significant reductions in branching [18]. These findings are consistent with the role for CK in promoting shoot branching. In this study, the RNA-Seq data revealed that CK signaling-related genes were down-regulated in SM at S2, including AHK4 and four ARRs genes (ARR2, 3, 8, and 9) in the CK response pathway. CKX functions in CK degradation and LOG functions in CK biosynthesis, both were not differentially expressed in this study. In a previous study, a mutant lacking type-A ARRs exhibited reduced shoot branching compared to wild-type plants [19].
SLs are a novel class of plant hormone that inhibit shoot branching. DAD2 may be involved in the SL signal transduction of shoot branching inhibition [28]. Mutations in petunia, rice and Arabidopsis elicit the typical highly branched SL-deficient phenotype [25,29,52]. DAD2 interacts with MAX2A in an SL-dependent manner [52]. Like DAD2, mutations in MAX2A are associated with a highly branched phenotype [27,53]. MAX2 is necessary for karrikins and SL responses, while KAI2 is not required for SL-mediated responses [30]. There may be MAX2-independent perception of SLs, possibly via KAI2 [25]. In this study, compared to SMP, DAD2 was up-regulated at the time of shoot branching (S2) in SM, and the SLs contents were lower in SM. At S2, DAD2 was down-regulated in the stem apex of SMP and MAX2A was down-regulated in the basal branches of SM. Additionally, SLs levels at S2 were lower than S1 (Figure 1). These findings are consistent with aforementioned descriptions of shoot branching numbers. MAX2 recognizes the branching-specific targets in willow [2]. It will be important for carrying out the interaction between DAD2, MAX2, and SLs in plant.
The role of GAs in shoot branching has been characterized, and its biosynthesis is crucial for shoot branching [8,41]. Mutants defective in GA biosynthesis displayed a highly branched shoot phenotype in rice [24,41]. Plants defective in GA biosynthesis showed typical GA-deficient phenotypes, such as dwarfism [42], which is often correlated with increased branching. The overexpression of GA catabolism genes (GA2ox1 and GA2ox8) reduces GA levels and leads to increased branching and dwarf phenotypes, respectively [54,55]. In this study in the stem apex, GA2ox1 and GA2ox8 were down-regulated at S2 (S2 vs. S1) in SM and SMP, and GA2ox1 was down-regulated in SM (SM vs. SMP) at S1. GA2oxs are responsible for reducing active GAs levels in plants [42,54,55]. GA4 was the dominant bioactive GA in willow [46]. Consistently, in this study, the GA4 contents were markedly higher than GA1 and GA3 ( Figure 1D). GA3 and GA4 possess similar abilities that stimulate branching from lateral buds of sweet cherry [56]. GA concentration measurements revealed that the GA3 and GA4 levels were reduced in the stem apex of SMP at S2. This result suggested that GAs inhibited branching. Conversely, GA was found to be a positive regulator of shoot branching in J. curcas due to overexpressing the GA biosynthesis gene, GA20ox1 [7]. In this study, GA20ox1 was up-regulated in the stem apex of SM (SM vs. SMP) at S2 and up-regulated at S1 (S2 vs. S1). Additionally, DELLA proteins are the main repressors of GA signaling. Recessive DELLA protein mutants exhibited reduced shoot branching and/or altered branching patterns [43]. We found that the DELLA protein GAI was down-regulated in SM at S2. These data indicate that branching in willow is considerably more sensitive to GAs at time of shoot branching (S2).
Based on the growth results, the hormone content and expression analysis of hormone-related genes, and results of previous studies on these genes in other plants [1,4,13,19,25,57], we propose a regulatory network model to explain the interaction between auxin, CKs, SLs, and GAs in the control of shoot branching in S. matsudana ( Figure 5). Collectively, auxin transported in the main stem through PAT could repress CK synthesis by inhibiting IPTs, which are dependent on AXR1 functions, and thereby inhibits bud activation [5,21,50]. With the exception of IPT, auxin regulates CK levels through the regulation CKX to control shoot branching. Additionally, CK treatment resulted in an increased accumulation of PIN3 via an unknown mechanism [19]. Moreover, AHK4 and B-type ARRs regulators are required for CK-stimulated PIN1 degradation [58]. In addition, auxin-SL interactions are critical for branching control. SLs are thought to act upstream of auxin by modulating PAT in the main shoot. Substantial auxin depletion results of increased branching are explained, at least in part by reducing SL synthesis [9]. The capacity of PATS is regulated by the MAX pathway, in which MAX1, MAX3, and MAX4 (and their homologs in pea, petunia, and rice) function in SL synthesis. SL regulates PIN1 levels through MAX2. Previously, MAX2 expression was repressed by GA and CK treatment in J. curcas [7]. SLs are also known to reduce auxin transport capacity, which is correlated with reduced PIN1 accumulation in stems. However, Ward et al. [2] found that SLs were only effective at inhibiting isolated willow buds in the presence of an apical auxin source. IAA negatively regulated shoot branching by maintaining high SL and low CK contents in intact plants. Furthermore, GA catabolism genes (GA2ox), GA biosynthesis gene (GA20ox), and DELLA protein mutants exhibited reduced or increased shoot branching. Clearly, GAs and CKs are corporately regulate shoot branching, this finding requires further investigation.

Conclusions
Our study showed that auxin, CKs, SLs, and GAs might play central roles in shoot branching control in willow. Most of DEGs were differentially expressed at the time of shoot branching happened in S. matsudana, which have more number of branching. This is well in accordance with our results of growth and hormone contents. The genes linked to auxin (PIN1, PIN3, and TIR1), CKs (IPT3, IPT5, AHK4, CKX, and ARRs), SLs (MAX2 and DAD2), and GAs (GA2ox, GA20ox, and GAI) were probably responsible for shoot branching. The genes with differential expression patterns were used to construct a regulatory network model to reveal the interaction between auxin, CKs, SLs, and GAs in the control of shoot branching in willow. These results will be helpful to the focus of our future research on the four hormones.