Genome-Wide Identification of the TCP Gene Family in Broussonetia papyrifera and Functional Analysis of BpTCP8, 14 and 19 in Shoot Branching

The plant-specific TCP family proteins play an important role in the processes of plant growth and development. Broussonetia papyrifera is a versatile perennial deciduous tree, and its genome data have been published. However, no comprehensive analysis of the TCP gene family in B. papyrifera has been undertaken. In this study, 20 BpTCP genes (BpTCPs) were identified in the B. papyrifera genome. Phylogenetic analysis divided BpTCPs into three subclades, the PCF subclade, the CIN subclade and the CYC/TB1 subclade. Gene structure analysis displayed that all BpTCPs except BpTCP19 contained one coding region. Conserved motif analysis showed that BpTCP proteins in the same subclade possessed similar motif structures. Segmental duplication was the primary driving force for the expansion of BpTCPs. Expression patterns showed that BpTCPs may play diverse biological functions in organ or tissue development. Transcriptional activation activity analysis of BpTCP8, BpTCP14 and BpTCP19 showed that they possessed transcriptional activation ability. The ectopic expression analysis in Arabidopsis wild-type and AtBRC1 ortholog mutant showed that BpTCP8, BpTCP14 and BpTCP19 could prevent rosette branch outgrowth. Collectively, our study not only established the first genome-wide analysis of the B. papyrifera TCP gene family, but also provided valuable information for understanding the function of BpTCPs in shoot branching.


Introduction
The TCP family is a plant-specific transcription factor family, which was first found in 1999 [1]. The name of TCP transcription factor family was derived from the four proteins originally discovered, TEOSINTE BRANCHED 1 (TB1) from Zea mays, CYCLOIDEA (CYC) from Antirrhinum majus and the PROLIFERATING CELL FACTORS 1 and 2 (PCF 1 and PCF2) from Oryza sativa [2][3][4]. The protein sequences of TCP transcription factors contain a conserved non-canonical basic helix-loop-helix (bHLH) motif of about 59 amino acids, which was called the TCP domain. According to the differences of TCP domains, the TCP family members can be divided into two subfamilies: class I (composed of the PCF subclade) and class II (composed of the CIN and CYC/TB1 subclades) [5,6]. The most noteworthy difference between these two subfamilies is that the basic region of TCP domain of class II family has four amino acids more than that of class I family. In addition, several members of class II have another conserved region outside the TCP domain named the R domain, which is an arginine-rich motif containing eighteen to twenty residues [6]. The R domain is hypothesized to form a coiled coil that may be involved in protein-protein interactions [7]. and expression patterns in different tissues. Meanwhile, the functions of BpTCP8, BpTCP14 and BpTCP19 in regulating shoot branching were validated by ectopic expression in Arabidopsis wild-type and AtBRC1 ortholog mutant.

Identification and Property Analysis of the TCP Genes in B. papyrifera
To identify the TCP genes in B. papyrifera, local BLASTp and HMM searches were performed against the B. papyrifera genome database using 24 known Arabidopsis TCP protein sequences. Subsequently, the non-redundant candidate BpTCP protein sequences of the two methods were submitted to the NCBI CD Search and InterPro database for further verification. Finally, a total of 20 BpTCP proteins were identified, which contained the conserved TCP domain. According to the B. papyrifera genomic information, a chromosome location map was constructed to illustrate the distribution of the BpTCP genes on each chromosome. The map showed that all BpTCP genes were unevenly mapped onto 10 out of 13 B. papyrifera chromosomes and the BpTCP genes were annotated as BpTCP1 to BpTCP20 in the light of their physical locations on the chromosomes ( Figure S1).
In order to understand the BpTCP proteins more comprehensively, the physical and chemical properties were analyzed ( Table 1). The length of BpTCP proteins varied from 165 (BpTCP11) to 566 (BpTCP12) amino acid residues. The molecular weight (MW) ranged from 18,140.32 (BpTCP11) to 59,824.02 (BpTCP12) Da. The protein isoelectric point (pI) distributed from 5.26 (BpTCP9) to 9.72 (BpTCP20) with a mean of 7.501. All the BpTCP proteins were predicted to localize in the nucleus. The prediction of the phosphorylation site displayed that the BpTCP proteins contained 1 (BpTCP15) to 9 (BpTCP8) phosphorylation sites, except BpTCP9 and BpTCP11 that had no phosphorylation site. Among them, 17 BpTCP proteins contained Ser site with the number ranging from 1 to 8, and 11 BpTCP proteins contained Thr site with the number ranging from 1 to 3. Only BpTCP15 and BpTCP18 contained Tyr site with the number of 1 and 3, respectively.

