Non-demented control (NDC), PSEN1A79V, PSEN2N141I, and APPV717I hiPSCs were differentiated into CD44−/CD184− neurons as previously described [11,22] (Fig. 1A; Supplementary Fig. 1A-B). RNA-seq and subsequent differential gene expression analysis identified a substantial number of differentially expressed genes (DEGs) relative to NDC in all three mutations (Fig. 1B), with 1339 common DEGs (Fig. 1C). Rank Rank Hypergeometric Overlap (RRHO) analysis revealed the strongest similarity overlap between the PSEN1A79V and PSEN2N141I mutations and the weakest similarity overlap between the APPV717I and PSEN2N141I mutations (Supplementary Fig. 1C). Next, we sought to determine the disease endotypes associated with these three mutations. To this end, we first carried out the fgsea and CERNO enrichment tests with the GO Biological Process and Hallmark ontology databases to identify common and distinct disease endotypes. This revealed positive enrichment of genesets related to activation of the cell cycle and dedifferentiation to non-ectoderm lineages (e.g., EMT), common across all mutations (Fig. 1D, Supplementary Fig. 1D). Interestingly, neuronal maturation and neuron function genesets (e.g., synaptic signaling) were negatively enriched in PSEN1A79V and APPV717I, these programs were modestly upregulated in PSEN2N141I. We next used curated gene lists for the key endotypes[11] observed here and performed tSNE clustering for pseudo-trajectory analysis of the three mutations studied here and two other mutations with an earlier AAO (PSEN1H163R, 42–47 years [74] and PSEN1A431E, 36–53 years [9]). This approach demonstrates that while PSEN1H163R has the most severe dysregulation across all endotypes, PSEN1A79V is similar to PSEN1A431E; in contrast, APPV717I and PSEN2N141I are comparatively less severe (Supplementary Fig. 2). To characterize the transcriptional regulation of these disrupted gene programs, we used ISMARA [39] (Fig. 1E) and DoRothEA [40] (Fig. 1F) to predict TF activity. ISMARA identified several regulators that are common with significant differential activities across all mutations associated with key endotypes, including early neuron lineage (NFIB, PAX6, NEUROD1, MEIS2) [75–78], axonal growth and synaptogenesis (CREB1) [79,80], mitochondrial energy and neuron function (NRF1) [81,82], non-ectoderm lineage (TEAD1/3) [83,84], pluripotency (GATA3 and TAF1), cell cycle (E2F1/7, FOXM1, ID4) [85], and inflammation (PRDM1) [86]. The neural differentiation repressor REST [87] was particularly activated in PSEN1A79V compared to APPV717I and PSEN2N141I. TF regulon analysis using DoRothEA revealed similar differential activity of TFs related to cell cycle (FOXM1, MYC, E2F1/4), neuron lineage and function (ASCL1, REST, EGR1) [88,89], neuron mitochondrial energy production (NRF1), and non-ectoderm lineage (TEAD2, TCF3, TFAP2C) [90]. Next, we performed co-expression module detection using the CEMItool [51] R package followed by module enrichment in each mutation relative to NDC using fgsea. Nine functional co-expression modules were detected with modules 1, 3 and 4 significantly enriched in all three mutations with a positive activity. In contrast, module 5 was enriched with a negative activity (Fig. 2A). Hypergeometric enrichment revealed an over-representation of genes associated with cell cycle, inflammation, non-ectoderm lineage, and early-stage neurogenesis (Fig. 2B-C). By integrating protein-protein interaction (PPI) and TF-target gene edges for module genes, neighboring genes, and key module TFs, key centroid genes for each module can be identified, such as the lineage regulators ASCL1 and ZIC2 for module 3 (Fig. 2D, Supplementary Fig. 3).
To determine whether endotype transcriptional changes are driven by modulation of chromatin topology, we performed ATAC-seq (Fig. 3A, Supplementary Fig. 4) and assessed differentially accessible regions (DARs) within promoter or PEREGRINE [61] enhancer regions for each mutation relative to NDC (Fig. 3B-C). Next, we sought to identify TFs with differential activity associated with chromatin accessibility by motif footprinting using HINT-ATAC and motif enrichment with GimmeMotifs maelstrom. HINT analysis using the CIS-BP motif database identified decreased footprinting activity of TFs controlling neuron differentiation (HEYL, PATZ1) [91,92], mitochondrial energy and neuron function (NRF1, GMEB1) [93], as well as synaptic plasticity (CREM) [94] and increased footprinting activity of early pro-neural TFs (ASCL1, NEUROG2, ARNT2) [95,96] across all three mutations (Fig. 3D-F). For GimmeMotifs maelstrom [97] motif enrichment using the SwissRegulon motif database, we categorized consensus ATAC-seq peaks into three categories: all, promoter-associated, and enhancer-associated. This revealed increased accessibility at TF motif sites related to pluripotency (NANOG, ZIC3), cell cycle (E2F8), non-ectoderm lineage (TEAD1, TBX3, EOMES), early neuron lineage (NFIC, ISL1, INSM1, NEUROD1), and neuronal repression (REST). On the other hand, we observed decreased accessibility at TF motif sites related to late-stage neuron lineage (PAX2, PBX3) mitochondrial energy and neuronal function (NRF1), as well as axonal growth and synaptogenesis (CREB1) (Fig. 4A). To uncover the functional programs associated with chromatin accessibility change, we performed differential peak enrichment using chipenrich with GOBP and Hallmark pathway databases. Promoter DARs with increased accessibility were commonly enriched for early neuron lineage, non-ectoderm lineage dedifferentiation, and repression of RNA metabolism. Promoter DARs with decreased accessibility were commonly enriched for cell cycle processes, processes modifying the chromatin state, and proteasome-controlled processes (e.g., mRNA translation and metabolic process) (Fig. 4B). Enhancer DARs were enriched for similar processes, particularly for genesets related to neuron differentiation, development, and non-ectoderm dedifferentiation (EMT, WNT β-Catenin Signaling) (Fig. 4C). chipenrich with the ENCODE-ChEA and ReMap TF-gene target databases revealed enrichment for targets of TFs associated with pluripotency (NANOG, SOX2), early neuron lineage (NFIC, PAX6, ASCL1), neuronal repression (REST), axonal growth and synaptogenesis (CREB1), and neuronal mitochondrial function (NRF1) in FAD mutations (Fig. 4D).
Previously, we demonstrated that chromatin accessibility change precedes and drives differential gene expression in PSEN1 mutant neurons[11]; therefore, we sought to determine the correlation between differential accessibility and gene expression in the FAD neurons studied here. We found the intersect of genes with a non-zero CONfident efFECT size (confect) in both gene expression and chromatin accessibility using the topconfects [50] R package, revealing a substantial number of genes in each mutation, particularly in PSEN1A79V (Fig. 5A-C). Most intersecting genes had DARs occurring in either promoter or enhancer regions, although some genes exhibited anti-correlated gene and accessibility change. To explore this further, we calculated the z-score correlation between gene expression and accessibility confect scores for all possible gene/peak pairs across the three mutations (Fig. 5D). This approach uncovered genes related to non-ectoderm dedifferentiation (SOX9, TEAD2, YAP1), early neuron lineage (ZIC2, OLIG2, ASCL1, PAX6), neuron differentiation (IRX2, MEIS2), and axonal growth and synaptogenesis (CREB1) with high correlation in at least one mutation. Next, we performed RNA-ATAC integration with diffTF to predict the differential activity of TFs using the CIS-BP motif database (Fig. 5E). By this approach, all three FAD mutations exhibited differential activity of factors involved in lineage development: increased activity of ZIC1/3 (activator) and decreased activity of IRX2 (activator). Further, we observed activation of regulators involved in early neuron lineage (ZIC1/2/3 [98], NFIA/C/X [99], PAX6) but deactivation of those controlling late stage neuron lineage (MEIS2) and mitochondrial energy and neuron function (NRF1, GMEB1) in both PSEN1A79V and PSEN2N141I mutations. Finally, we identified genesets with a strong correlation between both promoter and enhancer accessibility change with differential gene expression using intepareto [72], followed by rank-based CERNO enrichment (Fig. 5F-G). This approach identified the targets of key TFs related to chromatin modification (PCGF2), pluripotency (NANOG), neuronal differentiation (MYT1L) (Fig. 5H) and ontological geneset processes related to lineage commitment, dedifferentiation (e.g., EMT), and neuronal differentiation commonly enriched in all three mutations with strongest correlation in PSEN1A79V (Fig. 5I). In summary, this integration of RNA-seq and ATAC-seq demonstrates how the modulation of key disease endotypes, particularly reprogramming of non-ectoderm and neuronal lineages, are orchestrated via concerted chromatin and transcriptional changes. Indeed, the differential chromatin accessibility at both promoter and enhancer regions associated with transcriptional change for key endotype marker genes (CXCL12, DLX2) and transcriptional regulators with predicted differential activity (ZIC2, NEUROG2) highlight the correlation between chromatin accessibility, gene expression, and subsequent regulator activity change (Fig. 6A-D; Supplementary Fig. 5).