Molecular mechanism of lateral bud differentiation of Pinus massoniana based on high-throughput sequencing

Knot-free timber cultivation is an important goal of forest breeding, and lateral shoots affect yield and stem shape of tree. The purpose of this study was to analyze the molecular mechanism of lateral bud development by removing the apical dominance of Pinus massoniana young seedlings through transcriptome sequencing and identify key genes involved in lateral bud development. We analyzed hormone contents and transcriptome data for removal of apical dominant of lateral buds as well as apical and lateral buds of normal development ones. Data were analyzed using an comprehensive approach of pathway- and gene-set enrichment analysis, Mapman visualization tool, and gene expression analysis. Our results showed that the contents of auxin (IAA), Zea and strigolactone (SL) in lateral buds significantly increased after removal of apical dominance, while abscisic acid (ABA) decreased. Gibberellin (GA) metabolism, cytokinin (CK), jasmonic acid, zeatin pathway-related genes positively regulated lateral bud development, ABA metabolism-related genes basically negatively regulated lateral bud differentiation, auxin, ethylene, SLs were positive and negative regulation, while only A small number of genes of SA and BRASSINOSTEROID, such as TGA and TCH4, were involved in lateral bud development. In addition, it was speculated that transcription factors such as WRKY, TCP, MYB, HSP, AuxIAA, and AP2 played important roles in the development of lateral buds. In summary, our results provided a better understanding of lateral bud differentiation and lateral shoot formation of P. massoniana from transcriptome level. It provided a basis for molecular characteristics of side branch formation of other timber forests, and contributed to knot-free breeding of forest trees.

www.nature.com/scientificreports/ In addition, recently researches found that some transcription factors were involved in development of collateral. For example, AP2/ERF, WRKY, MYB, and TCP transcript factors were found to play important roles in the differentiation and elongation of lateral shoot [21][22][23][24] .
Woody plants, especially afforestation timber tree species, have different research purposes on lateral shoot development. Crops and ornamental plants with multiple branches are necessary to get high yield and beautiful appearance. Whereas the timer tree species are expected to have few or zero lateral shoot in order to cultivate large diameter and knot-free timber, increase output and utilization of timber, reduce processing cost, and enhance the appearance of solid wood products. Cultivation and physiology on lateral bud development have been studied by many researches 25 , Lei and Sumida 26 studied the changes in branch structure of spruce under different light conditions. Researches on lateral bud development of polar in molecular side has been studied. The role of DWARF14, BRANCHED1 and other genes 27,28 in SL signaling regulation were determined. Wang et al. 29 preliminarily revealed the relationship between activity of lateral buds and hormones in hybrids of Populus euramericana and poplar using transcriptomics. At present, the research on the regulation mechanism of collateral development of woody plants has just begun, and the theoretical research on the mechanism of hormonal regulation of lateral bud differentiation will become the key direction.
P. massoniana is a typical coniferous native tree species in southern China with strong adaptability, fast growth rate, which provides 16% forest growing stock in southern China and rosin production accounts for more than half of the world's turpentine 30,31 . P. massoniana has more thick lateral shoot than Pinus elliottii, Pinus taeda and Larix potaninii, which is not easy to pruning and unfavorable for knot-free timber cultivation. Many researches about density regulation and pruning have been carried out, however, there have not been any reports on the selection of superior varieties of knot-free timber. Therefore, the analysis of development mechanism of lateral shoot is the best way to improve quality of P. massoniana wood.
In our study, we observe the process of lateral bud differentiation by microtechnique, and analyzed hormone changes during lateral bud development as well as transcriptomes of terminal bud and lateral bud after decapitation and in normal growth plants by high-throughput sequencing technology (control). We identified differentially expressed genes (DEGs) by reference and de novo assembly methods. Functional enrichment analysis of these DEGs indicates that lateral bud differentiation is affected by hormones and regulated by related genes, which may be the putative regulatory center of lateral bud differentiation. Our results provide important genetic resources and theoretical basis for differentiation mechanism of P. massoniana, and provide a theoretical basis for knot-free timber cultivation.