Phylogenetic Relationship of BpTCP Proteins and Distribution of TCP Proteins in Plants
To comprehend the accurate classification of BpTCP proteins and the phylogenetic relationship with other known TCP proteins, we used the full length sequences of 65 TCP proteins from B. papyrifera, Arabidopsis, O. sativa and 3 TCP proteins (AmCIN, AmCYC and ZmTB1) with known functions to establish an unrooted phylogenetic tree through MEGA X software ( Figure 1). The results showed Plants 2020, 9, 1301 4 of 19 that the BpTCP proteins were obviously classified into two subfamilies: class I (namely PCF subclade) and class II (contained CIN and CYC/TB1 subclade). Moreover, the TCP proteins of Arabidopsis and O. sativa exhibited the same classification as previous studies [6,47], which verified the reliability of the phylogenetic tree. In accordance with the classification, there were 11 BpTCP proteins in the PCF subclade, 6 BpTCP proteins in the CIN subclade and 3 BpTCP proteins in the CYC/TB1 subclade.

Phylogenetic Relationship of BpTCP Proteins and Distribution of TCP Proteins in Plants
To comprehend the accurate classification of BpTCP proteins and the phylogenetic relationship with other known TCP proteins, we used the full length sequences of 65 TCP proteins from B. papyrifera, Arabidopsis, O. sativa and 3 TCP proteins (AmCIN, AmCYC and ZmTB1) with known functions to establish an unrooted phylogenetic tree through MEGA X software ( Figure 1). The results showed that the BpTCP proteins were obviously classified into two subfamilies: class I (namely PCF subclade) and class II (contained CIN and CYC/TB1 subclade). Moreover, the TCP proteins of Arabidopsis and O. sativa exhibited the same classification as previous studies [6,47], which verified the reliability of the phylogenetic tree. In accordance with the classification, there were 11 BpTCP proteins in the PCF subclade, 6 BpTCP proteins in the CIN subclade and 3 BpTCP proteins in the CYC/TB1 subclade. Meanwhile, in order to understand the evolutionary relationship of TCP family among different plants, we analyzed the composition of TCP family on 30 species (including angiosperms, Meanwhile, in order to understand the evolutionary relationship of TCP family among different plants, we analyzed the composition of TCP family on 30 species (including angiosperms, ferns, mosses and green algae) and the number of TCP proteins of each subclade in every species was annotated ( Figure 2, Table S1). As the figure shows, there are no TCP proteins in Volvox carteri, Dunaliella salina and Chlamydomonas reinhardtii. Only a few TCP proteins exist in Selaginella moellendorffii (6), Physcomitrella patens (7), Marchantia polymorpha (2), and they do not contain CYC/TB1 subclade. Besides, the basal angiosperm Amborella trichopoda, which has radially symmetrical flowers, contains 15 TCP proteins and also lacks CYC/TB1 subclade. Both monocots and eudicots have the CYC/TB1 subclade, the number of CYC/TB1 subclade members varies from 1 to 17 and the number of the whole TCP family members varies greatly among species. The range of TCP proteins in monocots is 17 to 21, with exceptions of 44 for Z. mays and 45 for Musa acuminate, and the number of CYC/TB1 subclade members is 17 in Z. mays and 3 in M. acuminate, which is a significant difference. The number of TCP proteins in eudicots ranges from 15 (Vitis vinifera) to 58 (Malus domestica), among which the corresponding number of CYC/TB1 subclade members is 3 and 4. B. papyrifera and Morus notabilis belong to Moraceae and the number of TCP proteins is similar, which is 20 and 22, respectively. However, the number of CYC/TB1 subclade is 6 in M. notabilis and 3 in B. papyrifera. In general, the number of TCP members in different species varies substantially. The CYC/TB1 subclade appeared after the emergence of radially symmetrical flowers and the proportion of CYC/TB1 subclade members in the whole TCP family was not linearly related to the total number. ferns, mosses and green algae) and the number of TCP proteins of each subclade in every species was annotated ( Figure 2, Table S1). As the figure shows, there are no TCP proteins in Volvox carteri, Dunaliella salina and Chlamydomonas reinhardtii. Only a few TCP proteins exist in Selaginella moellendorffii (6), Physcomitrella patens (7), Marchantia polymorpha (2), and they do not contain CYC/TB1 subclade. Besides, the basal angiosperm Amborella trichopoda, which has radially symmetrical flowers, contains 15 TCP proteins and also lacks CYC/TB1 subclade. Both monocots and eudicots have the CYC/TB1 subclade, the number of CYC/TB1 subclade members varies from 1 to 17 and the number of the whole TCP family members varies greatly among species. The range of TCP proteins in monocots is 17 to 21, with exceptions of 44 for Z. mays and 45 for Musa acuminate, and the number of CYC/TB1 subclade members is 17 in Z. mays and 3 in M. acuminate, which is a significant difference. The number of TCP proteins in eudicots ranges from 15 (Vitis vinifera) to 58 (Malus domestica), among which the corresponding number of CYC/TB1 subclade members is 3 and 4. B. papyrifera and Morus notabilis belong to Moraceae and the number of TCP proteins is similar, which is 20 and 22, respectively. However, the number of CYC/TB1 subclade is 6 in M. notabilis and 3 in B. papyrifera. In general, the number of TCP members in different species varies substantially. The CYC/TB1 subclade appeared after the emergence of radially symmetrical flowers and the proportion of CYC/TB1 subclade members in the whole TCP family was not linearly related to the total number.

