Phylogeny, Functional Annotation, and Protein Interaction Network Analyses of the Xenopus tropicalis Basic Helix-Loop-Helix Transcription Factors

The previous survey identified 70 basic helix-loop-helix (bHLH) proteins, but it was proved to be incomplete, and the functional information and regulatory networks of frog bHLH transcription factors were not fully known. Therefore, we conducted an updated genome-wide survey in the Xenopus tropicalis genome project databases and identified 105 bHLH sequences. Among the retrieved 105 sequences, phylogenetic analyses revealed that 103 bHLH proteins belonged to 43 families or subfamilies with 46, 26, 11, 3, 15, and 4 members in the corresponding supergroups. Next, gene ontology (GO) enrichment analyses showed 65 significant GO annotations of biological processes and molecular functions and KEGG pathways counted in frequency. To explore the functional pathways, regulatory gene networks, and/or related gene groups coding for Xenopus tropicalis bHLH proteins, the identified bHLH genes were put into the databases KOBAS and STRING to get the signaling information of pathways and protein interaction networks according to available public databases and known protein interactions. From the genome annotation and pathway analysis using KOBAS, we identified 16 pathways in the Xenopus tropicalis genome. From the STRING interaction analysis, 68 hub proteins were identified, and many hub proteins created a tight network or a functional module within the protein families.


Introduction
Transcription factors are usually classified into different families based on their sequence of functional DNA-binding or protein-binding domains, which are highly conserved among many species and include many members mediating cell fate allocation during animal and plant development [1][2][3][4][5][6][7][8][9][10][11]. The expression and activity of basic helix-loop-helix (bHLH) transcription factors can be regulated in response to cellcell signaling, leading to the transcription of specific sets of genes required for a cell to adopt particular fates. Due to their important functions found in various organisms, bHLH transcription factors have been the subject of many researches. The first report of bHLH transcription factors focused on the murine factors E12 and E47 [12]. Later, more and more bHLH proteins have been identified in living organisms. In 1997, Atchley and Fitch [1] proposed an organization for the classification of the bHLH proteins based on the phylogenetic analysis of the 122 bHLH domains combined with the presence or absence of another additional domain. Their analysis allowed for the defining of four different groups of bHLH protein families according to structural similarities [1]. This classification was performed using only the bHLH motif or domain, because the flanking regions for bHLH proteins are very divergent. Atchley and Fitch's classification led to the postulation of four distinct groups based on amino acid patterns and E-box-binding specificity [1]. In 2002, Ledent et al. [4] defined 44 orthologous families or sub-families and 6 supergroups based on the DNA-binding activities of bHLH transcription factors after large-scale phylogenetic analyses. After the revision of Simionato et al. [6] in 2007, animal bHLH proteins are reclassified into 45 families. Among these 6 supergroups, members of groups A and B are common bHLH proteins [1,[3][4][5][6]. Group A proteins bind to CACCTG or CAGCTG, while group B proteins bind to CACGTG or CATGTTG. The consensus DNA binding sequences for these bHLH proteins form the typical E boxes (CANNTG). Group C proteins are complex molecules with one or two PAS domains following the bHLH domain, being inclined to bind the core sequence ACGTG or GCGTG. They are mainly responsible for the regulation of midline and tracheal development, circadian rhythms, and gene transcription in response to environmental toxins. Group D proteins correspond to bHLH proteins that are unable to bind to DNA due to lack of a basic domain. Both, group D and group F, are proteins that lack basic parts and act as antagonist partners of group A proteins in the heterodimers. Particularly, group F are a kind of COE proteins characterized by the presence of an additional COE domain involved in both dimerization and DNA binding. Group E proteins are another type of special transcription factors. They usually contain two additional domains named "Orange" and "WRPW" peptides in their carboxyl termini and they bind preferentially to sequences typical of N boxes (CACGCG or CACGAG). Generally speaking, all of the bHLH transcription factors share a common bHLH motif or domain of approximately 60 amino acids, which contains a basic region and two helices separated by a loop (HLH) region of variable length [3][4][5]12]. The basic region is a DNA-binding domain, and the amphipathic -helices of two bHLH proteins can interact with each other. The HLH domain promotes dimerization, and interaction between the helix regions of two different bHLH proteins leads to the formation of homodimeric or heterodimeric complexes, while the basic region of each partner recognizes and binds to a core hexanucleotide DNA sequence [2][3][4]. In a couple of reports [13,14], Atchley et al. inferred a predictive motif for the bHLH domains based on 242 bHLH proteins, in which 19 conserved sites were found within the bHLH domains. It was found and proved that a sequence with no more than 9 mismatches could be a putative bHLH protein [15].
Recently, in many organisms whose genomes have been released and are available, more and more bHLH proteins have been identified and bHLH transcription factor families have been analyzed due to their important and pivotal regulatory functions displayed in various organisms . As well as Xenopus laevis the Xenopus tropicalis is a model organism for researches testing the developmental, behavioral, and neurological consequences of genetic variation [26][27][28]. The draft of Xenopus tropicalis genome assembly was submitted by American scientists at the Lawrence Berkeley National Laboratory in California [28], and the Xenopus tropicalis genome project is still underway. In previous work, the preliminary survey identified 70 bHLH transcription factors [16]. Recently, we found it was incomplete and the functional properties and regulatory networks of bHLH transcription factors were not fully analyzed. In this study, we used the criteria developed by Atchley et al. [13] and the 45 representative bHLH domains defined by Ledent et al. [4] and Simionato et al. [6] to do updated searches using BLAST search algorithms in the Xenopus tropicalis genomic database and identified 105 bHLH proteins. We next made large-scale phylogenetic analyses of the Xenopus tropicalis bHLH domains with the 118 human bHLH domains [6]; this allowed us to define the full set of bHLH orthologous genes and their related families. We further report the result of analyses of gene ontology (GO) annotations, functional pathways, and protein interaction networks based on the Xenopus tropicalis genomic databases.

