Research ArticleDISEASES AND DISORDERS

Dedifferentiation and neuronal repression define familial Alzheimer’s disease

See allHide authors and affiliations

Science Advances  13 Nov 2020:
Vol. 6, no. 46, eaba5933
DOI: 10.1126/sciadv.aba5933

Abstract

Identifying the systems-level mechanisms that lead to Alzheimer’s disease, an unmet need, is an essential step toward the development of therapeutics. In this work, we report that the key disease-causative mechanisms, including dedifferentiation and repression of neuronal identity, are triggered by changes in chromatin topology. Here, we generated human induced pluripotent stem cell (hiPSC)–derived neurons from donor patients with early-onset familial Alzheimer’s disease (EOFAD) and used a multiomics approach to mechanistically characterize the modulation of disease-associated gene regulatory programs. We demonstrate that EOFAD neurons dedifferentiate to a precursor-like state with signatures of ectoderm and nonectoderm lineages. RNA-seq, ATAC-seq, and ChIP-seq analysis reveals that transcriptional alterations in the cellular state are orchestrated by changes in histone methylation and chromatin topology. Furthermore, we demonstrate that these mechanisms are observed in EOFAD-patient brains, validating our hiPSC-derived neuron models. The mechanistic endotypes of Alzheimer’s disease uncovered here offer key insights for therapeutic interventions.

INTRODUCTION

Early-onset familial Alzheimer’s disease (AD) (EOFAD) is a dominantly inherited neurodegenerative disorder elicited by more than 300 mutations in the PSEN1, PSEN2, and APP genes (1). Hallmark pathological changes and symptoms observed, namely, the accumulation of misfolded β-amyloid (Aβ) in plaques and Tau aggregates in neurofibrillary tangles associated with progressive memory loss and cognitive decline, are understood to be temporally accelerated manifestations of the more common sporadic late-onset AD (LOAD). The complete penetrance of EOFAD-causing mutations has allowed for experimental models, which have proven integral to the overall understanding of AD (2). However, the failure of pathology-targeting therapeutic development suggests that the formation of plaques and tangles may be symptomatic and not describe the etiology of the disease (3, 4). In this work, we used an integrative, multiomics approach and systems-level analysis in human induced pluripotent stem cell (hiPSC)–derived neurons to generate a mechanistic disease model for EOFAD. Using patient-specific cells from donors harboring four unique mutations in PSEN1 differentiated into neurons, we characterized the disease-related gene expression and chromatin accessibility changes by RNA sequencing (RNA-seq), assay for transposase-accessible chromatin sequencing (ATAC-seq), and histone methylation chromatin immunoprecipitation sequencing (ChIP-seq). Here, we show that the defining disease-causing mechanism of EOFAD is dedifferentiation, causing neurons to traverse the lineage-defining chromatin landscape along an alternative axis to a mixed-lineage cell state with gene signature profiles indicative of less-defined ectoderm and nonectoderm lineages via REST (RE1-Silencing Transcription factor)–mediated repression of neuronal lineage specification gene programs and the activation of nonspecific germ layer precursor gene programs concomitant with modifications in chromatin accessibility. Further, a reanalysis of existing transcriptomic data from PSEN1-mutant patient brain samples demonstrates that the mechanisms identified in our experimental system recapitulate EOFAD in the human brain. Our results comprise a disease model which describes the mechanisms culminating in dedifferentiation that contribute to neurodegeneration.

RESULTS

Nondemented control (NDC), PSEN1M146L, PSEN1H163R, PSEN1A246E, and PSEN1A431E hiPSCs were generated and differentiated into neurons as previously described (510) (Fig. 1A and fig. S1A). Subsequent RNA-seq and differential expression analysis of CD44/CD184 purified neurons using the limma TREAT method (11, 12) with a threshold of |fold change (FC)| > 1.5 and adjusted P < 0.05 identified 11,187 differentially expressed genes (DEGs) in PSEN1M146L neurons, 5443 DEGs in PSEN1H163R neurons, 3069 DEGs in PSEN1A246E neurons, and 2508 DEGs PSEN1A431E neurons relative to NDC with substantial overlap between all four mutations at the level of shared DEGs and shared DEG ontology (Fig. 1, B and C, and fig. S1, B and C). Quantification of Aβ38, Aβ40, and Aβ42 levels illustrates that PSEN1 mutant neurons exhibit elevated Aβ40 and Aβ42, decreased Aβ38/Aβ42 ratios, and differential Aβ40/Aβ42 ratios relative to NDC (Fig. 1D). To characterize the functional programs modulated and to delineate the transcriptional regulators controlling them, we first carried out preranked gene set enrichment analysis (GSEA) using fgsea (13, 14) with the Gene Ontology (GO) database, revealing consistent enrichment of gene sets related to cell cycle, pluripotency, inflammation, and nonectoderm lineage dedifferentiation (mesendoderm lineage specification captured by cardiac system–related development gene sets) with positive gene regulation, combined with consistent enrichment of gene sets related to neuron lineage specification and neuron function with negative gene regulation across all four PSEN1-mutant neurons (Fig. 1E). Notably, the GO biological process (GOBP) synaptic signaling gene set, among other key neuronal function–related gene sets, was identified as significantly enriched by the CERNO statistical test (15) and the limma fry (16) rotational test with negative gene regulation in all PSEN1 mutant neurons (Fig. 1E and fig. S2A). Enrichment using the Hallmark database with fgsea and CERNO test confirmed the significant enrichment with positive gene regulation of cell cycle, inflammatory, and nonectoderm lineage [represented by the enrichment of epithelial-mesenchymal transition (EMT) gene set] programs observed using the GOBP databases, particularly for the PSEN1A431E and PSEN1A246E mutations (Fig. 1F). This test also revealed significant enrichment with positive gene regulation of gene sets related to Notch and transforming growth factor–β (TGFβ) signaling, which can contribute to the promotion of a nonectoderm lineage. As terminally differentiated neurons are generally thought to exist in the quiescent G0 phase, yet these four PSEN1 mutant neurons exhibit significant enrichment for cell cycle–related processes, we sought to determine the neuron cell cycle state. This revealed that PSEN1M146L and PSEN1A246E neuron populations have a significantly lower percentage of cells in the G0 and S phases and a higher percentage in the G2-M phases, suggesting that they have reentered the cell cycle (fig. S2B). To validate these results, we performed quantitative polymerase chain reaction (qPCR) in PSEN1M146L, PSEN1A246E, and PSEN1H163R hiPSC-derived neurons for key cell cycle–related and REST-repressed (GAD1) genes, recapitulating expression patterns observed in the corresponding RNA-seq data for each mutation (fig. S2C). To characterize the transcriptional regulation of these altered gene programs, we performed two complementary transcription-factor (TF) activity analysis approaches: motif-centric TF activity analysis by ISMARA (Integrated System for Motif Activity Response) (17) to characterize the TFs and microRNAs (miRNAs) whose motif activity most significantly changes between the PSEN1 mutant and NDC neurons; and TF-gene target-centric activity analysis by DoRothEA (Discriminate Regulon Expression Analysis) (18), which uses the viper activity inference algorithm with a curated TF target regulon generated from consensus ChIP-seq, expression signatures, and TF motif studies to predict TF activity based on the directional regulation of target genes by TFs. The top transcriptional regulators identified by ISMARA motif enrichment analysis coalesced into the same general categories of disease-causal endotypes identified by ontology enrichment analysis, with consistent activation of transcriptional regulators related to pluripotency (including SOX2 and SOX3, GATA2 and GATA3, and KLF4) (19, 20), cell cycle (E2F family members and FOXM1) (21), inflammation [interferon regulatory factor (IRF), signal transducers and activators of transcription (STAT), and nuclear factor κB (NFκB) transcriptional regulators] (22, 23), and dedifferentiation (including TEAD factors, SNAI2, and PRRX2) (Fig. 1G). A substantial number of regulators controlling neuron lineage specification and neuron function were found to have differential activity, including activation of the neuron lineage repressor REST and deactivation of neuronal and mitochondrial metabolic regulators such as NRF1 (Nuclear Respiratory Factor 1) (24, 25). The alternative DoRothEA approach identified significant activity change for TFs associated with the same endotypes identify by fgsea ontology enrichment and ISMARA analysis, with expanded activity characterization of regulators related to cell cycle activation, inflammation, dedifferentiation (SMAD and TEAD factors, and TCF7L2), and repression of neuronal lineage (deactivation of ASCL1 and SOX11) (Fig. 1H). Preranked GSEA by fgsea using the Encyclopedia of DNA Elements (ENCODE)–Chip-X Enrichment Analysis (ChEA) and ReMap TF databases showed analogous enrichment of cell cycle, pluripotency, nonectoderm lineage, and regulators identified by ISMARA and DoRothEA (fig. S2, D to E). Given the consistent enrichment of ontology and TF gene sets related to nonectoderm and neuronal lineage, we generated a modest custom TF database for transcriptional regulators controlling these programs from previously published ChIP-seq and motif data. Using the self-contained fry enrichment test, we observed significant enrichment of TFs related to dedifferentiation [YAP1 and its co-factor paralog TAZ (WWTR1), SNAI2, and NOTCH1(NICD)] (2628) and neuronal specification (MYT1L, NEUROD1, and SOX11) (Fig. 1I) (2931). As MYT1L has been previously described to function as a repressor of nonectoderm lineages and activator of neuron lineages and the target gene enrichment direction of MYT1L target genes were found to be both up-regulated (PSEN1H163R and PSEN1A431E) and down-regulated (PSEN1M146L and PSEN1A246E), we split the MYT1L gene set into MYT1L-activated and MYT1L-repressed and performed the fry enrichment test. All PSEN1 mutants, save for the PSEN1A431E condition, demonstrated significant enrichment with up-regulated gene expression for MYT1L-repressed targets (containing nonectoderm lineage genes) and significant enrich with down-regulated gene expression for MYT1L-activated targets (containing neuron lineage-specifying genes) (Fig. 1J). Further, two miRNAs previously been shown to promote neuronal lineage specification as well as to repress cell cycle and nonectoderm lineage processes, miR-124 and miR-9 (32, 33), were identified by ISMARA as having significant activity change and are known repressors of the most significantly activated TF identified, REST (34). To characterize the transcriptional regulation of these enriched gene sets, we used the miRTarBase (35) database to perform preranked GSEA miRNA target enrichment, revealing that the targets of functionally relevant miR-124 and miR-9 are significantly enriched with targets overall up-regulated in all four PSEN1 mutant conditions (fig. S2F). To determine whether this is due to loss of miR-124 and miR-9, we performed miRNA-qPCR, revealing that miR-124 and miR-9 levels are significantly decreased in PSEN1M146L and PSEN1A246E neurons relative to NDC (fig. S2G), which suggests that the increased expression of miR-124 and miR-9 targets is due to derepression caused by loss of miRNA levels.

