Comparative analysis of the TCP gene family in celery, coriander and carrot (family Apiaceae)

The Apiaceae family of plants includes more than 400 genera and 3,000 species^[1,2]. Several Apiaceae species such as carrot (Daucus carota), celery (Apium graveolens), and coriander (Coriandrum sativum) are cultivated as a vegetable or for medicinal purposes worldwide^[3,4].

Celery is an annual or biennial herbage species originating from the Mediterranean and Middle East^[5,6]. Besides being a vegetable, celery is also used as a medicinal plant^[7]. Coriander, which is also known as cilantro, is a popular herb and a major ingredient of curry powder^[8]. Carrot is an economically important vegetable with a high nutritional value^[9]. All three of these plants are diploid species although their chromosome number and genome size differ: celery and coriander each have 22 chromosomes (2n = 2x = 22) whereas carrot has 18 (2n = 2x = 18)^[3,10,11], and the assembled genome size of celery is 3.33 Gb, which is larger than coriander (2.11 Gb) and carrot (421.5 Mb)^[3,11,12].

The TCP gene family in plants is named after the first identified members—ie, TEOSINTE BRANCHED 1 (TB1) in maize (Zea mays)^[13], CYCLOIDEA (CYC) in snapdragon (Antirrhinum majus)^[14], and PROLIFERATING CELL FACTOR 1 (PCF1) and PCF2 in rice (Oryza sativa)^[15]. TCP genes regulate multiple processes in plant growth and development such as shoot branching^[16], seed germination^[17,18], gametophyte development^[19,20], leaf development^[21−25], leaf senescence^[26−28], mitochondrial biogenesis^[29], flower development^[30−32] and cell cycle^[29,33]. There are 24 TCP genes in Arabidopsis thaliana^[34], 22 in rice^[34], 24 in tomato^[35], 19 in plum^[26], 42 in switchgrass^[36], 36 in carrot^[37], and 32 in celery^[38]. However, to date, no TCP genes have been identified in coriander. High-quality genome sequences of celery, carrot and coriander were recently released^[3,11,12], which can facilitate comparative analyses of specific genes in Apiaceae.

In this study, we identified and characterized TCP genes in celery, coriander and carrot and performed comparisons with genes in Lactuca sativa (lettuce), Vitis vinifera (grape), and Arabidopsis. We mapped the TCP genes to chromosomes, identified orthologs and paralogs, detected collinearity and gene expansion or loss, and analyzed their expression patterns in plant tissues. Our results provide a basis for comparative studies on the function and evolution of TCP genes in plants.