Analysis of Gene Structure and Conserved Motifs
To get a better insight into the evolutionary relationships and explore the structure diversification of the BpTCP family in each subclade, we constructed a new phylogentic tree of

Analysis of Gene Structure and Conserved Motifs
To get a better insight into the evolutionary relationships and explore the structure diversification of the BpTCP family in each subclade, we constructed a new phylogentic tree of BpTCP proteins, and analyzed the gene structure and conserved motifs of BpTCP family. The phylogenetic tree of the BpTCP family was constructed with the conserved TCP domain sequences, and the result showed that it was also divided into three subclades ( Figure 3A). The gene structures of BpTCP genes were remarkably similar ( Figure 3B). There were 19 BpTCPs with one coding region and only BpTCP19 that belongs to the CYC/TB1 subclade contained one intron in the coding region. Conserved motifs occupy an important role in the characteristic analysis and classification of gene family. Using the Plants 2020, 9, 1301 6 of 19 MEME online program, 10 motifs were identified among the BpTCP proteins ( Figure 3C and Figure S2). As shown, BpTCP proteins in the same subfamilies possess similar motif composition. Motif 1 (TCP domain) exists in all BpTCP members. Motif 2 exists in all class I (PCF subclade) members (except BpTCP11). Motif 5 exists in all class II (CIN and CYC/TB1 subclade) members. In addition, motif 3 is distributed in 9 out of 20 BpTCP proteins, which are scattered throughout each subclade. Motif 4 and motif 10 are exclusively present in part of the PCF subclade members. Motif 6 is exclusively present in 2 CIN subclade members and all CYC/TB1 subclade members. Motif 7 exists in several PCF and CYC/TB1 subclade members. Motif 8 and motif 9 only exist in some PCF and CIN subclade members. These results suggest that the motifs exclusively existing in a certain subclade may be relative to the specific function.
BpTCP proteins, and analyzed the gene structure and conserved motifs of BpTCP family. The phylogenetic tree of the BpTCP family was constructed with the conserved TCP domain sequences, and the result showed that it was also divided into three subclades ( Figure 3A). The gene structures of BpTCP genes were remarkably similar ( Figure 3B). There were 19 BpTCPs with one coding region and only BpTCP19 that belongs to the CYC/TB1 subclade contained one intron in the coding region. Conserved motifs occupy an important role in the characteristic analysis and classification of gene family. Using the MEME online program, 10 motifs were identified among the BpTCP proteins ( Figures 3C and S2). As shown, BpTCP proteins in the same subfamilies possess similar motif composition. Motif 1 (TCP domain) exists in all BpTCP members. Motif 2 exists in all class I (PCF subclade) members (except BpTCP11). Motif 5 exists in all class II (CIN and CYC/TB1 subclade) members. In addition, motif 3 is distributed in 9 out of 20 BpTCP proteins, which are scattered throughout each subclade. Motif 4 and motif 10 are exclusively present in part of the PCF subclade members. Motif 6 is exclusively present in 2 CIN subclade members and all CYC/TB1 subclade members. Motif 7 exists in several PCF and CYC/TB1 subclade members. Motif 8 and motif 9 only exist in some PCF and CIN subclade members. These results suggest that the motifs exclusively existing in a certain subclade may be relative to the specific function.

Chromosome Location and Duplication Event Analysis
Gene duplications are considered to play a vital part in the expansion and evolution of the gene family, so the duplication event analysis of BpTCP genes was performed ( Figure 4). The results showed that there was no tandem duplication events observed, which suggested that segmental duplications may be the primary driving force for the expansion of BpTCP family. A total of four pairs of paralogous BpTCP genes were identified and they distributed on five chromosomes. Among these genes, four genes (BpTCP3, BpTCP5, BpTCP6 and BpTCP16) belong to the CIN subclade and three genes (BpTCP8, BpTCP14 and BpTCP19) belong to the CYC/TB1 subclade. The four segmental duplication events were BpTCP3/BpTCP5, BpTCP6/BpTCP16, BpTCP8/BpTCP14 and BpTCP8/BpTCP19, respectively.

Chromosome Location and Duplication Event Analysis
Gene duplications are considered to play a vital part in the expansion and evolution of the gene family, so the duplication event analysis of BpTCP genes was performed ( Figure 4). The results showed that there was no tandem duplication events observed, which suggested that segmental duplications may be the primary driving force for the expansion of BpTCP family. A total of four pairs of paralogous BpTCP genes were identified and they distributed on five chromosomes. Among these genes, four genes (BpTCP3, BpTCP5, BpTCP6 and BpTCP16) belong to the CIN subclade and three genes (BpTCP8, BpTCP14 and BpTCP19) belong to the CYC/TB1 subclade. The four segmental duplication events were BpTCP3/BpTCP5, BpTCP6/BpTCP16, BpTCP8/BpTCP14 and BpTCP8/BpTCP19, respectively.

