Resolving the three-dimensional interactome of Human Accelerated Regions during human and chimpanzee neurodevelopment

Abstract

Human Accelerated Regions (HARs) are highly conserved across species but exhibit a significant excess of human-specific sequence changes, suggesting they may have gained novel functions in human evolution. HARs include transcriptional enhancers with human-specific activity and have been implicated in the evolution of the human brain. However, our understanding of how HARs contributed to uniquely human features of the brain is hindered by a lack of insight into the genes and pathways that HARs regulate. It is unclear whether HARs acted by altering the expression of gene targets conserved between HARs and their chimpanzee orthologs or by gaining new gene targets in human, a mechanism termed enhancer hijacking. We generated a high-resolution map of chromatin interactions for 1,590 HARs and their orthologs in human and chimpanzee neural stem cells (NSCs) to comprehensively identify gene targets in both species. HARs and their chimpanzee orthologs targeted a conserved set of 2,963 genes enriched for neurodevelopmental processes including neurogenesis and synaptic transmission. Changes in HAR enhancer activity were correlated with changes in conserved gene target expression. Conserved targets were enriched among genes differentially expressed between human and chimpanzee NSCs or between human and non-human primate developing and adult brain. Species-specific HAR gene targets did not converge on known biological functions and were not significantly enriched among differentially expressed genes, suggesting that HARs did not alter gene expression via enhancer hijacking. HAR gene targets, including differentially expressed targets, also showed cell type-specific expression patterns in the developing human brain, including outer radial glia, which are hypothesized to contribute to human cortical expansion. Our findings support that HARs influenced human brain evolution by altering the expression of conserved gene targets and provide the means to functionally link HARs with novel human brain features.

Figure 1:. Overview of the study.

(A) Schematic illustrating how HARs and HGEs are defined. (B) Outline of the datasets and analyses used to identify and annotate HAR and HGE interactions as described in the Results. Schematics were generated using BioRender.com.

Figure 2.. Conserved and species-specific HAR and HGE interactions in human and chimpanzee NSCs.

(A) Left. Proportion of conserved interactions for HARs and their chimpanzee orthologs. The pie charts show the breakdown of HARs and chimpanzee orthologs by the percent of conserved interactions, in quartiles as shown in the legend. The violin plot shows the overall distribution of conserved interactions across HARs and their orthologs, and indicates the median value as well as the interquartile range. Right. Proportion of conserved interactions for HGEs and their orthologs, as shown for HARs on the left. (B) CHi-C interaction profiles of HACNS240 (top) and its chimpanzee ortholog (bottom). Left. Each point on the scatterplot represents a significant interaction, plotted by the degree of significance as shown on the y-axis. The relative distance from the HAR is given on the x-axis. A normalized density distribution of the interactions is shown below the scatterplot, and regions of high interaction density are marked by orange boxes (Methods). Right. Schematized UCSC Genome Browser view (in GRCh38 coordinates) for each interaction profile. The location of HACNS240 and its ortholog are shown in purple above the gene models, and the regions of high interaction density are shown in orange. Representative looping events are shown in gold. Protein-coding genes are shown in blue and lncRNA genes are shown in green.

Figure 3.. Functional annotation of HAR and HGE interactions in human and chimpanzee NSCs.

(A) Conserved and species-specific gene targets of HARs, HGEs and their chimpanzee orthologs based on the human GENCODEv43 annotation, as defined in the Results and Methods. The horizontal bar plot further subdivides human-specific, chimpanzee-specific, and conserved targets into protein-coding and long non-coding RNA (lncRNA) genes. (B) GO Biological Process enrichment analysis performed on the conserved set of protein-coding gene targets, using the human GENCODEv43 annotation. (C) Enrichment of conserved HAR or HGE gene targets within differentially expressed gene (DEG) sets called between hNSCs and cNSCs and their relative overrepresentation within the DEG set compared to species-specific gene targets (Methods). P-values were computed based on permutation using random sampling of the background followed by Bonferroni correction (n = 20,000 trials; Methods). (D) HAR or HGE gene target enrichments in gene co-expression modules previously identified using WGCNA in human fetal brain. Modules start with the identifier ‘HS’ as they are called in human samples. The circles are colored according to significantly enriched GO Biological Processes as shown in the legend (Methods). The enrichment adjusted P-value for each module was calculated based on a permutation test with n = 20,000 trials followed by Benjamini-Hochberg (BH) correction. (E) Doughnut plot showing the distribution of HAR and HGE interactions based on gene and phastCons constrained noncoding element annotations as described in the Results and Methods. (F) Bar plot showing the proportion of conserved and species-specific interactions for HARs, HGEs and their chimpanzee orthologs overlapping putatively functional elements based on H3K27ac, H3K27me3, CTCF, and RAD21 profiles generated in this study and ENCODE cCREs.

Figure 4.. HARs and HGEs are significantly associated with gene expression differences between human and chimpanzee.