Result
Effect of removal of apex on plant morphology. The removal of apex results in plant height decrease from 55-60 cm to 35-40 cm (Fig. 1a). There were a great number of lateral bud outgrowth after 8-10 days, which were coming from a cluster of needle leaves (Fig. 1b). After removal of apical dominance, lateral buds grew rapidly with the length of 0.3 cm and 4 cm in12th and 55th day, respective.
Morphological feature of buds. The bud tip of P. massoniana included meristematic zone, elongation zone and maturarion zone (Fig. 2a). In ball-shape meristematic zone, the promeristem was composed of tunica cells and corpus cells. Tunica cells have longitudinal division to enlarge superficial area while irregular corpus cells undergo longitudinal and transverse divisions, amplifying volume of stem tip (Fig. 2b). Vascular bundle was composed of primary phloem, fascicular cambium and primary xylem, and meristematic cell in fascicular cambium had strong activity, and grew out of periderm along direction of phloem rays and became lateral buds (Fig. 2c,d). As can be seen from Fig. 2e and f, axillary buds originated from parenchyma cells of vascular cambium. The flat and rectangular-shape fusiform cells, distributed in ascicular cambium, were the main bodies of vascular cambium. Axillary buds have formatted during normal growth and development, but apical dominance has not been stimulated.