Expression Pattern Analysis of BpTCP Genes
To provide more insight into the potential roles of BpTCP genes during growth and development, the expression patterns of BpTCP genes were investigated based on the RNA-Seq data of nine tissues (shoot apex, young leaf, developing leaf, mature leaf, immature stem, phloem of proximal stem, phloem of mature stem, phloem of root and root tip) ( Figure 5). As shown in the heat map, genes belonging to the same subclade did not always share similar expression patterns. Some genes possessed similar expression patterns in different tissues, while some genes exhibited significant tissue specificity. For example, among PCF subclade genes, BpTCP5, BpTCP6, BpTCP7, BpTCP12, BpTCP16, BpTCP17 and BpTCP20 were expressed in almost all tissues, whereas BpTCP11 was hardly expressed in all tissues. BpTCP2 and BpTCP4 were lowly expressed in all tissues. The expression level of BpTCP3 in young leaf was higher than other tissues. In the CIN subclade, BpTCP1, BpTCP10, BpTCP13 and BpTCP15 showed the similar expression patterns, which were highly expressed in shoot apex, young leaf, developing leaf and mature leaf. BpTCP9 and BpTCP18 were expressed at low level in all tissues. In the CYC/TB1 subclade, the expression levels of BpTCP14 and BpTCP19 were a little higher in shoot apex and immature stem than other seven tissues. BpTCP8 was slightly expressed in the shoot apex and extremely lowly expressed in other eight tissues.

Expression Pattern Analysis of BpTCP Genes
To provide more insight into the potential roles of BpTCP genes during growth and development, the expression patterns of BpTCP genes were investigated based on the RNA-Seq data of nine tissues (shoot apex, young leaf, developing leaf, mature leaf, immature stem, phloem of proximal stem, phloem of mature stem, phloem of root and root tip) ( Figure 5). As shown in the heat map, genes belonging to the same subclade did not always share similar expression patterns. Some genes possessed similar expression patterns in different tissues, while some genes exhibited significant tissue specificity. For example, among PCF subclade genes, BpTCP5, BpTCP6, BpTCP7, BpTCP12, BpTCP16, BpTCP17 and BpTCP20 were expressed in almost all tissues, whereas BpTCP11 was hardly expressed in all tissues. BpTCP2 and BpTCP4 were lowly expressed in all tissues. The expression level of BpTCP3 in young leaf was higher than other tissues. In the CIN subclade, BpTCP1, BpTCP10, BpTCP13 and BpTCP15 showed the similar expression patterns, which were highly expressed in shoot apex, young leaf, developing leaf and mature leaf. BpTCP9 and BpTCP18 were expressed at low level in all tissues. In the CYC/TB1 subclade, the expression levels of BpTCP14 and BpTCP19 were a little higher in shoot apex and immature stem than other seven tissues. BpTCP8 was slightly expressed in the shoot apex and extremely lowly expressed in other eight tissues.

Quantitative RT-PCR of BpTCP8, BpTCP14 and BpTCP19
AtBRC1 (AtTCP18) and AtBRC2 (AtTCP12), which belong to the CYC/TB1 subclade, were demonstrated to regulate the shoot branching by arresting bud development [11]. In addition, plant architecture is of great significance for the application of B. papyrifera. Therefore, we turned our focus to the CYC/TB1 subclade of the BpTCP family. Quantitative RT-PCR was carried out to compare expression level of BpTCP8, BpTCP14 and BpTCP19 in shoot apex, leaf axil, axillary bud, leaf, stem and root ( Figure 6). The results showed that BpTCP8 and BpTCP19 exhibited an extremely prominent expression level in leaf axil and axillary bud. BpTCP14 exhibited higher expression level in axillary bud, but much lower than BpTCP8 and BpTCP19.

Quantitative RT-PCR of BpTCP8, BpTCP14 and BpTCP19
AtBRC1 (AtTCP18) and AtBRC2 (AtTCP12), which belong to the CYC/TB1 subclade, were demonstrated to regulate the shoot branching by arresting bud development [11]. In addition, plant architecture is of great significance for the application of B. papyrifera. Therefore, we turned our focus to the CYC/TB1 subclade of the BpTCP family. Quantitative RT-PCR was carried out to compare expression level of BpTCP8, BpTCP14 and BpTCP19 in shoot apex, leaf axil, axillary bud, leaf, stem and root ( Figure 6). The results showed that BpTCP8 and BpTCP19 exhibited an extremely prominent expression level in leaf axil and axillary bud. BpTCP14 exhibited higher expression level in axillary bud, but much lower than BpTCP8 and BpTCP19.