Fig. 1 PSEN1 mutant hiPSC-derived neurons undergo dedifferentiation through the activation and repression of key endotypes.

(A) Patient-derived NDC, PSEN1M146L, PSEN1H163R, PSEN1A246E, and PSEN1A431E hiPSCs were differentiated into neurons and purified by CD44/CD184 selection. (B) RNA-seq volcano plots of DEGs (Ensembl ID) in the four PSEN1 mutations relative to NDC as determined by the limma TREAT method with a 1.5-fold change (FC) threshold and a false discovery rate (FDR)–adjusted P < 0.05. (C) Venn diagram overlap of DEGs across the four AD-causal PSEN1 mutant hiPSC-derived neurons. (D) Aβ38, Aβ40, and Aβ42 levels across NDC and the PSEN1 mutant neurons; numbers above indicate Aβ38/Aβ42 and Aβ40/Aβ42 ratios. (E and F) Ranked enrichment analysis of PSEN1M146L, PSEN1H163R, PSEN1A246E, and PSEN1A431E gene expression signatures using the (E) GOBP or (F) Hallmark database by the fgsea multilevel enrichment test (left) or tmod CERNO enrichment test (right); genes ranked by minimum significant difference (msd; CERNO) or signed msd (fgsea); *−log10 FDR–djusted P > 60. (G and H) Common TFs across the PSEN1 mutants with predicted significant activity change by (G) ISMARA motif analysis (based on z score, TF-gene Pearson correlation, and average gene target expression change) or (H) DoRothEA TF-gene target analysis (based on normalized enrichment score), curated into six key modulated disease-associated endotypes. ISMARA average target gene expression change indicated by up (increasing relative to NDC) or down (decreasing relative to NDC) arrows. (I) Common ranked enrichment analysis of PSEN1 mutant neuron gene expression signatures using a custom-defined list of TF-gene targets using the limma fry test; genes ranked by topconfects confect ranking. (J) Distribution of DEGs in PSEN1 mutant neurons ranked by topconfects for MYT1L-repressed (top) and MYT1L-activated (bottom) genes with corresponding enrichment statistic and direction of overall expression change (arrows) by the limma fry test. TF, transcription factor.

Next, we aimed to evaluate whether these observed transcriptional changes were due to modulation of the chromatin state. To this end, we performed ATAC-seq on NDC and all four PSEN1 mutant hiPSC-derived neurons, identifying a substantial number of regions of chromatin with differential accessibility occurring within promoter and transcription start site (TSS)–associated enhancer regions in PSEN1M146L, PSEN1H163R, and PSEN1A431E neurons relative to NDC (Fig. 2, A and B, and fig. S3, A and B). To identify the TFs mechanistically associated with these differential accessible regions of chromatin, we used two approaches: HINT (Hmm-based IdeNtification of Transcription factor footprints) TF footprinting analysis (36, 37) and motif enrichment analysis with GimmeMotifs (38). HINT analysis with the HOCOMOCOv11 (39) motif database identified a loss of TF footprinting activity of neuron mitochondrial regulators (NRF1, CREB1, and ATF1) (40) and KLF factors potentially involved in axon regeneration (41) across multiple PSEN1 mutant neurons (Fig. 2C). NRF1 TF footprinting accessibility, as characterized by Tn5 transposase insertions around the NRF1 motif by HINT, revealed a correlated loss of accessibility in PSEN1M146L, PSEN1H163R, and PSEN1A431E neurons (Fig. 2D). Next, we separated differentially accessible peaks into those occurring within the gene promoter region [−1000, +500 base pairs (bp) from TSS] or within a TSS-associated enhancer region (FANTOM5) (42) and performed GimmeMotifs motif enrichment with the Homer (43) algorithm and the HOCOMOCOv11 and SwissRegulon motif databases. This revealed significant enrichment of motifs related to cell cycle (E2F5), pluripotency (OCT4), chromatin modifications and interactions (ZBTB14 and CTCF), inflammation (IRF7), and cell damage and death associated with neurodegeneration (TFEB and BACH1) (44) within regions of increased accessibility with the promoter region and enrichment of neuron lineage and sensory neuron factors (RFX2 and ELK4) (45), neuron mitochondrial regulators (NRF1, CREB1, and ATF1), and chromatin modifiers [BAF170 (SMARCC2)] within regions of decreased accessibility (Fig. 2E and fig. S3, C and D). Peak enrichment using the chipenrich logistic regression enrichment test with the ENCODE-ChEA and ReMap TF databases revealed similar ontologies associated with the TFs identified, as nonectoderm lineage TFs (RUNX2, TBX2, MYOD1, and TEAD1) exhibited target gene enrichment within the increasing promoter peaks, while neuron-associated TFs (NRF1, CREB1, PML, and MAFB) (46, 47) exhibited target gene enrichment within decreasing promoter peaks (Fig. 2F). Expanding this approach to our custom TF database revealed enrichment of motifs for Notch signaling factors, which are substrates of the gamma secretase complex (NOTCH1/NICD and RBPJ) within regions of increasing promoter chromatin accessibility and motifs for neuronal lineage factors (RFX2, SOX11, OTX2, and NeuroD1) enriched within regions of decreasing promoter chromatin accessibility (fig. S3E). These ontological programs were also identified performing chipenrich using the GOBP and Hallmark gene set databases, with significant enrichment for nonectoderm lineage and cell cycle ontologies (Fig. 2G). Motif and chipenrich peak enrichment of enhancer-associated differential regions of accessibility identified significant enrichment of neuron lineage programs and regulatory TFs (ASCL1, NEUROD1, and REST), nonectoderm lineage programs and regulatory TFs (SNAI2, TWIST1, and RUNX2), and inflammatory processes (Fig. 2, H to J). Peak enrichment with our custom TF database revealed significant enrichment of targets of regulators involved in nonectoderm lineage and Notch signaling (YAP1, NOTCH1, and RBPJ) and both the repression (SOX9) and activation of neuronal lineage (MYT1L, SOX11, NEUROD1, and FOXP2) (fig. S3F). These results illustrate the differential accessibility occurring within neuron lineage and synaptic function genes, associated with later stages of lineage specification and function and controlled by regulators such as NRF1 and CREB1, primarily having loss of accessibility within the promoter regions, whereas genes associated with early stages of lineage definition and controlled by regulators such as REST and MYT1L-controlled genes primarily have differential accessibility within enhancer regions.

