Properties of STAT1 and IRF1 enhancers and the influence of SNPs
- Mohamed Abou El Hassan†1, 2, 3,
- Katherine Huang†1,
- Manoja B. K. Eswara1,
- Zhaodong Xu1,
- Tao Yu1,
- Arthur Aubry1,
- Zuyao Ni1, 6,
- Izzy Livne-bar1,
- Monika Sangwan1,
- Mohamad Ahmad1 and
- Rod Bremner1, 4, 5Email author
© The Author(s) 2017
Received: 20 December 2016
Accepted: 2 March 2017
Published: 9 March 2017
STAT1 and IRF1 collaborate to induce interferon-γ (IFNγ) stimulated genes (ISGs), but the extent to which they act alone or together is unclear. The effect of single nucleotide polymorphisms (SNPs) on in vivo binding is also largely unknown.
We show that IRF1 binds at proximal or distant ISG sites twice as often as STAT1, increasing to sixfold at the MHC class I locus. STAT1 almost always bound with IRF1, while most IRF1 binding events were isolated. Dual binding sites at remote or proximal enhancers distinguished ISGs that were responsive to IFNγ versus cell-specific resistant ISGs, which showed fewer and mainly single binding events. Surprisingly, inducibility in one cell type predicted ISG-responsiveness in other cells. Several dbSNPs overlapped with STAT1 and IRF1 binding motifs, and we developed methodology to rapidly assess their effects. We show that in silico prediction of SNP effects accurately reflects altered binding both in vitro and in vivo.
These data reveal broad cooperation between STAT1 and IRF1, explain cell type specific differences in ISG-responsiveness, and identify genetic variants that may participate in the pathogenesis of immune disorders.
IFNγ is a pleiotropic cytokine that plays essential roles in antiviral and anticancer immune responses (reviewed in [1, 2]). IFNγ binds to its receptor complex and activates receptor-associated JAK kinases, which phosphorylate a substantial fraction of cytoplasmic signal transducer and activator of transcription 1 (STAT1). Phosphorylated STAT1 forms homodimers that translocate to the nucleus and bind IFNγ activation sites (GAS). STAT1 recruits histone acetyltransferases (HATs) and other transcriptional co-activators to acetylate chromatin and facilitate transcription. Genomic studies showed that STAT1 binds at promoter proximal and distal sites, suggesting a role in remote gene regulation [3–6]. Indeed, IFNγ induces long range interactions between STAT1-bound enhancers and target promoters [7–9].
Interferon regulatory factor 1 (IRF1) is a primary target gene of STAT1. Like STAT1, IRF1 also acts as a transcription factor (TF), binding to IRF-E motifs and interferon-stimulated response elements (ISRE) [10, 11]. Access of both STAT1 and IRF1 to target enhancers requires the SWI/SNF chromatin remodeling complex to counter PRC2, which uses the histone methyl transferase EZH2 to deposit H3K27me3 and block the induction of many other cytokine and cytokine responsive loci [7, 12, 13]. IRF1 functions at the transcription initiation level by facilitating RNA Pol II recruitment to ISGs promoters [14, 15]. IRF1 also binds to remote enhancers of the CIITA locus that loop together to form a 3D interconnected hub with the promoter . Indeed, ChIP-chip and ChIP-seq studies show that IRF1 binds many remote enhancers [6, 16–18], and analysis of 128 transcription factors in K562 cells revealed that STAT1-IRF1 co-binding is a recurring pattern in IFNγ treated cells . Notably, STAT1 is essential but not sufficient for gene induction , and both STAT1 and IRF1 are required for the IFNγ-induced expression of CIITA, GBP1, and gp19 [14, 15, 20]. In addition, STAT1 complexes with IRF1 at the LMP2 promoter and maintains its constitutive expression .
Here, we studied the extent of STAT1 and IRF1 cooperation in HeLa cells within ISG-rich chromosomal segments encompassing ~10% of all known ISGs. Most of these loci responded to IFNγ in HeLa cells, leaving ~20% resistant ISGs. IRF1 binding sites outnumbered STAT1 sites 2 to 1. A large fraction of STAT1/IRF1 binding occurred at remote sites and looping studies confirmed the functional role of putative enhancers at the SOCS1 locus. Most STAT1 binding occurred at or near to IRF1 sites (dual binding), but IRF1 often bound isolated from STAT1. Dual STAT1 and IRF1 but not isolated IRF1 or STAT1 binding was linked to ISG responsiveness. Finally, several variants affecting STAT1/IRF1 motifs induce or impair binding.
Diverse gene responses to IFNγ
Signaling pathway target loci show cell-specific responsiveness, but the exact TF binding patterns that distinguish induction versus resistance in a specific cell type are unclear. Thus, we compared the pattern of STAT1 and IRF1 binding at different gene types. For this we compiled a database of ~ all known ISGs using our own and prior transcriptome data (Additional file 1: Table S2). As summarized in Fig. 1b, ISGs fell into 8 classes depending on whether IFNγ caused induction, no effect (resistant ISGs in HeLa cells), or repression, and whether induction/repression were early (detected at 6 h) or late (24 or 48 h), and strong (differential score ≥13, and ≥twofold change) or weak (differential score ≥13, <twofold). The microarray expression data was validated using reverse transcription and quantitative PCR (RT-qPCR), which confirmed 83% (20/24) of indISGs and 87% (21/24) of resISGs (Fig. 1c). Of all 95 known ISGs on the array, 31 (33%) genes were es-indISGs, 15 (16%) were ls-indISGs, 29 (31%) were ew-indISGs or lw-indISGs, and 20 were resISGs (Fig. 1b). Es-indISGs were distributed at an average density of 1.9/Mb within the studied regions (Additional file 1: Table S1). The highest density was observed at the IFIT and GBP clusters with an average of 4.0 es-indISGs/Mb.
No genes were repressed (IFNγ repressed genes, IRGs) at the early 6 h time point, while 19 and 23 were strongly or weakly repressed at later times, respectively (ls-IRGs and lw-IRGs; Fig. 1b), suggesting indirect regulation of IRGs (perhaps through activation of a repressor). The remaining genes that were not IFNγ-responsive either in this or any prior study were termed “Other Genes”. In summary, known ISGs fall into induced and resistant subclasses in HeLa cells, providing a useful system to define STAT1 and IRF1 binding patterns linked to responsiveness.
Validation of STAT1 and IRF1 ChIP-chip analyses
Basal TF binding
Unphosphorylated STAT1 has roles in regulating ISGs days after IFN treatment [23, 24], but its role in untreated cells is less clear, although STAT1 nuclear cytoplasmic shuttling occurs even in untreated cells [25–27]. Basal STAT1 binding is linked to the nuclear localization of unphosphorylated STAT1 and contributes to the constitutive expression of some targets [21, 28]. IRF1 is also expressed to low levels in unstimulated HeLa cells  and it cooperates with STAT1 to maintain low basal expression levels of LMP2 . In addition, there is also some STAT1 phosphorylation (below detectable levels) in untreated cells that contributes to basal activity, as shown elegantly by knockin studies in mice . We detected 2 STAT1 and 28 IRF1 binding sites in untreated cells, accounting for 2.2 and 14.3% of induced sites, respectively. Our data accords with another ChIP-chip analysis of STAT1 binding which reported that 6.5% of IFNγ-induced STAT1 sites in HeLa cells treated for 30 min with IFNγ (as opposed to 6 h in our case) are occupied in uninduced cells .
Further analysis suggests that basal TF binding detected here is physiologically relevant. Of the genes with basal STAT1 or IRF1 binding, we assessed 21 by microarray and/or RT-PCR and all were expressed in untreated cells (Additional file 1: Table S3). In contrast, of 26 randomly selected ISGs that lacked basal TF binding, only 13 were basally expressed. Indeed, constitutive expression of PSMB9 and TAP2 requires constitutive IRF1 binding [21, 30]. In addition, several loci with basal TF binding are in paralogous gene clusters suggesting conservation of high affinity binding sites during gene duplication (e.g. PSMB8 and PSMB9, GBP2 and GBP3, and IFIT1, IFIT2 and IFIT3). A high fraction (82%) of the 28 IRF1 basally occupied sites possessed IRF1 binding motifs. Thus, our data supports the notion that basal binding of STAT1 and IRF1 is physiologically relevant.
Remote IFNγ activated enhancers are common at ISGs
Unusual IRF1 distribution at MHC loci
As noted earlier, IRF1 exceeded STAT1 sites by ~twofold, but this varied at some regions, most notably at the MHC class I locus where the ratio was 3.5:1 (56 IRF1:16 STAT1 sites; Additional file 1: Table S4). The ratio was particularly skewed at the extended (6:1) versus classic (2.9:1) MHC class I region. 26 of all MHC class I IRF1 sites were within 5 kb of Known Gene starts and 17 within 5 kb of pseudogenes (Additional file 1: Tables S4, S5), giving a total of 77% (43/56) promoter proximal sites, which is higher than the 44% at all loci (Fig. 3a; Additional file 1: Table S4). However, whereas 66% (25/38) of IRF1 sites were promoter-proximal in the classical MHC class I region, this dropped to only 6% (1/18) at the extended MHC class I region, and was low even after including pseudogenes (4/18; Additional file 1: Tables S4, S5), leaving an unusually high fraction of remote IRF1 sites (78%). Thus, IRF1 seems to play a broader role than STAT1 at the MHC class I cluster, primarily at proximal elements in the classic region, but at remote elements in the extended region.
STAT1 and IRF1 induce CIITA, the master regulator of MHC class II expression (reviewed in ). The number of STAT1 and IRF1 sites was typically very low in the MHC class II region (Additional file 1: Table S4). Out of 13 MHC class II genes, 5 (DRB5, DQB1, DQB2, DQA2, and DOA) were resistant to IFNγ in HeLa cells, 5 (DOB, DRB1, DQA1, DPA1, DPB1) responded only after 24 h, a time of maximum production of CIITA , and 3 (DRA, DMB and DMA) were es-indISGS. With the exception of DOB, none of the resistant or late-induced genes exhibited STAT1 or IRF1 promoter binding. However, of the 3 es-indISGs, two had promoter proximal IRF1 binding while DMA had IRF1 binding fairly near (~8 kb) its promoter. Thus IRF1 may cooperate with CIITA at a subset of MHC class II promoters. Others reported CIITA-independent induction of MHC class II genes [36–39], which may, therefore, involve IRF1.
STAT1 and IRF1 binding is enriched at robustly induced ISGs
STAT1 and IRF1 binding sites were most highly associated with robustly induced IFNγ targets (es-indISGs; Fig. 6). This applied when STAT1 or IRF1 were considered together, separately, and at proximal or remote locations (Fig. 6b). Consistent with this finding, weakly induced genes (ew-indISGs) had fewer binding events and lower TERs (Fig. 6). Of 236 Other genes (never classified as an ISG in any study), a total of only 17 had 9 STAT1 and 16 IRF1 proximal peaks, mostly (13/17) located at the MHC and RT-qPCR confirmed no induction at 10/10 of these genes (Additional file 1: Table S6). IFNγ enhancers loop over large distances at CIITA  and SOCS1 (Fig. 5), so proximal and distal enhancers nearest to Other Genes may target neighboring ISGs. In summary, the data indicate a clear bias of STAT1 and IRF1 binding at rapidly and robustly induced ISGs, but not other gene classes.
Isolated or dual STAT1 and IRF1 recruitment is directed by binding motifs
JASPER analysis of IRF1 and STAT1 peak regions revealed that the cognate binding motif was observed at a statistically significant level relative to equal numbers of random peaks (Fig. 7b). 60% of isolated STAT1 peaks had a STAT1 motif, and only 30% had an IRF1 motif, while 70% of isolated IRF1 peaks possessed an IRF1 motif, but only 20% had a STAT1 motif. A strong correlation existed between STAT1/IRF1 binding and the presence of the corresponding motifs (Fig. 7c). Indeed 40% of dual STAT1/IRF1 sites had both binding motifs, whereas there were none at equal numbers of randomly generated sites (Fig. 7b). Dual sites which have only a STAT1 or IRF1 binding motif may reflect protein–protein interaction or DNA looping as seen at the SOCS1 and CIITA loci (Fig. 4c) [7, 8]. In summary, DNA sequence directs isolated or dual STAT1/IRF1 binding in IFNγ treated cells.
Dual STAT1 and IRF1 targeted enhancers distinguish responsive from resistant ISGs
Degree of TF binding and responsiveness in HeLa predicts ISG responsiveness in other cell types
We assessed the relationship between broad responsiveness, degree of induction, and STAT1/IRF1 binding. HeLa ChIP-chip data provided STAT1 and IRF1 binding information for 24 es-indISGs present in all 7 expression array datasets. Of these, 3/24 were induced exclusively in HeLa, 10/24 were induced in 2–4 lines and 11/24 were induced in 5–7 lines (Fig. 9b). Of note, genes induced exclusively in HeLa were up-regulated to a much lower extent than ubiquitously IFNγ-responsive targets (Fig. 9c). Greater induction of ubiquitously responsive loci was paralleled by a higher density of TF binding at promoter proximal sites (Fig. 9d). Thus, the level of induction is linked to the degree of STAT1 and IRF1 recruitment, and there is an unexpected link between the strength of ISG induction in one context (HeLa in this case) and competency to respond to IFNγ in other contexts.
SNPs modulate STAT1 and IRF1 binding in vitro
Defects in IFNγ signaling are linked to a wide range of disorders [40–44]. Several studies focused on the association between genetic variants and the risk of IFNγ related disorders, but at gene promoters or coding regions of ISGs rather than IFNγ responsive enhancers. Within the 16 Mb of DNA around ISGs studied here, there are a total of 7.1 × 105 dbSNPs [hg19; SNPs (141)]. Of these, 6648 dbSNPs lay within the 230 STAT1/IRF1 peaks. Only 7 of these 6648 dbSNPs were listed on the GWAS database. GWAS SNPs do not define all disease associated SNPs (DA-SNPs) because GWAS genotyping arrays provide low genomic coverage  and therefore the 6648 dbSNPs may encompass other DA-SNPs not mapped yet. None of the 7 DA-SNPs (GWAS database) overlapped with a STAT1/IRF1 motif, but 80 of the 6648 dbSNPs overlapped with 27 STAT1 and 47 IRF1 motifs (Additional file 1: Table S9).
We studied which of these 80 SNPs affect STAT1/IRF1 binding. First, we utilized the CisGenome “Known Motif Mapping” program to predict which of the variants may modulate STAT1/IRF1 binding motifs (see “Methods” section). CisGenome compares the position weight matrix (PWM) in the JASPAR CORE database and creates likelihood scores for the reference or variant allele. We calculated the fold change in likelihood scores (variant/reference allele) to assess the predicted relative effect. At a cutoff of 1.5-fold, the variant alleles of 34/80 dbSNPs were predicted to modulate the binding affinity of 24 IRF1 motifs and 10 STAT1 motifs (Additional file 1: Table S9).
rs9260102 affects IRF1 binding in vivo
STAT1 and IRF1 drive the induction of IFN induced genes, but the extent to which they act collectively is unclear. We report that most STAT1 binding (60%) occurs together with IRF1, but most IRF1 binding (72%) is isolated (Fig. 7a). Binding occurs where there are cognate binding motifs (Fig. 7), suggesting that most ChIP signals reflect direct recruitment. Both proximal and remote STAT1 and IRF1 binding is observed at robustly induced ISGs, but not at other loci (Fig. 6). In line with the importance of TF occupancy for responsiveness , every responsive locus exhibits a mixture of STAT1 and IRF1 bound enhancers (Fig. 8a). Moreover, dual bound enhancers distinguish induced vs resistant ISGs, whereas single bound enhancers are found with similar frequency at responsive or non-responsive ISGs (Fig. 8). This is not to say, however, that single bound enhancers are irrelevant. For example, while multiple remote SOCS1 enhancers recruit both TFs, the +50 kb element or promoter are targeted only by STAT1 or IRF1, respectively, yet both are involved in 3D looping (Figs. 4, 5). Similarly, while dual STAT1/IRF1 binding occurs at the active CIITA promoter, the −50 kb and +59 kb enhancers recruit only STAT1 or IRF1, respectively, yet contribute to 3D looping and are essential for responsiveness . Indeed, all the responsive genes we surveyed exhibit a mix of STAT1-only, IRF1-only, and STAT1/IRF1 dual enhancers (Fig. 8). Together, these results suggest that IFNγ-responsiveness requires cooperation between enhancers that bind both or either TF, but that STAT1- or IRF1-only enhancers are insufficient for gene induction. Irrespective, it is clear that responsive ISGs integrate information from both STAT1 and IRF1.
Previously, we showed that there is a pre-existing 3D structure at the silent CIITA locus, generated through looping between enhancers that subsequently recruit STAT1 and IRF1 upon IFNγ treatment . This was true even in the absence of BRG1, a chromatin remodeling enzyme that is critical to allow stable TF recruitment and thus IFNγ-responsiveness. Subsequent genome-wide analyses indicate that enhancer looping in the poised but silent state is common at inducible loci . We observed the same phenomenon at the IFNγ responsive SOCS1 locus (Fig. 5). Potentially, these contacts are mediated by pioneer factors that mark responsive enhancers, but their identity at ISGs is unknown. The data here and other studies show that STAT1 and IRF1 can bind some sites in the basal state [21, 28], so in theory, low/unstable binding (undetectable by ChIP) could poise ISG enhancers. It would thus be interesting to perform looping studies at ISGs in STAT1/IRF1 deficient cells. It is of note that the degree to which ISGs were induced in HeLa cells predicted whether they were likely to respond to IFNγ in other cells (Fig. 9). Thus, the chromatin at these genes is accessible in many contexts. The factors that mediate this broad poised, open state may also initiate the basal looping at ISGs.
Over 90% of the disease markers identified in GWAS studies lie within the non-protein-coding regions of the genome . These markers correlate with gene expression [49–52], and lie within gene regulatory regions [53–56]. There is thus considerable interest in identifying SNPs that influence TF binding and, therefore, gene regulation. We identified 80 SNPs within STAT1 or IRF1 motifs, and in silico assessment predicted that 34 may alter binding. In vitro quantification confirmed these predictions in 5/5 cases, arguing that most of/all the other predictions are accurate. The availability of a cell line heterozygous for one such SNP allowed us to test whether the prediction held up in vivo. Indeed, the T allele of rs9260102, which lies just upstream of the HLA-A locus, bound IRF1 whereas the G allele did not, as observed in silico and in vitro. These data serve as proof of principle that in silico prediction is a reliable tool to anticipate the effect of SNPs on STAT1 and IRF1 binding.
This study provides strong evidence for widespread cooperation between STAT1 and IRF1 at ISGs, and suggests that in silico predictions reliably predict the effect of nucleotide variants on binding in vivo.
Custom oligonucleotide ChIP Tiling array design
A custom oligonucleotide tiling array was designed to cover 11 genomic regions spanning a total of 16 Mb of human genomic DNA in 8 chromosomes (Additional file 1: Table S1). Regions covered from 1 to 5 Mb genomic sequences. Arrays consisted of 50 mers, in quadruplicate, with median probe spacing of 80 bp within non-repetitive DNA regions.
ChIP on tiled genome arrays (ChIP-chip) and ChIP–quantitative PCR
Details of primers and antibodies used in ChIP assays are in Additional file 1: Tables S11, S12. HeLa-ini1-11 cells (HeLa), were grown as described . Cells were left untreated or exposed to 300 units/ml of human IFN-γ for 6 h (PHC4834, BioSource International, Camarillo, CA, USA). Crosslinked chromatin was sonicated to an average size of about 500 base pairs and was incubated with STAT1 or IRF1 antibody. Bound fragments were purified by ChIP and amplified by ligation-mediated PCR, as described , then labeled and hybridized to the arrays. Hybridization intensities were normalized to internal standards and values from quadruplicate spots were averaged. Significantly different intensities between ChIP DNA and input DNA samples in three biological replicates (p < 0.0001) were determined with the Wilcoxon rank-sum test. Peaks representing significantly enriched DNA regions (p < 0.0001) where the ratio of ChIP to input DNA was 1.5-fold or more were visualized with the University of California at Santa Cruz Human (Homo sapiens) Genome Browser (Phast-Cons) and are plotted on a log2 scale. Peaks in a sliding window of 500 base pairs were merged with an in-house Perl script pipeline. ChIP—quantitative PCR was done as described , and in all cases, the low background signal obtained with a no-antibody control was subtracted.
Custom oligonucleotide ChIP tiling array data analysis
Raw intensities from three independent biological replicates, quality assessed by Nimbelgen SignalMap software, were quantile normalized , and averaged for each quadruplicate 50 mer. We developed a Wilcoxon Rank Sum test  based software to studying the difference between the intensities of the ChIP signal compared to the input DNA signal for each probe within a 500 bp sliding window. Genomic positions with statistically higher intensities (>1.5-fold, p < 10−4) from input DNA were merged to form a peak. ChIP-chip data were imported into UCSC genome browser (assembly hg17, NCBI build 35) as two sets of separate tracks for each antibody before and 6 h after IFNγ treatment (http://research.lunenfeld.ca/IFNy).
STAT1 and IRF1 motif analysis
We mapped STAT1 and IRF1 consensus motifs to the STAT1 and IRF1 ChIP-chip binding regions by using CisGenome “Known Motif Mapping” program. Motif occurrences were determined by using position frequency matrices (PFMs) of STAT1 (ID: MA0137.2) and IRF1 (ID: MA0050.1) from the JASPAR CORE database. The PFMs were converted to pseudo-count matrix for CisGenome’s input. A motif mapping location is selected by the cutoff of a likelihood ratio (LR) > 500. The LR is determined by comparing the motif’s PFM with a background model estimated from input ChIP-chip regions in the 3rd order. Seven sets of ChIP-chip peak regions were mapped and compared with background control regions: STAT1 peaks (92), IRF1 peaks (196), merged STAT1 and/or IRF1 peaks (231), dual STAT1 and IRF1 peaks (56), STAT1 peaks isolated from IRF1 peaks (isolated-STAT1, 56), IRF1 peaks isolated from STAT1 peaks (isolated-IRF1, 140), and basal IRF1 peaks. For each set of peaks, the same number and size of control background regions were randomly sampled from blank regions without any ChIP-chip bindings following the same frequency distribution as the real binding peaks within each cluster segments on the chromosome of ChIP-chip data. The frequencies of regions mapped with motifs were compared between ChIP-chip peaks and random control sites by using the two-sided probability test in R for each paired set of peaks.
RNA extraction, expression microarray analysis, and Reverse transcriptase-qPCR
RNA extraction and reverse transcription were done from HeLa cells left untreated or at 6, 24 or 48 h after IFNγ treatment as described previously . RNA quality was checked using both Nanodrop (Thermo Fischer Scientific; 260/280 ratio was ≥1.8) and Bioanalyzer (Agilent Inc.; RNA Integrity Number, RIN, ≥ 9.4, range 9.4–9.9). RNA samples were converted to cDNA, followed by a second strand synthesis, and cRNA was prepared using the Ambion kit (Applied Biosystems). The cRNA was column purified and quality was checked using Bioanalyzer (Agilent Inc.). A total of 1.5 µg of cRNA was hybridized to human whole-genome expression arrays (HumanRef-6 Expression BeadChip, Illumina, Inc.) using standard Illumina protocols. Slides were scanned on an Illumina Beadstation and analyzed using BeadStudio (Illumina, Inc). Genes induced by ≥twofold compared to controls and that achieved a differential score of ≥13 were classified as strongly induced ISGs. ISGs which achieved a differential score of ≥13 but fold induction less than 2 were considered weakly induced. Genes reduced by ≥twofold compared to control and achieved a differential score of ≤ −13 were considered strongly reduced. Genes that had a differential score of ≤ −13 but the fold reduction was less than 2 were considered weakly reduced. Three biological replicates were included for each treatment group.
RT-qPCR was performed much as described . Briefly, RNA was extracted from HeLa cells left untreated or at 6, 24 or 48 h after IFNγ treatment using Trizol (Invitrogen), and quality assessed by RIN and OD260/280 as above. cDNA was prepared from 1 mg RNA using random primers and SuperScript RT (Invitrogen). Amplification of cDNA was performed using gene specific primer pairs and SYBER Green Mix (ABI). PCR was ran on Applied Biosystems PRISM 7900HT. Primers were designed in the coding region of each gene (Additional file 1: Table S12). Human genomic DNA was used to prepare calibrators for the quantification of cDNAs. Dissociation curves were inspected to ensure a single product and all PCR products were also tested on a gel to confirm amplification specificity. In addition, no template controls (NTC) were included to ensure the absence of DNA contamination. Gene expression was normalized to multiple house-keeping/reference genes to control for the total amount of RNA. All experiments were done in triplicate.
Chromosome conformation capture (3C)
Assessment of TF binding distribution around different classes of ISGs
We compared the distribution of TFBS in the vicinity of es-indISGs and resISGs in STAT1 peaks, IRF1 peaks, and merged STAT1 and IRF1 peaks. We first aligned the TFBS within a range of 400 kb region around the TSS of each gene in a 1 kb resolution, which means that we divided each region into 1 kb windows with the window 0 centered at the TSS and others line up to the 200 kb end upstream and 200 kb end downstream, and then scored the frequency of TFBS at each window as the number of peaks whose center was within the window. If 400 kb extended beyond the ChIP-chip segments, the binding frequencies along the truncated regions were regarded as missing data. Then we plotted the average binding frequencies per 1 kb window versus the relative distance of each to the TSS. Missing values were discarded for averaging the frequencies.
Defining STAT1 and IRF1 functional SNPs
First, we queried dbSNPs located within the 230 ChIP-chip peaks (UCSC, Build hg19; Track, All SNPs(141); Table, snp141). We defined a total of 6648 dbSNPs. Next we defined SNPs that overlap with STAT1 or IRF1 binding motifs within the 230 ChIP-chip peaks. Then we selected a region of ±50 bp around the SNPs that overlapped with STAT1/IRF1 binding motifs and recovered the DNA sequence of these regions using CisGenome. Then we computationally evaluated the binding affinity of the reference or variant sequence using the likelihood scores obtained from the Cisgenome “known motif mapping” program with the sequences as input to map the STAT1 and IRF1 motif matrix. In some cases the introduction of the variant SNP renders the motif unidentifiable and in this case the sequence of the motif was indicated as “NULL” and the likelihood score was considered as zero (Additional file 1: Table S9). The cutoff value of affinity change was set at 1.5-fold.
ELISA-based DNA binding affinity assay and ChIP coupled with DNA sequencing
We designed 33-mers with either the reference or variant alleles (Additional file 1: Table S13). Control probes with wild type or dead mutant STAT1 or IRF1 motifs were also included. Probes were ordered biotinylated or biotin-free (competitors). Two pmol of biotinylated probes were immobilized per well of 96-well streptavidin-coated plates. Cell lysates where incubated with different concentrations of the competitor probes (probe-lysate mix) at 4 °C for 3 h to allow STAT1 or IRF1 binding. The probe-lysate mix was then added to streptavidin-coated plates with immobilized biotin probes and incubated overnight at 4 °C. To quantify bound TFs, wells were washed and probed with STAT1 or IRF1 primary antibodies, followed by IR-800 conjugated secondary antibodies. Excess antibodies were washed thoroughly and plates were scanned and quantified using Odyssey Infrared imaging system (LICOR). Signal from no-competitor well is considered as 100% and the % antibody signal was plotted against competitor probe concentration. IC50 values were calculated using Graphpad PRISM 5.2.
For in vivo studies, EBV-transformed lymphoblastic GM18857 cells, cultured as recommended by the supplier (Coriell Biorepositories), were treated with IFNγ for 6 h, fixed and harvested for ChIP analysis. Chromatin was immunoprecipitated using IRF1 antibody and isolated DNA was sequenced using Snapshot sequencing.
chromosome conformation capture
IFNγ activation site
interferon regulatory factor 1
interferon-stimulated response element
IFNγ stimulated genes
early strong induced ISGs
late strong induced ISGS
early weak induced ISGs
late weak induced ISGs
IFNγ repressed genes
late strong IRGs
late weak IRGs
single nucleotide polymorphisms
signal transducer and activator of transcription 1
transcription start site
MAEH generated and analyzed most of the data and KH performed most of the bioinformatics. MBKE and MAEH assessed TF binding in vitro, and were helped by AA. ZX helped with the bioinformatics. TY and ZN helped MAEH with ChIP. MS, MA, and IL helped generate reagents and analyze data. RB obtained funding, analyzed data, and wrote the paper with MAEH and KH. All authors read and approved the manuscript.
We thank Malcolm Macleod for designing the software to create TF binding maps.
Funding was provided by Canadian Institutes of Health Research (Grant No. MOP 111004); Canadian Cancer Society Research Institute (Grant No. 703079).
The authors declare that they have no competing interests.
Availability of data and materials
This work was funded by grants to R.B. from the Canadian Cancer Society Research Institute (CCSRI), the Canadian Health Research Institutes (CIHR), and the Krembil Foundation.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Vesely MD, Kershaw MH, Schreiber RD, Smyth MJ. Natural innate and adaptive immunity to cancer. Annu Rev Immunol. 2011;29:235–71.View ArticlePubMedGoogle Scholar
- Maher SG, Romero-Weaver AL, Scarzello AJ, Gamero AM. Interferon: cellular executioner or white knight? Curr Med Chem. 2007;14:1279–89.View ArticlePubMedGoogle Scholar
- Hartman SE, Bertone P, Nath AK, Royce TE, Gerstein M, Weissman S, et al. Global changes in STAT target selection and transcription regulation upon interferon treatments. Genes Dev. 2005;19:2953–68.View ArticlePubMedPubMed CentralGoogle Scholar
- Robertson G, Hirst M, Bainbridge M, Bilenky M, Zhao Y, Zeng T, et al. Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nat Methods. 2007;4:651–7.View ArticlePubMedGoogle Scholar
- Robertson AG, Bilenky M, Tam A, Zhao Y, Zeng T, Thiessen N, et al. Genome-wide relationship between histone H3 lysine 4 mono- and tri-methylation and transcription factor binding. Genome Res. 2008;18:1906–17.View ArticlePubMedPubMed CentralGoogle Scholar
- Au-Yeung N, Mandhana R, Horvath CM. Transcriptional regulation by STAT1 and STAT2 in the interferon JAK-STAT pathway. JAK-STAT. 2013;2:e23931.View ArticlePubMedPubMed CentralGoogle Scholar
- Ni Z, Abou El Hassan M, Xu Z, Yu T, Bremner R. The chromatin-remodeling enzyme BRG1 coordinates CIITA induction through many interdependent distal enhancers. Nat Immunol. 2008;9:785–93.View ArticlePubMedGoogle Scholar
- Abou El Hassan M, Bremner R. A rapid simple approach to quantify chromosome conformation capture. Nucleic Acids Res. 2009;37:e35.View ArticlePubMedPubMed CentralGoogle Scholar
- Harismendy O, Notani D, Song X, Rahim NG, Tanasa B, Heintzman N, et al. 9p21 DNA variants associated with coronary artery disease impair interferon-γ signalling response. Nature. 2011;470:264–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Ning S, Huye LE, Pagano JS. Regulation of the transcriptional activity of the IRF7 promoter by a pathway independent of interferon signaling. J Biol Chem. 2005;280:12262–70.View ArticlePubMedGoogle Scholar
- Schroder K, Hertzog PJ, Ravasi T, Hume DA. Interferon-gamma: an overview of signals, mechanisms and functions. J Leukoc Biol. 2004;75:163–89.View ArticlePubMedGoogle Scholar
- Pattenden SG, Klose R, Karaskov E, Bremner R. Interferon-gamma-induced chromatin remodeling at the CIITA locus is BRG1 dependent. EMBO J. 2002;21:1978–86.View ArticlePubMedPubMed CentralGoogle Scholar
- Abou El Hassan M, Yu T, Song L, Bremner R. Polycomb repressive complex 2 confers BRG1 dependency on the CIITA locus. J Immunol. 2015;194:5007.View ArticlePubMedGoogle Scholar
- Morris AC, Beresford GW, Mooney MR, Boss JM. Kinetics of a gamma interferon response: expression and assembly of CIITA promoter IV and inhibition by methylation. Mol Cell Biol. 2002;22:4781–91.View ArticlePubMedPubMed CentralGoogle Scholar
- Ramsauer K, Farlik M, Zupkovitz G, Seiser C, Kröger A, Hauser H, et al. Distinct modes of action applied by transcription factors STAT1 and IRF1 to initiate transcription of the IFN-gamma-inducible gbp2 gene. Proc Natl Acad Sci USA. 2007;104:2849–54.View ArticlePubMedPubMed CentralGoogle Scholar
- Shi L, Perin JC, Leipzig J, Zhang Z, Sullivan KE. Genome-wide analysis of interferon regulatory factor I binding in primary human monocytes. Gene. 2011;487:21–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Frontini M, Vijayakumar M, Garvin A, Clarke N. A ChIP-chip approach reveals a novel role for transcription factor IRF1 in the DNA damage response. Nucleic Acids Res. 2009;37:1073–85.View ArticlePubMedPubMed CentralGoogle Scholar
- Rettino A, Clarke NM. Genome-wide identification of IRF1 binding sites reveals extensive occupancy at cell death associated genes. J Carcinog Mutagen. 2013;(Spec Iss Apoptosis):S6–009.
- Xie D, Boyle AP, Wu L, Zhai J, Kawli T, Snyder M. Dynamic trans-acting factor colocalization in human cells. Cell. 2013;155:713–24.View ArticlePubMedPubMed CentralGoogle Scholar
- Kumatori A, Yang D, Suzuki S, Nakamura M. Cooperation of STAT-1 and IRF-1 in interferon-gamma-induced transcription of the gp91(phox) gene. J Biol Chem. 2002;277:9103–11.View ArticlePubMedGoogle Scholar
- Chatterjee-Kishore M, Wright KL, Ting JP, Stark GR. How Stat1 mediates constitutive gene expression: a complex of unphosphorylated Stat1 and IRF1 supports transcription of the LMP2 gene. EMBO J. 2000;19:4111–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 2009;459:108–12.View ArticlePubMedPubMed CentralGoogle Scholar
- Sung PS, Cheon H, Cho CH, Hong S-H, Park DY, Seo H-I, et al. Roles of unphosphorylated ISGF3 in HCV infection and interferon responsiveness. Proc Natl Acad Sci USA. 2015;112:10443–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Cheon H, Stark GR. Unphosphorylated STAT1 prolongs the expression of interferon-induced immune regulatory genes. Proc Natl Acad Sci USA. 2009;106:9373–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Reich NC, Liu L. Tracking STAT nuclear traffic. Nat Rev Immunol. 2006;6:602–12.View ArticlePubMedGoogle Scholar
- Lödige I, Marg A, Wiesner B, Malecová B, Oelgeschläger T, Vinkemeier U. Nuclear export determines the cytokine sensitivity of STAT transcription factors. J Biol Chem. 2005;280:43087–99.View ArticlePubMedGoogle Scholar
- Vinkemeier U. Getting the message across, STAT! Design principles of a molecular signaling circuit. J Cell Biol. 2004;167:197–201.View ArticlePubMedPubMed CentralGoogle Scholar
- Meyer T, Marg A, Lemke P, Wiesner B, Vinkemeier U. DNA binding controls inactivation and nuclear accumulation of the transcription factor Stat1. Genes Dev. 2003;17:1992–2005.View ArticlePubMedPubMed CentralGoogle Scholar
- Majoros A, Platanitis E, Szappanos D, Cheon H, Vogl C, Shukla P, et al. Response to interferons and antibacterial innate immunity in the absence of tyrosine-phosphorylated STAT1. EMBO Rep. 2016;17:367–82.View ArticlePubMedPubMed CentralGoogle Scholar
- White LC, Wright KL, Felix NJ, Ruffner H, Reis LF, Pine R, et al. Regulation of LMP2 and TAP1 genes by IRF-1 explains the paucity of CD8+ T cells in IRF-1−/− mice. Immunity. 1996;5:365–76.View ArticlePubMedGoogle Scholar
- Dimitriou ID, Clemenza L, Scotter AJ, Chen G, Guerra FM, Rottapel R. Putting out the fire: coordinated suppression of the innate and adaptive immune systems by SOCS1 and SOCS3 proteins. Immunol Rev. 2008;224:265–83.View ArticlePubMedGoogle Scholar
- Wu J, Ma C, Wang H, Wu S, Xue G, Shi X, et al. A MyD88-JAK1-STAT1 complex directly induces SOCS-1 expression in macrophages infected with Group A Streptococcus. Cell Mol Immunol. 2015;12:373–83.View ArticlePubMedGoogle Scholar
- Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007;39:311–8.View ArticlePubMedGoogle Scholar
- Wright KL, Ting JP. Epigenetic regulation of MHC-II and CIITA genes. Trends Immunol. 2006;27:405–12.View ArticlePubMedGoogle Scholar
- Ni Z, Karaskov E, Yu T, Callaghan SM, Der S, Park DS, et al. Apical role for BRG1 in cytokine-induced promoter assembly. Proc Natl Acad Sci USA. 2005;102:14611–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Chou SD, Khan AN, Magner WJ, Tomasi TB. Histone acetylation regulates the cell type specific CIITA promoters, MHC class II expression and antigen presentation in tumor cells. Int Immunol. 2005;17:1483–94.View ArticlePubMedGoogle Scholar
- Collinge M, Pardi R, Bender JR. Class II transactivator-independent endothelial cell MHC class II gene activation induced by lymphocyte adhesion. J Immunol. 1998;161:1589–93.PubMedGoogle Scholar
- Zhou H, Su HS, Zhang X. Douhan J rd, Glimcher LH. CIITA-dependent and -independent class II MHC expression revealed by a dominant negative mutant [In Process Citation]. J Immunol. 1997;158:4741–9.PubMedGoogle Scholar
- Magner WJ, Kazim AL, Stewart C, Romano MA, Catalano G, Grande C, et al. Activation of MHC class I, II, and CD40 gene expression by histone deacetylase inhibitors. J Immunol. 1950;2000(165):7017–24.Google Scholar
- Guerra SG, Vyse TJ. Cunninghame Graham DS. The genetics of lupus: a functional perspective. Arthritis Res Ther. 2012;14:211.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim H-J, Eom C-Y, Kwon J, Joo J, Lee S, Nah S-S, et al. Roles of interferon-gamma and its target genes in schizophrenia: proteomics-based reverse genetics from mouse to human. Proteomics. 2012;12:1815–29.View ArticlePubMedGoogle Scholar
- Kantarci OH, Hebrink DD, Schaefer-Klein J, Sun Y, Achenbach S, Atkinson EJ, et al. Interferon gamma allelic variants: sex-biased multiple sclerosis susceptibility and gene expression. Arch Neurol. 2008;65:349–57.PubMedGoogle Scholar
- Cooke GS, Campbell SJ, Sillah J, Gustafson P, Bah B, Sirugo G, et al. Polymorphism within the interferon-gamma/receptor complex is associated with pulmonary tuberculosis. Am J Respir Crit Care Med. 2006;174:339–43.View ArticlePubMedGoogle Scholar
- Thye T, Burchard GD, Nilius M, Müller-Myhsok B, Horstmann RD. Genomewide linkage analysis identifies polymorphism in the human interferon-gamma receptor affecting Helicobacter pylori infection. Am J Hum Genet. 2003;72:448–53.View ArticlePubMedPubMed CentralGoogle Scholar
- Grant SFA, Hakonarson H. Microarray technology and applications in the arena of genome-wide association. Clin Chem. 2008;54:1116–24.View ArticlePubMedGoogle Scholar
- Dogan N, Wu W, Morrissey CS, Chen K-B, Stonestrom A, Long M, et al. Occupancy by key transcription factors is a more accurate predictor of enhancer activity than histone modifications or chromatin accessibility. Epigenetics Chromatin. 2015;8:16.View ArticlePubMedPubMed CentralGoogle Scholar
- Jin F, Li Y, Dixon JR, Selvaraj S, Ye Z, Lee AY, et al. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature. 2013;503:290–4.PubMedPubMed CentralGoogle Scholar
- Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M, et al. An atlas of active enhancers across human cell types and tissues. Nature. 2014;507:455–61.View ArticlePubMedPubMed CentralGoogle Scholar
- Ellinghaus D, Ellinghaus E, Nair RP, Stuart PE, Esko T, Metspalu A, et al. Combined analysis of genome-wide association studies for Crohn disease and psoriasis identifies seven shared susceptibility loci. Am J Hum Genet. 2012;90:636–47.View ArticlePubMedPubMed CentralGoogle Scholar
- Nicolae DL, Gamazon E, Zhang W, Duan S, Dolan ME, Cox NJ. Trait-associated SNPs are more likely to be eQTLs: annotation to enhance discovery from GWAS. PLoS Genet. 2010;6:e1000888.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhong H, Beaulaurier J, Lum PY, Molony C, Yang X, Macneil DJ, et al. Liver and adipose expression associated SNPs are enriched for association to type 2 diabetes. PLoS Genet. 2010;6:e1000932.View ArticlePubMedPubMed CentralGoogle Scholar
- Franke A, McGovern DPB, Barrett JC, Wang K, Radford-Smith GL, Ahmad T, et al. Genome-wide meta-analysis increases to 71 the number of confirmed Crohn’s disease susceptibility loci. Nat Genet. 2010;42:1118–25.View ArticlePubMedPubMed CentralGoogle Scholar
- International Genetics of Ankylosing Spondylitis Consortium (IGAS), Cortes A, Hadler J, Pointon JP, Robinson PC, Karaderi T, et al. Identification of multiple risk variants for ankylosing spondylitis through high-density genotyping of immune-related loci. Nat Genet. 2013;45:730–8.View ArticleGoogle Scholar
- Ricaño-Ponce I, Wijmenga C. Mapping of immune-mediated disease genes. Annu Rev Genomics Hum Genet. 2013;14:325–53.View ArticlePubMedGoogle Scholar
- Maurano MT, Humbert R, Rynes E, Thurman RE, Haugen E, Wang H, et al. Systematic localization of common disease-associated variation in regulatory DNA. Science. 2012;337:1190–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Schaub MA, Boyle AP, Kundaje A, Batzoglou S, Snyder M. Linking disease associations with regulatory information in the human genome. Genome Res. 2012;22:1748–59.View ArticlePubMedPubMed CentralGoogle Scholar
- Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinform Oxf Engl. 2003;19:185–93.View ArticleGoogle Scholar
- Hollander M, Wolfe DA. Nonparametric statistical methods, solutions manual. 2nd ed. New York: Wiley. http://www.wiley.com/WileyCDA/WileyTitle/productCd-047132986X.html. Accessed 6 Dec 2016.
- Chen D, Pacal M, Wenzel P, Knoepfler PS, Leone G, Bremner R. Division and apoptosis of E2f-deficient retinal progenitors. Nature. 2009;462:925–9.View ArticlePubMedPubMed CentralGoogle Scholar