Quantitative RT-PCR of BpTCP8, BpTCP14 and BpTCP19
AtBRC1 (AtTCP18) and AtBRC2 (AtTCP12), which belong to the CYC/TB1 subclade, were demonstrated to regulate the shoot branching by arresting bud development [11]. In addition, plant architecture is of great significance for the application of B. papyrifera. Therefore, we turned our focus to the CYC/TB1 subclade of the BpTCP family. Quantitative RT-PCR was carried out to compare expression level of BpTCP8, BpTCP14 and BpTCP19 in shoot apex, leaf axil, axillary bud, leaf, stem and root ( Figure 6). The results showed that BpTCP8 and BpTCP19 exhibited an extremely prominent expression level in leaf axil and axillary bud. BpTCP14 exhibited higher expression level in axillary bud, but much lower than BpTCP8 and BpTCP19.

Transcriptional Activation Activity of BpTCP8, BpTCP14 and BpTCP19
The transcriptional activation activity of BpTCP8, BpTCP14 and BpTCP19 was tested by yeast one-hybrid assay. The degree of transcriptional activation activity can be measured by observing the ability of transformed yeast cells growing on the selective medium with 0-50 mM HIS3 protein competitive inhibitor 3-aminotriazole (3-AT). All transformants of BpTCP8, BpTCP14 and BpTCP19 could readily grow on the SD-Trp medium (Figure 7). In the SD-Trp-His medium without 3-AT, only the yeast cells containing BpTCPs could grow well. In the SD-Trp-His medium supplemented with different concentrations of 3-AT, the yeast cells containing BpTCPs exhibited different growth capacity. BpTCP19 exhibited the strongest activation activity, followed by BpTCP8 and finally BpTCP14.

Transcriptional Activation Activity of BpTCP8, BpTCP14 and BpTCP19
The transcriptional activation activity of BpTCP8, BpTCP14 and BpTCP19 was tested by yeast one-hybrid assay. The degree of transcriptional activation activity can be measured by observing the ability of transformed yeast cells growing on the selective medium with 0-50 mM HIS3 protein competitive inhibitor 3-aminotriazole (3-AT). All transformants of BpTCP8, BpTCP14 and BpTCP19 could readily grow on the SD-Trp medium (Figure 7). In the SD-Trp-His medium without 3-AT, only the yeast cells containing BpTCPs could grow well. In the SD-Trp-His medium supplemented with different concentrations of 3-AT, the yeast cells containing BpTCPs exhibited different growth capacity. BpTCP19 exhibited the strongest activation activity, followed by BpTCP8 and finally BpTCP14.

Roles of BpTCP8, BpTCP14 and BpTCP19 in Shoot Branching
To investigate whether the functions of BpTCP8, BpTCP14 and BpTCP19 are similar with AtBRC1 in regulating shoot branching, 35s::BpTCPs were transferred into Arabidopsis wild-type and AtBRC1 ortholog mutant brc1 by the way of ectopic expression. Two independent transgenic lines of each BpTCP gene were selected for phenotypic analysis (Figures 8A and 9A). The expression levels of BpTCP8, BpTCP14 and BpTCP19 were detected by RT-PCR and all transgenic lines exhibited elevated expression levels ( Figures 8B and 9B). Under the same growth conditions, the brc1 mutants exhibited an apparently luxuriant rosette branches than wild-type plants ( Figure 8A). The ectopic expression of BpTCP8, BpTCP14 or BpTCP19 in brc1 could reduce the number of primary rosette branches, which was a significant difference from brc1. Furthermore, BpTCP19 was sufficient to restore the number of rosette branches of brc1 to the wild-type ( Figure 8C). In wild-type plants, the ectopic expression of BpTCP19 could also lead to a decreasing number of primary rosette branches, which was apparently distinguishable from the wild-type, but BpTCP8 and BpTCP14 could not ( Figure 9C). On the other hand, all of the transgenic plants showed a similar number of primary cauline branches as the brc1 or the wild-type (Figures 8D and 9D). These results imply that BpTCP8, BpTCP14 and BpTCP19 prevent primary rosette branch outgrowth.

Roles of BpTCP8, BpTCP14 and BpTCP19 in Shoot Branching
To investigate whether the functions of BpTCP8, BpTCP14 and BpTCP19 are similar with AtBRC1 in regulating shoot branching, 35s::BpTCPs were transferred into Arabidopsis wild-type and AtBRC1 ortholog mutant brc1 by the way of ectopic expression. Two independent transgenic lines of each BpTCP gene were selected for phenotypic analysis (Figures 8A and 9A). The expression levels of BpTCP8, BpTCP14 and BpTCP19 were detected by RT-PCR and all transgenic lines exhibited elevated expression levels ( Figures 8B and 9B). Under the same growth conditions, the brc1 mutants exhibited an apparently luxuriant rosette branches than wild-type plants ( Figure 8A). The ectopic expression of BpTCP8, BpTCP14 or BpTCP19 in brc1 could reduce the number of primary rosette branches, which was a significant difference from brc1. Furthermore, BpTCP19 was sufficient to restore the number of rosette branches of brc1 to the wild-type ( Figure 8C). In wild-type plants, the ectopic expression of BpTCP19 could also lead to a decreasing number of primary rosette branches, which was apparently distinguishable from the wild-type, but BpTCP8 and BpTCP14 could not ( Figure 9C). On the other hand, all of the transgenic plants showed a similar number of primary cauline branches as the brc1 or the wild-type (Figures 8D and 9D). These results imply that BpTCP8, BpTCP14 and BpTCP19 prevent primary rosette branch outgrowth.