Fig. 2 Regions of differential chromatin accessibility are enriched for transcriptional regulators and pathways mirroring gene expression signatures.

(A) TSS and TSS-associated enhancer heatmap coverage plots of Tn5-accessible chromatin in NDC, PSEN1M146L, PSEN1H163R, PSEN1A246E, and PSEN1A431E hiPSC-derived neurons as determined by ATAC-seq. (B) Annotation and directionality of differential ATAC-seq peaks for the four PSEN1 mutant hiPSC-derived neurons relative to NDC. (C) HINT TF footprinting analysis in all differentially accessibly ATAC-seq regions using the HOCOMOCOv11 motif database to identify TFs with a change in footprinting activity. (D) Tn5 insertion density in each PSEN1 mutant relative to NDC around NRF1 motifs as determined by HINT. (E) TF motif enrichment of promoter-associated regions with increased accessibility (top) or decreased accessibility (bottom) using GimmeMotifs and the HOCOMOCOv11 motif database. (F and G) chipenrich enrichment analysis of differentially accessible promoter-associated regions with increased accessibility (top) or decreased accessibility (bottom) using (F) ENCODE-ChEA consensus and ReMap TF-gene target databases (E, ENCODE; C, ChEA; R, ReMap) or (G) GOBP and Hallmark ontology databases (H, Hallmark). (H) TF motif enrichment of differentially accessible enhancer–associated regions using GimmeMotifs and the HOCOMOCOv11 motif database. (I and J) chipenrich enrichment analysis of differentially accessible enhancer–associated regions using (I) ENCODE-ChEA consensus and ReMap TF-gene target databases or (J) GOBP and Hallmark ontology databases. All results presented P < 0.05 (GimmeMotifs) or FDR P < 0.05 (HINT, chipenrich).

Next, we investigated regions of promoters with common differential accessibility across all four PSEN1 mutations relative to NDC, focusing on genes related to neuron lineage and neuron function, two cellular programs significantly enriched in our RNA-seq and ATAC-seq analyses. This revealed chromatin accessibility loss in the promoters of three factors involved in controlling neuron lineage fate: DUSP4, POU3F1 (OCT6), and SOX11; notably, OCT6 and SOX11 were predicted by our RNA-seq analysis to have significantly decreased activity in all four PSEN1 conditions by ISMARA and DoRothEA TF activity determination approaches, respectively (Fig. 3A and fig. S4A) (48, 49). Further, key synaptic function genes, such as STXBP1, VGF, and SNAP25, display loss of promoter accessibility with corresponding and correlated loss of gene expression relative to NDC (Fig. 3B). This accessibility loss within the promoter region was similarly observed for additional genes involved in neuron lineage specification (BRINP1 and CUX2) (50, 51) and synaptic function (GABRG1, CNTN4, and BDNF) (fig. S4, A and B). We sought to expand this approach further, finding the intersection of genes with differential accessibility and directionally correlated gene expression change. Here, we found substantial overlap of DEGs with differential accessibility under the PSEN1M146L, PSEN1H163R, and PSEN1A431E conditions (Fig. 3C). Hypergeometric enrichment analysis of the DEGs with increased expression and accessibility or decreased expression or accessibility using the ENCODE-ChEA and ReMap TF databases demonstrates the mechanistic activation of genes regulated by TFs controlling dedifferentiation (TEAD1, TWIST1, SMAD3, SALL4, and TCF7L2), nonectoderm lineage (RUNX2, TBXT, and EOMES), and pluripotency (SOX2) processes and the repression of genes regulated by TFs controlling neuron lineage and function (REST), neuron energy production (NRF1, CREB1, and NFYB), and lineage-dependent chromatin modification (PCGF2 and SUZ12) (Fig. 3, D and E). The enrichment analysis using GOBP and Hallmark mirrors these results, revealing endotype-associated activation of dedifferentiation, TGFβ signaling, and pluripotency programs concomitant with repression of neuron lineage and synaptic function (Fig. 3F).

Fig. 3 Modulation of chromatin accessibility drives differential gene expression and dedifferentiation in PSEN1 mutant hiPSC-derived neurons.

(A and B) Left: Differential ATAC-seq peaks common across PSEN1 mutant hiPSC-derived neurons occurring in the promoter regions of (A) neuronal lineage regulators or (B) synaptic function genes. Right: Corresponding gene expression of each gene for each PSEN1 condition (PSEN1A246E, VGF: log2FC = −1.47, FDR P = 5.28 × 10−06). (C) Venn diagram of overlap of differential ATAC-seq peaks with corresponding directional differential gene expression change in PSEN1M146L, PSEN1H163R, PSEN1A246E, and PSEN1A431E hiPSC-derived neurons relative to NDC (dots and error bars, replicate expression values; *DEG). (D and E) Commonly enriched terms by hypergeometric (hg) enrichment analysis using the ENCODE-ChEA consensus and ReMap TF-gene target databases for genes with correlated (D) increased or (E) decreased gene expression and chromatin accessibility. (F) Commonly enriched terms by hypergeometric enrichment analysis using the GOBP and Hallmark databases for genes with correlated increased (top) or decreased (bottom) gene expression and chromatin accessibility.

As chromatin accessibility can be modulated through changes in histone methylation (52), and TF enrichment analysis revealed gene expression loss arising from differential chromatin accessibility in NRF1- and CREB1-controlled promoter regions and REST- and histone-modulating PRC2 (polycomb repressor complex 2)–controlled target gene enhancer regions, we performed ChIP-seq on NDC and PSEN1M146L and PSEN1A246E hiPSC-derived neurons for the activating mark histone H3 lysine K4 trimethylation (H3K4Me3) (primarily associated with promoters) and the repressive mark histone H3 lysine K27 trimethylation (H3K27Me3) (primarily associated with polycomb repressor–controlled regions) (53). This revealed a substantial number of differential H3K4Me3 and H3K27Me3 peaks in both PSEN1 conditions (Fig. 4, A and B). To investigate the potential transcriptional regulators associated with these methylation changes, we performed motif enrichment analysis of regions with an increase in activating H3KMe3 status (gain of H3KMe3 or loss of H3K27Me3) or an increase in repressive H3KMe3 (loss of H3KMe3 or gain of H3K27Me3) occurring within promoter and TSS-associated enhancer regions in PSEN1M146L and PSEN1A246E hiPSC-derived neurons. This identifying motifs for TFs related to cell cycle (E2F5), inflammation (IRF3), neuron lineage (NEUROD2 and PAX6), chromatin modification (EZH2 and CUX2), and neuron autophagy repression (ATF2) with an increasing in activating H3KMe3 marks and motifs for TFs related to neuron lineage (REST, RFX3, MAFB, and ELK1), neuron mitochondrial function (NRF1 and CREB1), and chromatin modification (BAF170) (Fig. 4C). These results were largely confirmed by chipenrich peak analysis using TF target and ontology databases, with observed increase in activating H3KMe3 status enriched for gene sets related to dedifferentiation (TWIST1, SMAD3, and SNAI2), nonectoderm lineage (RUNX2 and TBXT), pluripotency (NANOG and OCT4), and inflammation and an increase in repressive H3KMe3 status enriched for gene sets related to neuron lineage definition (REST, ASCL1, and ELF1) and neuron function (NRF1, CREB1, and FOXP2) (Fig. 4, D and E). Of these differentially methylated regions, 3429 and 1022 DEGs were found to have directionally correlated change in H3K4Me3 or H3K27Me3 status in PSEN1M146L and PSEN1A246E neurons, respectively (Fig. 4F). The gain of H3K4Me3 and loss of H3K27Me3 coverage or loss of H3K4Me3 and gain of H3K27Me3 coverage is apparent around genes with differential histone methylation and increased or decreased gene expression, respectively (Fig. 4G). Hypergeometric enrichment of these DEGs with differential methylation shares a similar trend with the enrichment analysis of differential methylation alone, save for the significant enrichment of cell cycle and FOXM1 target gene sets for genes with increased gene expression in the PSEN1M146L condition (Fig. 4, H and I). These results highlight that up-regulated gene expression changes related to dedifferentiation, cell cycle, and inflammation as well as down-regulated gene expression changes related to neuron lineage and function are associated with and potentially driven by change in H3K4Me3 and H3K27Me3 at promoter and enhancer regions.

