Large-scale repressive compartmentalization during hiPSC-to-neuron transition spares schizophrenia risk loci

Functionally mature induced dopaminergic (iDOPA), GABAergic (iGABA) and glutamatergic (iGLUT) neurons were generated from two control donors. Transient overexpression of ASCL1, LMX1B, and NR4A2 (also known as "NURR1") induced iDOPA neurons within 21 days, which were 92% positive for tyrosine hydroxylase, enriched in fetal midbrain dopaminergic neuron gene expression signatures, synthesized dopamine, and exhibited electrophysiologic hallmarks of midbrain dopaminergic neurons by day 35 [30]. Likewise, transient overexpression of ASCL1 and DLX2 induced iGABA neurons within 35 days that were 95-99% positive for GABA itself and glutamic acid decarboxylase 1/2, the majority of which (>75%) were of the SST+ subtype and possessed mature physiologic properties of inhibitory neurons by day 42 [31]. Finally, transient overexpression of NGN2 induced iGLUT neurons within 21 days that were >95% pure glutamatergic neurons, robustly expressed glutamatergic genes, released glutamate, exhibited spontaneous synaptic activity, and recapitulated the impact of psychiatric disease-associated genes [32,33,34,35,36,37,38].

Genome-wide Hi-C contact maps exhibited a greater proportion of intra-chromosomal versus inter-chromosomal contacts across all samples (Supplementary Figure 1), as expected, with the typical distance-dependent decay in contact frequencies (Supplementary Figure 2). The cis:trans ratios of chromosomal contacts were lower than expected (Supplementary Table 1), but neural cells may in fact have lower cis:trans ratios than other tissues [18], and this may be particularly so among in vitro induced neurons [13]. The stratum-adjusted correlation coefficient (SCC) [39] among samples was higher within specific cell types (median SCC of 0.942) than within donors (median SCC of 0.907; one-sided Wilcoxon rank sum test, p = 9.5 ×10) (Supplementary Figure 3), and sample SCC values clustered by cell type (Supplementary Figure 3). Hierarchical clustering with published datasets [12, 40, 41] revealed that hiPSC-derived neurons generally clustered most closely to brain cell types (Supplementary Figure 4A), and correlations to post-mortem NeuN+ neurons were significantly higher than to NeuN- glial cells (Supplementary Figure 4B). Given these findings and the general clustering of samples by cell type independent of donors, data from all donors were combined for each cell type to create higher coverage Hi-C maps, which were subsequently down-sampled to a consistent 185 million cis interactions for further analysis (Supplementary Dataset 1B).

The 3D genomes of hiPSCs and iDOPA, iGABA, and iGLUT neurons were segmented into active ("A") compartments and inactive ("B") compartments using principal component analysis (PCA) on the 100 kb Hi-C correlation matrix, orienting the first eigenvector scores such that positive values were associated with gene-dense regions of the genome as described by the 4D Nucleome Project [41] Chromosomal compartment architectures within a given cell type were reproducible across donors, as the fraction of concordant compartment calls within cell types was significantly greater than between cell types (Wilcoxon rank sum test, p = 9.24 ×10) (Supplementary Figure 5). In order to demonstrate cell-specific chromatin organization, we identified windows with highly variable eigenvector scores. Clustering was performed on the first eigenvector values, i.e., compartment scores, of the donor samples. hiPSCs separated from the three neuronal subtypes on principal component 1 ("PC1"), explaining 39.6% of the variance, and the three neuronal subtypes aligned on PC2, accounting for 18.4% of the variance, with the iGABA samples separating away from the iDOPA and iGLUT samples on this axis as well (Fig. 1A). An illustrative example region was identified around schizophrenia risk genes ANXA2 and RORA (Fig. 1F). Assessed at 250 kb resolution in the merged Hi-C maps, hiPSC-to-neuron differentiation was marked by large-scale A-to-B-compartmentalization (i.e., inactivation) (389-454 Mb, or 12.6-14.7% of total genome sequence, dependent on neuronal subtype; p = 1 x 10 and log odds ratio > 3.4 by Fisher's exact test). Only a minor portion (2.5-3.4% of the genome) converted from inactive B-compartment to active A-compartment status during hiPSC-to-neuron differentiation (Fig. 1B, upper panel; Supplementary Figure 6). In general, developmentally regulated compartment switching was similar across iDOPA, iGABA, and iGLUT neurons, with limited A/B compartment differences between the three neuronal subtypes (A-to-B: 4.6-8.0% and B-to-A: 5.5-7.0% of the genome) (Fig. 1B, lower panel). Compartment status has previously been shown to correlate with transcriptional states [42]. Utilizing RNA-sequencing data from the same cell types from a previous report [30], we observed dynamic changes in A/B compartmentalization during hiPSC-to-neuron differentiation that were broadly associated with transcriptomic changes (Fig. 1C). Regions of the genome undergoing compartment repression (A-to-B) more often encompassed repressed genes, while compartment activation (B-to-A) was more often associated with overexpressed genes; this directional consistency was significant across all neuronal subtypes and neurons as a whole compared to hiPSC (Fisher's exact test, p = 1.6 x 10 to 3.7 x 10).