Discussion
As a kind of vigorous pioneer woody plant, B. papyrifera is widely used in livestock feed, papermaking and ecological afforestation. Genome-wide identification and analysis of the BpTCP gene family are of great significance for understanding the development of leaf and branch of B. Different lowercase letters denote significant differences (p < 0.05).

Discussion
As a kind of vigorous pioneer woody plant, B. papyrifera is widely used in livestock feed, papermaking and ecological afforestation. Genome-wide identification and analysis of the BpTCP gene family are of great significance for understanding the development of leaf and branch of B. papyrifera. In this study, we carried out a multilevel analysis of BpTCP genes with the aim to understand the important and diverse roles involved in the growth and development of B. papyrifera.

Overview of BpTCP Gene Family
The analysis of evolution relationship exhibited that 20 BpTCP proteins were classified into three major subclades, which was similar to the previous studies of A. thaliana, O. sativa, V. vinifera, Z. mays and P. vulgaris [42][43][44]48]. In addition, the composition of conserved motifs in each subclade member showed its particular features. For instance, motif 2, motif 4 and motif 10 were unique to the PCF subclade. All CYC/TB1 subclade members (BpTCP8, BpTCP14 and BpTCP19) and two CIN subclade members (BpTCP10 and BpTCP15) contained motif 6 (that is, R domain), which is similar to that in VvTCPs, PmTCPs and FvTCPs [42,49,50]. R domain is hypothesized to form a hydrophilic α-helix or a coiled coil structure that may mediate protein-protein interactions [1]. The exon/intron organization of BpTCPs was most similar to FvTCPs, among which only one CYC/TB1 subclade gene had two exons in the coding region and the rest members had one exon [50]. In addition, the number of BpTCP family members was relatively conserved with A. thaliana (24), O. sativa (21), Sorghum bicolor (19), Medicago truncatula (21), Prunus persica (20) and M. notabilis (22), but was smaller than M. domestica (58), Glycine max (56), M. acuminate (45) and Z. mays (44) (Figure 2). Combined with genome size, it was found that the number of TCP family members was not related to the genome size. For example, the genome size of Phyllostachys heterocycla is 2075 Mb with only 17 TCP members, while the genome size of M. domestica is 742.3 Mb with 58 TCP members [51,52]. The diversity of the number of TCP family members in different species may be influenced by genome duplication events, such as whole genome duplication, segmental duplication or tandem duplication. The analysis of duplication events of BpTCP genes showed that there were no tandem duplication events, but there were four segment gene pairs. It suggested that segmental duplications may be the primary driving force for the expansion of the BpTCP gene family, as described previously, in V. vinifera and Gossypium raimondii [42,53].

Potential Functions of BpTCP Genes Inferred from the Expression Patterns
The RNA-Seq data of nine tissues (shoot apex, young leaf, developing leaf, mature leaf, immature stem, phloem of proximal stem, phloem of mature stem, phloem of root and root tip) were used to investigate the expression patterns of BpTCP genes. According to the expressive level of BpTCP genes in specific tissues, the development process in which they probably participate could be speculated. As previously reported, class I (PCF) genes usually promoted cell growth and proliferation [14,54]. In contrast, class II (CIN and CYC/TB1) members had an antagonistic effect on some biological processes compared with class I members. Class II members majorly played a vital role in preventing cell proliferation and tissue overgrowth [55]. CIN subclade genes mainly regulated leaf development.
The cin mutant showed a phenotype of negative leaf curvature because the regulation of cell division was disturbed [10]. In Arabidopsis, the tcp2 tcp4 mutant exhibited enlarged flat leaves and the quadruple mutant tcp2 tcp3 tcp4 tcp10 exhibited strongly crinkled leaves [9,29]. CYC/TB1 subclade genes mainly participate in lateral branch and floral symmetry. The cyc mutant showed semipeloric flowers and the tb1 mutant showed excessive axillary branches [2,56]. AtTCP1 is the ortholog gene of AmCYC, which could regulate the longitudinal elongation of the petioles, rosette leaves and inflorescent stems [57]. AtTCP18 and AtTCP12, two homologs of ZmTB1, were proved to suppress bud outgrowth [11,34]. In B. papyrifera, most BpTCP genes from the PCF subclade were relative highly expressed in almost all test samples and less tissue-specific expression patterns, which were similar with VvTCPs and MtTCPs [42,58]. It implied that BpTCP genes of the PCF subclade might be involved in multiple growth and development stages. Two-thirds of CIN subclade genes were abundantly expressed in shoot apex, young leaf, developing leaf and mature leaf. These findings indicated that BpTCP1, BpTCP10, BpTCP13 and BpTCP15 may play an important role in leaf development. In the CYC/TB1 subclade, combined the RNA-Seq and quantitative RT-PCR data, BpTCP8 and BpTCP19 were transcribed at relatively high levels in leaf axil and axillary bud. It suggested that BpTCP8 and BpTCP19 may regulate the axillary bud outgrowth. Taken together, BpTCP genes may play diverse biological functions in organ or tissue development and shoot branching.