Fig. 4 Dedifferentiation in PSEN1M146L and PSEN1A246E hiPSC-derived neurons is caused by changes in histone methylation.

(A) TSS heatmap coverage plots of H3K4Me3 and H3K27Me3 in NDC and PSEN1M146L and PSEN1A246E hiPSC-derived neurons as determined by ChIP-seq. (B) Annotation and directionality of differential ChIP-seq peaks in PSEN1M146L and PSEN1A246E hiPSC-derived neurons relative to NDC. (C) TF motif enrichment of regions with an increase in activating histone methylation status (increase in H3K4Me3 or decrease in H3K27Me3; top) or an increased in repressive histone methylation status (decrease in H3K4Me3 or increase in H3K27Me3; bottom) occurring within promoter or TSS-associated enhancer regions using GimmeMotifs and the HOCOMOCOv11 motif database. (D and E) chipenrich enrichment analysis of regions with an increase in activating histone methylation status (increase in H3K4Me3 or decrease in H3K27Me3; top) or an increased in repressive histone methylation status (decrease in H3K4Me3 or increase in H3K27Me3; bottom) occurring within promoter or TSS-associated enhancer regions using (D) ENCODE-ChEA consensus and ReMap TF-gene target databases or (E) GOBP and Hallmark ontology databases. (F) Venn-Euler diagram of overlap of DEGs with regions of differential H3K4Me3 or H3K27Me3 occurring within promoter or TSS-associated enhancer regions in PSEN1M146L (top) and PSEN1A246E (bottom) hiPSC-derived neurons. (G) Histone methylation density in PSEN1M146L around DEGs with corresponding gain or loss of H3K4Me3 or H3K27Me3 status within promoter and enhancer regions. (H and I) Commonly enriched terms by hypergeometric (hg) enrichment analysis using the (H) ENCODE-ChEA Consensus, ReMap, Transfac/JASPAR (T/J), and TRRUST (T) TF-gene target databases or (I) GOBP and Hallmark databases for genes with directionally correlated increased gene expression and gain of activating H3KMe3 status (top) or decreased gene expression and gain of repressive H3KMe3 status (bottom).

Next, we integrated the ATAC-seq data with the histone methylation ChIP-seq data, revealing that the intersection of DEGs with differential methylation concomitant with differential chromatin accessibility (Fig. 5A) showed a substantial overlap in the PSEN1M146L condition. Moreover, the genes with increased H3K4Me3 or decreased H3K27Me3 and increased chromatin accessibility and expression are enriched for gene sets related to dedifferentiation (SMAD3 and RUNX2), nonectoderm lineage (TBX2 and MYOD1), and inflammatory gene sets in addition to pluripotency lineage factors such as GATA2, OCT4, and PRDM14, while the genes with decreased H3K4Me3 or increased H3K27Me3 and decreased chromatin accessibility and expression are enriched for gene sets related to neuron lineage and synaptic function and transcriptional control by REST, PCGF2, and PRC2 complex factors (EZH2) (Fig. 4, B and C). The mechanistic switch from differentiated to dedifferentiated neuronal state via remodeling of the chromatin state is demonstrated by the activation of transcriptional regulators controlling dedifferentiation (TEAD2, SMAD3, and TCF7L2) (Fig. 5D) (54) and the repression of transcriptional regulators controlling neuron lineage (MYT1L and SOX11) and function (SCRT1) (Fig. 5E) (55). Furthermore, the locus around the two key neuron lineage miRNAs down-regulated in PSEN1M146L neurons, miR-9 and miR-124, show decreased chromatin accessibility and H3K4Me3 concomitant with increased H3K27Me3 (Fig. 5F). A heatmap corresponding to directionally correlated RNA expression, promoter accessibility, and H3KMe3 for an expanded list of genes belonging to the core gene programs identified here or previously identified as AD marker genes illustrates key gene expression changes occurring due to modifications of the chromatin landscape (fig. S4, C and D). The gene for the top TF with targets enriched for an increase in gene expression, chromatin accessibility, and permissive H3KMe3, the pluripotency and epigenetic regulator PRDM14, displays an increase in expression with a corresponding gain of chromatin accessibility and H3K4Me3 and loss of H3K27Me3. As REST is transcriptionally up-regulated and predicted to be activated as evidenced by ISMARA analysis and target gene repression, we carried out REST ChIP-seq in PSEN1M146L neurons. Here, we identified 1379 REST-bound target genes, a substantial portion of which are shared with other previously reported REST target lists (fig. S5, A and B) (56, 57). Furthermore, ontology enrichment analysis of these REST genes revealed significant enrichment for the canonical REST processes neuron differentiation and synapse development (fig. S5C). Enrichment analysis of these REST targets differentially expressed in PSEN1 mutant neurons identified similar neuron lineage and function gene sets as well as dedifferentiation endotype–associated processes enriched, such as EMT and blood vessel morphogenesis (fig. S5, D and E). Overall, this multiomics, integrative analysis illustrates that REST, PRC2, MYT1L, and SOX11 neuronal target genes are primarily transcriptionally repressed because of a change in chromatin accessibility driven by modification of H3K27Me3 and H3K4Me3 status, whereas NRF1 and CREB1 neuronal target genes are significantly repressed because of loss of chromatin accessibility in the promoter without change in the H3K4Me3 or H3K27Me3 methylation state.

Fig. 5 Repression of neuronal lineage and function and activation of nonectoderm lineage in PSEN1M146L hiPSC-derived neurons is orchestrated by modulation of chromatin accessibility via histone methylation.

(A) Venn diagram of overlap of DEGs, which have intersecting and directionally corresponding change in H3K4Me3 or H3K27Me3 status flanking regions of differential accessibility occurring within gene promoters or TSS-associated enhancers in PSEN1M146L hiPSC-derived neurons. (B and C) Hypergeometric enrichment analysis using the (B) ENCODE-ChEA consensus and ReMap TF-gene target databases or (C) GOBP and Hallmark databases for genes with increased gene expression, chromatin accessibility, and gain of activating H3KMe3 status (top; yellow) or decreased gene expression, chromatin accessibility, and gain of repressive H3KMe3 status (bottom; blue). (D to F) ATAC-seq, H3K4Me3, and H3K27Me3 read profiles around the promoter regions of (D) nonectoderm lineage–controlling TFs with increased gene expression, (E) neuronal lineage– and function-controlling TFs with decreased gene expression, or (F) lineage-defining miRNAs with decreased expression.

Patient-specific hiPSC-derived neurons are a powerful model system to investigate EOFAD-causing mutations in a mechanistic, systems-level fashion; to validate whether the mechanisms identified in PSEN1 mutant neurons faithfully capture the disease-related processes occurring in patients with AD, we reanalyzed a previous transcriptomic study (58) (GSE39420) of posterior cingulate cortex brain tissue from NDC and PSEN1 mutant (PSEN1M139T, n = 4; PSEN1V89L, n = 2; PSEN1E120G, n = 1) EOFAD donors (fig. S6A). Expression analysis revealed a substantial number of DEGs in PSEN1 mutant brain overlapping with PSEN1 mutant neurons, particularly in the case of the down-regulated genes (Fig. 6A and fig. S6B). This shared gene expression signature in PSEN1 mutant patient brain was similarly observed by transcriptional regulatory analysis approaches ISMARA (Fig. 6B) and DoRothEA (Fig. 6C), with significant activity change predicted for regulators related to neuronal differentiation (REST, PRC2 complex, and miR-124), nonectoderm lineage specification and dedifferentiation (TEAD1/3, SMAD1/3, PRRX2, and miR-124), pluripotency (GATA2), inflammation (IRF1/8 and STAT1/2/3), and cell cycle (E2F1/7, FOXM1, and miR-9) as observed in PSEN1 mutant neurons, with REST identified as the top enriched/activated TF in PSEN1 brain by both approaches (Fig. 6, B and C). Gene expression signature ranking of DEGs by topconfects for MYT1L-controlled genes in PSEN1 human brain resemble the trend of activation of MYT1L-repressed genes and activation of MYT1L-activated genes as observed in PSEN1 mutant neurons, suggesting that loss of MYT1L, an REST-repressed, and miR-124–promoted gene, contributes to loss of the neuronal state (Fig. 6D). Furthermore, preranked GSEA using the GOBP and Hallmark gene set databases revealed markedly similar enrichment of the six key disease-associated endotypes observed in PSEN1 neurons, as demonstrated by the up-regulation of cell cycle, inflammation, nonectoderm lineage, and pluripotency programs and down-regulation of neuron lineage and neuron function programs (Fig. 6E).

