Epigenetic signals that direct cell type-specific interferon beta response in mouse cells

Abstract

The antiviral response induced by type I interferon (IFN) via the JAK-STAT signaling cascade activates hundreds of IFN-stimulated genes (ISGs) across human and mouse tissues but varies between cell types. However, the links between the underlying epigenetic features and the ISG profile are not well understood. We mapped ISGs, binding sites of the STAT1 and STAT2 transcription factors, chromatin accessibility, and histone H3 lysine modification by acetylation (ac) and mono-/tri-methylation (me1, me3) in mouse embryonic stem cells and fibroblasts before and after IFNβ treatment. A large fraction of ISGs and STAT-binding sites was cell type specific with promoter binding of a STAT1/2 complex being a key driver of ISGs. Furthermore, STAT1/2 binding to putative enhancers induced ISGs as inferred from a chromatin co-accessibility analysis. STAT1/2 binding was dependent on the chromatin context and positively correlated with preexisting H3K4me1 and H3K27ac marks in an open chromatin state, whereas the presence of H3K27me3 had an inhibitory effect. Thus, chromatin features present before stimulation represent an additional regulatory layer for the cell type-specific antiviral response.

Graphical abstract

Figure 1.. ISG induction patterns in ESCs and MEFs.

(A) ESCs, NPCs differentiated in vitro from them and MEFs from the same mouse strain were studied to reveal the relation between cell type–specific chromatin features and IFNβ response. (B) Gene expression changes after IFNβ treatment. Red dots represent significant differentially expressed genes at P_adj < 0.05 and fold change ≥1.5 as computed with DESeq2. Four biological replicates for ESCs and two for MEFs were acquired for RNA-seq. Corresponding data for NPCs are shown in Fig S1A. (C) Overlap of all ISGs found at 1 h or 6 h of IFNβ treatment in ESCs and MEFs. (D) Uniform manifold approximation and projection for dimension reduction (UMAP) embedding of gene expression in ESCs (left) and MEFs (right) as two-dimensional representation of the transcriptome information. Each dot represents a single cell, and colors indicate IFNβ treatment duration. (E) Violin plots of scRNA-seq expression levels of the ISGs Ifit1 and Isg15 in single ESCs (top) and MEFs (bottom). For both genes, the number of transcripts detected largely increased in ESCs and MEFs from 1–6 h. (F) Expression levels of selected ISGs identified by bulk RNA-seq analysis from aggregated scRNA-seq in MEF clusters 0, 1, 2, and 3. Source data are available for this figure.

Figure S1.. Gene ontology (GO) terms and differential gene expression analysis of nascent RNA and single-cell data.

(A) Gene expression changes in NPCs after IFNβ treatment measured by RNA-seq. Red dots represent significant differentially expressed genes at P_adj < 0.05 and fold change ≥1.5 as computed with DESeq2. Four biological replicates were acquired. (B) Venn diagram of all ISGs found in NPCs either at 1 and 6 h of IFNβ treatment intersected with ISGs identified in ESCs and/or MEFs. (C) Enriched GO terms of all ISGs in ESCs, MEFs, and NPCs. Functional annotations of differentially expressed genes was done with the DAVID database (Huang et al, 2009) to identify GO terms. The visualization of the resulting terms was conducted with REVIGO (Supek et al, 2011). Orange terms were found in all three cell types, green terms in ESCs and NPCs, and violet terms were identified in the differentiated cell types MEFs and NPCs. Terms in grey were specific to a single cell type. (D) Analysis of intronic reads to identify differentially expressed genes on the nascent RNA level after 6 h of IFNβ treatment. Intronic reads were counted using HTSeq (Anders et al, 2015), and a modified GTF file containing only intronic sites was used for the differential gene expression analysis with DEseq2 (Love et al, 2014). Red dots represented differentially expressed genes with P_adj < 0.05 and fold change ≥1.5 as computed with DESeq2. A total of 82 (ESCs), 128 (NPCs), and 453 (MEFs) genes were up-regulated, whereas no down-regulated genes were detected. (E) Overlap of all (0 versus 6 h and 0 versus 1 h) ISGs detected by analysis of intronic reads detected in at least one differential gene expression analysis in ESCs, NPCs, and MEFs. (F) Single-cell embedding of gene expression in MEFs (UMAP). Each dot represented a single cell. Coloring depicts the score of principal component 2 (PC2) per cell, which was identified as main separator of the identified clusters. The plot shows that the clustering of MEFs into two groups was mainly driven by PC2. (G) Overrepresented KEGG pathways with positive (right, contribution >0.025) and negative contributors to PC2 (left; contribution < −0.025). Source data are available for this figure.