BpTCP8, BpTCP14 and BpTCP19 Prevent Branch Outgrowth
It was shown that the BRC1 gene was the focal points for multiple environmental and developmental signals that acted in the axillary buds to inhibit shoot branching. The homologs of AtBRC1 in tomato, pea, chrysanthemum, poplar and cucumber could inhibit lateral bud or branch outgrowth [59][60][61][62][63]. In chrysanthemum, DgBRC1s were mainly expressed in the nodes containing axillary buds and expression of DgBRC1s in Arabidopsis brc1 or wild-type could reduce the number of rosette branches [61]. In woody plant poplar, the loss-of-function mutant of PcBRC1 exhibited strongly enhanced bud outgrowth. The mutant of PcBRC2 showed an extreme bud outgrowth and had two ectopic leaves at each node, which was different with Arabidopsis. PcBRC2 may have retained or evolved the function of controlling leaf development [62]. In our work, we validated the function of BpTCP8, BpTCP14 and BpTCP19 on shoot branching in both Arabidopsis wild-type and brc1 mutant backgrounds. The functions of BpTCP8, BpTCP14 and BpTCP19 on inhibiting rosette branch outgrowth in Arabidopsis was similar to chrysanthemum [61]. In addition, the difference of our study from previous studies is that all BpTCP genes (BpTCP8, BpTCP14 and BpTCP19) in the CYC/TB1 subclade possessed the ability to suppress rosette branch outgrowth in Arabidopsis, not just the ortholog genes of AtBRC1. The analysis of duplication events showed that these three genes are involved in two segmental duplication events (BpTCP8/ BpTCP14 and BpTCP8/ BpTCP19), which may be the reason for the functional redundancy in regulating the branch outgrowth. The CYC/TB1 subclade genes are considered to regulate shoot branching and floral symmetry [64]. Therefore, BpTCP8, BpTCP14 and BpTCP19 could prevent lateral branch growth, same as shown in previous research, and may participate in leaf development like that of poplar.