The developmentally regulated A/B compartment map was superimposed with the most recent GWAS map for schizophrenia [3], comprising 291 common risk loci and 1111 PsychENCODE schizophrenia risk genes. Although the bulk of compartment changes upon neuronal differentiation were from A to B, schizophrenia-associated risk genes and SNPs were depleted in regions that were repressed (A-to-B) in each of the three neuronal subtypes (compared to randomly sampled, GC-content-matched windows, n = 1000, p = 1.0 × 10 to 3.0 x 10) (Fig. 1D; top panel, gray diamonds; "active in iPSCs"). Furthermore, credible schizophrenia risk-associated SNPs were significantly enriched in compartments that were active (A) in hiPSC and iGLUT (p = 4.9 x 10; Fig. 1D; top left panel, red diamond) and compartments repressed in hiPSC but activated in iGABA. We also found regions inactive in hiPSC but active in iGLUT were enriched for risk genes (p = 4.8 x 10) (Fig. 1D, top right panel, blue diamond). An example risk gene and its corresponding risk SNP with active compartments in neurons is shown for the region overlapping the gene DGK1 (Fig. 1G). Interestingly, neuronal subtype-specific A compartments (defined as A in one subtype but B in the other two) were enriched in risk genes (Fig. 1D; bottom right panel, blue diamonds). To assess disease specificity, we calculated Z-scores for overlaps of various GWAS catalog SNP sets [43] (≥100 SNPs) using the same bootstrapping procedure, then converted Z-scores to ranks (decreasing order). Switches showing significant overlap with the Psychiatric Genomics Consortium 3 (PGC3) credible set tended to have similar ranks in other schizophrenia-associated GWAS catalog SNP sets, suggesting a disease-specific enrichment over the unrelated sets (Supplementary Figure 7). Further supporting the biological relevance of compartment switching during in vitro differentiation, differentially expressed genes residing in regions that experienced B-to-A activation during differentiation were associated with biological processes related to neuronal development and functioning, while differentially expressed genes that underwent A-to-B inactivation were associated with processes related to the cell cycle, DNA replication, and organelle biogenesis (Fig. 1E), suggesting that during hiPSC-to-neuron differentiation, compaction of non-essential genes to repressed (B) compartments with associated gene repression is an important part of transcriptional regulation. Taken together, the broad inactivation of chromatin compartments upon neuronal differentiation largely spares schizophrenia risk loci and risk genes, which are instead localized to regions that remain active during neurodevelopment, with additional enrichment of risk loci in regions undergoing compartment activation in iGLUT and iGABA neurons.

