Transcriptome Analysis of a Multiple-Branches Mutant Terminal Buds in Betula platyphylla × B. pendula

: To investigate the molecular mechanism of the mutation of a multiple-branches birch mutant ( br ), we explored genes that were genome-wide di ﬀ erentially expressed in the main and lateral branches’ apical buds of br . The plant architecture not only has e ﬀ ects on the process of plant growth and development, but also a ﬀ ects the agronomic characters. In woody plants, branches determine the application value of timber. Therefore, analyzing genes that were di ﬀ erentially expressed in br apical buds will bring new insights to understand the molecular basis of plant architecture alteration. Wild type (WT) and the mutant br were used as materials to observe phenotype di ﬀ erences between br and the control lines (WT and OE2). The transcriptome sequencing of the main and lateral branches’ apical buds of br and controls were further performed to explore genes that were genome-wide di ﬀ erentially expressed in br . Compared to the control lines, br exhibited a multiple-branches and dwarf phenotype. In addition, biomass, rooting number, leaf area, internal diameter, and external diameter of stomata, and the size of terminal buds of br were less than that of WT and OE2. Transcriptome analysis results indicated that gene expression proﬁles of br were di ﬀ erent from the control lines. The genes that were di ﬀ erentially expressed in br apical buds were involved in multiple pathways, including organogenesis, fertility regulation, cell division and di ﬀ erentiation, plant hormone biosynthesis, and signal transduction. The multiple-branches, dwarf, and small leaves and buds of br might be due to the di ﬀ erentially expressed genes (DEGs) involved in organogenesis, cell division and di ﬀ erentiation, plant hormone biosynthesis and signal transduction.


Introduction
Plant tissues are derived from a group of cells called meristems that reside at the shoot tips and the root tips. The shoot apical meristem (SAM) is essential for the development of higher plants, and is responsible for the development of leaves, stems, and flowers. The SAM differentiates into leaf primordia at the vegetative growth stage, while it differentiates into inflorescence meristem, during the stage of transition from vegetative to reproductive growth. The axillary meristem (AM) can differentiate into the lateral bud primordia, which develops into the lateral buds and forms lateral branches of plants [1,2].
The growth and development of plants are determined by the structure and activity of SAM. Plant architecture and organ morphogenesis are influenced by both internal or environmental factors, WT, OE2, and br were cultured in a woody plant medium containing 0.2 mg/L NAA, for 30 days. Fifteen plants from each line were selected to examine the biomass and the growth traits. The stems and roots length were measured by using the scanner Image Scanner III GE-081 (GE Healthcare, Beijing, China). The 15 plants from each line were equally divided into three groups. The aboveground and underground parts of each group were separated to test the fresh weight and the dry weight, under a constant weight. The rooting index, water content (P), and dry to wet ratio (P w ) were calculated as follows: Rooting index = Roots number × Roots length × Rooting rate, Three-year-old WT, OE2, and br were used to measure the plant height, lateral branches number, and the leaf morphology. The bud morphology was observed using stereomicroscope SZX7 (Olympus, Tokyo, Japan). Each line was expanded to 15 plants. The lower epidermis of the fourth leaves of each line were torn off and used to observe the cells and stomas using an Olympus DP26 (Olympus, Tokyo, Japan). The stomatal inner diameter, outer diameter, and stomatal density were examined and each value of stomata characteristics was the average value of the three measurements from three individual plants.

Observation of Apical Buds Microstructure
The main and lateral branches' apical buds of three-year-old birch from WT, OE2, and br were fixed in FAA (formaldehyde, glacial acetic acid and 50% ethyl alcohol, V:V:V = 1:1:18) for 24 h, dehydrated in an ascending series of ethanol (70%, 80%, 85%, 95%, and 100%), cleared in xylene and embedded in paraffin. Sections of 10 µm thickness were made on an HM 340E microtome (Thermo, Waltham, MA, USA), and stained by safranine and fast green (Sigma, MO, USA). Apical bud microstructures of the samples were examined and photographed using Olympus DP26 (Olympus, Tokyo, Japan). The size of SAM was measured by using the CellSens Entry (Olympus, Tokyo, Japan). Each value of SAM size was the average of five measurements from three individual plants.