The genome sequences of coriander and celery were downloaded from the coriander genome database (http://cgdb.bio2db.com) and celery genome database (http://celerydb.bio2db.com), respectively^[3,10]. The Arabidopsis genome sequence was retrieved from The Arabidopsis Information Resource (TAIR10; http://www.arabidopsis.org). The sequences of carrot (version 2), lettuce (version 5), and grape (Genoscope.12X) were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html)^[11,39,40]. The Pfam database (https://pfam.sanger.ac.uk) was used to identify TCP family genes with the identifier PF03634 (E value < 1e−4)^[41]. The Simple Modular Architecture Research Tool (SMART v9.0) database and Conserved Domains Database (CDD) were used for domain validation^[42,43].

Phylogenetic analysis of TCP gene family

To analyze the evolutionary relationships of TCP genes in Apiaceae, multiple sequence alignment of the TCP amino acid sequences of Arabidopsis, grape, lettuce, carrot, coriander and celery was performed with ClustalW software (v2.0)^[44] and a phylogenetic tree was constructed with the neighborhood-joining method (bootstrap = 1,000) using MEGA X^[45]. The reconstructed TCP gene tree was compared to the actual species tree using Notung v2.9 software with default parameters^[46,47].

Chromosome location, gene structure and conserved motif analysis

The location of coriander, celery and carrot TCP genes on chromosomes was drawn using Tbtools software and the files were saved in general feature format (gff)^[48]. Gene structure was determined using Gene Structure Display Server 2.0 (https://gsds.cbi.pku.edu.cn)^[49]. Conserved motifs were analyzed using Multiple Expectation maximizations for Motif Elicitation suite 5.2.0 (http://meme-suite.org)^[50].

Analysis of orthologs and paralogs

Orthologous and paralogous TCP gene pairs in celery, coriander and carrot were analyzed using OrthoMCL software v2.0 (https://orthomcl.org/orthomcl)^[51]. The relationships of the genes among the three species was depicted using Circos software (v0.69)^[52].

Identification of collinear blocks and duplication types

MCScanX was used to identify collinear blocks and duplication types of the TCP genes^[53]. Whole-genome sequences were searched against themselves using BLASTp (E value < 1e−5). We extracted TCP genes located in the collinear blocks using Perl scripts. The duplication type of TCP genes was estimated using the subprogram duplicate_gene_classifier.

Calculation of Ka/Ks and estimation of divergence time

The nonsynonymous rate (Ka), synonymous rate (Ks), their ratio (Ka/Ks) and divergence time among orthologous gene pairs of the three species were calculated using KaKs_Calculator 2.0^[54]. The coding sequence of orthologous gene pairs were aligned using ClustalW (v2.0)^[44]; AXTconvertor software (v1.0) was then used to convert the alignment file to axt format. Lastly, the Ka value, Ks value and their ratio were calculated based on the Nei–Gojobori method^[54]. Ks was used to estimate the divergence time based on the formula T = Ks/2r, where r indicates neutral substitutions (5.2 × 10⁻⁹ for Apiaceae)^[12].

Analysis of selective pressures

We used the maximum likelihood method and codon substitution models to determine the likelihood ratio of positive selection. We analyzed each branch of the phylogenetic tree to infer ω (the ratio of nonsynonymous to synonymous distances) using CodeML implemented in PAML4.9^[47,55]. We adopted a complete deletion method for analyzing alignments with gaps and eliminated sequences with gaps in over 40% of their length. The likelihood ratio test between M0 and M1 and between M7 and M8 models were used to determine variation sites.

Identification of micro (mi)RNA target genes in the TCP gene family

Mature miRNA sequences of A. thaliana were downloaded from miRBase (release 22.1; http://www.mirbase.org)^[56]. TCP genes that are miRNA targets were predicted using psRNATarget Schema v2 (2017 release)^[57] with maximum expectation ≤ 3 and other default parameters. A miRNA–TCP gene network was constructed using Cytoscape v3.7.2 software^[58].

Expression of TCP genes

TCP gene expression data in celery and coriander were extracted from RNA sequencing (RNA-seq) datasets previously published by our group^[3,10] using Perl script; the values were normalized as reads per kilobase per million reads (RKPM). An expression heatmap was created using TBtools software (v1.0)^[48].

We identified 29 TCP genes in the genome of celery, 43 in coriander, and 50 in carrot (Supplemental Table S1 and S2). Additionally, 24 TCP genes were identified in Arabidopsis along with 20 in grape and 49 in lettuce (Supplemental Table S2). Thus, a total of 215 TCP genes were identified in the six species for further analysis.

Phylogenetic and functional analyses of TCP gene family

To classify the TCP gene family in plants, we constructed a phylogenetic tree of all 215 amino acid sequences from the six abovementioned species using MEGA X (Fig. 1). Consistent with the phylogenetic relationships described in Arabidopsis and grape, the phylogenetic tree had three groups according to the type of TCP protein domain including the PCF, CINCINNATA (CIN), and CYC/TB1 classes.

Figure 1.  Phylogenetic tree of TCP family genes in three Apiaceae species (carrot, celery and coriander) and lettuce, grape and Arabidopsis. The topology of the phylogenetic tree was determined using IQ-TREE with maximum likelihood (ML) based on the JTT+F+R8 model. The bootstrap was set to 1,000 replicates, and values > 40% are shown. The three classes were identified based on bootstrap values and phylogenetic topology.

In class PCF, there were ten AgTCP, 26 CsTCP, and 36 DcTCP genes; in class CYC/TB1, there were eight AgTCP, seven CsTCP, and seven DcTCP genes; and in class CIN, there were 11 AgTCP, ten CsTCP and seven DcTCP genes (Fig. 1 and Supplemental Table S2). Notably, in coriander and carrot there were more TCP genes in class PCF than in the other two classes.

The functions of most TCP family genes have been well studied in the model plant Arabidopsis. We inferred the function of homologous genes within the same taxonomic group in the phylogenetic tree in order to clarify the function of TCP genes in Apiaceae. For example, AT1G53230.1 (AtTCP3) is known to suppress the expression of CUP-SHAPED COTYLEDON (CUC), resulting in cotyledon fusion^[24]. We identified three TCP genes—namely, DcTCP34, AgTCP2, and CsTCP6—that clustered together with AtTCP3 (Fig. 1), suggesting that they are also related to cotyledon fusion. It may also be possible to deduce the function of other Apiaceae TCP genes based on the function of the homologous genes in Arabidopsis.

TCP gene family structure and conserved motifs

We carried out a gene structure analysis of TCP family genes to identify exons, introns and untranslated regions (Supplemental Fig. S1). Of the 122 TCP genes in Apiaceae, 94—including all 50 DcTCP genes—lacked introns. Most TCP genes had a single exon, although there were some exceptions. For example, CsTCP15 (class CIN) had four exons and DcTCP11 (class PCF) had three. In general, genes in the same class or subclass had similar gene structure and size.

As gene structure varied among genes, we performed a motif analysis to examine the structure in greater detail. We compared five motifs in Apiaceae species and found that motif 3 was found at the beginning of most genes, followed by motif 1 and motif 2 (Fig. 2). However, motif 2 was located at the start of the CsTCP7, DcTCP15, CsTCP23, CsTCP36, and CsTCP2 and motif 5 was present at the beginning of DcTCP11. Almost all TCP genes had motif 1 except for DcTCP and CsTCP in class PCF, indicating that this motif is highly conserved and plays an important role in Apiaceae. Motif 3 was also present in most TCP genes and is likely conserved in Apiaceae.

Figure 2.  Conversed motifs in TCP family genes of three Apiaceae species.

Most genes in classes CYC/TB1 and CIN lacked motif 2 except for DcTCP17 and CsTCP79. Only five genes in class CIN had motif 4, which was present in all class CYC/TB1 genes. Motif 5 was only found in class PCF and not in other classes. Interestingly, DcTCP14 (class PCF) did not have any of the 5 motifs, suggesting that they were lost during the course of the evolution of carrot. Thus, genes in the same class had similar motif composition, indicating that they are functionally similar.

Chromosomal distribution of TCP genes in Apiaceae species

In celery, 27/29 TCP genes mapped to nine chromosomes (Fig. 3a and Supplemental Table S1). Two celery TCP genes did not map to any chromosomes, and no TCP family gene was found on chromosomes Agr6 and Agr10. There were five genes that mapped to chromosomes Agr3, Agr5, and Agr11 while only one was located on chromosomes Agr1 and Agr2.

Figure 3.  Distribution of TCP transcription factors on each chromosome in three Apiaceae species. (a) Celery, (b) Coriander, (c) Carrot.

In coriander, 35/43 TCP genes mapped to nine chromosomes (Fig. 3b and Supplemental Table S1); none were found on chromosome Csa8. Chromosome Csa9 had the most TCP genes (10), followed by chromosomes Csa7 (6) and Csa1 (5). Chromosomes Csa2 and Csa3 each had just one TCP gene.

The 50 carrot TCP genes were unevenly distributed across nine chromosomes (Dca1–9) (Fig. 3c and Supplemental Table S1). Interestingly, TCP gene expansion was observed on chromosome 1, which had 16 genes. Additionally, 15 genes (DcTCP1–15) were clustered together, mainly through tandem and proximal duplication.

Identification of orthologous and paralogous TCP genes in Apiaceae

We examined orthologous and paralogous gene pairs in Apiaceae and found that there were 22 orthologous gene pairs between any two of celery, coriander and carrot (Fig. 4 and Supplemental Table S3), indicating a close phylogenetic relationship between these species. There were only three paralogous gene pairs in celery and coriander (Supplemental Fig. S2 and Supplemental Table S4) but 29 were found in carrot.

Figure 4.  Circle plot of orthologous TCP gene pairs among three Apiaceae species.

The Ks value was calculated to estimate the divergence time of orthologous TCP gene pairs among celery, coriander and carrot (Fig. 5 and Supplemental Table S5). The divergence time ranged from 14.03 to 116.15 million years between celery and coriander TCP genes, 24.94 to 81.44 million years between celery and carrot, and 22.59 to 94.72 million years between coriander and carrot orthologous TCP gene pairs. Therefore, the divergence time of most TCP genes was earlier than that of any two species (celery vs coriander, 11–13 Mya; carrot vs celery or coriander, 20–23 Mya)^[4].

Figure 5.  Ks and divergence time of TCP orthologs. Ks values (a) and divergence time estimation (b) of orthologous gene pairs between any two of three Apiaceae species.

Detection of duplication type for TCP family genes in Apiaceae

Various types of gene duplication can lead to the expansion of a gene family. We examined five types of gene duplication in celery, coriander and carrot (Fig. 6a, Table 1 and Supplemental Table S6)—namely, singleton, dispersed, proximal, tandem and whole-genome duplication (WGD). There were no singleton TCP gene in the three species of Apiaceae (Table 1). Dispersed and tandem duplications were the predominant types in celery and coriander. In celery, the percentage of genes showing dispersed duplication and WGD was 48.1% and 40.7%, respectively; in coriander, the percentages were 40.0% and 45.7%, respectively. WGD was also the predominant type in carrot (32.0%), which had a lower percentage of dispersed duplication (24.0%) than celery and coriander. Moreover, the percentage of the tandem type was higher in carrot (30.0%) than in celery (11.1%) and coriander (11.4%). These results demonstrate that WGD played an important role in TCP gene expansion in celery, coriander and carrot, which is supported by the previous suggestion that they underwent two WGD events since their divergence from lettuce^[3].

Figure 6.  Duplication and loss of TCP family genes. (a) Percentage of duplication types for TCP family genes and genes in the whole genome of three Apiaceae species. (b) Duplication (+) or loss (−) of TCP family genes in three Apiaceae species and three other representative species. Numbers after '+' and '−' are the number of genes. Agr, Apium graveolens (celery); Ath, Arabidopsis thaliana (Arabidopsis); Csa, Coriandrum sativum (coriander); Dca, Daucus carota (carrot); Lsa, Lactuca sativa (lettuce); Vvi, Vitis vinifera (grape).

TCP gene loss and duplication during the evolution of Apiaceae

We compared species and gene trees in celery, coriander and carrot to identify gene losses and duplications during the evolution of the TCP gene family. There were more gene losses in celery (20) than in coriander (7) and carrot (17) but more gene duplications in carrot (18) than in coriander (4) and celery (3). In the common ancestor of coriander, celery and carrot, there were 20 gene duplications and 14 gene losses (Fig. 6b).

Positive selection of TCP family genes in Apiaceae

We analyzed natural selection in the evolution of TCP genes in Apiaceae (Fig. 7). Strong positive selection was observed at the major nodes of the phylogenetic tree, which may have contributed to the functional divergence of Apiaceae species. We detected 35 positive selection sites overall; most were in class PCF (29), followed by class CYC/TB1 (5) and class CIN (1), indicating that TCP genes in class PCF underwent greater positive selection in the evolution of Apiaceae.

Figure 7.  Positive selection of TCP family genes in celery, coriander and carrot. Red stars represent branches in which positive selection occurred. The maximum likelihood (ML) phylogenetic tree was constructed using PhyML software.

miRNA target TCP genes in Apiaceae

We next sought to identify TCP family genes in Apiaceae that are regulated by miRNAs. We found nine miRNAs that regulated 20 TCP family genes including five genes in celery, nine in coriander and six in carrot (Supplemental Table S7 and Fig. 8). Of the nine miRNAs, miR-319 had the most target genes (11), followed by miR-172 (3) and miR-181 (3) (Fig. 8). Specifically, miR-319 regulated four TCP family genes in celery, four in coriander, and three in carrot. Our results are supported by other studies demonstrating that miR-319 regulates TCP family genes^[31,59−61]. We also found that four genes were regulated by more than one miRNA: DcTCP20 was regulated by miR-319 and miR-837, CsTCP33 was regulated by miR-319 and miR-8181, and AgTCP26 and AgTCP27 were regulated by miR-319 and miR-8181 (Supplemental Table S7 and Fig. 8).

Figure 8.  Interaction network of miRNAs and target TCP family genes in carrot, celery and coriander.

Expression of TCP family genes in celery and coriander

We analyzed the expression patterns of TCP genes in root, petiole and leaf tissues of celery and coriander. In celery, all TCP family genes in class CYC/TB1 had relatively low expression in the three tissues (Fig. 9a and Supplemental Table S8) whereas AgTCP22 had the highest expression level (RPKM > 60), suggesting a key role in celery growth and development. In coriander and celery, all TCP family genes in class CYC/TB1 were expressed at a relatively low level in the three tissues (Fig. 9b and Supplemental Table S8); meanwhile, several genes including CsTCP39, CsTCP14, CsTCP12, CsTCP31, and CsTCP40 had high expression. CsTCP6 and CsTCP30 were more highly expressed in leaf than in the other two tissues.

Figure 9.  Expression of TCP family genes in three replicates of plant tissues including root (R1, R2 and R3), petiole (P1, P2 and P3), and leaf (L1, L2 and L3). Hierarchical gene expression clustering of TCP genes in celery (a) and coriander (b). Expression levels were calculated based on RPKM.

Celery, coriander and carrot are typical members of the Apiaceae family. The draft genomes of these three species were recently released^[3,11,12] and there have been several studies on TCP family genes in carrot and celery based on the sequences^{[37,38,62,63]}. The latest versions of the celery, coriander and carrot genomes are of high quality with chromosomal-level assembly, allowing us to accurately and comprehensively analyze the TCP gene family in Apiaceae. To date there have been no reports on TCP genes in coriander. In this study, we identified 43 TCP genes in coriander as well as 29 in celery and 50 in carrot. Our results provide a resource for future studies on the TCP gene family in Apiaceae or related species.

TCP transcription factors have a 59-amino acid basic helix-loop-helix (bHLH) motif that is involved in DNA binding and protein–protein interaction^[64]. The bHLH-like domain of TCP differs from the canonical bHLH in its basic region^[34,65]. PCF1 and PCF2 interact with DNA-binding proteins that specifically bind to the PROLIFERATING CELL NUCLEAR ANTIGEN (PCNA) promoter^[15]. Our analyses of gene structure and motifs of TCP genes in Apiaceae revealed similarities within the same class or subclass.

Gene duplication is the main mechanism underlying the evolution of complex phenotypes^[66]. Many duplicated genes in plants were produced by WGD or whole-genome triplication^[67−70]. Most duplicated genes were functionally redundant and had one of four fates during the course of evolution namely: conservation, neofunctionalization, subfunctionalization and specialization^[67,71]. In conservation, the ancestral function was maintained in both copies, thus preserving gene dosage^[71]. In neofunctionalization, one copy retained the original function while the other acquired a novel function^[71,72]. In subfunctionalization, both copies were required to preserve the ancestral gene function^[71,73,74]. In specialization, subfunctionalization and neofunctionalization acted cooperatively, producing two gene copies that were functionally distinct from each other and from the ancestral gene^[71,75]. Functional redundancies of TCP genes have been reported in Arabidopsis^[76]. In carrot, 15 TCP genes were clustered on chromosome 1, and the number of paralogous gene pairs was greater in carrot (29) than in celery (3) and coriander (3). Although there were more gene losses than duplications in the evolution of celery, coriander and carrot, we found that WGD made a major contribution to TCP gene family expansion in Apiaceae, which is similar to what has been reported in most other gene families in higher plants^[47,77−80].

The broad range of functions of TCP family genes in plants can be attributed to the diverse structures of different members. Most TCP genes are highly expressed in meristematic tissues, suggesting that their main function is to promote plant proliferation and growth^[81]. However, some TCP genes, such as CIN and CYC/TB1, are known to negatively regulate plant proliferation and development^[82] (lateral organ development for CIN genes and flower and lateral shoot development for CYC/TB1 genes)^[26]. Our gene expression analysis showed that TCP gene expression in celery (AgTCP22) and coriander (CsTCP12) was nearly 2x higher in root and petiole than in leaf, suggesting roles in plant growth and development.

In conclusion, we identified and characterized TCP genes in three Apiaceae species. We described their chromosomal location, exon–intron structure, motifs, collinearity, positive selection and expression patterns in plant tissues. These results provide a basis for investigations on the molecular networks regulating growth and development in Apiaceae and other plants.