(A) Scatterplot of relative expression for all HAR and HGE gene targets (log2(human TPM/chimpanzee TPM)) against the relative H3K27ac level at the HAR or HGE (log2(human/chimpanzee)). Differentially expressed (DE) gene targets contacted by elements showing differential H3K27ac marking (DM) are plotted as green triangles. The polynomial of best fit (y ~ O(x³)) is plotted in blue and standard error is shaded in light blue; see Results and Methods for details. The inset shows the expected distribution of HARs and HGEs differentially marked by H3K27ac that also target a differentially expressed gene as obtained by permutation, compared to the observed value shown by the vertical dotted line. The enrichment P-value was computed based on random sampling of the background (n = 20,000 trials; Methods), * = P < 0.05. (B) Two examples of HAR gene targets whose species-biased differential expression is associated with a specific enrichment of H3K27ac (in green) or H3K27me3 (in blue) at the HAR and the target gene promoters as described in the Results. Fold enrichment and fold change (FC) ratios are shown on a log2 scale, and significant differences in marking and expression were identified using DESeq2 (* = BH-corrected P < 0.05). (C – D) Enrichment of HAR and HGE gene target sets within (C) human-specific DE genes and (D) human-specific differentially accessible region (CRE)-linked genes, respectively, as described in the Results. P-values were computed based on permutation using random sampling of the background (n = 20,000 trials; Methods), * = Bonferroni-corrected P < 0.05. The observed intersection of these gene sets with HAR and HGE targets in hNSCs is shown by the dotted line. (E) Enrichment of HAR and HGE gene target sets within differentially expressed genes identified in a comparison of human and rhesus macaque fetal brain as described in the Results. The DEGs were identified as differentially expressed in at least one brain region. P-values were computed as in (C – D) with * = Bonferroni-corrected P < 0.05, and display of observed intersection with hNSC gene targets is shown by the dotted line as in (C – D). HAR and HGE targets were combined for these analyses. Illustrations in these panels were generated using BioRender.com.

Figure 5:. Activation of repressed HAR and HGE interactions upon neuronal differentiation.

(A) The number of H3K27me3-marked interactions in hNSCs that maintain H3K27me3 marking in human neurons (black), interactions that switch to an H3K27ac-marked state in neurons (white) and interactions that lose H3K7me3 at either the HAR, HGE or the gene target (grey). (B) An example of a repressed interaction that is activated upon neuronal differentiation. HAR116 engages in an H3K27me3-marked interaction with its gene target TCF20 in hNSCs, which switches to an H3K27ac-marked state in neurons. The curved golden line indicates the HAR116-TCF20 interaction detected via CHi-C, and H3K27me3 (in blue) and H3K27ac (in green) profiles are shown for both hNSCs and iPSC-derived human neurons. Dashed black boxes highlight the location of the HAR and the TCF20 promoter relative to the H3K27me3 and H3K27ac peaks. (C) Comparison of of TCF20 and CDH23 expression in hNSCs and neurons (TPM values, P_BH < 0.05 by Mann-Whitney test). (D) Example of a repressed interaction in NSCs that is maintained in neurons, shown as in (B). (E) Whisker plots with pairwise comparisons of distributions of gene target expression for each category of H3K27me3-marked interactions across hNSCs (in teal) and neurons (in purple). Each box represents the 25^th to 75^th percentiles of the data, and the median of each distribution is marked with a solid line. The median increase in log2(TPM+1) between neurons and hNSCs for each of the three categories is denoted by Δ. P values for differences in expression between NSCs and neurons in each category were calculated using a Mann-Whitney U test (* = BH-corrected P < 0.05).

Figure 6.. HAR and HGE gene targets in NSCs show cell type-specific expression profiles in the fetal human brain.

(A) UMAP showing neural progenitors in scRNA-seq from eight embryonic and fetal human brain regions colored by regional identity with cortical outer radial glia (oRGs) and intermediate progenitor cells (IPCs) highlighted in yellow and silver, respectively. Heatmaps of clustered average scaled expression profiles across 29 neural progenitor subtypes, separated into gene sets that show cell-type specific expression in cortical oRGs and IPCs. CX – cortex; HC – hippocampus; HT – hypothalamus; CB – cerebellum; TH – thalamus; LGE – lateral ganglionic eminence; MGE – medial ganglionic eminence; and CGE – caudal ganglionic eminence. (B) Density plots for a subset of genes showing cell-type specific expression in progenitor cell types, visualized on the UMAP shown in A. (C) Visualization of the average expression profiles of 38 genes showing higher expression in hNSCs compared to cNSCs that share oRG-biased expression profiles.

Figure 7:. HAR and HGE gene targets in neurons show biased expression across multiple neuronal subtypes in the developing human brain.

(A) UMAP showing neurons in scRNA-seq from eight embryonic and fetal human brain regions colored by regional identity with rostral thalamic neurons (RTN) highlighted in yellow. Heatmaps of clustered average scaled expression profiles across 33 neuronal subtypes, separated into gene sets with biased expression in groups of brain regions (thalamus/hypothalamus/hippocampus and hippocampus/cortex) and specific neuronal subtypes (rostral thalamic neurons). CX – cortex; HC – hippocampus; HT – hypothalamus; CB – cerebellum; TH – thalamus, LGE – lateral ganglionic eminence; MGE – medial ganglionic eminence; and CGE – caudal ganglionic eminence. The neuron subtype abbreviations are provided in Table S3. (B) Density plots of expression in neuronal cell types (visualized on the UMAP shown in A) display the regional or cell type-specific biases of HAR and HGE gene target expression within each cluster.