RNA Extraction, Library Construction, and Sequencing
The main and lateral branches' apical buds of three-year-old WT, OE2, and br were harvested and frozen in liquid nitrogen. Total RNA was extracted from each sample, using the CTAB method [33]. RNA concentration was detected using NanoDrop 2000 (Thermo, Waltham, MA, USA). RNA integrity was assessed using Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperatures in the NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, CA, USA). Phusion high-fidelity DNA polymerase, universal PCR primers, and index (X) Primer were used to amplify the PCR products, after connecting the sequencing connector. Library quality was assessed on the Agilent Bioanalyzer 2100 system. After cluster generation, the library preparations were sequenced on the Illumina platform, with an average reading length of PE150. RNA-seq was performed by the Anoroad Technology Company (Beijing, China) with two independent replicates. The clean reads were obtained by removing the junction sequences and low-quality sequences. Clean reads were mapped to the birch reference genome, using TopHat2 [34]. Genes were aligned with the NCBI (National Center for Biotechnology Information) non-redundant (Nr), Swiss-Prot, Pfam, Kyoto Encyclopedia of Genes and Genomes (KEGG), and the Clusters of Orthologous Groups (COG) protein databases, using Blastx with an E value ≤10 −5 . Blast2GO was used to obtain gene ontology (GO) annotation of the genes, based on the Nr annotation. Putative metabolic pathways were assigned by performing Blastx searches against the KEGG pathway database, with an E value ≤10 −5 [35][36][37]. New genes were obtained using KOBAS 2.0 and were annotated after being compared with HMMER and Pfam [38,39].
To quantify the expression level of each gene, the specific formulae was as follows: Differentially expressed genes (DEGs) between each two-sample comparison were defined with fold change ≥2 and FDR (false discovery rate) <0.01 as a threshold, according to the statistical analysis performed by DESeq R package (1.10.1) [40]. Data analysis and functional annotation were carried out by the Biomarker Technology Company (Beijing, China).

Verification of qRT-PCR
Main and lateral branches' apical buds of three-year-old WT, OE2, and br were used to extract the total RNA using the RNA pure total RNA fast isolation kit (Bioteke, Beijing, China). cDNA was synthesized from 1 µg RNA of tested lines, using a ReverTra Ace ® qPCR RT Kit (Toyobo, Osaka, Japan), according to the manufacturers' instructions, respectively. Quantitative RT-PCR (qRT-PCR) was performed on a 7500 real-time PCR system (Applied Biosystems, Darmstadt, Germany), using SYBR ® Green PCR master mix (Toyobo, Osaka, Japan). Amplification system and procedure were carried out according to Huang [41]. The relative expression level of the genes was calculated using the 2 − Ct method [42]. Primer sequences are listed in Table S1. Bp18S rRNA was selected as an internal reference gene.

Statistical Analysis
Data analysis were performed using SPSS statistics software, version 16.0 (International Business Machines, Armonk, NY, USA). Differences between the means of each line on all traits were determined using one-way analysis of variance and the Duncan multiple comparison procedure. p value <0.05 was considered to be statistically significant and was labeled with different letters. The same letter indicates that the results were not significantly different.

Variations in Biomass
Thirty-day rooting seedlings of WT, OE2, and br were used to measure biomass. The results showed that the stem length, average root length, and root number of br were significantly lower. Dry weight and fresh weight of the aboveground and underground parts of br were significantly lower than the two control lines, except the dry weight of the aboveground parts (p < 0.05) ( Table 1). The rooting index and dry to wet ratio of br were obviously decreased but the water content was significantly increased, compared with WT and OE2 (Figure 1a-c). These results confirmed that the growth of br was delayed and was weaker, and the roots were underdeveloped at the seedling stage.

Variations in Growth
To further investigate whether the phenotypes of br would recover over time, the plant growth traits of three-year-old WT, OE2, and br were measured. There was a significant reduction of br height, which was 60.73% of WT and 66.25% of OE2, respectively. The number of primary branches of br was lower, while the number of secondary branches was significantly increased, compared with WT and OE2, which was 1.76 times of WT and 1.61 times of OE2 (Figure 2a, Figure S1). The leaf size and morphology of br also showed obvious changes. The leaf area was decreased by 56.12% and 60.62% of WT and OE2, respectively. The leaf length, leaf width, and the number of veins were also reduced. The leaf base angle of br was about 180 degrees, which was smaller than the other two lines. By contrast, the leaf apex angle of br was larger than the control lines ( Figure 2b). We then observed the lower epidermis cells and stomas, it was found that both br and OE2 had irregular shape lower epidermis cells, while the stomatal inner diameter, outer diameter, and stomatal density were significantly decreased in br ( Figure 2c, Table S2). These results indicated that br was dwarf, had multiple-branches, and small leaves. These variations in growth of br were not caused by environmental factors as it could not recover over time.

Variations in Growth
To further investigate whether the phenotypes of br would recover over time, the plant growth traits of three-year-old WT, OE2, and br were measured. There was a significant reduction of br height, which was 60.73% of WT and 66.25% of OE2, respectively. The number of primary branches of br was lower, while the number of secondary branches was significantly increased, compared with WT and OE2, which was 1.76 times of WT and 1.61 times of OE2 (Figure 2a, Figure S1). The leaf size and morphology of br also showed obvious changes. The leaf area was decreased by 56.12% and 60.62% of WT and OE2, respectively. The leaf length, leaf width, and the number of veins were also reduced. The leaf base angle of br was about 180 degrees, which was smaller than the other two lines. By contrast, the leaf apex angle of br was larger than the control lines ( Figure 2b). We then observed the lower epidermis cells and stomas, it was found that both br and OE2 had irregular shape lower epidermis cells, while the stomatal inner diameter, outer diameter, and stomatal density were significantly decreased in br (Figure 2c, Table S2). These results indicated that br was dwarf, had multiple-branches, and small leaves. These variations in growth of br were not caused by environmental factors as it could not recover over time.