Fig. 6 Mechanisms of neuronal dedifferentiation occurring in PSEN1 hiPSC-derived neurons are observed in EOFAD patient brains.

(A) nVenn quasi-proportional diagram of overlapping DEGs across the four AD-causal PSEN1 mutant hiPSC-derived neurons and EOFAD patient brains. (B to C) TFs with predicted significant activity change in PSEN1 mutant human brains by (B) ISMARA motif analysis or (C) DoRothEA TF-gene target analysis commonly enriched with PSEN1 mutant hiPSC-derived neurons. (D) Distribution of DEGs in PSEN1 mutant human brain samples ranked by topconfects for MYT1L-repressed (top) and MYT1L-activated (bottom) genes. (E) Enrichment analysis of PSEN1 mutation human brain gene expression signatures using the GOBP and Hallmark databases by fgsea multilevel enrichment. (F) Mechanistic integration of the six disease-associated endotypes and key transcriptional regulators leading to remodeling of the chromatin state and dedifferentiation in the EOFAD disease state. (G) A lineage specification landscape model. Repression of neuron specification and neuron lineage miRNAs concomitant with activation of cell cycle, pluripotency, inflammation, and dedifferentiation programs revert EOFAD neurons along an alternative axis into a less differentiated, precursor state with nonectoderm lineage characteristics.

DISCUSSION

In this work, we used a patient-specific hiPSC-derived neuron model system and an integrative multiomics approach to characterize the mechanistic changes occurring at the TF and chromatin dynamics level leading to transcriptional dysregulation in EOFAD due to PSEN1 mutations. Expanded enrichment analysis of PSEN1 mutant neuron transcriptomics, combined with chromatin dynamics analysis, revealed six modulated gene programs integral to the understanding of the overall disease mechanism: pluripotency, dedifferentiation, cell cycle reentry, inflammation, lineage miRNA, and neuronal specification encompassing lineage definition and synaptic function. Strikingly, all four distinct PSEN1 mutations demonstrated consistent dysregulation of the endotypes described herein by both RNA-seq and ATAC-seq analyses, reinforcing the consistency and pervasiveness of these mechanisms to remodel the chromatin and gene expression landscape in PSEN1-induced EOFAD. The integration of differential RNA-seq and ATAC-seq or histone methylation ChIP-seq genes offers mechanistic insight into the type of regulatory control of these programs; for example, while the cell cycle endotype is significantly enriched in all four PSEN1 mutants and particularly the PSEN1A431E condition, there is less observed directional overlap between the genes belonging to cell cycle processes at the intersection of RNA-seq and ATAC-seq. This may indicate more transient gene regulatory mechanisms of cell cycle genes alternative to chromatin modification and explains why there is less observed overlap between RNA-seq and ATAC-seq in the PSEN1A431E condition, as this mutation has the most severe dysregulation of the cell cycle endotype among the PSEN1 mutations studied. In contrast, the neuron lineage/synaptic function and dedifferentiation endotypes are particularly enriched at the intersection of RNA-seq and ATAC-seq approaches, suggesting that these gene expression changes are mechanistically driven by persistent changes in the chromatin landscape.

Further, the comparative enrichment for top transcriptional regulators and gene sets belonging to these endotypes across the EOFAD mutations may give insight into the documented severity of each; for example, there is evidence that the PSEN1M146L and PSEN1A431E mutants causes a particularly early onset (an average of 39 and 40 years, respectively) (59, 60) along with an accelerated progression of the disease, while PSEN1H163R and PSEN1A246E mutants may have a relatively delayed onset (an average of 45 and 53 years, respectively) (61). Our analysis revealed that PSEN1M146L mutants have the most severe repression of neuronal lineage and synaptic function, whereas PSEN1H163R and PSEN1A246E mutants have the most significant up-regulation of EMT-like dedifferentiation and nonectoderm lineage, and PSEN1A431E is particularly enriched for cell cycle reactivation. This systems-level analysis approach categorized into disease-associated endotypes highlights the utility of using multiple, complementary enrichment approaches for the analysis of omics data: Whereas ISMARA identified the PSEN1M146L condition as having the most significant reactivation of the neural repressor REST, DoRothEA identified the PSEN1A431E as the condition with most substantially increasing REST activation, due to the differences in TF activity determination but also the definition of REST targets via motif (ISMARA) or ChIP-seq validated by motif (DoRothEA) data. Perhaps, these key endotypes arising from our analysis of PSEN1 mutant hiPSC-derived neurons are analogously modulated and associated with neurodegeneration in EOFAD human patient brains from alternative mutations in PSEN1 (PSEN1M139T, PSEN1V89L, and PSEN1E120G), demonstrating that these mechanisms are not unique to our hiPSC-derived neuron model system and further validating their pan-PSEN1 mutant effect.

The integration of ATAC-seq and histone methylation ChIP-seq combined with TF target and binding site enrichment analysis of RNA-seq data in PSEN1M146L neurons revealed that transcriptional modulation of the core gene programs have distinct mechanisms underlying the chromatin accessibility changes; whereas up-regulation of nonectoderm lineage and dedifferentiation genes and REST-mediated repression of neuronal genes occur due to changes in H3K4Me3 and H3K27Me3 status leading to changes in accessibility, the down-regulation of NRF1 and CREB1-mediated neuronal maturation and mitochondrial function genes primarily occurs due to chromatin accessibility changes within promoter regions potentially driven by alternative histone modifications combined with the loss of the pioneering TF capabilities of NRF1 and CREB1 (62). Analysis of promoters with differential H3K4Me3 and H3K27Me3 status revealed significant enrichment of NRF1 and CREB1, suggesting that there are additional NRF1 and CREB1 genes with promoter regions primed for chromatin accessibility loss and subsequent gene expression loss.

However, these core gene programs and the transcriptional regulators orchestrating their modulation do not exist in isolation from one another and, in certain key aspects, are intertwined in how they reshape the chromatin landscape (Fig. 4, G and H). For example, commitment to a neuron lineage involves the coordinated repression of REST, repression of the REST-recruited, histone-modifying complex PRC2, and switching of the chromatin-modulating SWI/SNF complex components from an npBAF (neuronal precursor) to an nBAF (neuronal) state, which predominantly occurs through increased expression of the CREB1-activated miRNAs miR-124 and miR-9 (32, 34, 63). miR-124 and MYT1L have been shown to restrict commitment to an neuroectoderm lineage through their repression of nonectoderm lineage genes, while miR-9 plays an integral role in the cell cycle exit required for neuronal differentiation through its repression of key cell cycle regulators upon activation (33). The mechanistic link between these primary transcriptional regulators and more lineage-specific secondary transcriptional regulators identified by our RNA-seq, ATAC-seq, and histone ChIP-seq highlights the hierarchical interconnection of these lineage regulators: For example, the neural-promoting TF SOX11, found to have decreased promoter chromatin accessibility and gene expression in multiple PSEN1 mutants and enrichment with negative target gene expression by DoRothEA and fry enrichment tests with our custom TF database, is a REST-repressed target. Similarly, transcriptional regulators of dedifferentiation identified in our RNA-seq and ATAC-seq analyses, such as TEAD1, SNAI2, SOX9, and TWIST1, are repressed by miR-124 during neuronal differentiation (64). These multifaceted functions of REST, miR-124, miR-9, and MYT1L highlight their overarching identity as key neuron lineage determination factors; therefore, the combined reactivation of REST and down-regulation of miR-124, miR-9, and MYT1L cause EOFAD neurons to dedifferentiate along an alternative axis through the repression of neuronal differentiation concomitant with activation of dedifferentiation programs (Fig. 4I). Consequently, the resulting EOFAD neurons exist in a precursor-like state with transcriptional signatures of ectoderm and nonectoderm lineage alike.

