Brain-trait-associated variants impact cell-type-specific gene regulation during neurogenesis

Abstract

Interpretation of the function of non-coding risk loci for neuropsychiatric disorders and brain-relevant traits via gene expression and alternative splicing quantitative trait locus (e/sQTL) analyses is generally performed in bulk post-mortem adult tissue. However, genetic risk loci are enriched in regulatory elements active during neocortical differentiation, and regulatory effects of risk variants may be masked by heterogeneity in bulk tissue. Here, we map e/sQTLs, and allele-specific expression in cultured cells representing two major developmental stages, primary human neural progenitors (n = 85) and their sorted neuronal progeny (n = 74), identifying numerous loci not detected in either bulk developing cortical wall or adult cortex. Using colocalization and genetic imputation via transcriptome-wide association, we uncover cell-type-specific regulatory mechanisms underlying risk for brain-relevant traits that are active during neocortical differentiation. Specifically, we identified a progenitor-specific eQTL for CENPW co-localized with common variant associations for cortical surface area and educational attainment.

Keywords: cell-type specificity; common genetic variants; expression/splicing quantitative loci; genome-wide association study; neurogenesis; neuropsychiatric disorders; transcriptome-wide association study.

Figure 1

Study design and cell-type-specific expression (A) Study design illustrating the fetal brain tissue derived cell-type-specific system used to perform eQTL and sQTL analysis. (B) Immunofluorescence of the cells showed that undifferentiated progenitors were SOX2 (in red) and PAX6 (in green) positive, and 8-week differentiated neurons labeled with AAV2-hSyn1-EGFP were positive for EGFP (in green) (scale bar is 100 $\mu$m, DAPI in blue). (C) Principal component analysis of progenitor (purple) and neuron (green) transcriptomes from each donor indicates cell-type-specific clustering. (D) MA plot showing differentially expressed genes in progenitor versus neurons. log₂FC > 0 and adjusted p value < 0.05 indicates genes upregulated in neurons shown in green (neuron up), log₂FC < 0 and adjusted p value < 0.05 indicates genes upregulated in progenitors shown in purple (progenitor up) and genes not significantly differentially expressed between two cell types are shown in gray. Blue lines indicate |log₂FC| > 1.5.

Figure 2

Cell-type-specific eQTL analysis (A) Enrichment of progenitor eSNPs (left) and neuron eSNPs (right) within chromatin states in the fetal brain from chromHMM listed on the y axis. The x axis shows the effect size of enrichment with 95% upper and lower confidence interval and the plot is color-coded based on −log10(p value) value from enrichment analysis. Significant enrichments are shown with an asterisk. Enrichment was tested using eQTLs thresholded at the eigenMT-FDR p value. (B) Comparison of the effects of shared ASE sites and eQTLs in progenitors (left in purple) and neurons (right in green). Nonsignificant ASE sites are shown as darker colors for both cell types, and significant ASE sites are shown as lighter colors. Correlation coefficient (r) values are indicated in colors for each category and the red dashed line indicates y = x. (C) Overlap percentage of cell-type-specific eSNP-eGene pairs shared with fetal bulk eQTLs in progenitors and neurons at m-value > 0.9. Odds ratio (OR) test p values are shown. (D) The fraction of progenitor/neuron primary eGene-eSNP pairs that are true associations (π₁) in fetal bulk eQTLs. 95% upper and lower confidence interval are shown.

Figure 3

Cell-type-specific sQTL analysis (A) Differential splicing of two intron junctions within DLG4. Splice 1 (chr17:7,191,358–7,192,945) supports a previously validated nonsense-mediated decay transcript (ENST00000491753) with higher expression in progenitors, whereas splice 2 (chr17:7,191,358–7,191,893) has higher expression in neurons. (B) A schematic illustrating splicing QTL mapping. Association of variants locating within 200 kb distance from each end of intron junctions were tested. The T allele is associated with more frequent splicing of the shorter intron junction. (C) Two intron junctions supporting an alternative 3′ splicing site for TMEM216 regulated by variant rs11382548 located at the splice site. The regional association of variants to two introns is shown in the genomic tracks on the left colored by pairwise LD r² relative to variant rs11382548, association p values on the y axis, and genomic location of each variant on the x axis. Dashed line indicates significance threshold. Gene model of TMEM216 is shown in the upper right with the position of the variant rs11382548 (closest variant to the splice site), green box indicates the splice site. Boxplots in the lower right show quantile normalized PSI values for splice 1 (chr11:61,397,975–61,398,261) and splice 2 (chr11:61,397,975–61,398,270) at variant rs11382548. (D) Enrichment of cell-type-specific sSNPs within RNA-binding protein (RBP) binding sites based on a CLIP-seq dataset. The significantly enriched RBPs based on −log10(enrichment p value) are listed on the y axis, and the x axis shows the effect size from enrichment test with 95% upper and lower confidence interval, where data points colored by −log10(p value) from the enrichment test and cell-type-specific RBPs are colored with purple for progenitors at the left, and as green for neuron at the right. (E) Overlap percentage of cell-type-specific sSNP-intron junction pairs shared with fetal bulk sQTLs for progenitors and neurons at m-value > 0.9. Odds ratio (OR) test p values are shown.

