Recent comparative genomic studies have identified many human accelerated elements (HARs) with elevated substitution rates in the human lineage. However, it remains unknown to what extent transcription factor binding sites (TFBSs) are under accelerated evolution in humans and other primates. Here, we introduce two pooling-based phylogenetic methods with dramatically enhanced sensitivity to examine accelerated evolution in TFBSs. Using these new methods, we show that more than 6000 TFBSs annotated in the human genome have experienced accelerated evolution in Hominini, apes, and Old World monkeys. Although these TFBSs individually show relatively weak signals of accelerated evolution, they collectively are more abundant than HARs. Also, we show that accelerated evolution in Pol III binding sites may be driven by lineage-specific positive selection, whereas accelerated evolution in other TFBSs might be driven by nonadaptive evolutionary forces. Finally, the accelerated TFBSs are enriched around developmental genes, suggesting that accelerated evolution in TFBSs may drive the divergence of developmental processes between primates.
During the course of evolution, a subset of genes and regulatory elements may be subject to different pressures of natural selection in distinct species. These genomic elements often have varying substitution rates across species, which may be identified by phylogenetic models with lineage-specific substitution rates 1,2,3,4,5,6,7,8,9,10 . Notably, previous studies have revealed a few thousand human accelerated regions (HARs) with dramatically elevated substitution rates in the human lineage compared to other vertebrates 7,8,9,10,11,12,13,14,15 . A large proportion of HARs are neural enhancers 14,15,16,17,18,19,20 and frequently subject to strong positive selection in the human lineage 9,10,11,21 , suggesting that they may contribute to the adaptive evolution of human brain. Also, recent studies show that deleterious mutations in HARs may be associated with neurodevelopmental disorders 22,23,24,25,26 , highlighting the key role of HARs in maintaining the integrity of the central nervous system. Thus, characterizing genomic elements under accelerated evolution is of great importance for understanding the genomic basis of human evolution and disease.
While numerous studies have been conducted to examine accelerated evolution in humans and other species 7,8,9,10,11,12,13,14,15,27,28,29 , the existing studies may suffer from two critical limitations. First, most of the previous studies have focused on conserved noncoding elements under accelerated evolution. Because a large proportion of noncoding regulatory elements may be subject to frequent evolutionary turnover 30,31,32,33,34 , these studies may not be able to characterize accelerated evolution in non-conserved regulatory elements. Second, the previous studies have focused on identifying individual HARs with a genome-wide level of significance. Because of the small amount of alignment data in a single genomic element and the high burden of multiple testing correction associated with a genome-wide scan, these studies may have limited statistical power to detect weakly accelerated evolution driven by relaxed purifying selection or weak positive selection. Altogether, it remains unknown to what extent non-conserved genomic elements are subject to weakly accelerated evolution.
Substitutions in transcription factor binding sites (TFBSs) are a main driver of phenotype diversity between species 35,36,37 , implying that TFBSs may also be subject to accelerated evolution. However, to the best of our knowledge, accelerated evolution in TFBSs has not been systematically explored in previous studies, possibly because the majority of TFBSs are not highly conserved across vertebrates 30,31,32,33,38,39 . Also, TFBSs might be subject to weaker acceleration compared to conserved elements because the phenotypic effects of mutations are weaker in TFBSs than in conserved elements 40 . Therefore, previous phylogenetic methods dedicated to infer strong signals of accelerated evolution may be underpowered to detect TFBSs under weakly accelerated evolution.
Here, we introduce two novel phylogenetic methods for exploring TFBSs under accelerated evolution. Unlike previous methods that analyze individual elements separately 7,8,9,10,11,12,13,14,15 , our new approaches pool thousands of TFBSs with similar functions together to boost the statistical power to detect weak signals of accelerated evolution. These new methods allow us to rigorously test whether a group of TFBSs as a whole is significantly enriched with accelerated elements, despite that we may lack statistical power to identity individual TFBSs under accelerated evolution due to limited alignment data in a single TFBS. Using these methods, we show that TFBSs of numerous transcription factors are likely to be under accelerated evolution in Hominini, apes, and Old World monkeys. Compared to previously identified HARs, these TFBSs show weaker acceleration but are more abundant genome-wide. Among these accelerated TFBSs, binding sites of DNA-directed RNA polymerase III (Pol III) show the strongest signal of acceleration, which might be driven by strong lineage-specific positive selection on par with HARs. Taken together, accelerated evolution may be a common characteristic of TFBSs in Hominini, apes, and Old World monkeys.
In the current study, we introduce a novel software application, GroupAcc, which includes two pooling-based phylogenetic approaches with improved statistical power to infer weakly accelerated evolution. The key idea of GroupAcc is to group TFBSs by the bound transcription factor and then examine whether each TFBS group as a whole shows an elevated substitution rate in a lineage of interest. In this study, TFBSs refer to peaks in Chromatin immunoprecipitation-sequencing(ChIP-seq) 41 . By pooling alignment data from a large number of TFBSs, our new methods have significantly higher statistical power to detect weakly accelerated evolution at the group level even if the signals of acceleration are statistically insignificant at the level of individual TFBSs.
In the first method, we utilize a group-level likelihood ratio test (LRT) to infer whether a TFBS group as a whole shows an elevated substitution rate in a predefined foreground lineage compared to the other lineage (background lineage) (Fig. 1). To this end, we first fit a reference phylogenetic model to the concatenated alignment of all TFBSs, where we estimate the branch lengths of a phylogenetic tree, the gamma shape parameter for rate variation among nucleotide sites 42 , and the parameters of the general time reversible substitution model 43 . Assuming that the majority of TFBSs may not be under accelerated evolution, the reference phylogenetic model may represent the overall pattern of sequence evolution in TFBSs when accelerated evolution is absent. Given the reference phylogenetic tree, we then fit the group-level LRT to the concatenated alignment of a TFBS group, where we estimate two scaling factors, r1 and r2, for the foreground and the background branches, respectively. We interpret r1 and r2 as the relative substitution rates of the TFBS group in the foreground and background lineages throughout this study. We assume that r1 = r2 in the null model (H0), indicating that the TFBS group has evolved at a constant rate across all lineages. Conversely, we assume that r1 ≠ r2 in the alternative model (Ha), indicating that the TFBS group has evolved at different substitution rates between the foreground and background lineages. Because the TFBS group may consist of hundreds of TFBSs, we assume that the likelihood ratio statistic of the group-level LRT asymptotically follows a chi-square distribution with one degree of freedom. If the null model is rejected in the group-level LRT and r1 > r2, we consider that the TFBS group as a whole may be subject to accelerated evolution.
In the second method, we use a phylogenetics-based mixture model to estimate the proportion of accelerated TFBSs in a TFBS group (Fig. 1). To this end, we first perform an element-level LRT to infer evidence for accelerated evolution in individual TFBSs given that H0 is rejected in the group-level LRT. The element-level LRT is similar to the group-level LRT but is applied to the alignments of individual TFBSs rather than to the concatenated alignment of the TFBS group. Given the likelihood ratio statistics from the element-level LRT, we then calculate empirical p-values for individual TFBSs using parametric bootstrapping. Unlike the chi-square distribution in the group-level LRT, the parametric bootstrapping procedure provides accurate p-values even when there is a small amount of alignment data per test 9 . Finally, we estimate the proportion of accelerated TFBSs in the TFBS group by fitting a beta-uniform mixture model to the distribution of p-values 44 . The beta-uniform mixture model allows us to estimate an upper bound of the proportion of TFBSs generated from H0 ( \(<\hat<\pi >>_\) ). We consider \(1-<\hat<\pi >>_\) as a conservative estimate (lower bound) of the proportion of accelerated TFBSs.
To verify the power of GroupAcc to infer weakly accelerated evolution in TFBS groups, we conducted simulations under various lineage-specific evolutionary dynamics. In the first scenario, we assumed that all the binding sites in one group are under accelerated evolution in a specific lineage. The second scenario considered the heterogeneity of evolutionary patterns in each single binding site: only parts of each binding site (for example, motif) undergo accelerated evolution in a specific lineage. The third scenario considered the heterogeneity of evolutionary dynamics in groups of binding sites: only certain numbers of binding sites in one group undergo accelerated evolution in a specific lineage, while the other binding sites do not undergo accelerated evolution. Under each scenario, we verified the ability of the group-level LRT method to detect accelerated evolution in a specific lineage and estimate the fold of increase in substitution rate (r1/r2). We also compared the performance of the phylogenetics-based mixture model and traditional element-level LRT in estimating the number of elements under accelerated evolution in a given lineage.
In each scenario, eight cases were generated in which different lineages of primates were under accelerated evolution: (1) only human, (2) subtree of Hominini (human, chimp), (3) subtree of human, chimp, and gorilla, (4) subtree of Great apes (chimp, gorilla, orangutan) and human, (5) only chimp, (6) only gorilla, (7) only orangutan, (8) only macaque. For each case, we simulated alignments of 10,000 binding sites, each at the length of 200 bp based on the reference model plus those assumptions. We also simulated alignments of different numbers of binding sites (1000) and different lengths of each alignment (100 bp) to test the performance of the model in different settings. Both weak and strong accelerated evolution were taken into consideration: the fold of increase in substitution rate in foreground lineage (r1/r2) varied from 1.2 to 5.
Under the first scenario, all the 200 bp binding sites in one group were assumed to be under accelerated evolution in a defined lineage as each case (1–8) showed, for example, in case 1, all the 200 bp binding sites would be under accelerated evolution in only human. With foreground lineage matching with the accelerated lineage in each case, the group-level LRT method was able to tell the presence of accelerated evolution at the group level and accurately estimate the fold of increase in substitution rate in foreground lineages (r1/r2), even given weak accelerated evolution when the fold of increase in substitution rates is only slightly larger than 1 (Fig. 2a). The GroupAcc model performed better than element-level LRT in estimating the number of elements under accelerated evolution (Fig. 2b). We also tested if the model could detect accelerated evolution in a tip if a subtree containing the tip is under accelerated evolution (Fig. 3). In cases (1), (2), (3) and (4), when accelerated evolution happened in lineages such as human or subtrees containing human, taking human as foreground lineage, GroupAcc methods were able to identify the presence of accelerated evolution in human and estimate the number of elements under accelerated evolution in human with higher accuracy compared to traditional element-level LRT method (Fig. 3). In cases (5), (6), (7) and (8), when accelerated evolution occurred in lineages other than human, the GroupAcc methods were able to identify the fact that human is not undergoing accelerated evolution (Fig. 3).
We validated the ability of our methods to identify lineage-specific acceleration when only part of the TFBS is under accelerated evolution from simulation scenario 2. We generated 10,000 200 bp alignments standing for elements. Each alignment was composed of 200 × L bp generated with a scaled tree (with substitution rate increase) and 200(1 − L) bp generated from an unscaled tree (without substitution rate increase). Given that L = 0.1, 0.2, 0.5, 0.8. the group-level LRT was able to identify the presence of accelerated evolution, even under weak acceleration when the fold of substitution rate increase in foreground lineage was only 1.2 (Fig. 4a). The GroupAcc method outperformed the element-level LRT method in estimating the number of elements under accelerated evolution (Fig. 4b). Therefore, under situations with the heterogeneity of evolutionary patterns in each single binding site, our pooling based methods were able to identify lineage-specific accelerated evolution with a uniform scaling of substitution rates on the foreground lineages across the whole binding sites.
Under the third scenario, a specific proportion M of binding sites (M = 0.1, 0.2, 0.5, 0.8) in a group were under accelerated evolution. This scenario considered the heterogeneity of evolutionary dynamics in multiple binding sites of one transcription factor. We found group-level LRT method was able to tell the presence of accelerated evolution at the group level and estimate the fold of increase in the substitution rate of foreground lineages, even when the fold of increase in substitution rates of foreground lineage was slightly larger than 1 (Fig. 5). The GroupAcc model performed better than element-level LRT in estimating the number of elements under accelerated evolution (Fig. 5).
Using the group-level LRT, we examined accelerated evolution in 4,380,444 TFBSs of 161 transcription factors identified by ChIP-seq experiments in the ENCODE Project 45 . We tested whether each group of TFBSs bound by the same transcription factor had an elevated substitution rate in the human lineage. We used Multiz genome alignments of ten primate species 46 and defined the human lineage after the divergence of chimpanzees and humans as the foreground lineage. Unlike previous studies of HARs 7,8,9,10,11,12,13,14,15 , we did not include non-primate vertebrates to mitigate the impact of the evolutionary turnover of TFBSs on our analysis 30,31,32,33,38,39 . After Bonferroni correction, we observed that 15 TFBS groups had significantly different substitution rates between the foreground and background lineages (Supplementary Data 1), which all showed elevated substitution rates in humans compared to other primates (r1 > r2).
TFBS groups with elevated substitution rates in humans could be either directly under accelerated evolution or merely overlapping with other accelerated TFBS groups. To identify TFBS groups directly under accelerated evolution, we sought to partition the binding sites of the 15 TFBS groups with elevated substitution rates into non-overlapping, biologically interpretable TFBS groups. Because BDP1, BRF1, and POLR3G are components of the Pol III transcription machinery 47 , we defined a new TFBS group, Pol III binding, consisting of genomic regions bound by at least two of the three transcription factors. Similarly, since POU5F1 and NANOG can interact with each other to form a protein complex 48,49 , we defined another TFBS group, POU5F1-NANOG binding, consisting of genomic regions bound by both of the two transcription factors.
Then, we removed all binding sites overlapping more than one TFBS group, resulting in 17 non-overlapping TFBS groups. We applied the group-level LRT again to these non-overlapping TFBS groups. After Bonferroni correction, seven non-overlapping TFBS groups showed significantly elevated substitution rates in the human lineage (Fig. 6; Supplementary Table 1). These non-overlapping TFBS groups included Pol III binding, POU5F1-NANOG binding, BDP1, FOXP2, POU5F1, NANOG, and NRF1. Compared to previously identified HARs, the seven non-overlapping TFBS groups showed weaker acceleration as evidenced by their smaller increases in substitution rates in the human lineage (Fig. 6). We focused on the seven non-overlapping TFBS groups with evidence for weakly accelerated evolution in downstream analysis.
A recent study showed that many HARs may also undergo accelerated evolution in other apes 15 . Based on the simulation of scenario 1 (Fig. 3), we found that the GroupAcc methods were able to tell the presence of accelerated evolution in human lineage at the group level when accelerated evolution occurred in any subtrees containing human. To characterize when acceleration occurred during the evolution of TFBSs, we employed a model comparison approach to search for lineages with elevated substitution rates. Specifically, we evaluated the goodness-of-fit of seven phylogenetic models with different foreground lineages, denoted as M1 to M7 (Fig. 7a). All these models were based on Ha in the group-level LRT (Fig. 1), and the foreground lineages associated with these models corresponded to all the monophyletic clades that included the human lineage (Fig. 7a). These models effectively assumed that the change of substitution rate occurred at most once during the evolution of a TFBS group, which was designed to explore the most parsimonious explanations of accelerated evolution and to limit the number of tested foreground lineages. We used the Bayesian information criterion (BIC) as a measure of goodness-of-fit of these models.
Although the seven accelerated TFBS groups were originally detected using humans as the foreground lineage, our model comparison analysis showed that accelerated evolution may not be human specific (Fig. 7b and Supplementary Table 2). Specifically, for binding sites of Pol III, BDP1, and NRF1, a model with both apes and Old World monkeys as the foreground lineage (M5) showed the best goodness-of-fit. Similarly, for binding sites of FOXP2 and NANOG, a model with apes as the foreground (M4) showed the best goodness-of-fit. Moreover, for POU5F1-NANOG and POU5F1 binding sites, a model with Hominini as the foreground (M2) showed the best goodness-of-fit. Altogether, the acceleration of TFBS evolution might be driven by changes of selection pressure in Hominini, apes, and Old World monkey and, thus, might contribute to phenotypic differences between these species and other primates.
In this section, we sought to infer the total number of TFBSs under accelerated evolution. While the group-level LRT can examine whether a TFBS group as a whole was under accelerated evolution, it could not estimate the number of accelerated TFBSs in the TFBS group. Also, because the signal of acceleration might be weak in TFBSs (Fig. 6), previous phylogenetic models could not be used to estimate this number either 7,8,9,10,11,12,13,14,15 . To address this problem, we utilized the phylogenetics-based mixture method to estimate the proportion of accelerated TFBSs from the distribution of p-values associated with individual TFBSs in the same group (Fig. 1).
We observed that 78% of Pol III binding sites were under accelerated evolution (Table 1 and Supplementary Table 3), which translates to approximately 222 accelerated Pol III binding sites. Also, 20 and 25% of binding sites of BDP1 and NRF1 were under accelerated evolution in Old World monkeys and apes, which translates to approximately 90 and 466 accelerated TFBSs, respectively. Approximately 25% of TFBSs of FOXP2 and NANOG were under accelerated evolution in apes, suggesting about 5000 binding sites in these TFBS groups were accelerated elements. Furthermore, approximately 8% of TFBSs of POU5F1 and POU5F1-NANOG were under accelerated evolution in Hominini, indicating that about 300 binding sites of the two groups were accelerated in the clade consisting of humans and chimpanzees. In total, more than 6000 TFBSs spanning 1573kb were under accelerated evolution in Hominini, apes, and Old World monkeys (Table 1), which is more than the 3098 known HARs spanning 720 kb (see “Methods”).
To identify the major functions represented by the top accelerated binding sites in the seven TFBS groups, we utilized Genomic Regions Enrichment of Annotations Tool (GREAT) to first find the potential target genes by predicting both proximal and distal binding events, and then analyzed the functional significance of those top accelerated binding sites by applying GO enrichment test and pathway enrichment analysis to their potential target genes with background gene lists composed of all the genes associated with the whole TFBS group 56,57,58 .
We extracted the significant binding sites in each of the seven groups from the phylogenetics-based mixture model and defined them as the top accelerated binding sites. GREAT identified 2611 potential target genes for the top accelerated binding sites of FOXP2, 662 genes for the top accelerated binding sites of NANOG, 390 genes for the top accelerated binding sites of NRF1, 222 genes for the top accelerated binding sites of POU5F1, 163 genes for the top accelerated binding sites shared by POU5F1 and NANOG, 104 genes for the top accelerated binding sites of BDP1 and 143 genes for the top accelerated binding sites shared by Pol III TFs. Using default settings in GREAT, we built seven background gene lists for seven TFBS groups, respectively containing 9896 potential target genes for FOXP2 binding sites, 3745 genes for NANOG binding sites, 2931 genes for NRF1 binding sites, 1976 genes for POU5F1-NANOG binding sites and 478 potential target genes for POU5F1 binding sites.
After removing the redundant GO terms with high semantic similarity (0.7) and performing Bonferroni correction on the GO enrichment results, we found FOXP2 top accelerated TFBSs were associated with genes functioning in artery development and regulation of transforming growth factor signaling pathway. The concatenation of top accelerated binding sites in seven TFBS groups were associated with genes playing roles in development and cell proliferation processes (Fig. 9; Supplementary Data 2).
In the genes associated with other accelerated TFBS groups, no pathways or biological terms were found to be significant after correction. Benjamini–Hochberg correction has been applied to the GO enrichment results. After Benjamini–Hochberg correction, developmental process terms were also enriched for the genes nearby the top accelerated binding sites among the seven groups (Supplementary Data 2).
To investigate the accelerated evolution in primates, we applied the GroupAcc method to datasets that were not human-centric, including non-human primates’ ChIP-seq peaks. Vermunt et al. 59 identified histone H3 lysine 27 acetylation (H3K27ac) enriched regions in human, chimpanzee and rhesus macaque brain. The H3K27ac enriched regions were predicted to be active cis-regulatory elements(CREs), We applied the group-level LRT method to the predicted CREs in human, chimpanzee and rhesus macaque brain. Results revealed a slight increase in substitution rates of human and chimpanzee lineage in CREs of human and chimpanzee brain, compared to the fold of increase in substitution rate of rhesus macaque lineage in CREs of rhesus macaque brain.
Villar et al. 60 identified trimethylated lysine 4 of histone H3 (H3K4me3) enriched regions and H3K27ac enriched regions in liver of 20 mammals including human and rhesus macaque. The regions were classified into active gene promoters and enhancers. Enhancers were identified by regions only enriched for H3K27ac, while promoters defined as regions containing both H3K27ac and H3K4me3. We applied the group-LRT method to the promoters and enhancers in human and rhesus macaque. Results showed that enhancers tended to evolve faster than promoters in both species.
To identify accelerated evolution in tissue-specific genetic regulatory elements, we applied the GroupAcc method to the most abundant TFBS group: CTCF binding sites. We included 80074 CTCF binding sites across 29 tissues or cell types 61 . We found lower leg skin and tibial nerve CTCF binding sites undergo weak accelerated evolution in human.
In the current study, we present two pooling-based methods to infer genomic elements under accelerated evolution. Unlike previous methods that focus on analyzing individual elements 7,8,9,10,11,12,13,14,15 , our new methods group hundreds of genomic elements with similar biological functions to increase the sample size per test and reduce the multiple testing burden. Thus, our methods may have higher sensitivity to detect weak signals of accelerated evolution. To the best of our knowledge, our methods are the first statistical framework dedicated to inferring weakly accelerated evolution in non-coding regions.
Using the group-level LRT, we identify seven groups of non-overlapping TFBSs with significant evidence for accelerated evolution (Fig. 6). The model comparison analysis suggests that these TFBS groups may be under accelerated evolution not only in humans but also in other primate species (Fig. 7). In agreement with our finding, a recent study of HARs has shown that many HARs may also be subject to accelerated evolution in other ape species 15 . Therefore, accelerated evolution of regulatory elements may be a shared characteristic of primates rather than specific to the human lineage.
Among the seven groups of accelerated TFBSs, we show that Pol III binding sites may be subject to positive selection in the human lineage but find no evidence for positive selection in other accelerated TFBS groups (Fig. 8). In contrast, more than half of HARs may be subject to positive selection in the human lineage 21 , suggesting that positive selection may be the main driving force of accelerated evolution in HARs. Because previous studies of HARs have focused on identifying individual genomic elements with extremely high substitution rates in the human lineage, the higher frequency of detecting positive selection in HARs could partially reflect the lower power of previous methods in discovering weakly accelerated elements driven by evolutionary forces other than positive selection.
Although accelerated evolution is much weaker in Pol III binding sites than in HARs (Fig. 6), Pol III binding sites are subject to strong positive selection in the human lineage, on par with HARs (Fig. 8). Because HARs are highly conserved across species, they may have very low substitution rates in non-human primates, which in turn enhances the signals of accelerated evolution. In contrast, Pol III binding sites may not be highly conserved across species, resulting in a weaker signal of accelerated evolution despite strong positive selection in the foreground lineage. Taken together, weak signals of accelerated evolution may not always imply weak positive selection in the foreground lineage.
Other than lineage-specific positive selection, we find that nonadaptive evolutionary forces, such as relaxed purifying selection and GC-biased gene conversion, may drive the accelerated evolution of TFBSs 21,55 . GC-biased gene conversion has been found in NRF1 binding sites but not in other accelerated TFBS groups. The seven groups of accelerated TFBSs have reduced values of the fraction of sites under selection ρ in human comparing to the collection of 161 groups of TFBSs. Overall, the seven groups of TFBSs are under weaker selection constraints than other TFBSs. The widespread nonadaptive evolutionary forces do not indicate the lack of functional importance of those accelerated regions.
Notably, the seven groups of accelerated TFBSs may play key roles in developmental processes. First, recent studies suggest that disruptive mutations in subunits of Pol III, such as POLR3A, POLR3B and BRF1, may be associated with neurodevelopmental disorders 62,63,64 . Therefore, accelerated evolution in Pol III binding sites might be associated with the adaptive evolution of the central nervous system in apes and Old World monkeys (Fig. 7). Second, POU5F1 and NANOG are transcription factors necessary to the pluripotency and self-renewal of embryonic stem cells 65,66,67 . The colocalization of POU5F1 and NANOG in regulatory elements, referred to as POU5F1-NANOG binding in the current study, might trigger zygotic gene activation in vertebrates 48,68,69,70,71 . Third, FOXP2 is a highly conserved vertebrate protein with high expression in the central nervous systems during embryogenesis, and detrimental mutations in the FOXP2 gene may cause impaired speech development in humans 72,73,74 . Also, previous studies have shown that the protein sequence and expression of FOXP2 could be subject to accelerated evolution in humans 75,76,77 , echolocating bats 78 , and vocal learning birds 79 . Finally, NRF1 has been found to regulate the expression of GABRB1, a gene associated with neurological and neuropsychiatric disorders 80,81 . To summarize, the collection of seven TFBS groups may be functionally related to developmental processes. Specifically, when compared with a background gene list containing all the genes associated with the collection of seven TFBS groups, developmental process terms were enriched for the genes nearby the top accelerated binding sites among the seven groups (Fig. 9). Therefore, among the collection, the binding sites with strongest signals of accelerated evolution might be more crucial to the developmental processes. Together with the fact that a large proportion of HARs are neural enhancers and subject to accelerated evolution in humans and other primates 19,20,21 , we conclude that regulatory sequences of neurodevelopmental genes may be the main target of accelerated evolution in primates.
Due to the scarcity of ChIP-seq data in non-human primates, we have used human-based TFBS annotations to infer accelerated evolution. It may limit our ability to detect accelerated evolution present in non-human primates but not in humans. Thus, our estimate of the number of accelerated TFBSs is likely to be conservative (Table 1). In future studies, it is of great interest to investigate accelerated evolution in TFBSs identified in non-human primates, highlighting the urgent need to perform high-throughput functional genomic experiments in our close relatives.
Compared to conserved genomic elements explored in previous studies of HARs, TFBSs may have a higher evolutionary turnover rate 30,31,32,33,38,39 . To alleviate the impact of evolutionary turnover on our analysis, we have only included primate genomes in the current study and filtered out low-quality alignments. Nevertheless, a small proportion of TFBSs identified in the human genome may still be subject to evolutionary turnover in other primates 39 . We expect that the evolutionary turnover of TFBSs in non-human primates may not lead to false positive results in our analysis. Indeed, conditional on the presence of a TFBS in the human genome, the evolutionary turnover of the TFBS in non-human primates is more likely to increase the substitution rate in the background lineage and hence makes our analysis conservative. Moreover, conditional on the loss of an old binding site in the human genome, the sequences would not be annotated as TFBSs in the human genome. Given that we used human genome annotation, those regions functional in background lineages but not in humans were not included in our analysis. Once ChIP-seq data become available in multiple non-human primates in the future, the conservativeness of our analysis may be alleviated by including only species where the TFBS of interest is detected.
Our pooling-based methods have a potential to be extended in future studies. For instance, if multiple TFBSs overlap with each other, our current methods cannot distinguish between TFBSs directly under accelerated evolution from those overlapping other accelerated TFBSs. To address this problem, we have used a heuristic method to remove overlapping TFBSs in the current study, which may reduce the number of TFBSs in our analysis. In the future, it is of great interest to develop a rigorous method for inferring accelerated evolution in overlapping TFBSs. Motivated by the recent success of evolution-based regression models 34,82,83,84 , we propose that unifying our pooling-based methods and generalized linear models may be a promising direction to disentangle causal from correlational relationships in the analysis of accelerated evolution.
We obtained the Multiz alignment of 46 vertebrate genomes from the UCSC Genome Browser 46 . Then, we extracted a subset of alignments for ten primate species from the 46-way Multiz alignment. The ten primate species and their genome assemblies included Homo Sapiens (hg19), Pan troglodytes (panTro2), Gorilla gorilla (gorGor1), Pongo abelii (ponAbe2), Macaca mulatta (rheMac2), Papio hamadryas (papHam1), Callithrix jacchus (calJac1), Tarsius syrichta (tarSyr1), Microcebus murinus (micMur1), and Otolemur garnettii (otoGar1). Also, we downloaded 4,380,444 TFBSs for 161 transcription factors from the UCSC Genome Browser. These TFBSs were identified by ChIP-seq experiments in the ENCODE Project 45 . We extracted alignments of TFBSs across ten primate species using PHAST 85 . We removed the TFBSs overlapping with UTRs, CDSs, and previously identified HARs. To filter out low-quality alignments, we obtained informative alignment sites where unambiguous bases were found in at least five out of ten primate species in the Multiz alignment. We retained TFBSs with at least 50 informative alignment sites for downstream analysis.
We obtained a comprehensive list of previously defined HARs from https://docpollard.org/research/. We first combined the following genetic elements: Merged list of 2649 HARs(a set of HARs in noncoding regions built by Capra et al. 17 ), 284 human accelerated elements in mammal conserved regions with adjusted p-value p-value
We built a reference phylogenetic model using the alignment of ten primate genomes, assuming that the majority of TFBSs may not be subject to accelerated evolution. We first concatenated alignments of all TFBSs. We then fit a phylogenetic model to the concatenated alignment using the phangorn library in R 86 . In the phylogenetic model, we used the generalized time-reversible (GTR) substitution model to describe nucleotide sequence evolution and the discrete Gamma distribution with four rate categories to model substitution rate variation among nucleotide sites 42 . Also, we fixed the tree topology of the reference phylogenetic model to the one used from the UCSC Genome Browser (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/multiz46way/46way.nh). We estimated model parameters, including branch lengths, the shape parameter of the discrete Gamma distribution, and parameters of the GTR substitution model, using the optim.pml function in phangorn.
Given the reference phylogenetic model, we used a customized R program based on phangorn to perform the group-level LRT. First, we concatenated alignments for each TFBS group separately. Then, we fit two group-level phylogenetic models to the concatenated alignment of each TFBS group. In the null model (H0), we inferred a global scaling factor of branch lengths with maximum likelihood estimation and fixed all other model parameters to the ones in the reference phylogenetic model. We interpreted the estimated scaling factor as the relative substitution rate of TFBS sequences in both the foreground and the background lineage. In the alternative model (Ha), we estimated two scaling factors of branch lengths, r1 and r2, for the foreground and background lineages, respectively. The two scaling factors were interpreted as the relative substitution rates in the foreground and background lineages in the alternative model. For each TFBS group, we calculated a likelihood ratio statistic defined as the two-fold difference in the log likelihood between Ha and H0. Given the likelihood ratio statistic, we obtained a p-value for each TFBS group using a chi-square test with one degree of freedom. Finally, we calculated adjusted p-values using the Bonferroni correction.
From the group-level LRT, we found that TFBSs of 15 transcription factors showed elevated substitution rates in the human lineage. We further partitioned the TFBSs of the 15 transcription factors into 17 non-overlapping TFBS groups. These non-overlapping TFBS groups included genomic regions exclusively bound by one of the 15 transcription factors and two new TFBS groups: Pol III binding and POU5F1-NANOG binding. The Pol III binding group consisted of TFBSs bound by at least two of BDP1, BRF1, and POLR3G. Similarly, the group of POU5F1-NANOG binding consisted of TFBSs bound by both POU5F1 and NANOG. Then, we applied the group-level LRT again to the 17 non-overlapping TFBS groups and calculated adjusted p-values using the Bonferroni correction.
We utilized the R program for the group-level LRT to perform the element-level LRT. To this end, we applied the phangorn package in R language 86,87 to the alignment of each individual TFBS separately, after filtering out TFBSs with less than 50 informative alignment sites. Then, we performed parametric bootstrapping at the group level to calculate a p-value for each TFBS. Specifically, we first fit the H0 in the group-level LRT to the concatenated alignment of each TFBS group, which provided a global scaling factor to calibrate the branch lengths of the reference phylogenetic model. Second, we randomly sampled 10,000 TFBSs with replacement from each TFBS group and used the calibrated phylogenetic model to generate 10,000 simulated alignments of matched length. Third, we fit the element-level LRT to the simulated TFBS alignments from the same group, which provided an empirical null distribution of the likelihood ratio statistic for each TFBS group. Fourth, we compared the observed likelihood ratio statistic to the empirical null distribution to calculate a p-value for each TFBS. Finally, we fit a beta-uniform mixture model with probability density function (PDF)
$$f(x|a,\, \lambda )=\lambda+(1-\lambda )ato p-values from each TFBS group 44 . We considered a statistic, \(<\hat<\pi >>_>>>\) , from the beta-uniform mixture model as the upper bound of proportion or binding site without acceleration and, accordingly, \(1-<\hat<\pi >>_>>>\) as the lower bound of proportion of accelerated TFBSs.