While early-onset AD accounts for at most 5% of all cases, it undoubtedly shares the hallmark pathological signatures and symptoms with the late-onset, sporadic form of the disease. However, the similarities between these two forms of AD with respect to the mechanisms driving neurodegeneration occurring before the onset of pathology and cognitive decline are less well known. A previous analysis of more than 600 brains of sporadic AD patients and non-demented, age-controlled subjects identified EMT, a nonectoderm dedifferentiation process significantly enriched in our PSEN1 mutant neurons and the PSEN1 human brains, as the top key up-regulated gene program that facilitates the progression of healthy aging into neurodegenerative AD (65). miR-124 is the top enriched transcriptional regulator for the disease-contributing gene programs identified in this previous study (fig. S6B), and we found substantial gene overlap between the EMT genes dysregulated in patients with sporadic LOAD and PSEN1 mutant neurons presented here; all four PSEN1 mutant neurons showed significant enrichment with increased gene expression for the Alz module (65) of uniquely disease-associated genes (fig. S6C). Moreover, the National Institutes of Health Accelerating Medicines Partnership - Alzheimer’s Disease (AMP-AD) program’s Agora list of sporadic AD susceptibility genes identified by multiple omics approaches shows an overlap with EOFAD at both the gene and gene program level, with key drivers and markers of neuronal dedifferentiation, including REST, HDAC1, ELAVL4, GABRA4, SCN2A, and VGF similarly identified. This commonality of cellular programs altered in both the early-onset familial and late-onset sporadic forms offers a key insight into the underlying mechanistic basis for neurodegeneration in AD and provides a basis for innovative therapeutic intervention at an earlier state in the progression of the disease.

MATERIALS AND METHODS

hiPSC line generation

NDC, PSEN1M146L, PSEN1H163R, and PSEN1A431E fibroblasts were derived from skin biopsies in accordance with University of California (UC), San Diego Institutional Review Board (IRB) approval, whereas PSEN1A246E fibroblasts were obtained commercially (Coriell, catalog no. AG06840). The generation and characterization of the hiPSC lines were carried out as previously reported (5, 7) from fibroblasts by retroviral transduction using the reprogramming factors OCT4, SOX2, KLF4, and c-MYC. These five iPSC lines were used for all downstream neuron differentiation applications.

Human neuron preparation

The protocol used for neuron preparation was previously described (6, 9) with modifications. Briefly, PA6 cells were plated in a 10-cm dish and seeded with 100,000 iPSC next day. To enhance neural induction, cultures were treated with Noggin (500 ng/ml; R&D Systems) and 10 μM SB431542 (Tocris) for the first 6 days of differentiation. On day 12, neural stem cells (NSCs) were sorted using cell surface signature CD24+/CD184+/CD44/CD271. NSCs were expanded in NSC growth medium [DMEM/F-12, GlutaMAX (Thermo Fisher Scientific, catalog no. 10565018), 1× B-27 (Thermo Fisher Scientific, catalog no. 17504044), 1× N-2 (Thermo Fisher Scientific, catalog no. 17502001), 1× penicillin-streptomycin (Thermo Fisher Scientific, catalog no. 15070063), and human bFGF-2 (20 ng/ml; BioPioneer, catalog no. HRP-0011)]. At 80% confluence, the medium was changed to neuron differentiation medium (DMEM:F12 + GlutaMAX, 1× B-27, 1× N-2, and 1× penicillin-streptomycin) for 3 weeks of differentiation as previously described (8, 9). After differentiation, the cultures were dissociated with Accutase (Sigma-Aldrich, catalog no. A6964). Cells were resuspended in 200 μl of iMag buffer (1× neural differentiation medium, 0.5 μM EDTA, and 0.5% bovine serum albumin), followed by incubation with R-phycoerythrin (PE) mouse anti-human CD184 and CD44 antibodies (BD Biosciences, catalog nos. 561733 and 561858, respectively) for 15 min on ice in the dark. The mixture was washed with iMag buffer and subsequently incubated with anti-PE–conjugated magnetic beads (BD Biosciences) for 30 min at room temperature as previously described (9, 10). Magnetic bead separation was carried out for 8 min according to the manufacturer’s protocol (BD Biosciences). The supernatant containing purified CD184/CD44 neurons was then removed and spun down for downstream applications.

RNA-seq and data processing

For RNA-seq, total RNA from magnetically purified human NDC (replicates, n = 8), PSEN1M146L (replicates, n = 4), PSEN1H163R (replicates, n = 3), PSEN1A246E (replicates, n = 3), and PSEN1A431E (replicates, n = 4) hiPSC-derived neurons was prepared using the RNeasy Plus Micro Kit (QIAGEN, catalog no. 74034) according to the manufacturer’s protocol. On-column deoxyribonuclease digestion was subsequently performed on total RNA extracts according to the manufacturer’s recommendations to remove any genomic contamination (QIAGEN, catalog no. 79254). Libraries were prepared for RNA-seq using the TruSeq Stranded Total RNA Library Prep Kit (Illumina, catalog no. RS-122-2303) by the Ribo-Zero ribosomal RNA reduction method (Illumina, catalog no. MRZG12324). Samples were sequenced at the UC San Diego Institute for Genomics Medicine (IGM) sequencing core on an Illumina HiSeq 4000 generating paired-end, 100-bp reads with an average of 25 million reads per sample (Illumina, catalog no. FC-410-1001). We incorporated an additional hiPSC-derived neuron NDC line [Craig Venter B (CVB) line 1; replicates, n = 3] from a previously published study (66) using the same neuron differentiation method and processed the data in parallel with the data generated here.

RNA-seq data preprocessing was performed using the TrimGalore! package (67), removing sequencing adaptors and selecting for all paired-end reads above a quality score threshold (Phred Q > 20). Trimmed RNA-seq reads were mapped to the GRCh38.99 human transcriptome using kallisto v0.46.1 (68) with the options -bias, --rf-stranded, and -b 100, followed by transcript level summation to the gene level using the R package tximport v1.8.0 (69). Differential expression analysis was performed for all four PSEN1 mutant iPSC-derived neurons relative to NDC using the TREAT method in the R package limma (11, 12) with a null-hypothesis t test FC threshold of 1.5 and selection of DEGs with a false discovery rate (FDR)–adjusted P value of <0.05.

TF and miRNA activity analysis

TF activity analysis was performed using ISMARA (17) and DoRothEA (18). For ISMARA analysis, quality and adaptor trimmed fastq.gz RNA-seq files for all hiPSC-derived neurons were uploaded to ismara.unibas.ch for processing, followed by sample average. To determine a directional z-score for each enriched motif identified, the z-score for each given motif was multiplied by the sign of the Pearson correlation between each motif and its target genes and the direction of change in expression for said target genes (i.e., −1 for down-regulated genes and +1 for up-regulated genes). For motifs associated with miRNAs, qPCR expression (miR-9 and miR-124) or literature evidence (e.g., miR-218) was used to determine a positive of negative correlation with gene targets of the given motif. For DoRothEA analysis, the gene ranking for each PSEN1 mutation relative to NDC was generated by calculating absolute value of the minimum significant distance (msd), defined as the lower bound of the 95% confidence interval of the log2FC, using the tmodlimmaTest function in the tmod (15) R package and multiplying by the sign of the log2FC. TF activity analysis was performed using the msviper function in the viper (70) R package with the DoRothEA C regulon network. TF target analysis was also performed with additional TF gene set databases: Transfac/JASPAR and ENCODE-ChEA consensus from Enrichr (71, 72), incorporating the consensus genes targets identified for 104 TFs from the ENCODE project (73) and ChEA datasets, ReMap (74), and TRRUST (75); for miRNA target analysis, the miRTarBase v2017 (35) database was used. For TF target analysis of neuron-related TFs not well represented in the ISMARA database, DoRothEA regulons, or other TF databases, we generated a custom TF target database of endotype-associated TFs obtained from Gene Expression Omnibus (GEO) (data S1). ChIP-seq data for each TF were downloaded from GEO provided by the original authors, or raw sequencing data were uploaded to CRUNCH (76) for peak identification. Peaks were annotated with the Homer v4.10 function annotatePeaks, and each list of target genes filtered for HGNC (HUGO Gene Nomenclature Committee) protein coding genes. For the MYT1L target gene list, candidate MYT1L targets were identified using the MYT1L_200-623 primary motif from Mall et al. (29) and the Homer v4.10 (43) function scanMotifGenomeWide, subsequently annotated with annotatePeaks. This list of candidate genes was first filtered for HGNC protein coding genes and then further filtered for differential gene expression change upon MYT1L knockdown or overexpression previously characterized (29, 77), splitting into MYT1L-repressed and MYT1L-activated gene sets. For enrichment analysis using this custom neural TF database, gene expression signatures were ranked by confect value at a FDR generated by the topconfects (78) R package and directional enrichment performed using the fry ranked rotational test function in limma with confect value serving as gene weights. The self-contained fry test was selected over the competitive fgsea test used for enrichment analysis with the publicly available TF database above due to the smaller size of our custom neural TF database.