Figure 4

Colocalization of cell-type-specific eQTLs with GWAS for brain-related traits (A) Number of GWAS loci colocalized with progenitor (purple)- or neuron (green)-specific eQTLs or both cell types (orange). Each GWAS trait is listed on the y axis (SA, surface area; TH, thickness). (B) LD-based overlap of colocalized GWAS loci-gene pairs per trait combinations across progenitor, neuron, and fetal bulk eQTL colocalizations for the traits listed in (A). (C) Genomic track showing regional association of variants with educational attainment (EA), global surface area (GSA), and CENPW expression in progenitors and neurons, −log10 of association p values on the y axis, and genomic location of each variant on the x axis. Progenitor eSNP rs4897179 (3rd row) was coincident with index SNP (rs9388490) for both EA (1st row) and GSA GWAS (2nd row), and conditioning progenitor eSNP rs4897179 on rs9388490 showed colocalization of the two signals (5th row). Also, rs4897179 was colocalized with another variant (rs9388486) located in the chromatin accessibility peak at the promoter of CENPW (6th and 8th rows). Genomic tracks were color-coded based on LD r² relative to the variant rs9388486. Dashed line indicates significance threshold. (D) Plot showing the chromatin accessibility peak (chr6:126,339,531–126,340,960) in progenitors across different genotypes of rs938848. The C allele of rs9388486 disrupted binding motifs of transcription factors including CREM, ATF1, ATF2, and ATF4. (E) Boxplots showing chromatin accessibility across rs9388486 genotypes in progenitors (purple) and neurons (green) (top). Boxplots showing VST normalized CENPW expression across rs9388486 genotypes in progenitors (purple), neurons (green), and fetal bulk (blue) (bottom). (F) A schematic showing that one or more of the implicated transcription factors (TF) has decreased preference to bind at the C allele, which results in lower CENPW expression, increase in global surface area, and educational attainment.

Figure 5

Colocalization of cell-type-specific sQTLs with GWAS for brain-related traits (A) Number of GWAS loci colocalized with progenitor (purple)- or neuron (green)-specific sQTLs or both cell types (orange). Each GWAS trait is listed on the y axis (SA, surface area; TH, thickness). (B) LD-based overlap of colocalized GWAS loci-intron junction pairs per trait across progenitor, neuron, and fetal bulk sQTL colocalizations for the traits listed in (A). (C) Genomic tracks color-coded based on pairwise LD r² relative to the variant rs1222218 showing regional association of variants with SCZ and an unannotated alternative splicing event for ARL14EP in progenitors and neurons, association p values on the y axis, and genomic location of each variant on the x axis. A cryptic exon skipping splice site (chr11:30,323,202–30,332,866) was associated with progenitor sSNP (rs1222218) colocalized with SCZ GWAS index SNP (rs1765142). Dashed line indicates significance threshold. (D) Sashimi plots with the gene model of ARL14EP and the genomic position of the unannotated splice site (blue) overlapping with ARL14EP. Average INT normalized PSI values for the splice site are shown for each genotype group. Schizophrenia risk allele G increases the frequency of the exon skipping event in progenitors. (E) Boxplots showing INT normalized PSI values for splice across rs1222218 genotypes in progenitors and neurons.

Figure 6

Prediction of differential gene expression during human brain development via TWAS (A) Manhattan plots for schizophrenia, IQ, and neuroticism TWAS for progenitors (purple-gray, top) and neurons (green-gray, bottom) where the LD matrix used was based on a European population. Each dot shows −log10 (TWAS p value) for each gene on the y axis, gene names were color coded based on discovery also in colocalization analysis (orange), defined as the nearest gene to GWAS locus (dark pink), being in both these two categories (blue), and discovered only in TWAS analysis (black). Only joint independent genes are labeled (positively and negatively correlated genes represented by triangle and square, respectively, and red line used for TWAS significant threshold). (B) Manhattan plots for IQ TWAS, as described in (A). (C) Manhattan plots for neuroticism TWAS, as described in (A). (D) IQ TWAS results for B3GALNT2, regional association of variants to IQ trait shown at the top, and statistics from each TWAS study shown at the bottom (red line used for genome-wide significant threshold 5 × 10⁻⁸).