Observation of Morphology and the Microstructures of Apical Buds
The morphology and microstructures of the main and lateral branches' apical buds in three-year-old WT, OE2, and br were observed, as changes of growth traits in br might be related to the development of SAM. Either the main branches or the lateral branches' apical buds of br were found to be smaller, compared to the two control lines. The length and width of the br apical buds of the main branches were 62.08% and 74.43% of WT, and 67.93% and 68.04% of OE2, respectively. These indices of the br lateral branches' apical buds were 40.91% and 66.08% of WT, and 48.79 % and 67.66% of OE2, respectively (Figure 3a,b). The ratio of main branches SAM width/lateral branches SAM width of WT and OE2 were 1.51 and 1.29, respectively. However, that ratio of br was 0.95 (Figure 3c,d). As a result, the apical buds were smaller and the SAM size of the main branches and the lateral branches were similar in br.

Observation of Morphology and the Microstructures of Apical Buds
The morphology and microstructures of the main and lateral branches' apical buds in threeyear-old WT, OE2, and br were observed, as changes of growth traits in br might be related to the development of SAM. Either the main branches or the lateral branches' apical buds of br were found to be smaller, compared to the two control lines. The length and width of the br apical buds of the main branches were 62.08% and 74.43% of WT, and 67.93% and 68.04% of OE2, respectively. These indices of the br lateral branches' apical buds were 40.91% and 66.08% of WT, and 48.79 % and 67.66% of OE2, respectively (Figure 3a, b). The ratio of main branches SAM width/lateral branches SAM width of WT and OE2 were 1.51 and 1.29, respectively. However, that ratio of br was 0.95 (Figure 3c, d). As a result, the apical buds were smaller and the SAM size of the main branches and the lateral branches were similar in br.

Sequencing Quality Evaluation and qRT-PCR Verification
To analyze the changes of gene expression patterns, RNA-seq was performed using main and lateral branches' apical buds of three-year-old WT, OE2, and br. As a result, a total of 55.17 Gb clean data was obtained from 12 samples, and each sample reached 4.37 Gb. The Q30 percentages were all

Sequencing Quality Evaluation and qRT-PCR Verification
To analyze the changes of gene expression patterns, RNA-seq was performed using main and lateral branches' apical buds of three-year-old WT, OE2, and br. As a result, a total of 55.17 Gb clean data was obtained from 12 samples, and each sample reached 4.37 Gb. The Q30 percentages were all over 93.71%. The alignment of clean reads in the 12 samples, against the birch reference genome was from 93.59% to 94.62% (Table 2). Meanwhile, 1550 new genes were identified and annotated (Table 3). Samples-sequencing samples; clean reads-total pair-end reads in clean data; clean bases-total base number in clean data; GC content-percentage of G and C bases in clean bases; % ≥Q30-percentage of bases with clean data quality values greater than or equal to 30; total reads-number of clean reads, in single-end terms; mapped reads-the number of reads mapped to the reference genome and the percentage of clean reads. T1 and T2 represented the apical buds of the main branches of WT; T3 and T4 represented the lateral branches' apical buds of WT; T5 and T6 represented the apical buds of the main branches of OE2; T7 and T8 represented the lateral branches' apical buds of OE2; T9 and T10 represented the apical buds of the main branches of br; T11 and T12 represented the lateral branches' apical buds of br. To test the reliability of RNA-Seq, 12 genes were selected for verification, using qRT-PCR. The gene expression detected by qRT-PCR showed a similar pattern to that detected by the RNA-seq, which proved the accuracy of the transcriptome (Figure 4).

Transcriptomic Analyses of the Apical Buds of the Main Branches
Two transcriptome comparisons were performed, using the apical buds of the main branches, including WT vs. br and OE2 vs. br. The results showed that 859 genes (301 up-regulated and 558 down-regulated) were found to be differentially expressed in br, compared with WT. A total of 659 genes (262 up-regulated and 397 down-regulated) were differentially expressed in br, compared with OE2. Meanwhile, 406 genes were gathered in the overlap of DEGs between two transcriptome comparisons, including 149 up-regulated and 257 down-regulated genes ( Figure 5a).
DEGs which were expressed in the overlap of the two transcriptome comparisons, were assigned to the GO categories, according to the biological process, cellular component, and molecular function. In the biological process, DEGs were mainly concentrated in the cell process, the metabolic process, and single-organism process. For the cellular component, DEGs were mainly within the cell and cellular part. For molecular function, DEGs were mainly within the binding and catalytic activities (Figure 5b). There were 46 DEGs in biological regulation, 25 DEGs in development, 19 DEGs in signal transduction, and 4 DEGs in growth (Table S3).
Then, the DEGs were subjected to the KEGG pathway analysis and five pathways, including cellular process, environmental information processing, genetic information process, metabolism and organismal systems, were considered enriched (Figure 5c). Among them, one gene was involved in zeatin biosynthesis (ko00908, KEGG Orthology was K13495) ( Figure S2). Three genes were involved in diterpenoid biosynthesis (ko00904, KEGG Orthology were K04123, K04121, and K04123, respectively), which participated in the gibberellin biosynthesis ( Figure S3). Three genes were involved in the tryptophan metabolism (ko00380, KEGG Orthology were K11820, K11816, and K00128, respectively) and participated in the auxin biosynthesis ( Figure S4). Three genes were involved in the plant hormone signal transduction (ko04075, KEGG_Orthology were K14484, K14490, and K14505, respectively) including auxin response, cytokinin transduction, and cell division induced by brassinosteroids ( Figure S5). These results suggested that the specific changes of genes

