Introduction {#Sec1}
============

As sessile organism, plant growth and yield are strongly influenced by environmental stimuli such as cold and drought^[@CR1]^. To respond and adapt to these conditions, plants develop various mechanisms at both physiological and biochemical levels^[@CR2]^. It has been well established that many adaptation processes are regulated by stress-responsive gene expression^[@CR1],\ [@CR3]^. Transcription factors (TFs), which are a diverse family of regulatory proteins with DNA-binding domains, play a central role in mediating various aspects of cellular processes by regulating gene expression through interacting with cis-elements in the promoter regions of various downstream genes^[@CR4],\ [@CR5]^. Series studies previously have uncovered a group of TF genes, such as *AP*2/*ERF*, *MYB* and *bZIP*, which participate in various stress-induced physiological processes and regulatory networks in higher plants^[@CR6],\ [@CR7]^.

*TCP* genes encode plant-specific transcription factors, which are named after the first four functionally identified members: *TB1* (*TEOSINTE BRANCHED 1*) in *Zea mays*, *CYC* (*CYCLOIDEA*) in *Antirrhinum majus*, and *PCF1* and *PCF2* (*PROLIFERATING CELL FACTORS 1* and *2*) in *Oryza sativa* ^[@CR8]^. Typically, N-terminus of this class of transcription factors exhibits a highly conserved TCP domain, which contains a non-canonical basic-Helix-Loop-Helix (bHLH) structure involved in DNA binding, protein-protein interaction and protein nuclear localization^[@CR9],\ [@CR10]^. Based on the homology of the TCP domains, TCP proteins can be further divided into two major subfamilies, Class I (represented by the rice PCF proteins) and class II (represented by CYC and TB1)^[@CR8],\ [@CR11]^. The DNA-Binding site selection assays revealed that the consensus binding sequences for these two classes are slightly different, but overlapping with GGNCCC sequences. The DNA binding sequence for class I is GGNCCCAC while class II prefer to bind the DNA motif G(T/C)GGNCCC^[@CR10]^. Both class I and class II include TCPs that can function as transcriptional activators and repressors^[@CR12]^.

Increasing evidences have indicated that proteins of TCP family take part in the regulation of many biological processes during plant growth and development, including plant architecture^[@CR12],\ [@CR13]^, leaf morphogenesis^[@CR14]--[@CR16]^, hormone pathways^[@CR13],\ [@CR17]--[@CR19]^ and response to environmental stimuli among various species^[@CR20]--[@CR22]^. For example, studies in *Arabidopsis* suggest that *AtTCP14* appears to function in regulating embryonic growth of seeds^[@CR23]^. *AtTCP15*, along with its closest homolog *AtTCP14*, regulates cell proliferation in the developing leaf blade and specific floral tissues, and also modulates gibberellins and auxin responses^[@CR17],\ [@CR24]^. The expression of a repressor form of *AtTCP11* caused pleiotropic developmental alterations^[@CR25]^. *AtTCP2*, *AtTCP3*, *AtTCP4*, *AtTCP10* and *AtTCP24* have been earlier identified as targets of miRNA319, and are required for leaf development, petal growth, cell wall synthesis and JA synthesis^[@CR14],\ [@CR26],\ [@CR27]^. Creeping bentgrass (*Agrostis stolonifera*) plants overexpressing Osa-miR319a, in which four putative target genes, *AsPCF5*, *AsPCF6*, *AsPCF8*, and *AsTCP14* were down-regulated, displayed morphological changes and exhibited enhanced drought and salt tolerance associated with increased leaf wax content and water retention^[@CR20]^. In rice, overexpression of *OsTCP19* led to developmental abnormalities like fewer lateral root formation and contributed to better stress tolerance^[@CR21]^.

Widely cultivated in tropical area, cassava is considered as one of the most important economic crops worldwide, providing starch for food, feed and energy production^[@CR28]--[@CR30]^. Cassava can effectively utilize light, heat and water resource; it is resistant to drought but sensitive to low temperature^[@CR31]^. Cold and drought stress severely limit cassava plant growth and yield^[@CR32]^. Thus, uncovering the mechanisms underlying the resistance of cassava to these stresses may provide candidate genes for genetic improvement of stress resistance for cassava. Previously, we performed strand specific RNA sequencing for cassava TMS60444 shoot apices and young leaves under cold and PEG-induced drought stress, providing an excellent resource for analysis of stress responsive genes globally.

To date, numerous *TCP* family members have been identified in various species^[@CR33]--[@CR45]^. However, no evidence is available regarding the *TCP* family in cassava. Due to the critical role of *TCP* transcription factors in the control of plant development and abiotic stress responses, we performed for the first time the comprehensive analysis of the *MeTCP* transcription factor family in cassava. In the present study, a total of 36 non-redundant *MeTCP* transcription factor encoding genes were identified in cassava genome and were subsequently subjected to a systematic analysis, including phylogenetic relationships, gene structure, conserved motif and expression pattern among different tissues and hormone treatments. We also analyzed the expression of the *MeTCP* genes under normal growth conditions and various abiotic stressors. On the basis of the expression profiles of *MeTCP* genes and the phylogenetic analysis among the TCP domain proteins in *Arabidopsis*, rice and cassava, the functions of *MeTCPs* were predicted. Taken together, our genome-wide analysis of *MeTCP* gene family will contribute to future studies on the functional characterization of MeTCP proteins in cassava, as well as the identification and comprehensive analysis of the *TCP* transcription factor family in other species.

Results {#Sec2}
=======

Identification of *TCP* genes in cassava {#Sec3}
----------------------------------------

In order to extensively identify cassava *TCP* genes, two search strategies were used in this study: first, we used the protein sequences of all 24 *TCP* genes from *Arabidopsis* ^[@CR33]^ as BLASTp queries to perform multiple searches against the latest whole genome of cassava (<http://www.phytozome.net/cassava>); then, the TCP domain (PF03634)^[@CR46]^ was employed as query to perform a blast search against the same cassava genome database. After discarding redundant sequences, 36 candidate *TCP* genes were identified in cassava, and further conserved domain detection confirmed that all the identified *TCPs* harbor the conserved TCP domain that is the basic characteristics of this family. Due to the lack of standard annotation designated to the 36 cassava *TCP* genes, we named them *MeTCP2* to *MeTCP23* according to the *Arabidopsis* TCP proteins with highest sequence similarity. The length of the 36 newly identified TCP proteins varied from 73 to 562 amino acid residues and the relative molecular weight ranged from 8.3 to 58.1 kDa, with isoelectric points in the range of 5.53--10.08. The grand average of hydropathy (GRAVY) of this family genes showed that all of the proteins had a negative value, indicating that all the MeTCP proteins were hydrophilic (Table [1](#Tab1){ref-type="table"}).Table 1*TCP* genes identified in cassava genome.Gene IDNameLength(aa)MW(Da)PIGravyGenomic locusManes.01G187000.1MeTCP15a39242214.99.32−0.643Chromosome01:28455200..28457428 forwardManes.01G263300.1MeTCP9a36638916.338.94−0.283Chromosome01:33561357..33563642 forwardManes.01G020700.1MeTCP18a38744030.357.93−0.773Chromosome01:3502035..3504116 forwardManes.01G094200.1MeTCP13a35038533.918.37−0.667Chromosome01:21880606..21882675 forwardManes.02G055900.1MeTCP13b36139880.769.06−0.63Chromosome02:4209171..4211393 forwardManes.02G066400.1MeTCP1238342814.849.78−0.569Chromosome02:4925458..4927250 reverseManes.02G194200.1MeTCP23a42545005.567.98−0.572Chromosome02:15997522..15999455 forwardManes.04G016700.1MeTCP11d24026058.478.98−0.394Chromosome04:1959602..1961886 forwardManes.04G016800.1MeTCP11b15516551.849.45−0.415Chromosome04:1966346..1966813 forwardManes.04G088500.1MeTCP20a31233848.597.22−0.782Chromosome04:22580287..22581884 reverseManes.05G041000.1MeTCP9b34636898.797.8−0.325Chromosome05:2911279..2913198 forwardManes.05G100100.1MeTCP15b38841706.518.23−0.582Chromosome05:8461910..8463076 reverseManes.05G123700.1MeTCP20c738314.6410.08−1.034Chromosome05:14666983..14667397 forwardManes.05G119300.1MeTCP18b37242181.388.12−0.748Chromosome05:12207088..12208206 reverseManes.06G072800.1MeTCP18c41647024.779.48−0.879Chromosome06:18738590..18740594 forwardManes.06G083400.1MeTCP5a34037976.228.47−0.611Chromosome06:19670814..19672682 reverseManes.06G093900.1MeTCP15c39642673.97.37−0.73Chromosome06:20627278..20631189 forwardManes.06G141800.1MeTCP19358377475.53−0.478Chromosome06:24669381..24671688 reverseManes.07G022400.1MeTCP8a55056832.227.55−0.646Chromosome07:2093408..2096799 forwardManes.08G009200.1MeTCP11a18519802.376.97−0.497Chromosome08:675290..675847 reverseManes.09G051000.1MeTCP1640343569.257.58−0.11Chromosome09:6817330..6819454 forwardManes.10G120400.1MeTCP8b56258079.718.46−0.621Chromosome10:23214679..23218733 reverseManes.11G083000.1MeTCP20b30732859.718.99−0.671Chromosome11:11404469..11406607 forwardManes.11G108500.1MeTCP727327679.869.72−0.271Chromosome11:20090517..20092799 forwardManes.11G149000.1MeTCP11c19821355.389.43−0.335Chromosome11:26034972..26035568 reverseManes.12G007700.1MeTCP20d29932367.129.41−0.702Chromosome12:722180..724162 reverseManes.13G008300.1MeTCP20e28230554.129.52−0.656Chromosome13:815757..817107 reverseManes.13G138300.1MeTCP3b34337809.475.9−0.782Chromosome13:26570479..26573043 reverseManes.14G058400.1MeTCP3a33636987.575.87−0.699Chromosome14:4620482..4621572 reverseManes.14G077200.1MeTCP15d39842850.167.39−0.716Chromosome14:6245053..6246249 reverseManes.14G086500.1MeTCP5b38743276.037.21−0.713Chromosome14:6949867..6951866 forwardManes.14G097000.1MeTCP18d47452697.969.31−0.782Chromosome14:7827437..7828861 reverseManes.15G091000.1MeTCP442245753.276.17−0.636Chromosome15:6736704..6738513 reverseManes.15G123800.1MeTCP2a48152167.17.86−0.886Chromosome15:9374446..9382168 reverseManes.17G072800.1MeTCP2b48152472.37.06−0.926Chromosome17:21184852..21198102 reverseManes.18G103100.1MeTCP23b42545159.976.81−0.563Chromosome18:9144633..9147156 forward

Phylogenetic analysis of the *MeTCP* genes {#Sec4}
------------------------------------------

To characterize the evolutionary and phylogenetic relationships between cassava *TCP* genes and other known *TCPs*, an unrooted Neiboring-Joining (NJ) tree was constructed on the basis of multiple sequence alignment of 36 MeTCP complete protein sequences with 24 and 22 TCP protein sequences from *Arabidopsis* and rice, respectively (Fig. [1](#Fig1){ref-type="fig"}). At deep nodes, the phylogenetic relationship was unclear and the bootstrap values were low as a result of relatively large number of sequences. To verify the reliability of our phylogenetic tree, we also build the phylogenetic trees of TCP transcription family with Minimal Evolution methods (Fig. [S1](#MOESM1){ref-type="media"}). The tree topologies were robust within different tree-building methods, except at the deep nodes. Considering these results, the NJ tree was employed for further study.Figure 1Phylogenetic relationships of TCP transcription factors from cassava, *Arabidopsis* and rice. A total of 36 *MeTCPs* from cassava, 24 *AtTCPs* from *Arabidopsis* and 22 *OsTCPs* from rice were used to construct the Neighbor-Joining tree by MEGA 6.0 with 1000 bootstrap based on the full length sequences of *TCPs*. The eight subgroups are indicated with different colors.

Based on the bootstrap value of clade and topology of the tree, the MeTCP proteins could be distributed into 8 distinct groups, designated as Group A to Group H. In general, *TCPs* from cassava have closer relationships with the *TCPs* from dicot plant *Arabidopsis* than that from monocot plant rice, which is accord with the current understanding of plant evolutionary history. Additionally, the *TCP* genes showed an interspersed distribution in most clades, which is consistent with those found in previous analyses of *TCP* in *G. raimondii* and *Citrullus lanatus* ^[@CR35],\ [@CR43]^, indicating that the TCP family expanded before the divergence of the lineages. However, the *TCP* genes were not evenly distributed in some clades, such as the largest clade Group H has the maximum 7 members, whereas Group C contains only 3 MeTCP genes from cassava, suggesting the existence of a diversified *MeTCP* family in cassava with diverse functions (Fig. [1](#Fig1){ref-type="fig"}). Remarkably, many *Arabidopsis TCP* genes had more than three counterparts in cassava, such as *MeTCP11*, *MeTCP15*, *MeTCP18* and *MeTCP20*, indicating that *MeTCP* genes duplicated after the divergence of cassava and *Arabidopsis*. It also suggests that higher number of genes in cassava as compared to *Arabidopsis* is the result of more gene duplication events in cassava or higher frequency of retaining copies after duplication. Group C contained three cassava *TCPs*, two *Arabidopsis* members but there were no *TCP* from rice, implying this group was either acquired after the divergence of monocots and dicots or lost in rice. Remarkably, 5 of the group H members (*AtTCP2*, *AtTCP3*, *AtTCP4*, *AtTCP10*, and *AtTCP24*) are post-transcriptionally targeted by miRNA319 in *Arabidopsis*. The closest homologs of these *Arabidopsis* genes in cassava are these five genes: *MeTCP2a*, *MeTCP2b*, *MeTCP3a*, *MeTCP3b* and *MeTCP4*, all containing putative target site of miR319^[@CR47],\ [@CR48]^. This suggests that regulation of leaf development by a redundant set of miRNA-regulated homologous *TCP* genes occurs in cassava.

Gene structure and conserved motifs of cassava TCPs {#Sec5}
---------------------------------------------------

To further examine the structural features of cassava *TCP* genes, we investigated the exon/intron structures of individual *MeTCP* genes by alignment of cDNA sequences and corresponding genomic DNA sequences. Additionally, we also built an unrooted phylogenetic tree with MeTCP protein sequences (Fig. [2a](#Fig2){ref-type="fig"}), to determine whether the gene structure of MeTCPs is consistent with the phylogenetic subfamily. As illustrated in Fig. [2b](#Fig2){ref-type="fig"}, the number of introns of *MeTCP* genes varied from 0 to 4. For instance, 32 out of 36 *MeTCP* genes had no intron, while the other *MeTCP* genes possess 1--3 introns, with the exception of *MeTCP16* containing four introns. As expected, most of *MeTCP* genes in the same subfamily showed similar exon-intron distribution patterns in terms of exon length and intron number, which supports their close evolutionary relationship and the classification of subgroup. However, *MeTCP* genes in group H showed great variability in intron number and exon length.Figure 2The gene structure and conserved protein motifs of *MeTCP* genes according to the phylogenetic relationship. (**a**) The unrooted phylogenetic tree of all *TCP* genes in cassava was constructed using Neighbor-Joining method and the bootstrap test was performed with 1,000 iterations. (**b**) The gene structure with exon/intron organization of TCP genes of cassava. The orange boxes represent 5′-UTR or 3′-UTR, red boxes represent exons and blue lines indicate introns. (**c**) The conserved protein motifs in the TCP family were identified using MEME program. Each motif is indicated with a specific color.

To obtain more insights into the diversity of motif compositions among MeTCPs, conserved motifs were predicted by using MEME program^[@CR49]^, and annotated by ScanProsite program^[@CR50]^. A total of 20 conserved motifs in the MeTCP proteins, designated as motif 1 to motif 20, were captured by MEME (Fig. [2c](#Fig2){ref-type="fig"}, Fig. [S2](#MOESM1){ref-type="media"}). The results showed that the only motif that hit for the database was the conserved TCP domain (motif 1). TCP domain was found in all MeTCPs, except MeTCP20c, which contained a truncated TCP domain. In general, MeTCP proteins clustered in same subgroup share similar motif composition, while high divergence was observed among different subgroups. This observation indicated that the MeTCP members within the same subgroup may have redundant functions and that some motifs may contribute to the specific function of that subfamily, which is in agreement with the previous report^[@CR35],\ [@CR41]^. According to 36 MeTCP sequence features within the TCP domain, we determined that *MeTCPs* from Group A, B, C, D and E belong to class I subfamily while the other *MeTCPs* belong to class II subfamily, as for all species so far. As reported earlier, the group members belonged to class I subfamily have extended homology C-terminal from the TCP domain, while the class II subfamily has an extended basic region, and all groups have internally conserved, but distinct loop region sequences^[@CR25],\ [@CR51]^. The motif analysis also showed that sequence conservation outside the TCP domain was low and sequence length on both sides of the TCP domain varied greatly, resulting in proteins ranging from 73 (MeTCP20c) to 563 (MeTCP8b) amino acids. For example, MeTCP20c, the smallest predicted protein, is probably truncated by a frame shift mutation that cause premature termination, since sequence homology with *Arabidopsis* TCP20 extends well beyond the stop codon. This result is similar to SlTCP27, which encode the smallest TCP protein with 113 amino acids in tomato^[@CR36]^. Although MeTCP20c lacks the conserved C-terminal part of the TCP domain, which may alter its DNA binding ability, experimental evidences are required to establish the precise role of truncated TCP domain in the regulation of MeTCP20c activity.

Expression analysis of *MeTCP* genes in different tissues {#Sec6}
---------------------------------------------------------

To investigate the potential functions as well as to identify probable functional redundancy through similar expression patterns for the cassava *TCP* genes, the detection of their expression were carried out in different tissues including root, leaf, stem and shoot apex using qRT-PCR. As shown in Fig. [3](#Fig3){ref-type="fig"}, it is apparent that the expression levels in different tissues vary widely between the cassava *TCP* genes, as well as between different tissues for individual *TCP* genes, indicating functional specialization among *TCP* gene family members in cassava plant development. Of them, some genes were exclusively highly expressed in a specific tissue. For example, *MeTCP20e* and *MeTCP11d* genes exhibited higher transcriptional abundance in roots as compared to other organs; *MeTCP2a*, *MeTCP3a*, *MeTCP5b*, *MeTCP8a*, *MeTCP13a* and *MeTCP20d* showed specifically high expression in leaves, implying their specific roles in the corresponding tissues. Opposite to the tissue-specific expression pattern of *MeTCP* genes, many genes were more widely and less specifically expressed, such as *MeTCP9b*, *MeTCP13b*, *MeTCP20a*, *MeTCP20b*, and *MeTCP23a* genes, implying that these genes may play regulatory roles at multiple development stages. However, further studies are still needed to unravel the divergent roles of *MeTCP* genes.Figure 3Heatmap representation for expression profiles of 36 *MeTCP* genes across different tissues. The expression levels of *MeTCP* genes were obtain through quantitative real-time PCR. *MeACTIN* was used as the reference gene.

To address the conservation and specificity of *TCP* expression pattern beyond species, we compared the expression level of homologous *TCP* gene pairs in these four tissues between *Arabidopsis* and cassava. We found a subset of *TCP* genes were positively correlated with pearson correlation coefficient (PCC) higher than 0.3 between *Arabidopsis* and cassava (Fig. [S3](#MOESM1){ref-type="media"}--[4](#MOESM1){ref-type="media"}, such as *MeTCP2a/b* and *AtTCP2*, *MeTCP5a/b* and *AtTCP5*, *MeTCP8a/b* and *AtTCP8*, *MeTCP13a/b* and *AtTCP13*, *MeTCP15b/c* and *AtTCP15*, *MeTCP19* and *AtTCP19 et al*., indicating functional conserved expressional pattern of these genes. In contrast, some *TCP* genes show no correlation or negative correlation between *Arabidopsis* and cassava, such as *MeTCP9a/b* and *AtTCP9*, *MeTCP11a/b/c/d* and *AtTCP11*, suggesting their functions have been diversely changed in different species.

Expression patterns of *MeTCP* genes in response to hormone treatments {#Sec7}
----------------------------------------------------------------------

Multiple studies have been reported that TCP proteins regulate plant development and environmental stress adaption by mediating hormone biogenesis and response^[@CR17]--[@CR19],\ [@CR21]^. To study the total effect of plant hormones on *MeTCP* genes, the expression levels of 36 *MeTCP* genes were detected in response to abscisic acid (ABA), gibberellin (GA3), indole acetic acid (IAA), jasmonate acid (JA), zeatin (ZT) and 6-benzylaminopurine (6-BA) hormone treatments by quantitative RT-PCR (Fig. [4](#Fig4){ref-type="fig"}, Supplemental Table [S1](#MOESM2){ref-type="media"}). In general, hormone treatments resulted in a wide variety of *MeTCP* gene expression profiles. In JA treatment, 15 and 7 *MeTCP* genes were obviously induced and inhibited, respectively. Of them, the most up-regulated gene was *MeTCP20e*, and the most down-regulated gene was *MeTCP3a*. Similarly, 6-BA and ZT treatment led to 12 and 13 *MeTCP* genes were obviously induced, 11 and 5 *MeTCP* genes were inhibited, respectively. In GA treatment, most of *MeTCPs* were significantly induced and only two *MeTCPs* were inhibited. *MeTCP7* and *MeTCP23b* were found to be most up-regulated. As for IAA treatment, 10 and 5 *MeTCP* genes showed dramatic increase and decrease, respectively. *MeTCP7* and *MeTCP16* went through the largest increase and decrease, respectively. ABA plays a crucial role in the adaptive response of plants to abiotic stresses^[@CR52]^. We found 22 members showed strong sensitivity toward ABA, indicating that these genes may be regulated by ABA signal pathway. Among them, 19 genes had relative high levels of transcript abundance after ABA treatment. Notably, most of *MeTCP* genes responded to at least one treatment; Particularly, *MeTCP11b*, *MeTCP12* and *MeTCP21* responded to all hormone treatments, indicating these genes might play pivotal roles in the cross-talk of hormones, which would be candidates for further research in the field. However, we also found *MeTCP11d*, *MeTCP20a* and *MeTCP23a* were not able to respond to any treatments. Taken together, these results suggest the complicated regulatory mechanism of *MeTCP* genes in response to hormone treatments in cassava.Figure 4Heatmap representation for expression patterns of *MeTCP* genes under various hormone treatments. The expression profile data of *MeTCP* genes under JA, 6-BA, GA, ZT, ABA and IAA were obtain through quantitative real-time PCR. *MeACTIN* was used as the reference gene.

Expression profiles of *MeTCP* genes in response to cold and/or drought stress {#Sec8}
------------------------------------------------------------------------------

Plants are frequently challenged by abiotic stressors such as cold and drought. Recent studies have suggested that TCP proteins are widely involved in signaling and response to environmental stimuli^[@CR18],\ [@CR22]^. However, information on the involvement of TCP proteins in stress responses in cassava is limited. To investigate the potential roles of *MeTCP* genes in response to abiotic stresses, cassava seedlings of TMS60444 genotypes were subjected to cold (4 °C) and PEG (20% PEG 6000)-induced drought stress and then the leaves tissues were sampled to extract RNA for subsequent RNA-seq analysis. According to the transcriptome data, 18 (50%) and 24 (66.7%) *MeTCP* genes showed significantly change (fold change \>2) under cold and drought treatment, respectively. Among them, 7 (38.9%) and 11 (61.1%) *MeTCP* genes were up- and down-regulated by cold, respectively; 10 (41.7%) and 14 (58.3%) *MeTCP* genes were up- and down-regulated by drought, respectively (Fig. [5](#Fig5){ref-type="fig"}, Supplemental Table [S2](#MOESM2){ref-type="media"}). These results also showed that the number of *MeTCP* genes down-regulated by cold and drought was greater than that were up-regulated, suggesting the comprehensive response of *MeTCP* genes to cold and/or drought at transcriptional levels. There were 23 *MeTCP* genes in total differential expressed under both cold and drought treatments, which were categorized into 4 different classes: concordant response to cold and drought, discordant response, cold-specific and drought-specific. Among concordant response, *MeTCP20c/20e/11a* and *MeTCP18b/11c/12* were co-induced and co-repressed by two kinds of stresses, respectively. However, in discordant response class, three genes (*MeTCP15d*, *MeTCP15b* and *MeTCP16*) showed increased expression pattern under cold treatment, whereas down-regulated after drought treatment. By contrast, the expression levels of four genes (*MeTCP2b*, *MeTCP19*, *MeTCP13a* and *MeTCP13b*) were induced by drought treatment, but were repressed or unaltered after cold treatment. Meanwhile, *MeTCP8a/5a* and *MeTCP20a/20b/9b/18d/3b* were specifically response to cold and drought stress, respectively. Generally, the expression levels of *MeTCP* genes in response to cold and drought stresses were dramatically changed, implying their putative roles in stress tolerance.Figure 5Expression profiles of *MeTCP* genes in leaves and shoot apices after cold and drought treatment. The transcript data generated from three replicates of RNA-seq data. The relative expression values were log2 transformed. The bar represents relative expression values.

Expression analysis of *MeTCP* genes under various abiotic stress treatments {#Sec9}
----------------------------------------------------------------------------

To further assess the response of *MeTCP* genes to various abiotic stresses and related signaling pathway at transcriptional levels, 12 *MeTCP* genes (*MeTCP20c*, *MeTCP18a*, *MeTCP20e*, *MeTCP11a*, *MeTCP11b*, *MeTCP18b*, *MeTCP11c*, *MeTCP12*, *MeTCP18c*, *MeTCP4*, *MeTCP3a* and *MeTCP15c*) induced or repressed by cold and drought stresses based on RNA-seq data were chosen for further examination of their expression patterns after cold, drought and salt treatments (Fig. [6](#Fig6){ref-type="fig"}). We found most of the analyzed genes exhibited differential expression in response to at least one stress treatment, implying their putative roles in these stresses tolerance. Overall, in the cold, drought and salt conditions, expression of *MeTCP20c*, *MeTCP18a*, *MeTCP20e* and *MeTCP11a* genes were significantly up-regulated, and the largest expressional change of these genes were usually observed when responding to cold and/or drought treatments. This observation is well consistent with the RNA-seq data. 6 genes (*MeTCP11b*, *MeTCP18b*, *MeTCP11c*, *MeTCP12*, *MeTCP4* and *MeTCP3a*) were down-regulated following cold and drought treatments, among which, *MeTCP12* underwent the greatest change of mRNA expression levels. Interestingly, the expression levels of all these *MeTCP* genes showed obviously increase under salt treatment, suggesting they may play different roles in response to the three stresses. In addition, not too many changes were observed in *MeTCP18c* and *MeTCP15c* genes when any of three stresses were carried out. These data show the potential of some *MeTCP* genes for enhancing adversity resistant capacity in cassava.Figure 6Confirmation of the expression patterns of cold- and drought-responsive *MeTCPs* using qRT-PCR. The expression patterns of *MeTCP* under cold, drought and salt stress. The values shown are the means ± standard deviation of three replicates. *MeACTIN* was used as the reference gene.

Analysis of the regulatory cis-elements in the promoter of cold- and drought-responsive *MeTCPs* {#Sec10}
------------------------------------------------------------------------------------------------

The altered expression of 12 *MeTCPs* that were co-induced/repressed by cold and drought stresses indicates that they may be regulated by key stress regulatory genes. In order to elucidate the mechanism of transcriptional regulation of these genes, analysis of their promoter region was performed for the cis-elements. A 1.5 kb sequence upstream to the open reading frame of *MeTCP20c*, *MeTCP18a*, *MeTCP20e*, *MeTCP11a*, *MeTCP11b*, *MeTCP18b*, *MeTCP11c*, *MeTCP12*, *MeTCP18c*, *MeTCP4*, *MeTCP3a* and *MeTCP15c* was identified and subjected to MEME analysis. A number of common cis-acting elements were identified in the proximal promoters (Fig. [7](#Fig7){ref-type="fig"}, Fig. [S5](#MOESM1){ref-type="media"}). All the identified promoters had motif 1, and most of them also contained motif 3 and 10, except for *MeTCP20e*, *MeTCP11b*, *MeTCP11c*, *MeTCP12* and *MeTCP18c*, suggesting that these commonly present cis-acting elements may be involved in stress response. Notably, motifs such as 9, 13, 14 and 15 are more present in the promoters of *MeTCP* which are down-regulated by both cold and drought stresses, suggesting that these motifs may play key roles in the regulation of *MeTCP* members.Figure 7The conserved DNA sequence motifs analysis of cold- and drought-responsive *MeTCP* promoters. The conserved motifs in the TCP promoters were identified using MEME program. Each motif is indicated with a specific color.

Discussion {#Sec11}
==========

Adverse environmental conditions, such as cold and drought stress, impose severe effects on the plant growth, development and limit crop productivity and yield^[@CR53]^. TCP transcription factors are a class of plant-specific transcription that play very important roles during plant development and abiotic stress responses^[@CR22]^. To our knowledge, although a range of TCP family members have been described in various species^[@CR35],\ [@CR36],\ [@CR40]--[@CR45],\ [@CR54]^, no studies have been performed on *TCP* genes in cassava. Additionally, the mechanisms of cassava responds to abiotic stress are poorly understood. This background knowledge prompted us to identify the full complement and expression profile of this important gene family during development and under abiotic stresses in cassava.

In the present study, a comprehensive set of 36 non-redundant *TCP*-encoding genes were identified and characterized from the current version of the cassava genome. Previous studies have identified 24 TCP genes in *Arabidopsis* ^[@CR33]^, 22 in rice^[@CR34]^, 38 in cotton^[@CR35]^, 30 in tomato^[@CR36]^, and 27 in watermelon^[@CR43]^. The amplification of TCP gene family members in cassava can be explained by its larger genome size (\~760 Mb) compared to that of *Arabidopsis* (\~125 Mb) and gene duplication events in this family. Evolutionary analysis indicated that the cassava *TCPs* could be clustered into 8 subgroups, which is minimal different with previous evolutionary classification of TCPs in cotton and watermelon^[@CR35],\ [@CR43]^ (Fig. [1](#Fig1){ref-type="fig"}). The phylogenetic tree showed obvious differences in number of *MeTCPs* and ratio of *MeTCPs*/*AtTCPs* among subgroups. The change in the ratio of *MeTCPs*/*AtTCPs* suggested that the *MeTCP* family had undergone lineage-specific expansion and functional divergence during the course of evolution. Phylogenetic tree also showed that group C contained three cassava members and two *Arabidopsis* members, but there were no *TCP* from rice. On contrast, *TCPs* in subgroup A expanded in monocots but not in *Arabidopsis*, indicating that *TCP* genes of these subgroups expanded in a species-specific manner from common ancestral genes that were present prior to the diversification of the monocot and dicot lineages. The classification of TCP protein was further supported by gene structure analysis and conserved protein motif analysis. Gene structure analysis showed that the majority of *MeTCP* genes within the same subgroup exhibited very similar gene structure in terms of exon length and intron number. Furthermore, conserved protein motif analysis indicated that all the MeTCP proteins contained typical TCP domain, except MeTCP20c, which contains a truncated TCP domain. Similar to the gene structure, most subgroups of MeTCPs also exhibited conserved motif composition, with several motifs observed in some MeTCP subgroups, such as motif 4 and 17 for subgroup D, motif 7 for subgroup A (Fig. [2c](#Fig2){ref-type="fig"}). These unique motifs may contribute to the specific function of these subgroup members. In general, the majority of *MeTCP* genes in the same subfamilies are evolutionarily conserved, which supports their close evolutionary relationship and the classification of subgroups.

The expression pattern analysis of *MeTCP* genes helps us to assess their possible functions and provide a solid foundation for future functional studies. Generally, similar to the previous study, *MeTCP* genes exhibited greatly differential expression patterns across a variety of tissues, not only among subgroups but members within the same subgroups (*MeTCP5a* and *MeTCP5b*, *MeTCP13a* and *MeTCP13b*), suggesting that these *MeTCP* genes may function diversely in various tissues (Fig. [3](#Fig3){ref-type="fig"}). On the contrary, some *MeTCP* genes with extremely high sequence identity (*MeTCP2a* and *MeTCP2b*, M*eTCP18c* and *MeTCP18d*) showed conserved expression patterns (Fig. [3](#Fig3){ref-type="fig"}), implying they may play a redundant role in regulating plant growth.

To date, the role of plant hormones in regulating plant growth, development, and abiotic stress responses by modulating gene expression is well established^[@CR52],\ [@CR55]^. To our knowledge, although the relationship between TCP proteins and hormones has been widely known, such as cytokinins, JA and GA, the dynamic and spatially expression patterns of *MeTCP* genes response to various hormones was still obscure. Our current results revealed that the majority of *MeTCP* genes detected here displayed distinct changes under different hormone treatments. To expect, most of *MeTCP* were up-regulated or down-regulated by cytokinins (ZT and 6-BA). *MeTCP15b* and *MeTCP23b*, orthologs of *AtTCP15* and *AtTCP23*, respectively, had high expression levels under ZT and 6-BA treatments. In *Arabidopsis*, cytokinin treatments induce *TCP15* transcription and promote *TCP15* (and *TCP14*) protein activation by post-translational modification, which in turn promote cytokinin responses^[@CR56]^. *MeTCP4*, *MeTCP20b*, and *MeTCP20e*, orthologs of *AtTCP4* and *AtTCP20*, respectively, had altered expression patterns under JA. In young leaves of *Arabidopsis*, *AtTCP20* repressed the transcription of *LIPOXYGENASE2* (*AtLOX2*) gene, which is involved in JA synthesis and promotes leaf senescence, while this negative control is antagonized by *AtTCP4* as the leaf matures^[@CR19]^. It is noteworthy that 22 *MeTCPs* were response to ABA, indicating these genes might function as key mediators of stress responses through ABA signaling pathways. Taken together, these results suggested that MeTCPs play potential regulatory roles by modulating phytohormone signaling in plant development or in the responses to stresses. Therefore, it will be particular important to further investigate the potential function of *MeTCP* genes in hormone signaling in the future.

Previous reports have shown that plant *TCP* genes are involved in plant growth and development, as well as abiotic stress responses under normal and stressed growth conditions^[@CR22]^. In cassava, our data showed that over half of the *MeTCP* genes were significantly upregulated under cold and drought condition. Among them, 4 genes were up-regulated, and 8 genes were down-regulated by both cold and drought stress. Moreover, combined analysis of expression correlation and promoter content has revealed that most of these *MeTCP* genes exhibited differential expression in response to more than one stress treatments, suggesting the wide involvement of *MeTCP* genes in environmental adaptation. Previous reports have shown that knockdown of miR319-dependent *TCPs* (by constitutive miR319 overexpression) increases drought and salinity stress tolerance in bentgrass^[@CR20]^. Our data showed that *MeTCP3a* and *MeTCP4*, targets of miRNA319, had altered expression patterns under cold, drought and salt stress, suggesting that these genes might play important roles under abiotic stress conditions in cassava. These data indicated that MeTCP might function in resistance to abiotic stresses in cassava.

In conclusion, we identified 36 *TCP* transcription factor genes from cassava. Phylogenetic analysis of cassava, *Arabidopsis*, and rice indicated that these *MeTCP* genes could be divided into 8 groups, which is supported by further conserved protein motif, and gene structure analyses. Although nearly all the *MeTCP* genes were expressed in the examined tissues, some genes were up-regulated in one or several specific organs. mRNA accumulation was altered by a variety of hormone treatments (ABA, IAA, GA, JA, ZT and 6-BA), environmental conditions (drought, high salinity, and low temperature). These results suggested that MeTCP family proteins play critical roles in maintaining cassava normal growth under normal or stress conditions through complicated mechanisms. Thus, additional studies on the detailed functions of each gene are warranted in cassava.

Materials and Methods {#Sec12}
=====================

Identification and bioinformatics analysis of candidate genes {#Sec13}
-------------------------------------------------------------

To identify potential members of the cassava TCP protein family, The *Arabidopsis* TCP protein sequences were used as seed queries in BLASTp searches against the cassava database (Phytozome: <http://www.phytozome.net/cassava.php>)^[@CR57]^. The TCP domain (PF03634, Pfam; <http://pfam.sanger.ac.uk/>)^[@CR46]^ was also employed as query to perform a blast search against the same genome database. The identified MeTCP proteins were renamed as MeTCP2 to MeTCP23 according to the *Arabidopsis* TCP proteins with highest sequence similarity. Information on *MeTCP* genes, including exons and introns number, open reading frame (ORF) and amino acid (AA) lengths, was obtained from Phytozome database. The molecular weight, theoretical isoelectric point (PI) and grand average of hydropathy (GRAVY) of the MeTCP proteins were investigated using ExPASy online tools (<http://web.expasy.org/protparam/>).

Analysis of phylogenetic relationships and gene structure {#Sec14}
---------------------------------------------------------

Multiple sequence alignments were applied to confirm the conserved domains of predicted MeTCP proteins. The Clustal × 2.0^[@CR58]^ was employed to align the full-length MeTCP proteins from cassava, *Arabidopsis* and rice. Then, the bootstrap neighbor-joining evolutionary tree was created by MEGA 6.0 software^[@CR59]^ with 1000 bootstrap replicates based on the sequence alignments. The exon-intron organization of *MeTCP* genes was determined by comparing the coding DNA sequence (CDS) with its corresponding genomic sequences using the Gene Structure Display Server (GSDS) software (<http://gsds.cbi.pku.edu.cn/>)^[@CR60]^.

Identification of conserved motif of MeTCP proteins and promoters {#Sec15}
-----------------------------------------------------------------

By using the Multiple Expectation maximization for Motif Elicitation (MEME) program (<http://meme.nbcr.net/meme/cgi-bin/meme.cgi>), the conserved motifs in full-length cassava MeTCP protein sequences were identified with the following parameters: maximum number of motifs was 20 and the optimum width of motifs was set between 10 and 50^[@CR49]^. The identified protein motifs were further annotated with ScanProsite^[@CR50]^. *MeTCP* promoter sequences in cassava were submitted to online MEME program for identification of conserved motifs. The optimized MEME parameters were as follows: any number of repetitions and maximum number of motifs-15.

Plant materials and hormone/stress treatment {#Sec16}
--------------------------------------------

Cassava (*Manihot esculenta*) cultivar (TMS60444) was used in the present study. Segments cut from cassava stems were inserted into MS plates in a greenhouse at 26 ± 2 °C, with a photoperiod of 16 h light and 8 h dark. All hormone and environmental treatments were conducted when uniform-sized seedlings developed two fully opened trifoliate leaves (approximately two weeks after sowing). For hormone treatment, 14-day-old cassava seedlings were soaked in liquid MS medium with 100 μM indole acetic acid(IAA), 100 μM gibberellin (GA3), 100 μM Methyl jasmonate (MeJA), 100 μM abscisic acid (ABA),100 μM zeatin (ZT) and 100 μM, 6-benzylaminopurine (6-BA) for 3 h, respectively, and then the young leaves and shoot apex from at least ten separate seedlings/plants were harvested. Seedlings soaked in liquid MS medium without any hormone were used as control. For cold treatment, seedlings were placed at 4 °C for 24 h, and then the young leaves and shoot apex were collected for RNA isolation. For drought and salt stress treatment, cassava seedlings were treated with 20% PEG6000 and 100 mM NaCl, and harvested at 6 h after treatment, respectively. In all cases, parallel and untreated plants at the same stage were used as controls. All samples harvested were flash-frozen in liquid nitrogen, and stored at −80 °C until RNA isolation.

RNA isolation and expression analysis {#Sec17}
-------------------------------------

Total RNA was extracted from 0.1 g of tissue by using Plant RNA kit (OMEGA), following the manufacturer's instructions. Reverse transcription reactions were performed using 5 μg of RNA with PrimeScript RT reagent kit with gDNA Eraser (TIANGEN, Beijing, China). Quantitative reverse transcription PCR (qRT-PCR) was performed as described elsewhere^[@CR61]^, the PCR conditions were as follows: pre-incubation at 94 °C for 5 min, followed by 40 cycles at 94 °C for 10 s, 60 °C for 10 s, 72 °C for 30 s. After amplication was complete, a melting curve was obtained by holding at 95 °C for 5 s and then at 65 °C for 15 s, followed by heating slowly at 0.1 °C/s to 95 °C. Real-time PCR was performed with a Bio-Rad real-time thermal cycling system using SYBR® Premix Ex Taq™ II (TaKaRa, Japan) to assess gene expression levels. The relative expression levels of each gene were calculated by the 2^-ΔΔCt^ method. The cassava actin gene was used as internal control for normalization. The primers used are listed in Supplemental Table [S3](#MOESM2){ref-type="media"}. *MeTCP* genes that were up- or down-regulated by at least two-fold were considered as differentially expressed. The experiments were performed in triplicate.

Expression correlation analysis of level between cassava and *Arabidopsis* {#Sec18}
--------------------------------------------------------------------------

The expression levels of *AtTCPs* in different tissues were extracted from the microarray data^[@CR62]^. To compare the microarray gcRMA value of AtTCPs with qPCR relative value of MeTCPs, all values were normalized to z-score, and pearson correlation coefficient were performed for each homologous gene pairs.

Transcriptome analysis {#Sec19}
----------------------

Cassava shoot apices and youngest leaves of TMS60444 under normal conditions, cold and drought treatments were used to isolated total RNA for transcriptome analysis. As previously described^[@CR61]^, the total RNA isolation, whole transcriptome libraries preparation and deep sequencing were performed by the Annoroad Gene Technology Corporation (Beijing, PR China). A total of 140 gigabase in-depth sequencing of library was performed initially on a HiSeq. 2500 instrument that generated paired-end reads with 125 nucleotides. Data analysis was carried out by previously described^[@CR61]^. The generated transcriptomic data has been submitted to Sequence Read Archive (SRA) in NCBI with the accession number SRP101302.

Electronic supplementary material
=================================

 {#Sec20}

Supplementary Figures Supplementary Table

Ning Lei, Xiang Yu and Shuxia Li contributed equally to this work.

**Electronic supplementary material**

**Supplementary information** accompanies this paper at doi:10.1038/s41598-017-09398-5

**Publisher\'s note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This research was supported by the Program of Hainan Association for Science and Technology Plans to Youth R & D Innovation (HAST201627), the National Natural Science Foundation of China (31561143012), Central Public-interest Scientific Institution Basal Research Fund for Chinese Academy of Tropical Agricultural Sciences (1630052017022) and the Hainan Province Innovative Research Team Foundation (2016CXTD013).

S.X.L. and M.P. devised the study, N.L.; X.Y; S.X.L. and C.Y.Z. conducted the experiments and analyses, all authors contributed to data interpretation and writing of the manuscript.

Competing Interests {#FPar1}
===================

The authors declare that they have no competing interests.
