Abstract: The basic helix-loop-helix (bHLH) transcription factors are the products of main regulatory genes in metazoans. Regulatory genes are usually believed to play a large role in morphological and genetic diversifications and are often characterized by elevatedrates of evolution. In previous studies, it was showed there was much genetic diversification among the members of Hairy and Enhancer of split (Hairy/E(spl) or H/E(spl)) bHLH family. This promoted us to do further research on the mechanism of genetic diversification within the H/E(spl) bHLH family. With all the data available, we carried out an evolutionary analysis among sequence sites of the enlarged H/E(spl) bHLH family and reported that all the Hes gene lineages were under purifying selection, but the Hes2, Hes3, Hes5 lineages evolved faster than the other ones with relatively higher omega values. We found nine sites of Hes5 genes under significantly positive selection and concluded that the diversity of H/E(spl) genes has been established through relaxed selective constraint.
INTRODUCTION
Transcription Factor (TF) plays crucial roles in the regulation processes of gene expression of all living organisms at the transcriptional levels, such as development, growth and responses to environmental stimulus, etc. Some of these families are common to most eukaryotic organisms, while others are specific to particular taxonomic groups (Luscombe et al., 2000; Riechmann et al., 2000). The basic helix-loop-helix (bHLH) TFs constitute one of the largest families of functionally important proteins and were believed to play key roles in a wide range of developmental processes in metazoan and other organisms (Atchley and Fitch, 1997; Massari and Murre, 2000; Ledent and Vervoort, 2001; Jones, 2004). Murre et al. (1989) reported the first case of bHLH TFs E12 and E47 found in the murine organism. Many studies focused on the bHLH transcription factors and aimed at identification and characterization of their functions and classifications. Especially in recent years, the number of bHLH TFs characterized rapidly increased. More and more bHLH TF genes has been identified and bHLH TF families have been analyzed in many organisms whose genome drafts have been available (Morgenstern and Atchley, 1999; Atchley et al., 1999; Ledent et al., 2002; Toledo-Ortiz et al., 2003; Heim et al., 2003; Li et al., 2006a, b; Simionato et al., 2007, 2008; Stevens et al., 2008; Wang et al., 2009; Zheng et al., 2009; Pires and Dolan, 2010; Carretero-Paulet et al., 2010; Liu and Zhao, 2011). With the genome resources of interested organisms being available, it would be desirable to have a more refined classification scheme of various types of bHLH motifs, as well as a better understanding of their functions and evolutionary implications within or among species.
Previous studies found that there was much genetic diversification among the members of H/E(spl) bHLH family (Zheng et al., 2009; Liu and Zhao, 2011). However, the mechanism of genetic diversification within this gene family was still unclear. From the viewpoint of gene evolution, regulatory genes believed to play important roles in morphological and genetic diversifications are often characterized by the elevated rates of evolution. This promoted us to do further researches within the H/E(spl) bHLH family. In this study, with all of the data from previous studies and our unpublished data, we made phylogenetic analyses of the H/E(spl) bHLH transcription factors. We next carried out an evolutionary analysis among sequence sites of the enlarged H/E(spl) bHLH family and reported the results of selective pressure estimation among sites of H/E(spl) bHLH genes.
MATERIALS AND METHODS
Data collection, sequence alignment and selection: A number of closely related genes, known as the Hairy and Enhancer of split (Hairy/E(spl) or H/E(spl)) bHLH family, including Hes, Her, or ESR genes and enhancer of split related factors have been isolated from many vertebrates. Like the Drosophila E(spl) genes, many of their vertebrate homologues are expressed in response to Notch activity (Campos-Ortega, 1993, 1994, 1995) and the products of these genes are essential to implement many of the cell fate decisions mediated by Notch signaling, such as the selection of cells to become neural precursors (Artavanis-Tsakonas et al., 1995; Greenwald, 1998). More than 300 sequences were downloaded and finally a total number of 163 gene sequences from 15 different species were analyzed. The bHLH proteins alignment by ClustalX 2.0 (Thompson et al., 1997) was transformed into an aligned cds (coding sequence) fasta file using PAL2NAL (Suyama et al., 2006), which is a program to construct multiple codon alignments from matching amino acid sequences. Incomplete sequences and highly divergent regions or gaps resulting in uncertain alignments were excluded from the further analysis. The final data set included a total of 163 sequences from 15 different species (sequence data not shown), including Hey/Herp, Dec, Hes and Hes-like genes, i.e., Homo sapiens (17), Pan troglodytes (10), Macaca mulatta (1), Canis lupus (6), Sus scrofa (1), Bos taurus (16), Mus musculus (15), Rattus norvegicus (11), Danio rerio (21), Gallus gallus (6), Taeniopygia guttata (10), Xenopus tropicalis (13), Xenopus laevis (17), Apis mellifera (7) and Drosophila melanogaster (12).
Phylogenetic analyses of H/E(spl) bHLH factors: Phylogenetic analyses were conducted by MRBAYES 3.1.2 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) and PHYML 2.4.4 (Guindon and Gascuel, 2003), with the JTT substitution frequency matrix (Jones et al., 1992). The obtained bHLH sequences were used to construct phylogenetic trees of bayesian inference and maximum likelihood matching with the human bHLH domains. Maximum Likelihood (ML) analyses were performed using the frequencies of amino acids estimated from the data set and rate heterogeneity across sites modeled by one constant rate and eight γ-rates. Statistical support for the different internal branches was assessed by bootstrap resampling with 100 replicates in PHYML (Guindon and Gascuel, 2003). Phylogenetic analyses of Bayesian inference was performed with two independent Markov chains, each containing from 80 to 150 million Monte Carlo steps until the standard deviation of split frequencies was below 0.01 with sample frequency saved every 1000 generations. Finally, the trees obtained in the two runs of Markov chains were meshed and the first 25% of the trees were discarded as ‘burnin’ and the 50% majority consensus trees were edited and displayed.
Evolutionary analysis and selective pressure estimation: To determine whether the H/E (spl) split bHLH family was subject to positive (diversifying) or negative (purifying) selection, the substitution rates at synonymous (dS) and nonsynonymous (dN) were calculated and compared on the basis of nucleotide coding regions of Hes and Hes-like genes. The dN should be significantly smaller than the dS (ω=dN/dS<1) for a protein-coding gene under functional constraint or purifying selection, while the dN would be significantly larger than the dS (i.e., ω=dN/dS>1) under positive selection, which drives the rapid evolution of genes (Yang and Nielsen, 2002; Yang et al., 2005; Yang, 2007; Anisimova and Yang, 2007). Different regions with different functions may experience different selection pressure and positive selection would likely affect a few amino acids at particular genes. Thus, variable selective pressure among sites should be detected carefully.
In total, a number of 163 gene sequences from 15 different species were analyzed. Numbers of dS and dN nucleotide substitutions per site between orthologous gene sequences and likelihood ratios were computed using the Maximum Likelihood approach with alternative models using the CODEML program implemented in the PAML 4.4 package (Yang, 2007). In the maximum likelihood approach, site-specific models M1 (neutral), M2 (selection), M7 (beta), M8 (beta and ω; Nielsen and Yang, 1998; Yang et al., 2000) and M8a (beta and ω=1; Swanson et al., 2003) are used with CODEML program in the PAML 4.4 package (Yang, 2007). Model M1 (neutral) allows two classes of sites with ω0=0 and ω1=1 in proportions p0 and p1=1-p0, respectively. Model M2 (selection) has an additional class with ω2, which takes on any nonnegative value and applies to a proportion p2 of sites, now with the constraint p0+p1+p2=1. Test for positive selection should be done by comparing twice the log-likelihood difference between M1 and M2 with a χ2 distribution in the LRT (Yang et al., 2000). Model M7 (beta) assumes a β-distribution for 0≤ω≤1. Model M8 (beta and ω) adds to M7 an extra category, with proportion p1 of sites with ω1, while the rest of sites (at frequency p0=1-p1) have ω from the β-distribution between 0 and 1. Here we compare twice the log-likelihood difference between M7 and M8 with a χ2 distribution to test for positive selection (Yang et al., 2000; Anisimova et al., 2001). Model M8a is similar to model M8 except that the category ω1 is fixed at ω1=1.
We estimated the selective pressures acting on the coding regions by taking the phylogenetic-based Maximum Likelihood (ML) analyses. These analyses were conducted under different competing evolutionary hypotheses. We first investigated whether the distribution of selective constraints acting on the each gene fluctuated across lineages by comparing the fit to data of the one ratio model (M0), which assumes a constant selective pressure across branches, with the free ratios model (FR). We next examined other evolutionary scenarios to detect sites under positive selection in Hes gene lineages, using three codon-based ML substitution models that are site-specific (i.e., models with different ω ratios among sites) (Yang et al., 2000, 2005; Yang and Nielsen, 2002) but assume the same selection pattern for a site in all lineages. The likelihood Ratio Test (LRT) was used to compare the fit to data of two nested models, assuming that twice the log likelihood difference between the two models (2ΔL) follows a χ2 distribution with a number of degrees of freedom equal to the different numbers of free parameters.
RESULTS AND DISCUSSION
Phylogenetic analyses of H/E(spl) bHLH factors in different animal species: Determining the phylogenetic relationships of the bHLH proteins is an important step for elucidating the evolutionary and functional divergence of this gene family as well. Meanwile, with genome sequence data for more and more species becoming available, it is now feasible to compare the bHLH gene family among different animal species at the genome-wide level. Herein, phylogenetic analyses of Bayesian Inference (BI) and maximum likelihood estimate (ML) were used to identify putative orthologous relationships of the H/E (spl) split bHLH members in different phylogenetic trees. The H/E (spl) split bHLH family, including diverse hairy/enhancer of split related genes (mainly Hes proteins and Hes-like factors). Previously, we showed most invertebrate species have an invariable number of either 11 or 12, while the vertebrate species have 6-15 members in the H/E(spl) bHLH family (Liu and Zhao, 2011). The phylogenetic tree of hairy/enhancer of split like orthologous genes from human, mouse, rat, zebrafish, Xenopus tropicalis and chicken was explored by a maximum likelihood method with bHLH protein sequences and the zebrafish HEYL being used as out-group (Fig. 1). It was found that Hes genes from human, mouse, rat, zebrafish, Xenopus tropicalis and chicken form clear monophyletic groups (except Hes4), indicating that each Hes lineage has its own ancestral sequence (Fig. 1).
Detecting variable selective pressure among sites of the H/E(spl) bHLH family: The analysis of nonsynonymous to synonymous substitution rate ratio was used to detect selective pressure on the enlarged H/E(spl) gene sequences. We estimated ω as an average over all sites and branches from the Hes and Hes-like paralogus sequences and the ratio was substantially smaller than 1 (one ratio model, ω=0.2742) that indicated that purifying selection had been the predominant force acting on the evolution of these genes. However, selective constrains, are unevenly distributed across the phylogeny (data not shown). Models of variable ω ratios among sequence sites were used to test for the presence of sites under positive selection (with ω=1) and to identify them. Positive selection would likely affect a few amino acids at specific lineages on the phylogeny and models estimating ω ratios averaged by codons or amid acids are highly conservative. Variations in ω among sites could occur along with particular lineage heterogeneity, which can be tested by fitting the data to a model comprising different site classes. The test results of site models for each lineage were shown in Table 1.
The M3 and M0 LRT was significant, indicating that one category of ω was not fit data well to describe the variability in selection pressure across amino acid sites. The estimates of ω for Hes1, Hes2, Hes3, Hes4, Hes5, Hes6, Hes7, Hes-like, Hey/Herp and Dec gene phylogenies were 0.0439, 0.2149, 0.2150, 0.0985, 0.2464, 0.1297, 0.1303, 0.1280, 0.0850 and 0.1300, respectively (Table 1). The tests contrasting the models M1a against M2a resulted in the p values less than 1 for all the lineages or groups suggesting a lack of power and the amino acid changes within each cluster were assumed neutral or under negative selection. M1a, the parameter estimates for the least parameter rich model describes that most sites with low ω estimates (indicative of selective constraints). The estimates of ω for Hes1, Hes2, Hes3, Hes4, Hes5, Hes6, Hes7, Hes-like, Hey/Herp and Dec gene phylogenies were 0.1518, 0.4312, 0.4392, 0.1226, 0.3011, 0.378, 0.142, 0.3691, 0.1264 and 0.2801, respectively (Table 1).
| Fig. 1: | Phylogenetic tree of the H/E(spl) bHLH family members from human, mouse, rat, zebrafish, Xenopus laevis, and chicken, with the zebrafish HEYL as outgroup. Figures around the node are the transferred bootstrap values of maximum likelihood estimates (MLEs) for corresponding branches. The phylogenetic tree of Hes factor motifs revealed that Hes1, Hes2, Hes3, Hes5, Hes6, and Hes7 had their own common ancestor sequences, respectively |
| Table 1: | Selection pressures detection and model comparison among sequence sites for the H/E(spl) gene lineages |
| Codon sites with ambiguity data were removed from the alignment before analysis. Marker ** Means sites inferred under selection at the 99% level (p>99%) and marker * means those at the 95% level (p>95%) by Bayes Empirical Bayes (BEB) analysis (Yang et al., 2005) | |
The test using M7 against M8, which allows for beta-distributed site-specific ω ratio, detected no gene lineage or group under positive selection. The estimates of ω for Hes1, Hes2, Hes3, Hes4, Hes5, Hes6, Hes7, Hes-like, Hey/Herp and Dec gene phylogenies were 0.0467, 0.2250, 0.2378, 0.0952, 0.2422, 0.1319, 0.1325, 0.1325, 0.0834 and 0.1319, respectively (Table 1). But for Hes5 genes, Model M8 found nine sites were significantly under possible positive selection with Bayes Empirical Bayes (BEB) analysis (six at the 0.05 level and three at the 0.001 level).
The Naïve Empirical Bayes (NEB) (Nielsen and Yang, 1998; Yang et al., 2000) and the Bayes Empirical Bayes (BEB) (Yang et al., 2005) methods are implemented since PAML version 3.14 for calculating the posterior probabilities for site classes and identifying sites under positive selection if the likelihood ratio test is significant. NEB uses the MLEs (Maximum likelihood Estimates) of parameters but do not account for their sampling errors, while BEB deals with the sampling errors by applying a Bayesian prior. BEB is available and used under models M2a and M8 only. From these three pair-wise nested comparisons, we found all the 163 Hes gene lineages were under purifying selection and/or strong selective constraints, but there were different evolutionary rates among these lineages. From the selection and neutral model (M2a and M1a) analyses, Hes2, Hes3, Hes5 lineages obviously evolved faster than the other ones with relatively higher ω values (i.e., 0.4312, 0.4392, 0.3011, Table 1). From the site model (M8 and M7) analyses, all the Hes gene sequences, except Hes5 genes, have no site under significant or strong positive selective pressure. No positively selected sites were found in the Hes3 gene lineages. Similar to Hes5 genes, some sites from Hes1, Hes2, Hes4, Hes6 and Hes7 lineages appeared to be under putative positive selection, but the p values detected with BEB analysis were not significant (p>0.05).
Hairy and Enhancer of split (H/E(spl) in brief), including several Hes or Her or ESR genes in vertebrates, is one of the Drosophila melanogaster neurogenic loci required for normal segregation of neural and epidermal cell progenitors (Campos-Ortega, 1993, 1994, 1995). The E(spl) locus was found to be composed of a complex of several genes clustered together, called E(sp1)-C (enhancer of split complex), required to allow epidermal development of neuroectodermal cells (Campos-Ortega, 1993, 1994, 1995). These genes and others formed the so-called H/E(spl) gene family. Apart from their presence in epidermoblasts, H/E(spl) transcripts have occasionally been detected in the neuroblasts after separation of the lineages. Molecular data suggested that at least seven transcription units are members of the gene complexes, they were found to be expressed in nearly identical patterns (Knust et al., 1987, 1992). These seven genes of enhancer of split complex in Drosophila melanogaster encode basic-helix-loop-helix transcription factors (Md, Mb, Mg, M3, M5, M7 and M8) which are components of the Notch signaling pathway. Their expression are in response to Notch activation and mediate some effects of the Notch signaling pathway by regulating the expression of target genes (Artavanis-Tsakonas et al., 1995; Greenwald, 1998).
Notch function is essential for many cell fate decisions in Drosophila that involve the singling out of one cell from a group of competent cells (Muskavitch, 1994). In these processes, the Notch protein appears to function as a receptor in a cell-cell signaling pathway whose other components include two ligands ‘Delta’ and ‘Serrate’, the intracellular transducer suppressor of hairless and the nuclear proteins encoded by the enhancer of split complex (Artavanis-Tsakonas et al., 1995). The most immediate transcriptional target genes of Notch activation in Drosophila melanogaster encode seven bHLH proteins (Md, Mb, Mg, M3, M5, M7 and M8) which are clustered in the E(spl)-C (Jennings et al., 1999). The expression of E(spl) genes are required downstream of Notch since the accumulation of E(spl) proteins depends on Notch activation and they mediate the effects of constitutively active forms of Notch during neurogenesis (Artavanis-Tsakonas et al., 1995; Greenwald, 1998). Therefore, mutations of E(spl) genes are deleterious during neurogenesis of embryo. Detailed genetic analyses suggested that there is some overlap in the functions of these proteins (Celis de et al., 1996). Consistent with the proposed functional redundancy, the patterns of expression of E(spl) genes are also very similar during embryo neurogenesis (Knust et al., 1987, 1992; Celis de et al., 1996). This could explain why there are many E(spl) genes in different vertebrates and whether there is possible positive selection driving the genetic diversity of the E(spl) proteins.
Moreover, all the seven E(spl) proteins have two additional helix motifs and terminate with the same tetrapeptide ‘WRPW’, besides the HLH domain. And the degree of conservation is variable in the two additional regions (helix III/IV, or helix III and IV). Alignment showed a region of 41 amino acids present in seven proteins of the E(spl)-C. Besides the great similarity in the basic domain of the E(sp1) proteins, including the proline residue (Fig. 1), the other domains are also conserved, such as the HLH domain and the C-terminal tetrapeptide ‘WRPW’. Although there is some amino acid sequence conservation in the hypothetical helix III/IV region of the E(sp1) proteins, the ability to form an a helix is only apparent in helix IV and less obvious in helix III. It is obvious that the two helix motifs are more variable than the other regions and drive the diversification of E(sp1) proteins. This was verified by our analysis of selective pressure on the nine sites of Hes5 genes. Among Hes5 genes, the nine sites under putative positive selection were 76T**, 158S**, 169S*, 171K**, 172G**, 174P*, 176A*, 177T** and 179K** (corresponding sites before deleting alignment gaps; the markers, ** and *, denote sites under selection at the 99 and 95% level, respectively). It is found that only one site (76T) lies in the HLH domain, while the other eight sites locate in the rest part following the HLH domain including the additional Hairy Orange domain.
Determining selective pressure among sites of the H/E(spl) bHLH domains: The analysis of nonsynonymous to synonymous substitution rate ratio was further used to detect possible selective pressure on the H/E(spl) bHLH and/or HLH domains. To determine whether there has been putative positive selection, we also did three analyses (M3 and M0, M2a and M1a, M8 and M7) and estimated ω values over all sites and branches from the Hes and Hes-like paralogus sequences. From these pairwise nested comparisons, we found the 163 Hes bHLH and/or HLH domains were under purifying selection or strong function constraints, although there were different evolutionary rates among different gene lineages. From the selection and neutral model (M2a and M1a) analyses, Hes3, Hes5, Hes7 lineages obviously evolved faster than the other ones with relatively higher ω values (i.e., 0.2605, 0.1941, 0.1202, Table 2). From the site model (M8 and M7) analyses, all the Hes gene sequences have no site, except one site (37N, corresponding site before deleting alignment gaps) in Hes2 genes, under putative positive selective pressure (p>0.05) and Hes2, Hes3, Hes5 lineages obviously evolved faster than the other ones with relatively higher ω values (i.e., 0.1697, 0.1156, 0.1175, Table 2). However, all the estimated ratios were all substantially smaller than 1 (Table 2), indicating that there might be predominant purifying selection affecting on the evolution of bHLH or HLH domains.
| Table 2: | Selection pressures detection and model comparison among sequence sites for the H/E(spl) bHLH domains |
| Marker ** means sites inferred under selection at the 99% level (p>99%) and marker * Means those at the 95% level (p>95%) by Bayes Empirical Bayes (BEB) analysis (Yang et al., 2005) | |
CONCLUSIONS
The basic helix-loop-helix (bHLH) transcription factors are key regulators in cell and organ and organism, especially animal’s muscle development and the adaptation to hypoxia (Massari and Murre, 2000; Velleman et al., 2010; Liu et al., 2010; Fakheri et al., 2010; Gonchar and Mankovska, 2010; Song et al., 2010). Present study showed that a considerable number of bHLH genes were found to have a multi-member distribution pattern in human, mouse, rat, zebrafish, chicken and Zebra Finch bHLH gene families (Liu and Zhao, 2011). This case suggests that they should arise through gene duplication at least before the divergence of vertebrates from invertebrates. In the research, phylogenetic analysis of the H/E (spl) family revealed that more than four gene duplication events should have occurred at an early date. The primary objective of this study was also to determine whether accelerated rates of evolution in the H/E (spl) members or bHLH domains are due to increased positive selection or decreased constraint. Present results provide little support for the positive selection hypothesis. By modeling the variable selective pressure among sequence sites of the enlarged H/E(spl) bHLH family, we found that all the 163 Hes gene lineages were under purifying selection and/or strong selective constraints, but the Hes2, Hes3, Hes5 lineages evolved faster than the other ones with relatively ω values. We found that nine sites of the Hes5 gene lineage had been under significantly positive selection. We further explored and described the evolutionary analysis for bHLH and/or HLH domains and found no site under positive selection. More generally, the failure to detect positive selection or adaptive substitution in E(spl) proteins and other transcription factors suggests that selective function constraint is relaxed and most substitutions are neutral in variable domains of these genes (Streisfeld and Rausher, 2007). This provides us an evolutionary scenario in which the diversity of E(spl) genes has been established through relaxed selective constraint. The present study deepens our knowledge of the animal Hairy and Enhancer of split bHLH transcription factors.
ACKNOWLEDGMENTS
We are grateful to the anonymous reviewers for their constructive comments and suggestions. This study was supported by Research Funds to LWY (2011PTFY04ZD, KJ2012A216).