Transcriptomic Analyses of the Apical Buds of the Main Branches
Two transcriptome comparisons were performed, using the apical buds of the main branches, including WT vs. br and OE2 vs. br. The results showed that 859 genes (301 up-regulated and 558 down-regulated) were found to be differentially expressed in br, compared with WT. A total of 659 genes (262 up-regulated and 397 down-regulated) were differentially expressed in br, compared with OE2. Meanwhile, 406 genes were gathered in the overlap of DEGs between two transcriptome comparisons, including 149 up-regulated and 257 down-regulated genes ( Figure 5a).
DEGs which were expressed in the overlap of the two transcriptome comparisons, were assigned to the GO categories, according to the biological process, cellular component, and molecular function. In the biological process, DEGs were mainly concentrated in the cell process, the metabolic process, and single-organism process. For the cellular component, DEGs were mainly within the cell and cellular part. For molecular function, DEGs were mainly within the binding and catalytic activities (Figure 5b). There were 46 DEGs in biological regulation, 25 DEGs in development, 19 DEGs in signal transduction, and 4 DEGs in growth (Table S3).
Then, the DEGs were subjected to the KEGG pathway analysis and five pathways, including cellular process, environmental information processing, genetic information process, metabolism and organismal systems, were considered enriched (Figure 5c). Among them, one gene was involved in zeatin biosynthesis (ko00908, KEGG Orthology was K13495) ( Figure S2). Three genes were involved in diterpenoid biosynthesis (ko00904, KEGG Orthology were K04123, K04121, and K04123, respectively), which participated in the gibberellin biosynthesis ( Figure S3). Three genes were involved in the tryptophan metabolism (ko00380, KEGG Orthology were K11820, K11816, and K00128, respectively) and participated in the auxin biosynthesis ( Figure S4). Three genes were involved in the plant hormone signal transduction (ko04075, KEGG_Orthology were K14484, K14490, and K14505, respectively) including auxin response, cytokinin transduction, and cell division induced by brassinosteroids  Figure S5). These results suggested that the specific changes of genes involved in the development, growth, plant hormones biosynthesis, and signal transduction in the apical buds of the main branches were important to the unique phenotype of br.

Transcriptomic Analyses of the Lateral Branches' Apical Buds
Two transcriptome comparisons were carried-out using the lateral branches' apical buds, including WT vs. br and OE2 vs. br. The results showed that 807 genes (241 up-regulated and 566 down-regulated) were found to be differentially expressed in the br lateral branches' apical buds, compared with the WT and 727 genes (284 up-regulated and 443 down-regulated) were differentially expressed, compared with the OE2. Furthermore, 396 DEGs (125 up-regulated and 271 downregulated) were gathered between the two sets ( Figure 6a).
Based on gene ontology, DEGs were classified into cell process, metabolic process, and singleorganism process in the biological process. For the cellular component, there were many DEGs involved in the cell and cell parts. For the molecular functions, DEGs were involved in binding and in the catalytic activity (Figure 6b). Compared with WT and OE2, 45 DEGs participated in the biological regulation, 18 DEGs participated in the development, 2 DEGs participated in the growth and 12 DEGs participated in the signal transduction in br (Table S4).
The KEGG pathway analysis revealed that three genes were enriched in diterpenoid biosynthesis, and three genes were enriched in the tryptophan metabolism, which were consistent

