Epigenomic priming of immune genes implicates oligodendroglia in multiple sclerosis susceptibility

Abstract

Multiple sclerosis (MS) is characterized by a targeted attack on oligodendroglia (OLG) and myelin by immune cells, which are thought to be the main drivers of MS susceptibility. We found that immune genes exhibit a primed chromatin state in single mouse and human OLG in a non-disease context, compatible with transitions to immune-competent states in MS. We identified BACH1 and STAT1 as transcription factors involved in immune gene regulation in oligodendrocyte precursor cells (OPCs). A subset of immune genes presents bivalency of H3K4me3/H3K27me3 in OPCs, with Polycomb inhibition leading to their increased activation upon interferon gamma (IFN-γ) treatment. Some MS susceptibility single-nucleotide polymorphisms (SNPs) overlap with these regulatory regions in mouse and human OLG. Treatment of mouse OPCs with IFN-γ leads to chromatin architecture remodeling at these loci and altered expression of interacting genes. Thus, the susceptibility for MS may involve OLG, which therefore constitutes novel targets for immunological-based therapies for MS.

Keywords: Polycomb; chromatin; genome-wide association studies; histone modifications; major histocompatibility complex; multiple sclerosis; myelin; neuroimmunology; oligodendrocyte; single-nucleotide polymorphisms.

Figure 1.. scATAC-seq reveals primed and increased CA at immune gene loci in OLG from the EAE mouse model of MS

(A) Experimental setup for scATAC-seq of OLG from the EAE mouse model of MS. (B) Uniform manifold approximation and projection (UMAP) clustering based on CA of 10x Genomics chromium scATAC-seq. (C) Origin of individual cells, projected on top of UMAP clustering. (D) Label transfer from matched scRNA-seq data (Falcão et al., 2019) projected on top of CA UMAP. (E) GSEA of the nearest genes to enriched CA peaks for immune Hallmarks categories. (F) Integrative genomics viewer (IGV)-merged normalized tracks of CA with 100 randomly selected cells for each cluster, with MiGl clusters grouped together. scATAC co-accessibility connections are shown. Highlighted with gray boxes are regions with differential CA in specific clusters or promoter priming and connections between regulatory regions. Genomic coordinates are shown. OPC, oligodendrocyte precursor cell; VLMC, vascular leptomeningeal cell; PLC, pericyte-like cell; MiGl, microglia; NFOL, newly formed oligodendrocytes; MOL, mature oligodendrocyte.

Figure 2.. Primed CA of immune genes in single OPCs and MOLs

(A) Genes in OPCs (left) or MOL1/2 (right) are clustered based on gene expression differences between EAE versus Ctr and correlation with CA activity score (CA over 500-bp promoter region). Top GO terms for Type1 genes (with increased expression and CA in EAE; [Type1a] low and [Type1b] high CA in Ctr-OPCs), Type2 (with increased expression in EAE, but no change in CA, [Type2a] high and [Type2b] low CA in Ctr-OPCs), and Type3 (with reduced expression in EAE, but no change in CA). Type4 (with reduced expression and CA in EAE) had no GO terms. (B) UMAP based on CA (right) and RNA-seq (left) of 10x Genomics multi-ome (simultaneous scATAC and RNA-seq) of Sox10-GFP cells sorted from the spinal cord of Ctr (2) and EAE mice (2, at disease peak). Label transfer from matched scRNA-seq data (Falcão et al., 2019). (C) (Left) IGV tracks of CA for each selected cluster, with MiGl clusters grouped together. scATAC co-accessibility connections are shown. Highlighted with gray boxes are regions with differential CA in specific clusters or promoter priming and connections between regulatory regions. Genomic coordinates are shown. (right) Violin plots depicting the expression of Ifit2 and Nlrc5 in each cluster. Cell acronyms as in Figure 1.

Figure 3.. Primed CA at immune gene loci in primary mouse OPCs

(A) Volcano plots showing differential gene expression in RNA-seq (left) and CA at promoter regions in ATAC-seq (right) between Ctr-OPCs and OPCs treated with 100 ng/mL IFN-γ for 48 h. Genes with adj. p value < 0.05 and log₂ fold change >1.5 are shown in red. (B) Genes in IFN-γ-treated OPCs and Ctr-OPCs are clustered based on their chromatin activity score (CA over 500-bp promoter region) and gene expression correlation. Top GO terms are shown for Type1–Type4 as defined in Figure 2A. (C) IGV tracks are shown for ATAC-seq and RNA-seq in IFN-γ-treated OPCs and Ctr-OPCs for selected genes. Highlighted with gray boxes and arrows are regions with differential CA in IFN-γ-treated OPCs or promoter priming; black arrows, promoter regions. Merged tracks of 3 biological replicates are shown for ATAC-seq and 4 biological replicates for RNA-seq.

Figure 4.. STAT1 and BACH1 have increased motif accessibility (MA) in OLG from EAE and are involved in IFN-γ-mediated regulation of immune genes in OPCs

(A) ChromVAR clustering of TF motif variability from scATAC-seq. Each row presents a TF motif, whereas each column represents a single cell. Scale, blue (low TF MA) to yellow (high). (B) TF motif variability projected on top of UMAP of scATAC-seq. (C) Volcano plots showing differential expressed genes in RNA-seq upon transfection of primary OPCs with siRNAs targeting Bach1 before treating with IFN-γ for 6 h. 4 biological replicates. Genes with adj. p value < 0.05 and log₂ fold change >1.5 are shown in red. Hmox1 coordinate values are shown (since beyond y-axis range). (D) IGV tracks for RNA-seq in IFN-γ-treated OPCs and Ctr-OPCs after transfection with siRNAs targeting Bach1 for selected genes. (E) qRT-PCR targeting selected genes upon transfection of primary OPCs with siRNAs targeting Stat1 before treating with IFN-γ for 6 h. 3–4 biological replicates. Error bars, mean ± SEM. (F) IGV tracks for RNA-seq in IFN-γ-treated OPCs and Ctr-OPCs after transfection with siRNAs targeting Stat1 for selected genes. (G) IGV tracks showing STAT1 binding in OPCs upon IFN-γ treatment and in Ctr-OPCs, assessed with CUT&Tag.3 biological replicates. (H) Heatmaps depicting STAT1 binding in OPCs treated with IFN-γ at genes that were up- or downregulated by STAT1 knockdown (as assessed by adjusted p value from bulk RNA-seq), centered at peaks.