In peripheral cells, stripes represent sequentially ordered contacts generated at domain boundaries that are anchored at cohesion docking sites, mark regulatory hubs such as super-enhancers [44], and drive coordinated expression of developmentally programmed gene expression and cell-type identity [45]. The role of stripes in mediating neuronal gene expression has yet to be explored. Stripenn [46] was used to identify donor-merged neuronal stripes from the Hi-C data and revealed 217 iDOPA, 575 iGABA, 193 iGLUT, and 623 hiPSC stripes (Supplementary Dataset 3A). These stripes were collapsed across subtypes into a consolidated set of 1071 which were then scored across each of the donor samples to determine contact enrichments (total observed/expected, O/E, contact frequencies). We performed PCA on the total O/E values for each donor at each stripe, which distinguished neurons from hiPSCs (Fig. 2A). Neuron- and hiPSC-specific stripes were identified by comparing all neuronal subtypes vs. hiPSC. Significantly different stripes were defined as having a nominal p-value less than 0.1 (Student's t-test). A prominent example of a neuronal stripe was observed between the PTPRU risk gene (Fig. 2B). We aggregated brain, U87, and astrocyte-associated super-enhancers from a previous study which broadly screened 86 cell and tissue types [47] and defined this as a "brain super-enhancer" set. The neuron-specific stripe near PTPRU formed a boundary adjacent to two such brain super-enhancers (Fig. 2B).

Neuron-specific stripes (Supplementary Dataset 3B) were enriched in constitutively active A compartments conserved across the hiPSC-to-neuron transition (Fisher's exact test, p = 9.5 x 10) (Fig. 2C). Conversely, hiPSC-specific stripes were enriched in neuron-repressed compartments (hiPSC(A)/neurons(B): p = 1.0 x 10) and constitutively inactive compartments (hiPSC(B)/neurons(B): p = 2.5 x 10) (Fig. 2C). Schizophrenia GWAS SNPs [3] were also overrepresented in neuron-specific stripes (Fisher's exact test, LD > 0.1: p = 4.3 x 10; LD > 0.6: p = 1.7 × 10) (Fig. 2D). To explore their potential functional roles, we examined the 155 hiPSC-specific and 199 neuron-specific stripes for consistency between their contact intensities and genes expressed within them. Unsurprisingly, O/E frequencies were higher in neuronal subtypes for neuron-specific stripes (Fig. 2E). Neuron-specific, differentially expressed genes were found to be significantly associated with neuron-specific stripes, supporting a positive correlation between stripe interaction frequency and gene expression (Fig. 2E, Supplementary Dataset 3C). Likewise, genes upregulated in hiPSCs were more often associated with hiPSC-specific stripes (Fig. 2E, Supplementary Dataset 3D). The directional consistency of differentially expressed genes and differential stripes in which they lied was significant (Fisher's exact test, p = 7.8 x 10). Stripes specific to neurons or hiPSCs contained genes that demonstrated enrichment in biological processes that were relevant to the respective cell type (Fig. 2F). Additionally, we evaluated the extent of enrichment of our neuron-specific stripes and each of the 86 cell and tissue super-enhancers previously identified by Hnisz, et al (2013) [47] (Fisher's exact test) and found the brain-associated set (brain, U87, and astrocytes) to be significantly more enriched than other super-enhancers (Wilcoxon Rank Sum Test, p = 4.11 x 10) (Fig. 2G). No such enrichment was observed in hiPSC-specific stripes.

Finally, we explored the genomic proximity of each of our hiPSC-specific and neuron-specific stripes and found an hiPSC-specific stripe downstream of a schizophrenia risk locus (rs2022265) overlapping a risk gene, SNAP91, whose transcription start site was found to be interacting significantly more than expected with distal upstream regions that varied between neurons and hiPSCs (see Methods, "Significant interactions with specific bin(s)") (Fig. 2H). Taken together, we provide evidence that neuron-specific stripes and their loop anchors are regulatory domains for developmentally relevant gene expression at sites of brain-specific super-enhancers. These stripes emerge predominantly in constitutively active A compartments upon neuronal differentiation and may serve as 3D structural hubs connecting risk loci to distal regulatory targets through neuron-specific loops.

We aimed to identify distal point-to-point chromosomal contacts over the linear genome using a previously developed machine learning algorithm, Peakachu [48]. Using 10 kb resolution to make calls on each of the donor-merged Hi-C contact matrices, we found a wide range of loop calls, with noticeably fewer in neurons, from 625 in iGLUT to 22,211 in hiPSCs, consistent with prior observations [49]. After collapsing overlapping coordinates, we detected 39,963 unique loops (Supplementary Dataset 4A). In order to assess the reproducibility and biological variation in looping structure between cell types and given the low number of calls found in iGLUT, we compared the top 500 loops from each type and found 232 loops shared among all four cell types (Fig. 3A). Additionally, more than half of these top peaks were observed in at least one other cell type-66.8%, 69.2%, 68.8%, and 58.2% for iGABA, iDOPA, iGLUT, and hiPSC, respectively. Using donor-level O/E frequencies at these top 500 loops to perform PCA, iGABA and hiPSC samples clustered by cell type while iDOPA and iGLUT showed less distinct separation (Fig. 3B). We suspect our ability to detect nuanced neuronal subtype-specific loops was diminished given the lower number of cis interactions (185 million). When evaluating our full set of loop calls for cell-type specificity, we found a significant enrichment of neuronal subtype-specific loops within neuronal subtype-specific stripes (binomial test; iDOPA-specific, p = 3.4 x 10; iGABA-specific, p = 1.5 x 10; iGLUT-specific, p = 4.2 x 10; neuron-specific, p = 2.2 x 10) (Fig. 3C).

We first investigated the genes contained within loops overlapping a schizophrenia risk locus on either side of the loop for enrichment in biological process gene ontologies (GO:BP) and found that loops across the three neuronal subtypes targeted genes associated with cell-cell adhesion (Fig. 3D, top panel). Genes contained within loops that did not overlap a risk locus were involved in distinct brain-related terms such as axonogenesis, signal transduction, and synapse organization (Fig. 3D, bottom panel). To uncover potential novel risk genes not identified by linear proximity-based methods and better understand 3D genome regulatory mechanisms of established risk genes, we distinguished two broad sets of schizophrenia GWAS loci-associated loops. The first set, "Risk Gene-Connected Loops," is defined by loops with one anchor positioned at or near an established risk gene proximal to a risk locus and another anchor at a distal, non-coding sequence; we hypothesized that these loops may represent potentially novel gene-regulatory loops impacting previously established schizophrenia risk genes. The second set, "Risk Locus-Connected Gene Loops," is defined by loops with one anchor overlapping a schizophrenia risk locus and another anchor positioned at or near a distal target gene; we hypothesized that these loops may reflect long-range gene-regulatory structures that bring risk loci into 3D spatial proximity to potentially novel risk genes not yet identified by traditional, proximity-based approaches [13]. Genes within Risk Locus-Connected Gene Loops represented new potential risk genes for schizophrenia and were enriched in processes related to homophilic cell adhesion, largely driven by genes from the protocadherin gene family, while risk genes previously identified by linear proximity-based methods were enriched in a variety of processes related to neuronal functioning (Fig. 3E; Supplementary Dataset 4, see Methods). Importantly, this pattern of biological process enrichment was similar to the enrichment seen for genes forming distal chromatin interactions identified from two prior studies [12, 13] (see Methods).

To assess the association of loops with regions of functional interest, we compared the intersection of loops anchored in schizophrenia risk loci and brain enhancers with genes that were and were not differentially expressed from our group's prior RNA-sequencing datasets on the same cell types [30]. Loops specific to iDOPAs (Fisher's exact test, p = 6.1 x 10), iGABAs (p = 4.7 x 10), iGLUTs (p = 4.5 x 10) and neurons overall (p = 6.7 x 10) (versus hiPSCs) overlapped often with fetal brain-specific enhancers previously inferred from segmentation of the genome using chromHMM [50]. (Fig. 3F, top panel). Additionally, iGABA-specific loops were significantly overrepresented in schizophrenia risk genes (p = 5.0 ×10) (Fig. 3F, bottom panel). We did not observe significant enrichment of loops linking schizophrenia risk loci (PGC3, LD < 0.1 or LD < 0.6) with differentially expressed genes, but we suspect our low loop calling sensitivity may be limiting. A representative region is illustrated in Fig. 3G showing two differentially expressed genes, SRR and HIC1, that border the schizophrenia risk locus marked by rs12943566 and rs3752827, where unique combinations of neuronal subtype-specific loops extend across the risk locus (Fig. 3G).