Transcriptomic Analyses of the Lateral Branches' Apical Buds
Two transcriptome comparisons were carried-out using the lateral branches' apical buds, including WT vs. br and OE2 vs. br. The results showed that 807 genes (241 up-regulated and 566 down-regulated) were found to be differentially expressed in the br lateral branches' apical buds, compared with the WT and 727 genes (284 up-regulated and 443 down-regulated) were differentially expressed, compared with the OE2. Furthermore, 396 DEGs (125 up-regulated and 271 down-regulated) were gathered between the two sets ( Figure 6a).
Based on gene ontology, DEGs were classified into cell process, metabolic process, and single-organism process in the biological process. For the cellular component, there were many DEGs involved in the cell and cell parts. For the molecular functions, DEGs were involved in binding and in the catalytic activity (Figure 6b). Compared with WT and OE2, 45 DEGs participated in the biological regulation, 18 DEGs participated in the development, 2 DEGs participated in the growth and 12 DEGs participated in the signal transduction in br (Table S4).
The KEGG pathway analysis revealed that three genes were enriched in diterpenoid biosynthesis, and three genes were enriched in the tryptophan metabolism, which were consistent with those found in the apical buds of the main branches ( Figures S3 and S4). Four genes were enriched in the plant hormone signal transduction pathway (KEGG_Orthology were K13449, K14484, K13463, and K14505, respectively.). These genes were involved in the auxin response, salicylic acid response, jasmonic acid response, and brassinolide-induced cell division (Figure 6c, Figure S5). These results suggested that the specific changes of genes related to development, growth, plant hormones biosynthesis, and signal transduction in the lateral branches' apical buds, also had effects on the unique phenotype of br. with those found in the apical buds of the main branches ( Figure S3, S4). Four genes were enriched in the plant hormone signal transduction pathway (KEGG_Orthology were K13449, K14484, K13463, and K14505, respectively.). These genes were involved in the auxin response, salicylic acid response, jasmonic acid response, and brassinolide-induced cell division (Figure 6c, Figure S5). These results suggested that the specific changes of genes related to development, growth, plant hormones biosynthesis, and signal transduction in the lateral branches' apical buds, also had effects on the unique phenotype of br.

Transcriptomic Analyses of the Apical Buds
Apical buds of the main branches transcriptome data and lateral branches' apical buds transcriptome data were combined for further analysis. The results showed that 1170 DEGs (403 upregulated and 767 down-regulated) and 991 DEGs (417 up-regulated and 574 down-regulated) were found in br, compared to the WT and the OE2, respectively. Additionally, 653 DEGs (253 up-regulated and 400 down-regulated) were common to two sets (Figure 7a).
Most of the DEGs identified were included in the cellular process, the metabolic process, and the single-organism process for biological processes based on gene ontology. There were many DEGs involved in the cell and cell components, while others were involved in the binding and catalytic activity for molecular function. Moreover, 58 DEGs related to biological regulation, 34 DEGs related to development, 4 DEGs related to growth and 23 DEGs related to signal transduction, expressed specifically in br, compared to WT and OE2 (Figure 7b, Table S5).
The DEGs were further subjected to the KEGG pathway analysis. The results showed that three genes were enriched in the diterpenoid biosynthesis ( Figure S3) and three genes were enriched in the tryptophan metabolism ( Figure S4). There were nine genes that participated in plant hormone signal transduction pathway (KEGG_Orthology were K14484, K14490, K13449, K14488, K13946, and K14505, respectively), which were involved in the auxin biosynthesis and response, salicylic acid

Transcriptomic Analyses of the Apical Buds
Apical buds of the main branches transcriptome data and lateral branches' apical buds transcriptome data were combined for further analysis. The results showed that 1170 DEGs (403 up-regulated and 767 down-regulated) and 991 DEGs (417 up-regulated and 574 down-regulated) were found in br, compared to the WT and the OE2, respectively. Additionally, 653 DEGs (253 up-regulated and 400 down-regulated) were common to two sets (Figure 7a).
Most of the DEGs identified were included in the cellular process, the metabolic process, and the single-organism process for biological processes based on gene ontology. There were many DEGs involved in the cell and cell components, while others were involved in the binding and catalytic activity for molecular function. Moreover, 58 DEGs related to biological regulation, 34 DEGs related to development, 4 DEGs related to growth and 23 DEGs related to signal transduction, expressed specifically in br, compared to WT and OE2 (Figure 7b, Table S5). The DEGs were further subjected to the KEGG pathway analysis. The results showed that three genes were enriched in the diterpenoid biosynthesis ( Figure S3) and three genes were enriched in the tryptophan metabolism ( Figure S4). There were nine genes that participated in plant hormone signal transduction pathway (KEGG_Orthology were K14484, K14490, K13449, K14488, K13946, and K14505, respectively), which were involved in the auxin biosynthesis and response, salicylic acid response, cytokinin transduction, and brassinolide-induced cell division (Figure 7c, Figure S5). These results suggested that there were many changes of genes involved in the development, growth plant hormones biosynthesis, and signal transduction, in both the apical buds of the main branches and the lateral branches' apical buds, and these changes were crucial to the variations of br. response, cytokinin transduction, and brassinolide-induced cell division (Figure 7c, FigureS5). These results suggested that there were many changes of genes involved in the development, growth plant hormones biosynthesis, and signal transduction, in both the apical buds of the main branches and the lateral branches' apical buds, and these changes were crucial to the variations of br.