Figure 5.. H3K27ac, CTCF binding, and enhancer-promoter contacts at immune genes in mouse OPCs are altered upon IFN-γ treatment

(A) Volcano plots showing differential H3K27ac (left) and CTCF binding (right) between IFN-γ-treated OPCs and Ctr-OPCs, assessed with Cut&Run.3 biological replicates. Genes with adj. p value < 0.05 and log₂ fold change >1.5 are shown in red. (B) Number of genes (y axis) in Ctr-OPCs (blue) and IFN-γ-treated OPCs (green) with n predicted interactions (x axis). (C) IGV tracks showing CTCF binding and H3K27ac occupancy, assessed with Cut&Run, ATAC-seq in IFN-γ-treated OPCs, and Ctr-OPCs for MHC-I and MHC-II loci. Predicted enhancer/promoter contacts computed by the activity-by-contact (ABC) model (Fulco et al., 2019) based on CA and H3K27ac-HiChIP. Highlighted with gray boxes are regions with increased H3K27ac, CTCF binding, CA, and/or predicted interactions in IFN-γ-treated OPCs. Highlighted with a red arrow is an enhancer region interacting with multiple genes in the MHC-I and MHC-II loci. Merged tracks for 3 biological replicates per condition.

Figure 6.. H3K27me3 modulation is involved in IFN-γ-mediated immune gene activation in OPCs

(A) Volcano plots for H3K27me3 and H3K4me3 in IFN-γ-treated OPCs versus Ctr-OPCs assessed with Cut&Run. Two replicates. Genes with adj. p value < 0.05 and log₂ fold change >1.5 are shown in red. (B) Cut&Run IGV tracks for selected MHC-I and MHC-II genes with increased H3K4me3 or decreased H3K27me3 in OPCs upon IFN-γ treatment. Merged tracks of two replicates. (C) Upset plots of Cut&Run Ctr-OPCs and IFN-γ-treated peak intersection in Type2 genes. Top barplot shows the number of intersecting peaks per combination, left barplot shows the size of each peak dataset, and the matrix shows the Cut&Run peak sets (dots) and shared (connecting line) in each combination. (D) Venn diagram showing the number of genes enriched in OPCs upon treatment with 1.5 μM EZH2 inhibitor EPZ011989 (EZH2i) for 4 days, with and without subsequent co-treatment with 100 ng/mL IFN-γ for last 6 h (and the genes enriched in both) and the top gene ontology biological terms for the genes in each category. (E) RNA-seq IGV tracks for MHC-I, MHC-II, and cytokine genes with increased expression upon EZH2i in IFN-γ-spiked OPCs. Merged tracks of three biological replicates.

Figure 7.. Primed CA at immune gene loci in human neural cells

(A) UMAP based on CA (right) and RNA-seq (left) of 10x Genomics multi-ome of brain gray matter from two healthy individuals. (B) IGV tracks of CA for each selected cluster (hg38). Violin plots depicting the expression of selected genes in each individual cluster. (C) Overlap with MS-associated GWAS variants. For each peak set, expected (x axis) versus observed (y axis) number of SNP hits overlapping the human healthy individuals hg19 multi-ome scATAC-seq peaks, scATAC-seq cell-type-specific peaks from Ctr and EAE mice, and Ctr-OPC and IFN-γ-treated primary OPC ATAC-seq peaks. Dot size scaled to adjusted p value and adjusted p values < 0.01 in red. ENDO, endothelial cells; ASTRO, astrocytes; EXCNEU, excitatory neurons; INHNEU, inhibitory neurons; MIGL, microglia; OLIGO or MOL, mature oligodendrocyte; OPC, oligodendrocyte precursor cell; VLMC, vascular leptomeningeal cell; PLC, pericyte-like cell.

Figure 8.. Overlap between CA peaks and MS susceptibility SNPs in OLG and activation of associated genes in mouse OPCs by IFN-γ

(A) Dot plot depicting overlap of CA peaks in specific cell types from human healthy individuals hg19 single-cell multi-ome and scATAC-seq cell types from Ctr and EAE mice, with individual MS susceptibility MHC, non-MHC, and suggestive SNPs (International Multiple Sclerosis Genetics Consortium, 2019) and outside variants (Factor et al., 2020). (B and C) SNP coordinates for six MS SNPs in the hg38 human genome reference are shown with corresponding CA regions derived from merged scATAC-seq populations from the adult brain single-cell multi-ome. (B) Corresponding locations are shown in the mouse mm10 genome reference with IGV tracks of CA in 50 randomly selected individual cells from scATAC-seq from EAE and CFA-Ctr mice. Red boxes show scATAC-seq peaks from mouse overlapping with SNP location. (C) Corresponding locations are shown in the mouse mm10 genome reference with IGV tracks for bulk ATAC-seq, RNA-seq, ABC model, and Cut&Run with antibodies against H3K27me3, H3K27ac, CTCF, H3K4me3 in IFN-γ-treated and Ctr-OPCs. Red boxes show CA peaks from mouse overlapping with SNP location and their ABC connections. Cell acronyms as in Figure 7.