Contents of hormones in different parts.
The order of IAA content was apical buds, lateral buds and kryptoblasts. The content of IAA was significantly different between apical buds and kryptoblasts (P < 0.05). The content of IAA was also the highest in needle leaves around apical buds, but there were no significant differences compared to other materials (Fig. 2b). The highest and the lowest amount of ABA was in apical buds and lateral buds, respectively, and there were significant differences among other ones (P < 0.05). The needle leaves also showed the same phenomenon (Fig. 2b). Content of Zea was the highest in the lateral buds and leaves around lateral buds, and there is significant difference compared with other materials (Fig. 2b). The lowest content and the highest content of SLs were in shoot apex and leaves around shoot apex, axillary buds, respectively. Although SL contents in lateral buds were closed to those of axillary buds, there were significant differences (P ≤ 0.05. Transcriptome data assembling and annotation. We sequenced a total of 9 c-DNA libraries from 3 samples (in duplicates). We obtained a total of 232,736,473 high-quality reads containing 69,512,197,035 bases with GC contents of 45-48% (Supplementary Table S1).
All data were assembled by the Trinity software, resulting in 756, 676 Transcripts and 604, 122 unigenes with N50 of 801 and 395, respectively. The average length of unigene was more than 200 bp, and the length of 200-500 bp unigenes occupied 85.4% (Fig. S1, Table S2). 355, 314 unigenes were higher similar to known proteins from seven above-mentioned protein databases (Table S3). To further identify the interactions of all the annotated unigenes, we carried out GO functional enrichment and KEGG pathway analysis. A total of 174, 942 were assign to the three main GO categories (Fig. S2). Also, and total of 1219, 4910, 4628 from L v T, L v K and K v T were assign to KEGG pathways (Fig. S2 Table S5). Gene ontology analysis was performed by mapping each unigene to the GO database. 12.37% unigenes were annotated in Mapman and found to be associated with 34 metabolism pathways (Figs. 3c, Fig. S4, Table S6). 90% DEGs were found to be similar to known protein from Nr database to provide basis for metabolism analysis (Fig. S2, Table S7).
Differential expression of transcription factor-related genes. We used Mapman software to understand the involvement of transcription factors during the development of lateral buds after top treatment. The  DAL21 (MADS family) showed up-regulated expression, and all others were down-regulated. HSP family TFs all showed up-regulated expression, which showed positive regulation of lateral bud development. Only Glycine-rich RNA-binding protein 2 (GRP2, GRF family) showed negative regulation of lateral bud development. The AuxIAA family TFs was expressed as a positive regulation of lateral bud development, and ARF2 and ARF19 (Auxin response factor) (ARF family) appeared as negative regulatory genes. The ethylene-responsive transcription factor LEP (ERE-LEP) was a positive regulatory gene in AP2-EREBP family TFs, and ERE3, ERE23, ERE38, and ERE12 negatively regulated lateral bud development. Overall, our study showed that a large number of transcription factors were involved in development of lateral buds, some families were all positively regulated, some families had positive and negative regulatory genes, and key genes in several families TFs have been identified.
Differential expression of genes related to secondary metabolism. We used Mapman software to understand the involvement of transcription factors in the development of lateral buds after cutting, and the results showed that 701 genes involved in 19 secondary metabolic pathways including carotenoids, chalcones, lignin and lignans, MVA pathway, Non MVA pathway, and terpenoids were involved lateral bud development (Figs. 6, S6, Table S10). In L v D sample, we found that carotenoids was up regulation; while in L v K sample, except for cytochrome P450 family monooxygenase, other genes were down regulation. Some genes negatively regulated, including acetyl-CoA C-acetyltransferase, acetyl-CoA synthetase-like protein, alcohol dehydroge- www.nature.com/scientificreports/ nase, caffeic acid ortho-methyltransferase, cinnamoyl-CoA reductase, GroES-like protein, NADP-dependent alcohol dehydrogenase from lignin and lignans pathway, e acetyl-CoA C-acetyltransferas from MVA pathway, other ones showed positively regulated, such as hydroxycinnamoyl CoA shikimate, phenylalanoyl CoA ligase, geranylgeranyl transferase and 1-seoxy-d-xylulose 5-phosphate synthase genes from non MVA pathway, isopimaradiene synthase, limonene synthase, longifolene synthase and other genes from terpenoids metabolic pathway.

qRT-PCR.
To test the reliability of RNA-Seq data, we selected 53 DEGs to be examined by real-time quantitative (qRT-PCR). These genes included genes related to AUXIN resistant, Phytohormone auxin, carotenoid cleavage dioxygenases, More axillary growth, alpha/beta hydrolase, Sigaling repressor, Decreased apical dominance, which played key roles in lateral bud development, as well as genes from WRKY, TCP, MYB, RAV. The qRT-PCR results were consistent with RNA-Seq (Fig. 7), indicating that the RNA-Seq sequencing results were accurate and effective.

Discussion
At present, there are few researches on the growth, development and stress response mechanisms of P. massoniana 32,33 . Development of lateral shoot is an important factor affecting the growth, product and quality of P. massoniana, but studies on response mechanism and molecular response mechanism are rarely reported.
The sequencing technology provides us a direct insight into the genes response and allowi us to systematically, rapidly and comprehensively elucidate the molecular mechanisms for development of lateral shoot. The study used high-throughput sequencing technology to obtain a total of 756, 676 Transcript and 604, 122 Unigene. In addition, through the comparative analysis of materials, a large number of DEGs related to the lateral bud differentiation of P. massoniana were identified. The annotated DEGs covered plant growth and development, signal transduction, transcriptional genes, and secondary biomass biological processes. It provides basic for the molecular mechanism of lateral shoot differentiation of P. massoniana.   www.nature.com/scientificreports/ Figure 6. Differential expression of genes related to secondary metabolism. The pathway was generated using MapMAN software (Version 3. 6.0RC1). Totally 701 genes were involved in different pathway. Each square means a gene. Blue color means up-regulation and red color means down-regulation。Data was from differential expression of L v K. www.nature.com/scientificreports/ The growth regulation of lateral buds is complex, and plants need to respond to multiple exogenous and endogenous signals simultaneously 34 . The growth of lateral shoot is regulated by a variety of endogenous hormones, including auxin and SLs, which inhibit the growth of lateral buds, and CKs, which promote the growth of lateral buds. The studies of a large number of model plants aim to explore interaction and effects of three hormones in axillary buds. The interaction between these hormones forms a balanced network to regulate lateral shoot development 2,35 .
After decapitation of apex, the reduce of auxin synthesis results in decrease in auxin content that transported downward, thereby inducing lateral bud outgrowth. In our study, the content of IAA in K sample was high, and the IAA content in the L sample also increased after removal of apical dominance, which indicated that inhibition of the lateral bud was removed after topping 9 . SLs are considered to be new hormones that directly regulate plant branching 12,14 . Our study found that zea also played a role in lateral bud development, and Zea was a prerequisite for CK synthesis 36 .
The interaction between auxin and SLs was found to be the most important mechanism in plant lateral shootl differentiation. Auxin polar transport (PATS) and the second messenger model were accepted to explain regulation 37,38 .
The discovery of multiple mutants accelerated the study of auxin-SLs pathway-associated genes and overall regulatory networks 2 . The auxin signal mainly relies on PIN, PGP, AUX1/LAX for transport and conduction, and regulates the synthesis process of SLs by regulating AUXIN RESISTANCE PROTEIN 1-AUXIN SIGNALLING F-BOX PROTEIN (AXR1-AFB) gene. CAROTENOID CLEAVAGE DIOXYGENASE (carotenoid cleavage double oxygenase 7, CCD7), CCD8 (carotenoid cleavage double oxygenase 8), D27 (a ferritin), MAX1 (Cytochrome P450) were involved in the biosynthesis process of SLs. Carotenoids synthesized active SL 5-deaza-dual-gold alcohol, and then synthesized a variety of sole gold. Lactones by a series of biochemical processes to synthesize the 5,12,15 . Signaling transduction genes of SLs included D14 (encoding α-β folding hydrolase), MORE AXILLARY GROWTH 2 (MAX2, encoding SCF complex F-box protein, also known as Leu-rich repeat F-box protein), TCP transcription factor, etc.
These literatures verified our research data. In our study, a large number of signaling pathway genes were activated and exhibited different expression patterns, indicating that they had various functions in regulation of collateral differentiation. Among them, D alpha/beta hydrolase (DWARF) 39 , D14 40 , carotenoid cleavage dioxygenases (CCD) 24 , MORE AXILLARY GROWTH (MAX) 41 of SL signaling pathway have been shown to be importance in collateral differentiation.
In addition, we also found that some hormone signaling pathway genes except for auxin, SL, and ABA were involved in the lateral shooting of P. massoniana, such as GID from GA metabolism, ERF1 from ethylene pathway, TGA from SA and TCH4 frome Brassinosteroid pathway. Although IAA, GA and ABA were involved in collateral germination and regulation 42,43 , regulatory mechanism was still unclear.
Transcription factors were another important regulatory genes that regulated gene expression. A series of physiological, biochemical, metabolic and defense mechanisms have been formed, which played important roles in plant growth and development and in adapting to stresses. There were 75 WRKY transcription factors regulating plant growth, senescence, hormonal signaling pathways in Arabidopsis, and these genes responsed to biotic and abiotic stresses 44 . Although we have a comprehensive understanding of hormone involvement in lateral bud differentiation, little is known about the role of transcription factors in collateral differentiation 2 .
Studies have shown that WRKY71 regulates EXCESSIVE BRANCHES1 (EXB1) gene expression through auxin pathway and plays an important role in lateral bud differentiation, and also positively regulates lateral bud differentiation through H 2 O 2 , ABA signaling pathway in Arabidopsis. WRKY71 increased the number of branches through up-regulated expression of RAX2 and IAA efflux-related PINI2 and confirmed the direct interaction with RAX2 22,23,45 . Studies have confirmed that WRKY70 is an important member of salicylic acid and jasmonic acidregulated defense signals. WRKY70 was activated by salicylic acid and repressed by JA, and regulates SA and JA to participate in defense 26 . At the same time, many genes play negative regulatory roles in the defense process 46 , and WRKY70 negatively regulated growth of lateral roots and root hairs in root development 47 . Our study showed that not all genes in the WRKY transcription factor family are up-regulated, such as WRKY50, WRKY51, and confirmed that WRKY1 gene regulated insect-resistant defense of P. massoniana (data not shown) through jasmonic acid, GA and other hormonal signaling pathways. We hypothesized that WRKY transcription factors may also participate directly or indirectly in lateral bud development through one or several hormonal signaling pathways.
The BRANCHED1 (BRC1) from TCP family played an important role in Arabidopsis lateral shoot differentiation 48 , and PsBRC1 regulates collateral differentiation of pears via the SL and CK pathways 49 . Our study showed that the TCP transcription factor family gene showed up-regulated expression during lateral bud differentiation, indicating that TCP promotes collateral development, but whether it is regulated by SL or CK pathway remains to be further studied.
The tomatoTrifoliate (Tf) gene from tomato MYB transcription factor family plays an important role in the development of axillary buds, which was similar to LATERAL ORGAN FUSION1 (LOF1) and LOF2 genes in Arabidopsis 50,51 . ERF BUD ENHANCER (EBE) gene from AP2/ERF transcription factor directly regulates the differentiation of Arabidopsis collaterals 52 , CmERF053 plays an important role in lateral shoot differentiation of chrysanthemum 21 . Loss of genes such as TB1, TCP18, TCP12, and BRC2 leads to an increase in collaterals 2,53 . We initially identified 55 transcription factors involved in the lateral bud differentiation of P. massoniana, the further verification of synergistic reaction in the lateral bud differentiation of P. massoniana is needed. After breeding under outdoor condition for five months, seedlings were placed in an artificial climate chamber for cultivation (light/dark: 8/12 h, light intensity 1200 lx, air humidity 75%) with consistent management as previous one. After 15 days, we excised the shoot tip from shoot axis of six-m-old seedlings to form lateral buds (L), while terminal buds (T) from normal seedlings and kryptoblasts (K) on stems were chosen as control groups. Lateral bud samples were removed when the length was 0.05-0.1 cm (Fig. 1). Three biological replicates were used for each sample, and 10 seedlings were treated per sample. Samples for transcript analysis and fluorescence quantitative expression analysis were stored in liquid nitrogen. Microscopic observation samples were stored in FAA fixative (38% formaldehyde 5 mL, glacial acetic acid 5 mL, 70% alcohol 90 mL). Fresh samples were used for hormone analysis.

Material and method
Micro-observation. The fresh tissues were fixed in Formalin fix liquid-acetic acid-ethanol (FAA) for more than 24 h. The tissues were soften by acid solution, dehydrated by graded alcohol series with dehydrator, soaked wax, and placed them in embedding machine. After trimming paraffin blocks into appropriate size, we cut 4 ~ 10 μm sections from each sample. The paraffin ribbons were placed in water bath at 50 °C, flatterned and mounted onto slides, incubated in 60 °C and stored at room temperature. SafraninO-fast green was used for staining. Finally, examination was carried out by microscope (LEICA, DM3000, Germany), image acquisition was performed by software Leica Application Suite V4 and analyzed by image analysis system DMC4500+. All of the tools are from LEICA company, Germany. De novo assembly and functional annotation. The clean data were assembled into unigenes using Trinity software 54 . BlastX alignment between unigenes and NR、Swiss-Prot, GO, COG, KOG, KEGG were performed. Amino acids of unigenes were aligned to the Pfam database 55 using HMMER 56 to obtain annotation. Coding sequences of unigenes were predicted by TransDecoder.

Measurement
Differential expression analysis. Aligement between Reads and unigenes was performed using Bowtie 57 and gene expression levels were estimated by RSEM. Unigene expression abundance was calculated using the FPKM method. Pearson's Correlation Coefficient (r) was used to evaluate the index of biological repeat correlation 58 . Differential expression analysis between sample groups was performed using DESeq 59 . The p-values were adjusted by Benjamini and Hochberg's approach to minimize false discovery rate. The differentially expressed genes (DEGs) were identified according to criteria: an adjusted p-value ≤ 0.05, a FDR < 0.01, a Fold Change ≥ 2.
Functional annotation and GO enrichment analysis. We imported Nr results into Blast2GO program and get all the selected genes with annotated GO terms. WEGO software was used to classify the GO terms. GO enrichment was performed by topGO. We assigned the assembled sequences by KEGG pathway.
Quantitative real-time polymerase chain reaction (qRT-PCR). Primers were designed for 52 isoform sequences (Table S11). actin1 gene was used as an internal control. ACT1 was amplified through the following primers as standard control: 5′ -CAG TGT CTG GAT TGG AGG TTC-3′ as primer 1 and 5′ -TCT GTG GAC GAT GGA AGG AC-3′ as primer 2 60 . Reserve transcription of mRNA was done according to manufacturer's instructions of M-MLV. qRT-PCR reactions were performered using SYBR Premix Ex Taq II (TaKaRa) on a Light Cycle 480 II system. The relative expression level was calculated by means of 2 -ΔΔCt method 61 .