DEGs Related to Flower Development
As seven-year-old br was different from the WT and OE2, it did not form infructescence ( Figure  S6). We hypothesized that br might have had flower development defects which caused fertility deficiency and delayed the formation of infructescence. Therefore, 25 DEGs related to flower development were selected to analyze the expression changes in the main and lateral branches' apical buds. It was found that gene 10077 (

DEGs Related to Flower Development
As seven-year-old br was different from the WT and OE2, it did not form infructescence ( Figure S6). We hypothesized that br might have had flower development defects which caused fertility deficiency and delayed the formation of infructescence. Therefore, 25 DEGs related to flower development were selected to analyze the expression changes in the main and lateral branches' apical buds.  Table S6). These results showed that changes of genes involved in floral development, flowering time, and plant fertility might have caused the fertility deficiency of br.  Table S6). These results showed that changes of genes involved in floral development, flowering time, and plant fertility might have caused the fertility deficiency of br.

DEGs Related to Cell Division, Differentiation, and SAM Activity
As br was dwarf, weak, slow-growing, multiple-branches, and its apical buds were smaller than the control lines, 39 genes related to cell division and differentiation, SAM activity, and organ formation were selected and analyzed. It was found that gene 15444, gene 17510, gene 22973, gene 27269, gene 28818, gene 29755, gene 5923, gene 701, and gene 7873 were up-regulated in br, all these genes were involved in the meristem initiation, vegetative to reproductive phase transition of meristem and negatively regulated growth, based on Go annotation. In addition, there were 30 genes which were involved in regulating the cell cycle, cell division and differentiation, meristem activity, and organ formation, which were significantly down-regulated. (Figure 9a, b, Table S7). These results showed that the changes of these genes were related to the dwarf, multiple-branches, and small buds of br.

DEGs Related to Cell Division, Differentiation, and SAM Activity
As br was dwarf, weak, slow-growing, multiple-branches, and its apical buds were smaller than the control lines, 39 genes related to cell division and differentiation, SAM activity, and organ formation were selected and analyzed. It was found that gene 15444, gene 17510, gene 22973, gene 27269, gene 28818, gene 29755, gene 5923, gene 701, and gene 7873 were up-regulated in br, all these genes were involved in the meristem initiation, vegetative to reproductive phase transition of meristem and negatively regulated growth, based on Go annotation. In addition, there were 30 genes which were involved in regulating the cell cycle, cell division and differentiation, meristem activity, and organ formation, which were significantly down-regulated. (Figure 9a,b, Table S7). These results showed that the changes of these genes were related to the dwarf, multiple-branches, and small buds of br.

DEGs Related to Plant Hormones
Based on the KEGG and GO analysis, 28 genes involved in auxin biosynthesis and transduction, cytokinin biosynthesis, gibberellin biosynthesis, and response, were selected for analyses in br. The genes related to auxin biosynthesis, including gene15480, gene16284, and gene21708, were downregulated. Similarly, the expression of seven auxin-responsive genes was down-regulated. While, the expression of gene4940 and gene701, which participated in auxin transduction and gene10815, Figure 9. Expression pattern of genes involved in cell division, differentiation and SAM activity in three-year-old WT, OE2, and br. (a) Thermal map of the genes involved in cell division, differentiation, and SAM activity, in apical buds of the main branches; and (b) thermal map of the mentioned genes in the lateral branches' apical buds. Analysis was performed using the normalization for 39 genes identified by the RNA-seq.

DEGs Related to Plant Hormones
Based on the KEGG and GO analysis, 28 genes involved in auxin biosynthesis and transduction, cytokinin biosynthesis, gibberellin biosynthesis, and response, were selected for analyses in br. The genes related to auxin biosynthesis, including gene15480, gene16284, and gene21708, were down-regulated. Similarly, the expression of seven auxin-responsive genes was down-regulated. While, the expression of gene4940 and gene701, which participated in auxin transduction and gene10815, gene19098, gene22059, and gene24144, which participated in auxin-responses, were up-regulated. The expression of gene11684, related to cytokinin biosynthesis, decreased in br. There were four genes involved in the gibberellin biosynthesis, in which gene18832 and gene9935 were down-regulated and gene14231 and gene19590 were up-regulated. The expression of genes participating in the gibberellin transduction reduced in br. (Figure 10a,b, Table S8). These results suggested that the changes of genes related to phytohormone of br led to abnormality in biosynthesis, distribution, and transduction of plant hormones. gene19098, gene22059, and gene24144, which participated in auxin-responses, were up-regulated. The expression of gene11684, related to cytokinin biosynthesis, decreased in br. There were four genes involved in the gibberellin biosynthesis, in which gene18832 and gene9935 were down-regulated and gene14231 and gene19590 were up-regulated. The expression of genes participating in the gibberellin transduction reduced in br. (Figure 10a, b, Table S8). These results suggested that the changes of genes related to phytohormone of br led to abnormality in biosynthesis, distribution, and transduction of plant hormones. (a) Thermal map of genes involved in biosynthesis, transduction, and response of plant hormones in the apical buds of the main branches; and (b) thermal map of the mentioned genes in the lateral branches' apical buds. Analysis was performed using normalization for the 28 genes identified by RNA-seq. Genes involved in auxin biosynthesis, transduction, and response, were from gene10815 to gene701; genes involved in cytokinin biosynthesis and response, were from gene11396 to gene23259; gibberellin-related genes were from gene121 to gene9935.