Gene set enrichment and pathway analysis

Gene set enrichment and pathway analysis for RNA-seq was performed primarily by two complementary approaches: (i) GOBP (79) and Hallmark Pathway (80) enrichment of differential gene expression signatures using the fgseamultilevel function in the fgsea (14) R package and (ii) the tmodCERNOtest function in the tmod (15) R package. For the weighted, directional fgsea statistical enrichment test, genes were ranked by the directional msd statistic as described for DoRothEA analysis, while for the CERNO test, genes were ranked by msd. For miRNA target analysis, the miRTarBase v2017 database was used with both fgsea and a hypergeometric test of up-regulated and down-regulated DEGs using the tmodHG function in tmod. For enrichment using fry rotational test and GOBP gene sets, genes were ranked by confect value with topconfects as in enrichment analysis using our custom TF database.

ATAC-seq and data processing

ATAC-seq transposition experiments were performed as previously reported (81) on 50,000 NDC (replicates, n = 5), PSEN1M146L (replicates, n = 2), PSEN1H163R (replicates, n = 3), PSEN1A246E (replicates, n = 2), and PSEN1A431E (replicates, n = 3) hiPSC-derived neurons, using the Illumina Nextera DNA Sample Preparation Kit (Illumina, catalog no. 15028523) and the QIAGEN MinElute PCR Purification Kit (QIAGEN, catalog no. 28004). ATAC-seq libraries were generated from transposed DNA using the Kapa Biosystems Real-Time Library Amplification Kit (Kapa Biosystems, catalog no. 07959028001) according to the manufacturer’s recommendations, monitoring amplification by qPCR and stopping the reaction when all samples reached a fluorescence amplification intensity between standards 1 and 3. ATAC-seq libraries were further purified using the QIAGEN MinElute PCR Purification Kit and sequenced at the UC San Diego IGM sequencing core on an Illumina HiSeq 4000 platform generating paired-end, 50-bp reads with an average of 25 million reads per sample.

ATAC-seq data preprocessing was performed with TrimGalore!, removing sequencing adaptors and selecting for all paired-end reads above a quality score threshold (Phred Q > 20). Trimmed reads were aligned to the GRCh38 human genome (GCA_000001405.15 with no alternative analysis) using BBMap v37.95 in the BBTools suite with the options maxindel = 20 ambig = random, followed by sorting and indexing of bam files using SAMtools v1.3 (82), and annotation of PCR duplicates using the Picard v2.3.0 MarkDuplicates function with the option VALIDATION_STRINGENCY = LENIENT. All duplicates and mitochondrial, chromosome X, chromosome Y, and Epstein-Barr virus (EBV) reads were removed using SAMtools v1.3 command view with the options -b -h -f 3 -F 4 -F 8 -F 256 -F 1024 -F 2048. One PSEN1A246E ATAC-seq sample was found to be unsuitable for further analysis, so was removed, and the other PSEN1A246E ATAC-seq sample was randomly subsampled into two using the SAMtools view function. To determine regions of open chromatin, i.e., those accessible by the Tn5 transposase, HMMRATAC v1.2.5 (83) was used to call peaks on the ATAC-seq data and determine regions of open chromatin with the options –m 50,200,400,600 --score all. The specifically called open regions of chromatin were then passed to the R package Diffbind v2.8.0 (84) to determine regions of differential accessibility between each NDC and PSEN1 mutant conditions using the CONDITION_CONSENSUS peak selection and edgeR-block differential peak method. Differentially accessible regions of DNA were annotated using the R packages ChIPseeker (85), defining the promoter region −1000 to 500 bp from the TSS. Enhancer-associated ATAC-seq regions were defined differential peaks occurring within the region list of TSS-associated enhancers generated by the FANTOM5 project (42), finding the intersecting DNA regions with the join_overlap_inner function in the plyranges R package (86).

ATAC-seq enrichment analysis

To determine the TF footprints associated with the gain or loss of chromatin accessibility and motifs enriched for these regions, we used two approaches: HINT-ATAC for TF footprinting and GimmeMotifs for motif enrichment. HINT v0.12.3 (36, 37) was used with the following parameters: rgt-hint function footprinting with the options --atac-seq --paired-end --organism = hg38 to identify TF footprints in each sample; rgt-motifanalysis function matching to match footprints to known TFs in the HOCOMOCO v11 (39) human pwm database; and rgt-hint function differential with the option --window-size 200. Differential HINT footprinting was performed for all differentially accessible regions over the entire HMMRATAC-called open chromatin region. For TF motif enrichment with GimmeMotifs (38), differential ATAC-seq regions were split into peaks occurring within promoters or within enhancers. The gimme motifs function was used with both the HOCOMOCOv11 (39) and SwissRegulon (87) pwm databases and the Homer (43) motif-finding algorithm. For enrichment analysis of differential ATAC-seq regions, we used the logistic regression model test function chipenrich in the chipenrich R package (88) with the locus definition nearest_tss for enrichment of promoter-located peaks and 1kb_outside for enhancer-located peaks using the ENCODE-ChEA and ReMAP TF-gene target databases and GOBP and Hallmark ontology databases. Genes with differential chromatin expression were defined as those called as DEGs (|FC| > 1.5 and FDR P < 0.05) as described above and had a differentially accessible peak (passing the threshold of |FC| > 1.5 and FDR P < 0.05) occurring in a promoter or enhancer region in the corresponding direction (e.g., decreased accessibility with decreased gene expression). For TF-gene target and ontology enrichment, we performed a hypergeometric test using the tmodHG function in the tmod (15) R package with the ENCODE-ChEA and ReMap TF databases and the GOBP and Hallmark databases. Promoter and TSS-associated enhancer heatmaps were generated with the deepTools (89) functions computeMatrix and plotHeatmap. ATAC-seq coverage plots were generated in IGV v2.8.6 (90).

Histograms of ATAC-seq read density were generated with the Homer v4.10 function annotatePeaks. For ranking of differential accessibility within the promoter or enhancer regions of target genes of different TFs, the confect ranking [i.e., the log2FC value at which each ATAC-seq region achieved significant (FDR P < 0.05) differential accessibility] as determined by the function edger_confects in the topconfects R package was used, and distribution figures were made with the barcodeplot function in the limma R package.

Integration of RNA-seq and ATAC-seq was performed by finding all overlapping Ensembl genes with differential gene expression and differential chromatin accessibility occurring in the same direction with the threshold of |FC| > 1.5 and a FDR P < 0.05 for both RNA-seq and ATAC-seq. For enrichment analysis, the hypergeometric test in the tmod R package was used with the ENCODE-ChEA, ReMap, and our custom TF databases, as well as the GOBP and Hallmark gene set databases with a defined background list of all genes with a called ATAC-seq peak in each condition.

ChIP-seq and data processing

For H3K4Me3 and H3K27Me3 ChIP-seq experiments, NDC, PSEN1M146L, and PSEN1A246E (n = 1) hiPSC-derived neurons were processed using the Active Motif ChIP-IT High Sensitivity Kit according to the manufacturer’s recommendations (Active Motif, catalog no. 53040). Briefly, cells were washed with cold 1× phosphate-buffered saline (PBS) and collected in a 1.5-ml tube in PBS, followed by formaldehyde fixation for 10 min at room temperature. After addition of stop solution, cells were incubated at room temperature for 5 min and subsequently washed with 1× PBS wash buffer and resuspended in chromatin prep buffer. Following incubation on ice for 10 min, cells were lysed by performing 30 plunges in a chilled dounce homogenizer. Samples were then resuspended in Active Motif ChIP Buffer and chromatin sonication performed using a probe tip sonicator on ice for a total of 10 min per sample (30 × 20-s intervals). Optimal sonicated chromatin fragment lengths were confirmed by agarose gel and used for immunoprecipitation (IP) reactions. Five hundred nanograms of purified, sonicated chromatin was incubated with 5 μl of H3K4Me3 (Active Motif, catalog no. 39159) or H3K27Me3 (Active Motif, catalog no. 39155) rabbit polyclonal antibody overnight at 4°C. After decross-linking, ChIP DNA was purified using the QIAGEN MinElute PCR Purification Kit (QIAGEN, catalog no. 28004); libraries for H3K4Me3, H3K27Me3, and 1% input were generated for NDC and PSEN1M146L conditions at the UC San Diego IGM sequencing core using Illumina TruSeq LT adaptors (Illumina, catalog no. FC-121-2001). Sequencing was performed on an Illumina HiSeq 4000 generating single-read, 75-bp reads with an average of 20 million reads per sample.

