Accelerated Evolution of the Regulatory Sequences of Brain Development in the Human Genome

Genetic modifications in noncoding regulatory regions are likely critical to human evolution. Human-accelerated noncoding elements are highly conserved noncoding regions among vertebrates but have large differences across humans, which implies human-specific regulatory potential. In this study, we found that human-accelerated noncoding elements were frequently coupled with DNase I hypersensitive sites (DHSs), together with monomethylated and trimethylated histone H3 lysine 4, which are active regulatory markers. This coupling was particularly pronounced in fetal brains relative to adult brains, non-brain fetal tissues, and embryonic stem cells. However, fetal brain DHSs were also specifically enriched in deeply conserved sequences, implying coexistence of universal maintenance and human-specific fitness in human brain development. We assessed whether this coexisting pattern was a general one by quantitatively measuring evolutionary rates of DHSs. As a result, fetal brain DHSs showed a mixed but distinct signature of regional conservation and outlier point acceleration as compared to other DHSs. This finding suggests that brain developmental sequences are selectively constrained in general, whereas specific nucleotides are under positive selection or constraint relaxation simultaneously. Hence, we hypothesize that human- or primate-specific changes to universally conserved regulatory codes of brain development may drive the accelerated, and most likely adaptive, evolution of the regulatory network of the human brain.