Figure S2.. Gene expression thresholds and expression of IFN signaling genes.

(A) Normalized gene expression levels (log) at 0 h IFNβ in ESCs (left), MEFs (middle), and NPCs (right). Frequency curves of gene expression were shown in blue. The data were fitted to two Gaussian distributions to represent actively expressed (green) and repressed (red) genes. Their intersection points were marked by a dotted line and define the thresholds used to distinguish actively expressed and repressed genes. (B) Normalized gene expression levels (TPM) of factors involved in IFN signaling. Gene expression levels of interferon receptors (Ifnar1, Ifnar2, Ifngr1, and Ifngr2), JAK-STAT cascade kinases (Jak1, Jak2, Tyk2, and Cdk8), and associated transcription factors (Stat1, Stat2, and Irf9) were shown for unstimulated (0 h) ESCs (green, n = 4), MEFs (purple, n = 2), and NPCs (magenta, n = 4). Significance was assessed with a Wilcoxon test at the P_adj < 0.05 level. The red line represents the calculated threshold in ESCs to distinguish actively expressed and repressed genes. Source data are available for this figure.

Figure 2.. Cell type–specific ISG induction and protein expression.

(A) Normalized gene expression levels of selected ISGs from bulk RNA-seq in ESCs (top, n = 4) and MEFs (bottom, n = 2). Gene expression is given as transcripts per kilobase million (TPM). (B) Western blots of IFNβ-stimulated ESCs and MEFs at 0, 1, and 6 h time points. The top row shows total levels of STAT1 (left) and STAT2 (right). The lower row shows phosphorylation of STAT1 at position 701 (left) and 727 (right). GAPDH was used as a housekeeping gene control. (C) Normalized gene expression levels from bulk RNA-seq of selected cell type–specific ISGs in ESCs (n = 4) and MEFs (n = 2). The red line represents a cell type–specific threshold to distinguish active and repressed genes. Top: expression of ISGs Ccnd2, Ifi27, and Nsg2 was only induced in ESCs. Bottom: expression of ISG Ccl2, Gbp6, and Ifit1bl1 was induced in MEFs. Source data are available for this figure.

Figure 3.. Binding of STAT1 and STAT2 in ESCs and MEFs.

(A) ChIP-seq of STAT1_p701 and STAT2 upstream of Stat1 (top) and at the Irf9 promoter (bottom). Tracks show one replicate for each condition. (B) STAT1_p701 and STAT2 peaks in ESCs and MEFs. The STAT1/2-binding sites were defined by the overlap of STAT1_p701 and STAT2 peaks from the combined list of peaks detected at 1 and 6 h of IFNβ treatment. Sample numbers are given in Table S1. (C) Enrichment of transcription factor binding motifs in STAT1p701, STAT1/2, and STAT2 peak sets identified in ESCs and MEFs. Motif color scheme: STAT-family (STAT1, STAT3, STAT3+IL21, STAT4, and STAT5), red; IRF-family (IRF1, IRF2, IRF3, IRF8, and ISRE [IRF9]), blue; other, black. Four biological replicates for ESCs and two for MEFs were analyzed. (D) Distribution of STAT1p701, STAT1/2, and STAT2 peaks at all annotated promoters, exons, introns, and intergenic regions annotated from the Ensembl database. (E) Overlap of STAT-binding sites between ESCs and MEFs for STAT1p701, STAT1/2, and STAT2. Source data are available for this figure.

Figure S3.. Cell type–specific binding of STAT1 and STAT2.