ChIP-seq data preprocessing was performed using the BBDuk package (v37.95) in the BBTools suite, removing sequencing adaptors and selecting for all paired-end reads above a quality score threshold (Phred Q > 10). Trimmed ChIP-seq reads were aligned to the GRCh38 human genome (GCA_000001405.15 with no alternative analysis) using BBMap v37.95 in the BBTools suite with the options maxindel = 20 ambig = random, followed by sorting and indexing of bam files using SAMtools v1.3, annotation of PCR duplicates using Picard v2.3.0 MarkDuplicates with the option VALIDATION_STRINGENCY = LENIENT, and removal of PCR duplicates using SAMtools v1.3 (82) command view with the options -b -h -F 1024. Differential ChIP-seq peaks between the NDC and PSEN1M146L or NDC and PSEN1A246E conditions were identified with the sicer_df function in the sicer2 package (91) with the filter threshold |FC| > 1.5 and a FDR P < 0.05. Promoter and TSS-associated enhancer annotation, motif enrichment using GimmeMotifs, TF target and ontology enrichment using chipenrich, and hypergeometric enrichment analysis of integrated ChIP-seq peaks and DEGs were performed as described above for ATAC-seq with the additional TRRUST and Transfac/JASPAR gene set databases. The Homer function annotatePeaks was used to generate histograms of ChIP read density at histone methylation peak centers for genes with differential expression and differential histone methylation with the options -size 5000 -hist 10.

Integration of RNA-seq and ChIP-seq was performed by finding all overlapping Ensembl genes with differential gene expression and differential H3K4Me3 or H3K27Me3 occurring in the same direction as the activating or repressing function of each histone methylation mark (e.g., increased gene expression with gain of H3K4Me3 or loss of H3K27Me3) with the threshold of |FC| > 1.5 and a FDR P < 0.05 for both RNA-seq and ChIP-seq. For enrichment analysis, the hypergeometric test in the tmod R package was used with the ENCODE-ChEA, ReMap, and our custom TF databases, as well as the GOBP and Hallmark gene set databases with a defined background list of all genes with a called ChIP-seq peak in each condition.

Human brain transcriptomic analysis

Microarray data corresponding to NDC and PSEN1 mutation EOFAD human brain samples from a previously published study (58) were retrieved from the National Center for Biotechnology Information (NCBI) GEO database (GSE39420) as raw CEL intensity files. Thirteen samples from this dataset, all obtained from the posterior cingulate cortex region, were selected for analysis: six NDC subject brains samples and seven EOFAD brain samples from patients harboring PSEN1M139T (replicates, n = 4), PSEN1V89L (replicates, n = 2), and PSEN1E120G (replicates, n = 1) mutations. Samples were background-corrected and log2-normalized using the “RMA” function (9294) in the oligo package v1.44.0 (95) in R. DEGs were determined using the limma v3.36.2 (11) package, selecting for genes with a FDR-adjusted P (q) < 0.05. Venn diagram overlaps of EOFAD brain DEGs with PSEN1 mutation hiPSC-derived neurons were generated with the nVennR v0.2.2 package (96) in R. ISMARA, DoRothEA, and fgsea enrichment analyses were performed as described for RNA-seq data.

Quantification of Aβ levels

Aβ peptide levels were assayed as previously described (7). Briefly, a Meso Scale human (4G8) Aβ 3-plex ultrasensitive kit (Meso Scale Discovery, catalog no. K15141E-2) and a custom human total Aβ kit (Meso Scale Discovery, catalog no. N45CA-1) were used to assay Aβ peptides in media from sorted, purified neurons following 5 days from plating. The NDC, PSEN1M146L, and PSEN1A246E Aβ quantification data presented here are from Liu et al. (7), whereas PSEN1H163R and PSEN1A431E Aβ quantification data were generated for this work.

miRNA qPCR assay

Total RNA from magnetically purified NDC, PSEN1M146L, and PSEN1A246E (replicates, n = 3) human iPSC-derived neurons was prepared using the miRNeasy Micro Kit (QIAGEN, catalog no. 217084), on the basis of the manufacturer’s procedures. complementary DNA (cDNA) was produced using the TaqMan Advanced miRNA cDNA Synthesis Kit (Life Technologies, catalog no. A28007) according to the manufacturer’s instruction. Briefly, miRNA was polyadenylated at 37C for 45 min using 20 ng of total RNA in a 5-μl reaction. After heat stop of the polyadenylation reaction at 60°C for 10 min, the reaction was held on ice for 2 min. The adaptor was then ligated to the miRNAs at 16°C for 1 hour using RNA ligase in a 15-μl reaction. Next, cDNA was synthesized from the adapter-ligated miRNAs. Quantitative reverse transcription (qRT)–PCR for miRNAs was carried out using advanced TaqMan Advanced miRNA assay (Life Technologies, catalog no. A25576), with TaqMan Fast Advanced Master Mix (Life Technologies, catalog no. 4444557) and TaqMan prime probe sets for miR-124 (Hsa-miR-124-3P and 477879_mir) and miR-9 (Hsa-miR-9-5P and 478214_mir). The qRT-PCR program was as follows: 95°C for 10 min, leading to 40 cycles at 95°C for 10 min, followed by 60°C for 20 s. The log2FC for each gene was calculated by the ΔΔCt method.

Cell cycle assay

Magnetically purified NDC, PSEN1M146L, and PSEN1A246E (replicates, n = 3) human iPSC-derived neurons were plated in 96-well plates. Cells were fixed in 70% ethanol and stored at 4°C until the day of assay. Cells were stained with FxCycle Far Red Stain (Thermo Fisher Scientific, catalog no. F10348) according to the manufacturer’s protocol. Stained cells were analyzed by flow cytometry using an Accuri C6 flow cytometer (BD Biosciences).

Statistical analysis

GraphPad Prism v7.0e software was used for statistical analysis to calculate either student’s t test or one-way analysis of variance (ANOVA) with posttest analysis (Dunnett’s multiple comparison test) for cell cycle and qPCR assays.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/46/eaba5933/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank K. Jepsen and the UC San Diego IGM for assistance with the sequencing performed in this research and M. Pachkov at the University of Basel for assistance with ISMARA data processing. Funding: This study was supported by the Alzheimer’s Association New Investigator Research Award (NIRG-14-322164) (to S.H.Y.); NIH grants P50 AGO5131 (to D.R.G.), U01 NS 074501-05 (to S.L.W.), R01 LM012595 (to S.S.), U01 CA198941 (to S.S.), U01 DK097430 (to S.S.), R01 HD084633 (to S.S.), and R01 HL106579-07 (to S.S.); NSF grant STC CCF-0939370 (to S.S.); Veterans Affairs RR&D 1I01RX002259 (to S.L.W.); and Cure Alzheimer’s Fund (CAF) grants (to S.L.W.). Author contributions: A.B.C. performed the RNA-seq, ChIP-seq, and ATAC-seq experiments, carried out the data processing and systems analysis of the RNA-seq, ChIP-seq, and ATAC-seq experiments, developed the mechanistic framework, and wrote the manuscript with contributions to Materials and Methods from other authors. A.B.C. and S.S. were involved in the characterization of the mechanisms and networks associated with EOFAD. S.H.Y. was involved in the design of hiPSC differentiation into neurons and the validation experiments along with S.L.W. and S.S. Q.L., supervised by S.H.Y., carried out the hiPSC-neuron preparation, RNA isolation, qPCR validation, Aβ peptide quantitation, and cell cycle assay experiments. G.P.S. contributed to the RNA-seq, ChIP-seq, and ATAC-seq. D.R.G. was involved in part of the study design. S.S. was involved in the overall study design, analysis, and revisions of the manuscript. Competing interests: The authors declare that they have no competing interests. The contents do not represent the views of the U.S. Department of Veterans Affairs or the U.S. government. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All sequencing presented in this work is available at the NCBI GEO under the SuperSeries accession GSE123998, encompassing RNA-seq (GSE123891), ATAC-seq (GSE123997), and ChIP-seq (GSE123884) datasets. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article