Identification of Putative BpTCP Genes in B. papyrifera
To identify the TCP family in B. papyrifera genome, two approaches were used. Local BLASTp searches were performed against the B. papyrifera genome database using the 24 known Arabidopsis TCP protein sequences, setting the threshold of the E-value to 1 × 10 −5 for the initial identification. Additionally, we downloaded the hidden Markov model (HMM) seed file of TCP domain (PF03634) from the Pfam database (http://pfam.xfam.org/), and used the HMMER software to carry out the HMM searches against the local B. papyrifera protein sequence database, setting the threshold of the E-value to 0.01. Summarizing the results of both methods and removing the redundant sequences, the remaining sequences were the candidate TCP protein sequences of B. papyrifera. Subsequently, all candidate BpTCP protein sequences were submitted to the NCBI Conserved Domain Search Service (CD Search) (https: //www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) [65] and InterPro (http://www.ebi.ac.uk/interpro/) database for the further verification. The Arabidopsis TCP proteins were downloaded from The Arabidopsis Information Resource (TAIR) database (http://www.arabidopsis.org). The molecular weight (MW) and isoelectric point (pI) of each BpTCP protein was predicted by the ExPASy program (https://web.expasy.org/protparam/). The subcellular localization of each BpTCP protein was predicted by the CELLO server (http://cello.life.nctu.edu.tw/). The online tool P3DB (http://www.p3db.org/) was used to carry out the phosphorylation analysis.

Phylogenetic Analysis
For phylogenetic analysis, the full-length amino acid sequences of TCP proteins from Arabidopsis, O. sativa, B. papyrifera, Z. mays (ZmTB1) and A. majus (AmCYC, AmCIN) were aligned by MEGA X software with default parameters [66]. The phylogenetic tree was constructed by the maximum likelihood (ML) method with 1000 bootstrap replicates. The JTT model was selected as the optimal model. The O. sativa TCP protein sequences were retrieved from the PlantTFDB (http://planttfdb.cbi. pku.edu.cn/). The sequences of ZmTB1, AmCYC and AmCIN were downloaded from the Phytozome database v12.1 (https://phytozome.jgi.doe.gov/pz/portal.html).

Gene Structure and Conserved Motif Analysis
The gene structure of BpTCP genes was analyzed by the online website GSDS 2.0 (http://gsds.cbi. pku.edu.cn). The conserved motif was predicted by the Multiple Em for Motif Elicitation (MEME) server (http://meme-suite.org/), and the program was performed with the following settings: site distribution = any number of repetitions; 0-order model of sequences for the background motif; motif width = 15-60 and maximum number of motifs = 10.

Chromosomal Localization and Duplication Event Analysis
The chromosomal location information of all BpTCP genes was obtained from the B. papyrifera genome database. The visualized diagram of chromosomal localization and the duplication events in the B. papyrifera were acquired from the TBtools software [67].

Expression Analysis of BpTCP Genes in Different Tissues
Expression level of 20 BpTCPs in different tissues was detected by an Illumina HiSeq 2500 platform. The shoot apex, young leaf, developing leaf, mature leaf, immature stem, phloem of proximal stem, phloem of mature stem, phloem of root and root tip were collected from a five-year-old female B. papyrifera and immediately frozen in liquid nitrogen. The total RNA was isolated and purified to construct the RNA-Seq libraries, then sequenced. Each tissue had three biological replicates. Gene expression was measured as fragments per kilobase of transcript per million fragments mapped (FPKM) using Cuffquan and CuffnormGene in Cufflinks. The expression values of each gene in different tissues were averaged and presented as a log value. The heat map was exhibited using the TBtools software [67]. The FPKM value expressed in different tissues are listed in Table S2.

Plant Materials and Growth Conditions
Arabidopsis ecotype Columbia-0 (Col-0) was used as the wild-type. The mutant brc1 was described previously [11]. Murashige and Skoog medium with or without 25 mg/mL hygromycin B was used for screening or growing the plant seeds. The plates with seeds were put at 4 • C for two days for synchronization before cultivation at 22 • C under long-day conditions with a photoperiod of 16 h/8 h (light/dark), photosynthetic photon flux density of 200 µmol m −2 s −1 and 60% relative humidity. The culture matrix was a mixture of vegetative soil and vermiculite with a ratio of 2:1.
The plantlets of B. papyrifera were cultured on the MS medium at 26 • C with a photoperiod of 14 h/10 h (light/dark) and photosynthetic photon flux density of 80 µmol m −2 s −1 . A month later, the shoot apex, leaf axil, axillary bud, leaf, stem and root were collected for quantitative RT-PCR.

Quantitative RT-PCR
Total RNA was extracted from different tissues using the TransZol Plant Mini kit (TransGen, Beijing, China) according to the manufacturer's instructions. The first-strand cDNA was synthesized by using a PrimeScript TM RT reagent kit with gDNA eraser (Takara, Dalian, China). The quantitative RT-PCR reactions were carried out on a StepOne TM real-time PCR system (Applied Biosystems, Forster City, USA) using the SYBR ® Premix Ex Taq TM II kit (Takara, Dalian, China). Each reaction was operated in a 20 µL volume, which contained 2 µL diluted cDNA, 0.4 µL forward/reverse primer (10 µmol/L), 6.8 µL sterilized ddH 2 O, 0.4 µL ROX Reference Dye II and 10 µL SYBR ® Premix Ex Taq II. Three technical replicates were performed for each biological sample. The Glyceraldehyde-3-phosphate Dehydrogenase (GAPDH) gene of B. papyrifera was used as an internal control. The 2 −∆∆Ct method was used to calculate the data. The primer sequences of each gene involved in quantitative RT-PCR are listed in Table S3.

Transcriptional Activation Activity Analysis
The coding region of the BpTCP8, BpTCP14 and BpTCP19 was amplified and cloned into the pBridge vector to generate a fusion protein with the GAL4 DNA binding domain, respectively. The fusion constructs and control vector were transferred into the AH109 yeast strain. Then the transformed yeast cells were dropped on SD-Trp-His medium with 0-50 mM 3-AT to detect their activities. The plates were cultured at 30 • C for 3-5 days before photographing.

Vector Construction and Plant Transformation
The coding region of the BpTCP8, BpTCP14 and BpTCP19 was inserted into the pCAMBIA1300 vector containing a cauliflower mosaic virus (CaMV) 35S promoter. The constructs were then transformed into the Agrobacterium tumefaciens strain EHA105. Both wild-type and brc1 plants were transformed by the floral dip method. In addition, transgenic plants were selected on the MS medium with 25 mg/mL hygromycin B for one week. Homozygous plants were used for the branching phenotype analysis in this study. The eukaryotic initiation factor 4A (eIF4A) gene of Arabidopsis was used as a normalization control. Primers used for semi-quantitative are listed in Table S3.

Conclusions
In this study, we identified 20 BpTCP genes, which distributed on 10 chromosomes. These BpTCPs were divided into three subclades according to the phylogenetic relationships and structural properties. Segmental duplication was the predominant duplication event, which induced the expansion of BpTCP genes. Expression patterns in different tissues suggested that BpTCP genes may play important roles in B. papyrifera growth and development. BpTCP8, BpTCP14 and BpTCP19 were verified that they possessed transcriptional activation ability, indicating that they were functional transcription factors. In addition, the function of BpTCP8, BpTCP14 and BpTCP19 in the regulation of shoot branching was confirmed. Ultimately, these findings will provide a solid foundation for the further functional investigation of BpTCP genes and the improvement of new cultivars via genetic engineering. Accession Number: Sequence data from this article can be found in the NCBI database under Sequence Read Archive (SRA) ID SRP136442.

Conflicts of Interest:
The authors declare no conflict of interest.