Discussion
RNA-seq is an important method that has been widely used to analyze gene expression profiles and functions in different plants. Mu   (a) Thermal map of genes involved in biosynthesis, transduction, and response of plant hormones in the apical buds of the main branches; and (b) thermal map of the mentioned genes in the lateral branches' apical buds. Analysis was performed using normalization for the 28 genes identified by RNA-seq. Genes involved in auxin biosynthesis, transduction, and response, were from gene10815 to gene701; genes involved in cytokinin biosynthesis and response, were from gene11396 to gene23259; gibberellin-related genes were from gene121 to gene9935.

Discussion
RNA-seq is an important method that has been widely used to analyze gene expression profiles and functions in different plants. Mu et al. and Lin et al. performed transcriptome analysis of B. pendula 'Dalecarlica' and B. pendula 'Purple Rain', to explore the mechanism of the leaf shape and leaf color, respectively. Their results laid a molecular foundation for the functional analysis of genes involved in plant growth and development [43,44]. In this study, we identified a multiple-branches mutant br from the BpCCR1 overexpression lines. Compared with WT and OE2, br exhibited phenotypes of dwarfing, multiple-branches, and small buds (Figure 2a, Figure 3). To investigate the molecular mechanism of the unique phenotype of br in birch, transcriptome sequencing of WT, OE2, and br were further performed. Genes involved in the SAM activity, cell division and differentiation, organ morphogenesis, plant hormone biosynthesis, and signal transduction were differentially expressed in br.
Most of the 39 genes involved in cell division and differentiation, SAM morphology, and activity were down-regulated in br ( Figure 9). For example, E2F (gene13790) encoded a cell proliferation protein controlled by ubiquitin-mediated degradation, which affected the cell division in SAM. Lateral organ boundaries (LOB, gene 2378) was specifically expressed at organ boundaries, and played an important role in maintaining SAM integrity and organ differentiation. Teosinte branched1/Cycloidea/proliferating cell nuclear antigen factor1 (TCP13, gene9677) promotes the mitosis of petal cells and affects the leaf morphological development [51][52][53][54]. In addition, genes of the cyclin B family (gene16589 and gene21913) and genes which affect the maintenance and initiation of meristem (gene11092 and gene8839) were also down-regulated in br. Our results indicated that the reduction of these genes might have resulted in defects of the SAM structure, the distinct plant architecture, and in the growth traits in br.
Genes involved in biosynthesis, transduction, and response of the plant hormones were also analyzed. Auxin is one of the important hormones which regulates many biological processes throughout the life cycle of plants, such as stem and root development, leaf morphological development, vascular bundle formation, and differentiation [55,56]. There are four major tryptophan-dependent auxin biosynthesis pathways, including indole-3-acetamide (IAM), indole-3-pyruvic acid (IPA), tryptamine (TAM), and indole-3-acetaldoxime (IAOX) pathways [57][58][59]. YUCCA proteins are crucial enzymes for auxin biosynthesis, and mutants of YUCCA can lead to the reduction of auxin content [60]. Auxin signal transduction begins with the specific binding between auxin and the auxin receptor, which then activates the AUX/IAA and promotes the expression of the response genes. In the br main and lateral branches' apical buds, YUCCA7 (gene16284) was down-regulated while auxin-response genes, such as IAA4 (gene19098), IAA12 (gene22059), IAA18 (gene10815), and SUR1-like (gene12158) were up-regulated. Our results indicated that the decrease of YUCCA7 might affect the auxin biosynthesis and reduce the IAA content in the apical buds of br. However, overexpression of the auxin-responsive genes suggests that there were no defects in signal transduction (Figure 11a).
Cytokinin controls plant cell division and regulates leaf senescence, chloroplast formation, and seed dormancy. In addition, cytokinin can also cooperate with other hormones to affect the lateral bud development, thereby, changing the plant architecture [61][62][63]. Cytokinin exists mainly in the form of isoprenoids in plants and there are two pathways (breakdown of tRNA and de novo synthesis) for zeatin biosynthesis [64,65]. Zeatin is classified into ubiquitous trans-zeatin and less active cis-zeatin. Cis-zeatin might work at a specific developmental stage of the plant, and it might be related to male sterility [66,67]. cisZOG1 encodes a cis-zeatin-specific O-glucosyltransferase which regulates the biosynthesis of cis-zeatin [68]. Our results indicated that the reduction of cisZOG1 (gene11684) in br has effects on cell division and differentiation in SAM, plant architecture, and fertility (Figure 11b).
Gibberellin is a diterpenoid plant hormone involved in the control of many important developmental processes, including plant growth, seed growth, anther growth, and flowering time [69][70][71]. It can also interact with auxin, cytokinin, abscisic acid, ethylene, and jasmonic acid to affect plant growth, development, and resistance [72][73][74][75][76][77][78]. The biosynthesis of gibberellin is divided into three steps. The first step is to form ent-kaurene by cyclization of the GGDP precursors twice. In the second step, ent-kaurene is oxidized into GA12. In the third step, GA12 forms different intermediates and Gas [70]. In br, gene18832 participating in ent-kaurene biosynthesis was down-regulated, while gene14231 and gene19590 encoding Kao were up-regulated. We inferred that the changes of these gene expression might have decreased the GA content and have had effects on the plant growth, floral development, and fertility in br (Figure 11c).
In summary, we found that there were many changes of gene expression patterns in the br main and lateral apical buds, compared with WT and OE2. These genes were involved in the process of cell division and differentiation, establishment of flowers and leaves, biosynthesis and signal transduction of plant hormones. However, we still do not know if the phenotypic mutation of br is caused by the changes of one gene or multiple genes, which need to be verified by subsequent experiments. (gene11684) in br has effects on cell division and differentiation in SAM, plant architecture, and fertility ( Figure 11b). Gibberellin is a diterpenoid plant hormone involved in the control of many important developmental processes, including plant growth, seed growth, anther growth, and flowering time [69][70][71]. It can also interact with auxin, cytokinin, abscisic acid, ethylene, and jasmonic acid to affect plant growth, development, and resistance [72][73][74][75][76][77][78]. The biosynthesis of gibberellin is divided into three steps. The first step is to form ent-kaurene by cyclization of the GGDP precursors twice. In the second step, ent-kaurene is oxidized into GA12. In the third step, GA12 forms different intermediates and Gas [70]. In br, gene18832 participating in ent-kaurene biosynthesis was down-regulated, while gene14231 and gene19590 encoding Kao were up-regulated. We inferred that the changes of these gene expression might have decreased the GA content and have had effects on the plant growth, floral development, and fertility in br (Figure 11c).
In summary, we found that there were many changes of gene expression patterns in the br main and lateral apical buds, compared with WT and OE2. These genes were involved in the process of cell division and differentiation, establishment of flowers and leaves, biosynthesis and signal transduction of plant hormones. However, we still do not know if the phenotypic mutation of br is caused by the changes of one gene or multiple genes, which need to be verified by subsequent experiments.