(A) Position weight matrices (PWMs) of the top five STAT- and IRF-family motifs identified in STAT1p701 and STAT2 ChIP-seq peaks based on HOMER annotation (known motifs). PWMs show the probability of each nucleotide on the y-axis and the position within the motif on the x-axis. For IRF1, IRF2, and IRF8 the reverse complementary sequences were similar to the ISRE and IRF9 annotated motif. The source of the motifs is indicated as HOMER (H) or JASPAR (J) database. (B) Top de novo identified motif by HOMER for each subset (STAT1, STAT1/2, and STAT2) of the STAT1p701 and STAT2 ChIP-seq peak overlaps. The de novo motifs were presented as PWMs with the probability of each nucleotide on the y-axis and the position within the motif on the x-axis. The log P-value, the percentage of the motifs within the target or in background, and best match and similarity scores were calculated with the HOMER de novo motif annotation. (C) Visualization of the HOMER similarity scores of the top 10 most similar motifs to the identified de novo motifs from (B). (D) Scatter plot of normalized gene expression (log₁₀ of TPM) before (0 h, x-axis) and after 6 h (y-axis) IFNβ stimulation. Genes with STAT_p701-, STAT2-, or STAT1/2-binding sites at promoters in ESCs (top) and MEFs (bottom) were shown. Red dots indicate ISGs identified by differential gene expression analysis of RNA-seq data, whereas black dots showed STAT-bound genes at promoters with no significant changes of expression. Source data are available for this figure.

Figure S4.. Co-accessible STAT1/2-binding sites and distal target genes.

(A) Overlap between ISGs (blue) and the nearest gene to a STAT1/2 bound site (grey) in ESCs (left) and MEFs (right). The nearest gene list was calculated with GREAT (McLean et al, 2010). Many STAT1/2 sites were nearest to a gene that is not induced by IFNβ. (B) Chromatin accessibility at STAT1/2 bound sites at 0 h IFNβ- (black), 1 h IFNβ- (red), and 6 h IFNβ (blue)-treated ESCs (top) and MEFs (bottom) in bulk (left) and scATAC-seq data (right). A general increase of accessibility at STAT1/2 bound sites upon IFNβ treatment was apparent except for MEF IFNβ 6 h scATAC-seq data. (C) Venn diagram of STAT1/2-mediated ISG regulation mechanisms in different cell types. (D) ISGs with a co-accessible link to another distal STAT1/2 bound ISG promoter and/or an exonic, intronic, or intergenic STAT1/2 bound bona fide enhancer in ESCs (left), epithelial-like MEFs (mid), and mesenchymal-like MEFs (right). (E) ISG expression (TPM) before induction computed from the bulk RNA-seq data for ISGs with different mechanisms of regulation by STAT1/2 binding in the indicated cell types. Promoter bound STAT1/2 (independent of the presence of additional links to distal sites), blue; gained co-accessible link to a distal STAT1/2 peak, green; ISGs that lost a co-accessible link to a distal STAT1/2 peak after IFNβ treatment, red; other ISGs, grey. Source data are available for this figure.

Figure 4.. Regulation of ISG expression by distal STAT1/2 binding.

(A) UMAP embedding of chromatin accessibility in ESCs (left) and MEFs (right). Each dot represents one cell and is colored according to treatment. (B) Same as panel (A) for MEFs with single cell coloring according to k-nearest neighbor clusters. (C) Same as panel (B) with single cell coloring according to MEF subtypes derived from scRNA-seq data integration by gene activities. (D) Co-accessibility before and after 6 h of IFNβ induction of ESCs in a region around the Uba7 ISG. Top: browser tracks of aggregated pseudo-bulk chromatin accessibility from single cells. Middle: co-accessible links between the indicated intronic STAT1/2 bound site 371 (differential STAT1/2 peak after 1 h of IFNβ treatment in ESCs) and other genomic loci. Experimentally identified ISG promoters (blue) and sites with bound STAT1/2 after 1 and/or 6 h (green) are marked. Bottom: gene expression levels from scRNA-seq. Transcription from Inka1 and Rnf123 was not detected. (E) Same as panel (D) but for three intergenic STAT1/2 bound sites 125, 126, and 127 (differential STAT1/2 peaks after 1 and 6 h of IFNβ treatment in MEFs) in the Ly6 ISG cluster in MEFs. (F) ISG regulation mechanisms according to STAT1/2 binding after IFNβ treatment. Promoter bound STAT1/2 (independent of the presence of additional links to distal sites), blue; gained co-accessible link to a distal STAT1/2 peak, green; ISGs that lost a co-accessible link to a distal STAT1/2 peak after IFNβ treatment, red; other ISGs, grey. (G) Expression changes of ISGs for the different STAT1/2-dependent regulation types shown in panel (F) from bulk RNA-seq data. P-values from a Wilcoxon rank-sum test are indicated as *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. Source data are available for this figure.

