Wnt activity reveals context-specific genetic effects on gene regulation in neural progenitors

Gene regulatory effects in bulk-post mortem brain tissues are undetected at many non-coding brain trait-associated loci. We hypothesized that context-specific genetic variant function during stimulation of a developmental signaling pathway would explain additional regulatory mechanisms. We measured chromatin accessibility and gene expression following activation of the canonical Wnt pathway in primary human neural progenitors from 82 donors. TCF/LEF motifs, brain structure-, and neuropsychiatric disorder-associated variants were enriched within Wnt-responsive regulatory elements (REs). Genetically influenced REs were enriched in genomic regions under positive selection along the human lineage. Stimulation of the Wnt pathway increased the detection of genetically influenced REs/genes by 66.2%/52.7%, and led to the identification of 397 REs primed for effects on gene expression. Context-specific molecular quantitative trait loci increased brain-trait colocalizations by up to 70%, suggesting that genetic variant effects during early neurodevelopmental patterning lead to differences in adult brain and behavioral traits.

Supplementary Figure 2| Effects of CHIR and WNT3A on proliferation and Wnt activity Cellular assays show that WNT3A and CHIR stimulation increase canonical Wnt pathway signaling and hNPC proliferation. Diagram of luciferase reporter assay measuring β-catenin mediated Wnt pathway activation by CHIR or WNT3A (A). Effects of 48h CHIR (B) or WNT3A (C) exposure on Wnt pathway activation, reported as the log of luciferase luminescence (activated by Wnt stimulation) normalized by renilla luminescence (from a constitutively active reporter transgene). Representative flow cytometry scatter plots from proliferation assays depict newly synthesized DNA from a 2h hr EdU pulse vs total DNA content following 48h vehicle (left) or 2.5uM CHIR (right) exposure (D). Percentage of cells in S-phase (%EdU+) exposed to vehicle or increasing doses of CHIR for 48h as measured by flow cytometry in a subset of cell-lines (E), and in all cell-lines used for ATAC-seq and RNA-seq in this study (F). Representative immunocytochemistry images (100x magnification) from proliferation assay of hNPCs following 48h vehicle (left) and 5nM WNT3A (right) exposure (G). Green (GFP) labels cells in S-phase during the EdU pulse, and blue (DAPI) stains all nuclei. Percentage of cells in S-phase (%EdU+) exposed to vehicle or increasing doses of WNT3A for 48h as measured by flow cytometry in a subset of cell-lines (H). CHIR and WNT3A concentrations (boxed) that maximize Wnt activation and proliferative responses were used in this study for ca/eQTL mapping.
Supplementary Figure  Percentage of increased (shrunken log2FC > 0.5 in dark blue, shrunken log2FC > 0 in blue) or decreased (shrunken log2FC < -0.5 in red, shrunken log2FC < 0 in pink) accessibility or expression are shown for each category. Nonsignificant peaks and genes (FDR-adjusted P >=0.1) are colored gray. Greater changes were observed in both chromatin accessibility and gene expression for CHIR, a potent Wnt activator, as compared to WNT3A, an endogenous Wnt ligand. Simultaneous WNT3A activation and inhibition via XAV as compared to WNT3A activation yielded few differentially expressed genes and chromatin elements, likely due to the low sample size. We also note that 93,853 peaks and 9,890 genes were filtered out (adjusted P value were set to NA) by DESeq2 2,3 as those have low mean read counts.
Supplementary Figure 7| Wnt-pathway related transcription factor binding sites are enriched in upregulated WREs The canonical Wnt-pathway downstream effectors, ß-catenin, LEF1, and TCF7L2 binding sites previously identified by ChIP-seq experiments in HEK cells 1,4 are more often overlapped with WREs opening due to WNT3A (A) or CHIR (B) than closing WREs. The figures show the percentage of overlaps with the indicated TF binding sites in either the 2,000 most upregulated peaks or the 2,000 most downregulated peaks based on shrunken logFC. Statistical significance was estimated using Fisher's exact test on the number of overlaps.
Supplementary Figure 8| Stimulus-specific regulatory element-gene correlation We identified stimulus-specific regulatory elements which regulate gene expression by correlating expression with chromatin accessibility. The number in the plot indicates significant and positively correlated peak-gene pairs. In total, 12,643 peak-gene pairs are detected only under the stimulus condition (highlighted in red), demonstrating new stimulus-specific regulatory elements.
Supplementary Figure  Supplementary Figure 10| ca/eQTL effect sizes compared to previous data We compared the effect size of ca-QTL (A) and eQTL (B) in vehicle condition to previous studies using the same hNPC cell lines with slightly larger sample size 5,6 . Strongly correlated effect sizes for unstimulated ca/eQTLs were found between the two datasets. In part A, the current study beta values are derived using RASQUAL as a scaled allelic ratio ranging from -0.5 to 0.5, which is not on the same scale as the beta values from Liang et al., which are calculated using a linear mixed effect model. One SNP showed opposite effect directionality across the caQTL datasets between the current and previous datasets (rs2076179). We suspect this was due to an error in allele assignment for that SNP in the previous analysis caused by sampling differences incorrectly changing the labeled minor allele for a SNP with allele frequency close to 0.50. We corrected this error for the one SNP prior to plotting.
Supplementary Figure 11| eQTL-caQTL overlaps Overlaps between caQTLs and eQTLs were called where a caSNP was within 1Mb of an eGene and at least 1 significant SNP in both datasets were in LD r 2 ≥ 0.8. The percentage of eQTLs (A) and caQTLs (B) overlapped within vehicle and WNT stimulus conditions are shown. The effect size of caQTL-eQTL sites which shared the same SNP position were compared within vehicle (C) and WNT stimulus conditions (WNT3A and CHIR shown in D and E respectively). SNPs selected to influence chromatin accessibility had relatively little overlap with those influencing gene expression likely because genetic variation affecting chromatin accessibility often does not lead to changes in gene expression (about <2% of caQTLs had a shared eQTL based on LD overlap (r 2 >= 0.8). However, SNPs selected to influence gene expression more often also influence chromatin accessibility in that~27% of eQTLs have a shared caQTL within condition, a comparable number observed in our previous study (34.9% caQTL-eQTLs overlap 5 ).

Supplementary Figure 12| Effect size differences across conditions
We compared effect sizes of index caQTLs (A-B) / eQTLs (C-D) in stimulated conditions versus vehicle conditions (WNT3A or CHIR in (A,C) and (B,D), respectively). Dots in red indicate significant interaction effects were observed (r-QTLs). Error bars are standard errors of beta. We observed lower correlation in r-QTLs as compared to non-r-QTLs. We note that some SNPs are not tested for a baseline model due to low expression thus could not be plotted because no beta value was calculated.

Supplementary Figure 13| Response and Non−Response caPeaks enrichment within chromHMM Annotations
Enrichment of response and non-response caPeaks within chromHMM states defined within fetal brain male (E081) 7 . *, **, *** indicate enrichments of overlaps relative to the entire genome, evaluated with a binomial test, similar to the GREAT test 8 , with P-values < 0.05, .01 and .001, respectively. We observed a significant enrichment of both response and non-response caPeaks within promoters, enhancers, and bivalent regions. TssA, TssFlnk, and Enh states showed a significant difference in overlap counts between response and non-response caPeaks (P = 3.9x10 -5 ; 2.8x10 -4 ; 7.4x10 -7 , respectively), as evaluated with Fisher's exact test. There was significantly less enrichment of response caPeaks in active TSS and enhancers as compared to non-response caPeaks, perhaps because these response caPeaks flag novel condition specific enhancers not annotated in post-mortem fetal brain tissue.
Supplementary Figure  Supplementary Figure 19| Stimulus-specific ENO4 eQTL colocalizing with regional cortical surface area GWAS (A) Chromatin accessibility (chr10:116973711-116975130; left) and ENO4 gene expression regulated by rs11197861. eQTL of ENO4 was detected only in the stimulated condition (A,  B). The eQTL is colocalized with several brain-related traits including surface area of Insula (B, C). (B) P values from Insula surface area GWAS, eQTL, caQTL for vehicle, CHIR condition are shown respectively. ETV2:FOXl1 motif is disrupted by the rs11197861-T located in the peak, which may result in decreased chromatin accessibility and ENO4 gene expression. (C) Z scores from regional cortical surface area 9 . In addition to Insula, surface area of precentral and precuneus are also associated with this SNP (indicated by asterisks; p < 5 x 10 -8 ).
Supplementary Figure 20| Optimizing control for known and unknown technical confounding in eQTL We tested different numbers of PC variables to determine the number of PCs to include for correcting expression values prior to running limix_qtl. (A) The plot shows the number of discovered eGenes with respect to the number of PCs at FDR-adjusted P < 0.1. and the number of overlapping found eGenes across models (B-D). To identify independent eQTLs, we repeated eQTL mapping by including the index eSNP in the association model until no SNP passed the genome-wide threshold (FDR-adjusted P < 0.1). In table (E), we show the number of SNPs discovered in each round of conditional associations.