Conclusions
In this study, we identified br which exhibited dwarf, multiple-branches, small leaves, and apical buds. About 406 DEGs, 396 DEGs, and 653 DEGs were obtained in the br apical buds of the main branches, lateral branches' apical buds, and apical buds from the overlap of DEGs between WT vs. br and OE2 vs. br, respectively. These genes involved in the SAM activity, organogenesis, cell division and differentiation, plant hormone biosynthesis, and signal transduction, based on annotation, were important to the unique phenotype of br.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/10/5/374/s1, Figure S1: Growth performance of three-year-old WT, OE2, and br, Figure S2: Ko00908 pathway, in which the red marker gene (gene 11684) was contained in DEGs and significantly down-regulated, Figure S3: Ko00904 pathway, in which the red marker genes (gene14231, gene18832, and gene19590) were contained in DEGs, Figure S4: Ko00308 pathway, in which the red marker genes (gene15480, gene16284, and gene21708) were contained in DEGs and down-regulated, Figure S5: Plant hormone signal transduction, three red marker genes (gene 22059, gene 23259 and gene 7394) were contained in DEGs of main branches apical buds and up-regulated, four green line genes (gene12948, gene7394, gene22059, and gene24822) were found in DEGs of lateral branches' apical buds, nine genes (gene12948, gene7394, gene10815, gene19098, gene22059, gene4940, gene24119, gene11396, and gene3405) marked using blue lines were found in apical buds (including main and lateral branches apical buds), Figure S6: The infructescence of seven-year-old WT, OE2, and br, white arrows represent the infructescence of birch. Table S1: The primer sequences used in qRT-PCR, Table S2: Multiple comparison of the stomatal morphological index in the tested lines, Table S3: DEGs of apical buds of the main branches involved in biological regulation, development, growth, and signal transduction based on gene ontology, Table S4: DEGs of lateral branches' apical buds involved in biological regulation, development, growth, and signal transduction based on gene ontology, Table S5: DEGs of main and lateral branches' apical buds involved in biological regulation, development, growth, and signal transduction based on gene ontology, Table S6: The gene ID, sequences, and annotations of the 25 flowering-related genes, Table S7: The gene ID, sequences, and annotations of the 39 genes involved in cell division, differentiation and the SAM activity, Table S8: The gene ID, sequences, and annotations of the 28 genes involved in biosynthesis, transduction, and response of plant hormones.