Figure S5.. Clustering of chromatin features at STAT1/2-binding sites.

(A) STAT1/2 binding sites used in the chromatin context analysis. A combined set of 392 STAT1/2-binding sites was obtained from the STAT ChIP-seq analysis in ESCs and MEFs at 1 and 6 h of IFNβ treatment. Chromatin features in a genomic region of ±1 kb around the binding sites center were analyzed. STAT1/2-binding sites at promoters were annotated according to transcription start sites in the Ensembl database. (B) Silhouette score calculated as the mean Silhouette coefficient over all samples plotted against the cluster number. When varying the cluster number between 2 and 20 clusters, it was seen that the selected number of five clusters was appropriate for this data set. (C) Heatmap with unsupervised k-means clustering of histone modifications (H3K4me1, H3K4me3, H3K27ac, H3K9ac, H3K27me3, and H3K9me3) and chromatin accessibility (ATAC) data from unstimulated ESCs and MEFs at 392 STAT1/2-binding sites. Five biologically relevant clusters with distinct signatures were identified and annotated as “Active Promoter,” “Active Enhancer,” “Bivalent,” “Poised,” and “Repressed” chromatin states. Source data are available for this figure.

Figure 5.. Contribution of chromatin features to STAT1/2 binding.

(A) Genomic regions around the ISGs Ifi27, Usp18, and Gbp6 in ESCs (top) and MEFs (bottom) with the different sequencing readouts and the promoter regions marked by boxes. Gene annotation was based on Ensembl, and the positions of the DNA binding motif IRSE were extracted from the HOMER database. Each browser track shows one representative biological replicate. (B) Heatmap of unsupervised k-means clustering of histone modifications and ATAC data at 392 STAT1/2-binding sites. The indicated five main chromatin states were identified. Data from unstimulated ESCs and MEFs and ESCs, at 1 and 6 h IFNβ treatment were used. Sample numbers are given in Table S1. (C) Chromatin state comparison between untreated ESCs and MEFs at STAT1/2-binding sites based on the data in panel (B) and corresponding coloring of the five different chromatin states. The lines link the same binding sites between conditions and do not represent a differentiation path between ESCs and MEFs. (D) Absolute numbers of STAT1/2-binding sites according to chromatin states in unstimulated ESCs and MEFs. (E) Distribution of 116 ESC-specific and 184 MEF-specific STAT1/2-binding sites according to the chromatin state. Source data are available for this figure.

Figure 6.. Correlation of STAT1/2 binding with preexisting chromatin features.

(A) Correlation between STAT1/2 binding after 1 h of IFNβ treatment and preexisting chromatin features before IFNβ treatment. The STAT1/2-binding signal was computed as the average signal of STAT1 and STAT2 after 1 h IFNβ treatment in ESCs (top) and MEFs (bottom). The chromatin features were quantified by counting the normalized read counts at the STAT1/2-binding sites before induction. ESC-specific STAT1/2-binding sites are shown in black and MEF-specific ones in red. Ellipses indicate the area, in which 75% of all data points are located. Density distributions are shown along the x- and y-axis. The blue line shows the linear regression of the combined set of ESC- and MEF-specific STAT1/2-binding sites. Sample numbers are given in Table S1. (B) Correlation between STAT1/2 binding and chromatin features determined for the data in panel (A) in ESCs (black) and MEFs (red). (C) Scheme of ISG induction via chromatin context-dependent STAT1/2 binding. Binding to ISREs at promoters or enhancers can be facilitated or repressed via preexisting chromatin states marked by the indicated chromatin features that can differ between cell types. ISGs that lack an ISRE and STAT1/2 binding at the promoter can also be activated by STAT1/2 binding to distal enhancers. Source data are available for this figure.