INTRODUCTION
In their seminal work almost 40 years ago, King and Wilson (1975) proposed a key role for regulatory modifications of noncoding DNA in shaping the evolution of our species. Indeed, the human genome was recently discovered to contain noncoding DNA segments that are conserved in other species and some of the sections showed evidence of lineage-specific accelerated changes (Bird et al., 2007;Bush and Lahn, 2008;Lindblad-Toh et al., 2011;Pollard et al., 2006;. These DNA segments, known as human-accelerated elements (HAEs), are genomic regions that are highly conserved throughout vertebrate evolution but are strikingly different across humans, suggesting a regulatory contribution to human-specific traits. Most HAEs are located within noncoding DNA, which perform transcriptional regulatory functions in a specific manner such as promoters, enhancers, insulators, or silencers. The change of DNA sequence in this noncoding element affects the regulatory landscape of gene expression by a loss of function or a gain of function (Spielmann and Mundlos, 2016). Therefore, HAEs in noncoding DNA elements have significant implications for regulatory evolution for human specific traits (McLean et al., 2011).
However, empirical evidence supporting regulatory func-tions has been obtained for only a few of these elements (Kamm et al., 2013;Prabhakar et al., 2008). In contrast to conserved sequences (Bejerano et al., 2004;Pennacchio et al., 2006), it is difficult to systematically test the function of human-acquired sequences in model organisms. Furthermore, whether these changes are the outcomes of directional selection or neutral fixation processes such as GC-biased gene conversion remains uncertain (Duret and Galtier, 2009;Katzman et al., 2010;Kostka et al., 2012;Lindblad-Toh et al., 2011;Pollard et al., 2006;Prabhakar et al., 2009). Epigenomic signatures can mark the location of functional elements and provide systematic information on the spatiotemporal specificity of their regulatory activities (Maurano et al., 2012). DNase I hypersensitive sites (DHSs), which are major epigenetic signatures, are genomic regions that represent loose chromatin and thereby function as regulatory regions by recruiting transcription factors (TFs) . The other crucial epigenetic signatures are histone modifications, in particular monomethylated and trimethylated histone H3 lysine 4 (H3K4me1 and H3K4me3, respectively), which imply transcriptional regulatory functions as promoters and enhancers, respectively (Heintzman et al., 2009). Longrange chromatin interaction controls gene expression within the noncoding fraction of the genome so that it represents a large reservoir of gene regulatory networks (Hnisz et al., 2013;Sandhu et al., 2012). A recent approach revealed that human-accelerated regions show regulatory activity during neural development in an autism spectra disorder model (Doan et al., 2016). However, evaluation of this regulatory function was not based on experimental data as described above but on TF motifs with chromatin interactions only.
In this study, we used a set of epigenetic data comprising transcriptional regulatory functions to validate the regulatory potential of HAEs. In contrast to previous studies which focused on individual experimental validation or sequential evolutionary implications for the functionality of HAEs, our approach was to perform a comprehensive analysis covering whole HAEs to unravel their regulatory function at the level of transcription.

Collection and processing of histone modification data
Eleven datasets for histone modifications in brain samples with fetal and adult origins (fetal brain, germinal matrix, neurosphere ganglionic eminence-derived, neurosphere cortex-derived, substantia nigra, mid frontal lobe, inferior temporal lobe, hippocampus middle, cingulate gyrus, anterior caudate, and angular gyrus) were obtained from the NIH Roadmap Epigenomics data portal. Seven different histone modifications, i.e., H3K4me1, H3K4me3, H3K27me3, H3K-9ac, H3K36me3, H3K9me3, and H3K27ac were included. The Roadmap data for 16 histone modifications in H1-ESCs, H1-ESC-derived mesenchymal stem cells, and H1-ESC-derived neuroprogenitor cells were also obtained. HOMER was run with the "-style histone" option to identify histone modification peaks. The chromosomal coordinates of the histone modification peaks were used for subsequent analyses.

Calculating the frequency of HAE or ultra-conserved element (UCE) overlapping
We overlapped the HAEs and UCEs with the DHS peaks and histone modification peaks using the intersectBed command in BedTools. The number of DHS peaks or histone modification peaks that overlapped with HAEs or UCEs was divided by the total number of DHSs or histone modifications peaks in the given sample, respectively. If the small number of overlapped peaks was too small as compared to the total number of peaks, the overlapping ratio was multiplied by 10,000 to adjust for "overlaps per 10 4 DHSs" or "overlaps per histone modification". By using the shuffle command in BedTools, the same number of DNA segments with the same size distribution as the HAEs or UCEs was randomly captured from each chromosome to generate a false set of HAEs and UCEs. The frequency of DHS overlapping was obtained repeatedly for a set of 1,000 random HAEs and UCEs in the same manner.

Lineage-specific acceleration of DHS sequences
Lineage-specific acceleration was estimated from the DNA sequences of the DHSs. First, phyloP 19 was run to assess the evolutionary significance of base substitutions in the human lineage by using the subtree option and the neutral tree model for the primate subset. The Multiz alignment file and primate neutral model were obtained from the UCSC Genome Browser. The number of significantly (P < 5.0 × 10 -4 ) accelerated nucleotides was divided by the total number of nucleotides contained in all DHSs in a given set (one of the five cell-type groups) and then multiplied by one million, leading to an acceleration estimate as the number of significant base substitutions per Mb. The primate subtree with the mammalian neutral model and the subtree of placental mammals with the vertebrate neutral model were used to estimate primate-specific and mammalian-specific acceleration, respectively. The Multiz alignment file and the neutral tree models were obtained from the UCSC Genome Browser.

Mixed conservation and acceleration of DHS sequences
For each DHS, regional conservation was measured as the average phastCons 24 score for the bases contained in an entire DHS region. The minimum (negative) phyloP 19 score was used as an estimate of outlier point acceleration for each DHS region. The estimates of regional conservation and outlier acceleration were obtained by using the pre-calculated baseby-base scores for the primate clade, which are available in the UCSC Genome Browser. Relative point acceleration was computed by multiplying the regional conservation score and the outlier point acceleration measure for each DHS. We used phastCons and phyloP for the above purposes respectively, because phastCons estimates the probability that each nucleotide belongs to a conserved element (ranging between 0 and 1) whereas phyloP measures conservation (positive scores) and acceleration (negative scores) separately at individual nucleotides and ignores the effects of their neighbors.

Analysis of gene expression in developing brains
Normalized reads per kilobase per million (RPKM) expression values from the RNA-seq data of 578 developing brain samples spanning 13 developmental stages were obtained from http://www.brainspan.org (RNA-Seq Gencode v3c summarized to genes). The expression values for samples from subjects of the same age and from the same area of the brain were averaged. Mean centering of the gene vector was performed before plotting the expression levels of thyroid hormone receptor beta (THRB) and snail family transcriptional repressor 2 (SNAI2).

Identification of TF motifs
A total of 379 TF motifs from the transcription factor database (TRANSFAC) and 26 additional motifs from the JASPAR database were obtained, leading to 405 unique motifs in Homo sapiens. The FIMO search tool (http://meme-suite.org/ doc/fimo.html?man_type=web) was used to search for within-DHS motifs at the threshold P value of 10 -5 .

Significant association of HAEs with epigenomic signatures
To obtain a comprehensive landscape of epigenomic signatures for unraveling regulatory potential of whole HAEs, we compiled 2,745 HAEs and overlapped them with DHSs with different developmental tissues from different origins. Sixty-six percent of the HAEs coincided with the examined DHSs ( Supplementary Fig. S1). Remarkably, the mean frequency of overlaps between HAEs and DHSs across fetal brains and spinal cords was higher than the mean frequency of overlaps across adult brains, non-brain fetal tissues, ESCs, and non-brain adult cells ( Table S1). The observed variations among samples with the same tissue origins (Fig. 1A) may have resulted from experimental errors or may reflect variable biological features of the sample donors. The higher number of fetal brain DHSs that overlapped with HAEs compared to adult brain DHSs indicates that human-specific evolutionary pressure would have led to the evolution of functional regulatory circuitry during human brain development. In addition, the HAEs were highly enriched in the histone markers H3K4me1 and H3K4me3 across different brain samples (Fig.  1B, Supplementary Table S2), which implies they function as active enhancers and active promoters, respectively. Specifically, the mean overlap frequency of HAEs with histone markers of H3K4me1 and H3K4me3 in fetal brains was higher than those in adult brains ( Supplementary Fig. S4). These results imply that HAEs have functional potential, which is explicitly associated with fetal brains.

Human-specific regulatory evolution of fetal brain DHSs
The HAEs were analyzed to determine whether they are under specific evolutionary pressure in the context of their functional implications. To provide a unified and quantitative measure of human-specific regulatory evolution, we computed the sequence acceleration because the divergence of humans and chimpanzees is based on a primate neutral model . We then examined how much each DHS deviated evolutionarily from chimpanzees to humans by calculating the frequency of significantly accelerated nucleotides for the relevant DHS. The calculated frequency was highest in fetal brain DHSs among DHSs with other tissue origins ( Fig. 2A, left) and was particularly high in DNA segments accessible only in fetal brains ( Fig. 2A, right). A similar magnitude of primate-specific acceleration and no considerable evidence supporting acceleration in the mammalian clade were found ( Fig. 2A, right), implying that most changes to the brain developmental sequences occurred after the split of primates and humans from other mammalian lineages. One of the biggest differences between primates and mammals is the size of the brain (DeCasien et al., 2017). Primates including humans have relatively large brain sizes (neocortex and cerebellum), as compared to other mammals, which results in the higher intelligence of primates (Barton and Venditti, 2014). Also, the higher relative cortex volume and neuron packing density of primates allow for more cortical neurons than other mammals with the same brain size (Roth, 2015). The accelerated nucleotides which we examined are reported to be related to the function of transcriptional enhancers during nervous system development and to genes associated with unique human features, such as complex language (Caporale et al., 2019;Kamm et al., 2013). From this perspective, our analyses and the supporting studies indicate that fetal brain DHSs are under human-specific evolutionary pressure in the context of brain development.

Coexistence of conservation and human-specific evolution in regulating brain development
We also hypothesized that regulatory regions of fetal brain could have evolutionarily conserved signatures due to their universal importance in brain development (Lu et al., 2019).   H e a rt − 9 6 H e a rt − 1 0 5 H e a rt − 1 1 0 H e a rt − 1 4 7 K id n e y − 8 5 K id n e y − 9 1 K id n e y − 1 0 5 K id n e y − 1 1 7 K id n e y − 1 2 2 K id n e y − 1 4 7 L u n g − 8 5 L u n g − 9 6 L u n g − 1 0 1 L u n g − 1 1 7 L u n g − 1 2 2 M u s c le H e a rt − 9 6 H e a rt − 1 0 5 H e a rt − 1 1 0 H e a rt − 1 4 7 K id n e y − 8 5 K id n e y − 9 1 K id n e y − 1 0 5 K id n e y − 1 1 7 K id n e y − 1 2 2 K id n e y − 1 4 7 L u n g − 8 5 L u n g − 9 6 L u n g − 1 0 1 L u n g − 1 1 7 L u n g − 1 2 2 M u s c le   indicate highly conserved genomics regions across multiple species, by filtering out overlapping HAEs and mapping them to DHSs with different tissue origins. Paradoxically, fetal brain DHSs were highly enriched not only in HAEs but also in UCEs (Fig. 2B, Supplementary Fig. S5). For example, we discovered a fetal brain DHS containing an HAE and a UCE together (Fig. 3). Two TF-binding motifs, one for THRB (P value: 7.0 × 10 -7 ) in the HAE and the other for SNAI2 (P value: 4.0 × 10 -6 ) in the UCE, were identified as the most significantly enriched motifs. This region was accessible in all fetal brain samples as well as in samples from 2-and 3-year-old medulloblastoma patients but not in any adult brain samples. Although SNAI2 is expressed throughout the prenatal period, THRB is specifically expressed in the mid-and late-prenatal stages and during early infancy (Fig. 4), in agreement with the role of thyroid hormone during early brain development (Morreale de Escobar et al., 2004). Fetal brain DHSs significantly interacted with the promoter of the nearest gene FAXC, which was located 163 kb away (Fig. 3A). Mutations in the Drosophila homolog of FAXC cause failed axon connections (Hill et al., 1995). The proximal DHS of FAXC is also accessible in fetal brains and is in linkage disequilibrium with an SNP (rs2132683) that is a strong genetic determinant of brain structure (Stein et al., 2010). Therefore, the constitutive and developmental stage-specific regulation of FAXC by SNAI2 and THRB may contribute to conserved and human-specific axon connections, respectively, during brain structural development. Finding that HAEs and UCEs coexist on the same DHSs indicates the high variability of sequence conservation across a developing brain DHS.

High relative point acceleration of fetal brains
To assess whether this was a general pattern, we measured evolutionary conservation of the DNA sequences of all DHSs based on multiple alignments in primates. The degree of selective constraint over an entire DHS region, or the regional conservation, was estimated as the average phastCons score (Siepel et al., 2005) for all nucleotides in a given region. On the contrary, the magnitude of primate acceleration was measured based on the phyloP score  of the most accelerated nucleotide in a given DHS region and was termed "outlier point acceleration." Fetal brain DHSs showed the highest levels of regional conservation (Fig. 5A) and outlier point acceleration (Fig. 5B) simultaneously, suggesting that DNA sequences that regulate brain development were selectively constrained in general but that particular nucleotides also underwent positive selection or constraint relaxation at the same time. Regulatory sequences in the adult brain also showed accelerated signatures (Fig. 5B), but their conservation levels were low (Fig. 5A). We estimated the relative functional importance of accelerated substitutions by considering the conservation level of background sequences for each DHS region (see Materials and Methods section). This measure, termed "relative point acceleration," was highest for fetal brains and lowest for adult brains (Fig. 5C). In other words, sequence substitutions in adult brain DHSs generally occurred in less-constrained regions, suggesting that relaxed constraints rather than positive selection likely resulted in point acceleration. Meanwhile, the DNA sequences of non-brain fetal DHSs were relatively well conserved (Fig.  5A); however, the relative point acceleration levels were high (Fig. 5C), which likely explains the relatively high acceleration of kidney, lung, and muscle DHSs (Figs. 1A and 2A). Taken together, the regulatory sequences in developing brains are subject to strong evolutionary constraint with outlier point acceleration, which highlight the functional potential of these sequence changes.

DISCUSSION
The regulatory implications of human-accelerated noncoding elements have remained largely unexplored, despite their potential to reveal molecular mechanisms underlying human-acquired traits. The main focus of previous studies has been to identify sequence features that can distinguish between positive selection and biased gene conversion, which promotes

Ultraconserved element
Human-accelerated element SNAI2 THRB Figure 3 A B Fetal brain DHS Fig. 3. Transcriptional regulation of coexistence of UCE and HAE in fetal brain DHS region. (A) Chromatin-interaction landscape between a fetal brain DHS containing an HAE and UCE and the promoter of FAXC in dorsolateral prefrontal cortex tissue. Hi-C interaction frequency maps were plotted using the 3DIV database, available at http://kobic.kr/3div/ (Yang et al., 2018). The location of the fetal brain DHS is indicated by the red dot. The green line indicates the cut-off for the distance-normalized interaction frequency. (B) A brain DHS that contains an HAE and UCE simultaneously and that is active only during prenatal periods and early infancy. Shown are the locations of the HAE and UCE, primate phastCons (sky blue) and phyloP (blue for conservation and red for acceleration) scores, SNAI2 and THRB binding motifs, and the frequency of each base among primates with the human reference sequences at the bottom. For the two motifs, THRB was identified as a single motif with statistical significance (P value: 7.0 × 10 -7 ) and SNAI2 showed the highest statistical significance (P value: 4.0 × 10 -6 ) among the discovered motifs with biological relevance to brain development.     the fixation of neutral or weakly deleterious mutations at recombination hotspots (Duret and Galtier, 2009;Katzman et al., 2010;Kostka et al., 2012;Lindblad-Toh et al., 2011;Pollard et al., 2006;Prabhakar et al., 2009). In this work, we used a different approach to examine epigenetic regulatory signatures and discovered biased evolutionary acceleration in the genomic regions that become accessible only in specific tissues and at certain developmental stages, lending support to the idea that directional selection, rather than neutral processes, likely play a role in shaping the landscape of human-accelerated evolution. Additionally, we found strong evolutionary constraints in the regulatory sequences of developing brains, indicating that failure of the normal development of the central nervous system may be universally detrimental to the fitness of organisms. Acceleration within deeply conserved regions highlights the functional importance of these sequence changes. In conclusion, our results elucidate the regulatory implications of human lineage-specific genetic alterations and will likely facilitate further related studies.
Note: Supplementary information is available on the Molecules and Cells website (www.molcells.org).