BLAST Searches and Retrieval of bHLH Domains.
At first, we followed the criteria developed by Atchley et al. [1,13] to define a bHLH protein [13]. These searches initially yielded a few bHLH transcription factors (up to 20 protein sequences). The deduced predictive protein consensus motif of Atchley et al. [13] is "+ + X (3)(4)(5)(6) (2) A XY X (2) L" where + = K, R; = I, L, V; Φ = F, I, L; = I, V, T; E, R, K, A, and Y are as defined; X = any residue; X ( ) = any residues; and X ( -) = to of any residues. We also used the 45 representative bHLH domains from the tables provided by Ledent et al. [4] and Simionato et al. [6] to make multiple TBLASTN and BLASTP searches of bHLH domains against the Xenopus Genome Resources built by NCBI (http://www.ncbi.nlm.nih.gov/genome/guide/ frog/) and Xenbase (http://www.xenbase.org/) for all putative bHLH proteins. Then, PSI-BLAST searches were conducted against the nonredundant database of Xenopus genomes at NCBI using the representative bHLH domain sequences. All of the TBLASTN, BLASTP, and PSI-BLAST searches were conducted with the methods and similar parameter settingups in the previous works [7,16]. With these BLAST searches above, we obtained all of the putative bHLH proteins with no more than 9 mismatches among the 19 amino acids residues [15]. Moreover, we also did TBLASTN searches of frog EST data against the Xenopus Genome EST databases with a stringency set as ≤ 0.0001 and an identity of 90% or higher as candidacy. The obtained EST data were translated into protein sequences using online analysis tools (http://www .genoscope.cns.fr/agc/tools/) to verify the putative bHLH sequences found.

Manual Improvement and Sequence Alignment.
Protein sequence accession numbers and genomic contig numbers were finally obtained by BLASTP and TBLASTN searches against the Xenopus tropicalis protein databases and genome sequence assembly (reference only) with the amino acid sequence of each identified bHLH domain. All of the obtained sequences were aligned using ClustalX 2.0 [29]. Redundant sequences of candidates were removed according to their corresponding serial numbers of the scaffold or clone or genomic contig, gene ID, protein ID, coding region, and alignment information. The, finaly, aligned bHLH domains were shaded using GeneDoc 2.6.02 [30] and copied into an RTF file for further annotation.

Analyses of Gene Ontology (GO) Annotations and Pathways.
A functional annotation analysis of Xenopus tropicalis bHLH transcription factor genes was conducted. Gene ontology (GO) function enrichment was analyzed using DAVID Functional Annotation Bioinformatics Tools [31,32], which use the ontology hierarchy tree and calculates and report statistical significance for GO term categories with a hypergeometric value and enrichment scores. This approach directly scores predefined gene sets and/or pathways based on given gene lists.
All of the bHLH transcription factor genes were also subjected to KOBAS analysis (http://kobas.cbi.pku.edu.cn/ home.do), and significant pathways were retrieved at the default value ≤ 0.5. We applied KOBAS vocabulary to first annotate all genes with corresponding KO and then identify both, the most frequent the statistically significantly enriched pathways. With rather strict cutoff of FDR ≤ 0.05, KOBAS found statistically significantly enriched pathways, as shown in Table 3.
We could thus identify and select significantly enriched gene ontology terms and pathways using bioinformatics databases DAVID [31,32] and KOBAS [33][34][35], respectively. We selected the functional categories that were more likely to be biologically meaningful by statistical significance of each functional category in the input set of genes versus all annotated genes in the Xenopus tropicalis genome.

Protein Interaction Network Analysis.
To investigate possible interactions between the gene lists from our updated surveys, the STRING search tool was used for the creation of protein interaction network (PIN) files as previously described [36][37][38]. To increase the completeness of our results, this search was set to include full information extracted from the STRING biological interaction databases. The created networks were explored and compared based on their topological characteristics and gene products (proteins) by default with a confidence of score 0.15 [38].

Phylogenetic
Analyses. The putative Xenopus tropicalis bHLH protein sequences, together with the human bHLH domains, were used to construct phylogenetic trees based on bayesian inference (BI) by MRBAYES 3.1.2 [37,38] and maximum likelihood estimation (MLE) by PHYML 2.4.4 [44] with the JTT substitution frequency matrix [45], respectively. Phylogenetic analyses by BI and MLE were performed with the methods and similar parameter setting-ups in the previous works [7,16]. Briefly, the BI analysis was performed with two independent Markov chains, each containing from 800 to 1100 million Monte Carlo steps until the standard deviation of split frequencies was below 0.01, with sample frequency saved every 1000 generations. Finally, all of the obtained trees were edited and displayed by means of the software package MEGA 4.0 [46].

Data Retrieval and Identification of bHLH Transcription
Factors. The names and related information of the putative Xenopus tropicalis bHLH proteins are listed in Table 1. All of the bHLH domains obtained had more than 10 conserved amino acids [15]. The putative bHLH proteins were named according to their phylogenetic relationship with its corresponding human orthologs and paralogs. If a human bHLH sequence had two or more Xenopus tropicalis orthologous genes, we used "a, " "b, " and "c" or "1, " "2, " and "3" and so on, to number them. In the present work, 34 frog hypothetical and/or predicted proteins belonged to novel bHLH members and were reannotated in this study, that is, NP 001096226. In total, 105 putative Xenopus tropicalis bHLH protein sequences were identified with the BLASTP, TBLASTN, and PSI-BLAST searches and manual examination of the 19 conserved amino acid sites (Table 1, Figure 1). Among these putative bHLH protein sequences, most of these hypothetical proteins were newly produced in the Xenopus tropicalis genome project. We further identified and verified these hypothetical proteins with corresponding EST sequences obtained by TBLASTN searches against the expressed sequence database (data not shown).
In summary, two proteins identified belonging to none of these groups were classified as "orphans, " while the other 103 bHLH members belonged to 43 families with 46, 26, 11, 3, 15, and 4 bHLH members in the corresponding high-order groups A, B, C, D, E, and F, respectively. Figure 1 showed the domain alignment of 105 Xenopus tropicalis bHLH proteins. In addition, the members of Delilah and Mist families were not found in this research.

Phylogenetic Analyses and Identification of Putative bHLH
Proteins. Phylogenetic trees of MLE and BI showed the diversity of the frog bHLH transcription factor family. All of the data of phylogenetic trees for Xenopus tropicalis bHLH proteins are available upon request. The topologies of these two inference methods agreed well with each other (Table 1). It was found that both human and frog proteomes have a number of lineage-specific bHLH families and their members. For example, in the Xenopus tropicalis proteomes, no orthologous genes for human TF12, Hath1, Hath4a, Hath4b, Hath5, and Id1 could be found in the present research. However, the Xenopus tropicalis proteomes also have multiple orthologous genes corresponding to one human gene, such as SREBP1a, SREBP1b, and SREBP1c (orthologous genes of human SREBP1); Hes1a and Hes1b (orthologous genes of human Hes1); Hes6a and Hes6b (orthologous genes of human Hes6); Hes5a, Hes5b, Hes5c, Hes5d, Hes5e, and Esr9 (orthologous genes of human Hes5).

Enriched Functional GO Annotations.
Gene ontology (GO) annotations including biological process (BP), molecular function (MF), and cellular component (CC) were downloaded and investigated from the gene ontology database (http://www.geneontology.org/), and the genes were grouped according to their GO hierarchy annotations. To explore functional properties and identify groups of genes coding for proteins with similar function or with participation in common regulatory pathways, all of the retrieved putative bHLH genes were grouped and functionally classified and enriched according to available GO annotations, information from curated pathways, and known protein interactions. In the present work, the 105 frog bHLH genes were grouped into 7 supergroups according to Ledent et al. [4] and Simionato et al. [6] to get available GO annotations and their enrichment by categories (cutoff of ≤ 0.05). With gene accessions, protein accessions, and the other eligible sequence information in DAVID Bioinformatics Database [32] for Xenopus tropicalis bHLH transcription factors, we retrieved all of the significant GO annotations (cutoff of ≤ 0.05). There were 96 genes fitting the record of DAVID Bioinformatics Database [31,32] and these genes obtained significant GO annotations, while the other nine genes did not get significant GO annotation and were discarded (mainly group D, F, and Orphans; Table 2).
Among the genes, more than half were annotated as exhibiting "transcription regulator activity" and/or "regulation of transcription" or similar terms related to DNAdependent regulation of transcription, DNA binding, or regulation of RNA metabolic process in the BP and MF categories. There were three significant KEEG pathways, that is, circadian rhythm (KEGG ID: 480089074, value 0.0039), TGF-beta signaling pathway (KEGG Id: 480089058, value 0.0024), and Notch signaling pathway (KEGG Id: 480089056, value 0.046), and few significant GO terms for bHLH genes identified in the CC category for Xenopus tropicalis bHLH proteins in DAVID Bioinformatics Database [32]. In the BP category, a total of 47.83% of the significant GO annotations were annotated as transcription and transcription (factor) activity and/or regulation of transcription, while 28.26% of GO annotations were connected to muscle cell development or differentiation and 26.09% of GO annotations were related to negative regulation of cellular biosynthetic or macromolecular metabolic processes. Several genes in the BP category were associated with neural tube development, floor plate development, sensory organ development, chordate embryonic development, hormone receptor binding, and so forth. In the MF category, 56.25% of GO annotations were connected to transcription factor binding or transcription regulator activity, while 3 out of 16 of the GO annotations were related to DNA binding.
DNA binding, protein dimerization, and transcription coactivator activity are important functional activities of bHLH domains. The DNA binding activity of bHLH proteins is mainly determined by the basic region [2]. Site-directed mutagenesis experiments and the crystal structure studies of bHLH proteins showed that the Glu-9/Arg-12 pair forms the CANNTG recognition motif, the critical Glu-9 contacts the first CA in the DNA-binding motif, and the role of Arg-12 is to fix and stabilize the position of the Glu-9 [35][36][37][38]47]. To further understand the functions of Xenopus tropicalis bHLH genes as a whole, we collected GO enrichment data on the 105 Xenopus tropicalis bHLH genes with significant hypergeometric values. Among all of the GO terms, 65 significant GO terms ( ≤ 0.05) were identified showing key cellular components, molecular functions, biological processes, and KEGG pathways for the 105 Xenopus tropicalis bHLH genes (Table 2). Muscle organ development, embryonic development ending in birth or egg hatching, chordate embryonic development, sensory organ development, neural tube development, camera-type eye development and eye development, floor plate development, and muscle fiber and tissue development have high frequencies when taking no account of the frequent GO term categories of transcriptional factors such as (negative) regulation of transcription and regulation of metabolism and biosynthetic processes. It has been well known that the bHLH genes in various groups have special recognition motifs of DNA-binding sites such as E-box and G-box. So, how about the gene functions of each group? To explore these issues, we calculated the hypergeometric distribution enrichment score of gene molecular functions from group A to group F based on GO annotations of GO term categories including biological process, molecular function, cellular component, KEGG pathways, and other key words. However, only significant enriched annotations (cut off ≤ 0.05) in deeper layers (sublayers) are shown in Table 2. GO statistics analyzed with a brief summary of subtypes describing each subgroup are also listed in Table 2.
Our analysis focused on significant GO terms for all of the whole Xenopus tropicalis bHLH gene family and for each subgroup ( Table 2). We found that each subgroup (except for D and F with few members identified) of bHLH transcription factors has its own specific GO term categories (Table 2), when common GO terms of transcription such as transcription regulator activity, regulation of transcription, and DNA binding and protein dimerization activity are discounted. Group A is characterized with muscle organ development such as (striated) muscle cell differentiation and development, (skeletal) muscle fiber development, (extraocular) skeletal muscle tissue development, and striated muscle and pharyngeal muscle development. In addition, digestive system development, pharynx development, and sensory organ development are also included in this group ( Table 2). The functions of bHLH members of group B and group C are mainly composed of transcription, transcription regulator activity, and regulation of transcription. However, group B is different from group C with some GO terms such as transcription coactivator activity, transcription cofactor activity, and (nuclear) hormone receptor binding (Table 2). Group E is composed of some functionally diversified transcription regulators whose GO terms are enriched in many aspects of transcription, such as transcription regulator activity, (negative) regulation of transcription, (negative) regulation of RNA metabolic process, (negative) regulation of transcription from RNA polymerase II promoter, (negative) regulation of nucleobase, nucleoside, and nucleotide and nucleic acid metabolic process, (negative) regulation of biosynthetic process, DNA binding, and protein heterodimerization activity. There are some special GO terms in group E, such as chordate embryonic development, floor plate development, neural tube development, anterior/posterior pattern formation, and (negative) regulation of muscle development (Table 2). KEGG terms, like TGF-beta signaling pathway and Notch signaling pathway, also provide key annotations and insights for bHLH members in group E.    Fig  Fig  92  100  Xenopus tropicalis bHLH genes were named according to their human orthologous genes' names (or common abbreviations) and the referred nomenclature was mainly from the tables and additional tables provided by Ledent et al. [4] and Simionato et al. [6]. Bootstrap values were converted from phylogenetic analyses with human bHLH sequences using BI and MLE algorithm, respectively. MLE bootstrap value a refers to the result from maximum likelihood estimate in phylogenetic analysis, and BI posterior probability b refers to the result from BI in phylogenetic analysis. The numbers in the phylogenetic trees are converted into percentages. c The accession numbers were retrieved from the following resources; this sequence was verified by many EST TBLASTN search hits, such as EG651417.1 and CX503003.2 (EST accession number). These numbered as "NP" were from the RefSeq protein database and those numbered as "XP" were from the Build protein database. Notes in the brackets are also gene symbols according to records in NCBI and Xenbase. All of the bHLH genes are organized in the order of bHLH families manifested in Table 1 of Ledent et al. [4]. The question mark means no matching; mark n/m * means no monophyletic group with single particular orthologous gene sequences, but formed a monophyletic group with two or more orthologous gene sequences of the family; mark n/m denotes the case of lower bootstrap value estimated less than 50%. e The accession numbers were retrieved from the ab initio protein database.

Pathways Analysis.
We could identify and select significantly enriched gene ontology terms and pathways using DAVID [31,32] and KOBAS [33][34][35] in the present study. We selected functional categories that were more likely to be biologically meaningful by calculating the statistical significance of each functional category in the input set of genes versus all annotated genes in the Xenopus tropicalis genome. After the GO annotations of Xenopus tropicalis bHLH transcription factors with the DAVID Bioinformatics Tools, all of the bHLH transcription factor genes were also subjected to KOBAS analysis (http://kobas.cbi.pku.edu.cn/home.do) and significant pathways were retrieved at the default values. We applied KOBAS to first annotate all of the genes with KO and to then identify both the most frequent and the statistically significantly enriched pathways. With the strict cutoff of FDR ≤ 0.05, KOBAS found statistically significantly enriched pathways in public databases, such as KEEG, Reactome, and PANTHER, as shown in Table 3. Using this threshold, we identified 16 pathways as induced in the Xenopus tropicalis genomic gene samples (Table 3). Among these pathways, 11 pathways were from KEEG database, while six pathways were at the significant level of ≤ 0.05. Interestingly, four of the main central cell signaling systems, that is, Notch signaling pathway, Wnt signaling pathway, TGF-beta signaling pathway, and MAPK signaling pathway, were identified. There were two most significant components related to Notch signaling pathway (corrected value 0.0024084 and 0.0150668) and circadian clock and/or circadian rhythm regulation (corrected value 0.0001219 and 0.0398896), respectively. The Jak-STAT signaling pathway, which is regarded as one of [ [39][40][41][42][43] and bHLH domains were shaded using GeneDoc. Family and bHLH protein names and high-order groups were organized according to Table 1 in the paper of Ledent et al. [4]. Highly conserved sites are shaded in black and indicated with asterisks on the top.  the central cell signaling systemS for muscle development, was identified too. It was the same case that many bHLH proteins were enriched in TGF-beta signaling pathway and Notch signaling pathway as annotated using DAVID Bioinformatics Resources. Furthermore, many interesting pathways were also identified as significantly, such as ErbB signaling pathway, Fanconi anemia pathway, and herpes simplex infection.

Concluding Remarks
In this research, we have identified 105 bHLH domains and their protein sequences in the Xenopus tropicalis genome databases by TBLASTN, BLASTP, and PSI-BLAST searches with the 45 representative bHLH domains as query sequences. Among these bHLH members, 34 hypothetical proteins, such as LOC100124777, were newly annotated by computational analysis and verified by EST searching in this research. These uncharacterized putative bHLH proteins may be novel transcription factors, which need further validation. The prediction of Xenopus tropicalis bHLH transcriptional factors will be very useful for the experiment identifying novel bHLH transcription factors and the construction of transcriptional regulatory network of Xenopus tropicalis. Through phylogenetic analyses of the Xenopus tropicalis bHLH protein domains with human bHLH orthologous protein sequences, we assigned the 105 Xenopus tropicalis bHLH genes to 43 families and two orphan genes according to the 45 defined bHLH families [3,11]. Two families, for example, Mist and Delilah, were not found in the study.
Further analysis of the Xenopus tropicalis bHLH transcription factors and their functional properties showed that 96 out of 105 bHLH genes could be annotated and only four supergroups' GO enrichment by categories were available [4]. GO enrichment statistics showed 65 significant GO annotations of biological processes and molecular functions counted in frequency. Besides common GO term categories of bHLH transcriptional factors, a large number of Xenopus tropicalis bHLH genes play significant role in muscle and organ development, chordate and neural development, floor plate and eye development, and so forth [39][40][41][42][43][48][49][50][51][52][53][54][55]. Moreover, as the group analysis results described, different groups of proteins have their special gene functions when taking no account of the common GO term categories. The trends of the gene function enrichment may be led by their DNA-binding specificity [51][52][53][54][55]. Therefore, the biology function of the uncharacterized genes or proteins can be predicted through the function GO annotation of the group analysis. To explore the functional pathways, regulatory gene networks and/or related gene groups coding for Xenopus tropicalis bHLH proteins, the identified bHLH genes were put into the databases KOBAS and STRING to get the signaling information of pathways and protein interaction networks according to available public databases and known protein interactions. From the KOBAS genomic annotation and pathway analysis, we identified 16 pathways in the Xenopus tropicalis genome. From the STRING interaction analysis, 68 hub proteins were identified and many hub proteins created a tight network or a functional module within their protein families.
The present research deepens our knowledge of frog bHLH transcription factors and provides a solid framework for further research on the functional and evolutionary aspects of Xenopus tropicalis bHLH transcription factors.