We next sought to experimentally target chromatin loops to generate support for gene-regulatory relationships between the identified loci and distal 3D genomic contacts. We first considered Risk Gene-Connected Loops, loops connecting a schizophrenia risk gene proximal to a GWAS locus to a distal element. We identified a significant interaction unique to iGABA neurons between the promoter region of SNAP91 and a distal non-coding region approximately 150 kb further upstream, which we selected because SNAP91 has a robust brain eQTL [51] and an established role in vesicle-mediated neurotransmitter reuptake at presynaptic terminals [52,53,54,55]. This loop was targeted in iGLUT neurons using dCas9 effectors fused to an abscisic acid (ABA)-inducible dimerization system [56] (SunTag-enhanced chromatin loop reorganization using CRISPR dCas9, "SunCLOUD9") (Fig. 4A, B; see Methods and Supplementary Figure 8A-J for an overview of the technique). Upon presumed loop formation, SNAP91 mRNA expression decreased by ~25% at both DIV14 (Student's t-test, p = 0.026) and DIV21 (Student's t-test, p = 0.029) (Fig. 4C). In the absence of ABA, SunCLOUD9 effectors did not alter SNAP91 expression (Supplementary Figure 10, DMSO -gRNA versus DMSO +gRNA, one-way ANOVA, adjusted p = 0.46). We found suggestive evidence that SunCLOUD9 targeting of the SNAP91 loop decreased the total neurite length per neuron (p = 0.078, r (effect size) = 0.24) (Fig. 4D) without affecting cell number (Supplementary Figure 11A, Student's t-test, p = 1.0). Additionally, SunCLOUD9 targeting of the loop decreased the average number of branch points per unit of neurite length (p = 0.0070, r = 0.35). Moreover, maturation-dependent changes in neuronal activity were disrupted by presumed creation of the SNAP91 loop in iGLUT neurons, as we observed blunted increases in population-wide weighted mean firing rate (WMFR) (β = -0.13 (-0.24 - -0.01), p = 0.032) and network synchrony (β = -0.21 (-0.37 - -0.05), p = 0.010) (Fig. 4E). Finally, network bursting activity was disturbed, as evidenced by more frequent (β = 0.86 (0.28 - 1.45), p = 0.004) but shorter (β = -0.71 (-1.11 - -0.31), p = 0.001) network bursts in the SunCLOUD9 iGLUT neurons (Fig. 4E).

To probe the role of a Risk Locus-Connected Gene Loop involving a long-range interaction between a non-coding schizophrenia GWAS locus and a potentially novel risk gene, we prioritized an iGABA-specific loop that connected a risk locus at rs10957321 to the promoter region of the gene BHLHE22 (Fig. 5A, B). The presence of this loop was associated with increased BHLHE22 expression in iGABA neurons (log(CPM) = 0.49 in iGABA neurons versus log(CPM) = -1.00 in iGLUT neurons, log(CPM) = -5.04 in iDOPA neurons, and log(CPM) = -5.53 in hiPSCs) (Fig. 5B); therefore, we hypothesized that the region containing the risk locus formed an enhancer loop with BHLHE22. BHLHE22 has been implicated in cortical development [57], axon guidance [58, 59], and genetic risk for depression [60], and thus served as an attractive target for validation. Furthermore, the loop anchor overlapping rs10957321 is also enriched in brain tissue enhancer markers [61]. Consistent with a potential role of the loop as an enhancer, presumptive SunCLOUD9-mediated loop-formation in iGLUT neurons resulted in an approximate 80% increase in BHLHE22 expression compared to a scramble-gRNA control (Fig. 5C, Wilcoxon-Mann-Whitney test, p = 0.0059, r = 0.59). SunCLOUD9-iGLUTs had increased neurites per neuron (Wilcoxon-Mann-Whitney test, p = 0.054, r = 0.20) but decreased neurite length (Wilcoxon-Mann-Whitney test, p = 0.00098, r = 0.34) and branch points (Wilcoxon-Mann-Whitney test, p = 0.00022, r = 0.38) (Fig. 5D), with no change in cell number (Supplementary Figure 11B, Student's t-test, p = 0.55). The percentage of electrodes with bursting activity was consistently decreased in SunCLOUD9-iGLUTs (Fig. 5E, β = -0.31 (-0.61 - 0.00), p = 0.049). Moreover, there was a significant interaction between DIV and treatment condition such that SunCLOUD9-iGLUTs in which the BHLHE22 loop was targeted developed increased synchrony at a faster rate (Fig. 5F, β = 0.03 (0.00 - 0.05), p = 0.022).

Altogether, targeting the loops, whether one connecting distal elements to a proximal schizophrenia risk gene (Risk Gene-Connected Loop) or one involving a potentially novel risk gene and a distal risk locus (Risk Locus-Connected Gene Loop), altered gene expression, morphology, and neuronal activity in iGLUT neurons.