High Dose IFN-β Activates GAF to Enhance Expression of ISGF3 Target Genes in MLE12 Epithelial Cells

Interferon β (IFN-β) signaling activates the transcription factor complex ISGF3 to induce gene expression programs critical for antiviral defense and host immune responses. It has also been observed that IFN-β activates a second transcription factor complex, γ-activated factor (GAF), but the significance of this coordinated activation is unclear. We report that in murine lung epithelial cells (MLE12) high doses of IFN-β indeed activate both ISGF3 and GAF, which bind to distinct genomic locations defined by their respective DNA sequence motifs. In contrast, low doses of IFN-β preferentially activate ISGF3 but not GAF. Surprisingly, in MLE12 cells GAF binding does not induce nearby gene expression even when strongly bound to the promoter. Yet expression of interferon stimulated genes is enhanced when GAF and ISGF3 are both active compared to ISGF3 alone. We propose that GAF may function as a dose-sensitive amplifier of ISG expression to enhance antiviral immunity and establish pro-inflammatory states.


INTRODUCTION
The type I interferon (IFN) response is essential to host antiviral immunity. Respiratory epithelial cells, particularly alveolar type II cells, are the primary target of respiratory pathogens and are capable of both producing and responding to IFNs (1)(2)(3)(4)(5)(6)(7). Type I IFN signaling induces expression of IFN stimulated genes (ISGs), which code for antiviral effector molecules as well as chemokines, cytokines, and other immune effector molecules (5,8). Consequently, defects in type I IFN signaling result in susceptibility to various infections, while overexuberant IFN responses are associated with heightened and pathological inflammation (9)(10)(11)(12)(13)(14). Therefore, the regulation of type I IFN signaling is critically important, as its mis-regulation in either direction can lead to disease.
Type I IFN-b signal through the type I IFN receptor (IFNAR) composed of IFNAR1 and IFNAR2 subunits (15). IFNAR is associated with the Janus tyrosine kinases JAK1 and TYK2. Upon ligand binding, the IFNAR bound JAK1/TYK2 kinases cross-phosphorylate and activate signal transducer and activator of transcription (STAT) proteins STAT1 and STAT2. Phosphorylated STAT1 and STAT2 interact with IRF9 to form the trimeric transcription factor (TF) ISGF3, which binds to interferon-stimulated response elements (ISREs) defined by the sequence AGTTTCn 2 TTTC (16).
STAT1 may also form a homodimer (17) named g-activated factor (GAF), which is typically induced by stimulation with type II IFN, IFN-g (18). GAF binds to gamma-activated sequences (GAS), defined by the motif TTCn 2-4 GAA (19), and IFN-g has a well-described role in immune cells such as macrophages (20). Notably, however, type I IFN signaling can also activate GAF in both immune and non-immune cells (1,18,21,22). Although this has been observed for decades, the function of type I IFN induced GAF is unclear. It is possible that GAF and ISGF3 have redundant roles in gene regulation, as type I and type II IFNs were reported to induce highly overlapping gene expression programs (23,24). To date, however, few studies have systematically examined the role of type I IFN induced GAF using genome-wide approaches.
Here we report that in a murine lung epithelial cell line IFN-b activates both ISGF3 and GAF but in differing amounts depending on IFN-b dose. We demonstrate that at high-doses of IFN-b, GAF enhances the ISGF3-responsive expression of ISGs genome-wide, but that binding of GAF alone is insufficient for gene expression.

Cell Culture
MLE-12 cells obtained from ATCC (Manassas, VA) were cultured in DMEM (Corning Inc., Corning, NY) supplemented with 10% FBS, 100 IU/ml penicillin, 100 mg/ml streptomycin, 2 mM L-glutamine at 37 0 C at 5% CO 2 . Cells were verified to be free of mycoplasma contamination. Recombinant murine IFN-b (PBL Assay Science, Piscataway, NJ) and IFN-g (BD Biosciences, San Jose, CA) were used for stimulation. For cycloheximide experiments, cells were incubated for 20 minutes with media containing 10 mg/ml cycloheximide prior to stimulation with IFN.

ChIP-Seq Libraries
MLE-12 cells grown to 80% confluence in 15-cm plates were double crosslinked with 2 mM disuccinimidyl glutarate in PBS for 30 minutes and 1% methanol-free formaldehyde in PBS for 10 minutes at room temperature (RT). Crosslinking reactions were quenched with 125 mM of glycine for 5 minutes at RT, followed by washing of the cells twice with cold PBS supplemented with 1 mM EDTA. Crosslinked cells were then scraped in PBS with EDTA, collected by centrifugation, and frozen at -80°C. Cell pellets were thawed in Lysis Buffer 1 (50 mM HEPES-KOH pH 7.6, 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100, and 1x protease inhibitor cocktail (Thermo Fisher Scientific, Waltham, MA)) on ice and incubated at 4°C for 7 min. Cells were then sonicated with a Diagenode Bioruptor 300 sonication system at medium intensity for 15 seconds ON/30 seconds OFF for 3 cycles in 1.5mL TPX microtubes, and centrifuged to isolate nuclei. Nuclei were washed in Lysis Buffer 2 (10 mM Tris-HCl pH 8.0, 200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, and 1x protease inhibitor cocktail) at RT for 10 mins. After centrifugation, nuclei were resuspended in Lysis Buffer 3 (10 mM Tris-HCl pH 8.0, 100 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 0.1% Na Deoxycholate, 0.5% N-lauroylsarcosine, and 1x protease inhibitor cocktail) and sonicated in TPX microtubes at low intensity for 30 sec ON/30 sec OFF for a total of 12 cycles while inverting and pulsespinning after every 4 cycles. After centrifuging the samples at max speed for 10 mins at 4°C, supernatant containing fragmented chromatin was diluted with 6 volumes of Dilution Buffer (10 mM Tris-HCl pH 8.0, 160 mM NaCl, 1 mM EDTA, 0.01% SDS, 1.2% Triton X-100, 1x protease inhibitor cocktail) and incubated with Protein-G Dynabeads (Thermo Fisher) for 2 hours at 4°C for pre-clearing. The chromatin was then incubated with 5 mg of STAT1 antibody (Cell Signaling Technologies, Danvers, MA; CST-9172) for every ½ plate of cells at 4°C overnight.
T h e a n t i b o d y -c h r o m a t i n c o m p l e x e s w e r e t h e n immunoprecipitated by incubating with protein-G Dynabeads for 5 hours at 4°C. Beads were collected and washed twice with the following buffers: Low Salt buffer (50 mM HEPES-KOH pH 7.6, 140 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% Na Deoxycholate, 0.1% SDS), High Salt buffer (50 mM HEPES-KOH pH 7.6, 500 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% Na Deoxycholate, 0.1% SDS), LiCl buffer (20 mM Tris-HCl pH 8.0, 250 mM LiCl, 1 mM EDTA, 0.5% Na Deoxycholate, 0.5% NP-40), and TE buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA). Immunoprecipitated chromatin complexes were treated with RNAse A (Thermo Fisher) at 37°C for one hour, and crosslinks were reversed with 10% SDS and 0.6 mg/ ml Proteinase K (New England Biolabs, Ipswich, MA) overnight with shaking at 65°C. Immunoprecipitated DNA fragments were purified with AMPure XP SPRI beads (Beckman Coulter, Brea, CA) at a 0.95 volume ratio according to manufacturer's instructions.
DNA was quantified with Qubit dsDNA HS (Life Technologies, Carlsbad, CA) using a Qubit 2.0 fluorometer. Libraries were prepared for sequencing using NEBNext Ultra II DNA Library Prep Kit (New England Biolabs) with NEBNext Multiplex Oligos (New England Biolabs). Each library was prepared from 1-5 ng of starting DNA. Input samples were pooled from samples from the same experiment. Final libraries were checked for quality by agarose gel, quantified with Qubit, and multiplexed with a maximum of 24 samples per sequencing reaction. Libraries were sequenced on an Illumina HiSeq 2500 with single end 50bp reads at the UCLA Broad Stem Cell Research Center.

BMDM Data
Raw ChIP-seq and RNA-seq data (21) were downloaded from the NCBI Sequence Read Archive under accession numbers SRP149943 (STAT1, STAT2, and IRF9 ChIP-seq) and SRP149944 (RNA-seq). Data was downloaded as fastq files and processed in the same manner as data generated in our lab.

ChIP-Seq Analysis
The low quality 3'ends of reads were trimmed (cutoff q=30), and remaining adapter sequences were removed using cutadapt (25). Reads were aligned using bowtie2 (26) using --non-deterministic and --very-sensitive parameters. Peaks were called in each sample against input control by MACS2 peak caller (27) with FDR threshold of 0.01 with the following parameters: keep-dup 1 --call-summits --nomodel --extsize 200. Peaks per sample were merged to generate a reference peak list. Each bam file was normalized based on its sequencing depth with deepTools (28). Sequencing-depth normalized reads were counted within every peak by multiBigWigSummary command from deepTools to generate a count table. Inducible peaks were identified using a threshold of two-fold induction from the basal condition in at least two timepoints.
One replicate of MLE12 IFN-b stimulation at 0.5 hour had particularly high signal to noise ratio, resulting in 10-fold higher number of peaks called compared to other samples. We included this sample in the identification of possible STAT1 peaks but excluded it for finding inducible peaks, thereby maximizing the number of peaks considered but stringently defining the inducible peaks. K-means clustering analysis was performed on scaled log 2 normalized ChIP-seq tags. Silhouette and elbow plots were examined to identify the best k of 2 from the data. De novo motif analysis was performed within peak regions using the findMotifsGenome function in the HOMER suite (29). Heatmaps were generated using the pheatmap package from R (30). The April 2019 version of the Ensembl database was used to extract transcription start site information and to annotate peaks relative to genes. Bigwig files were generated with the bamCoverage function in deeptools, and tracks were visualized in IGV (31).

ChIP-Seq Motif Clustering
Categorization of STAT1 peaks based on motifs was done by searching the HOMER database (29) for position weighted matrix files of STAT/GAS and ISRE/IRF motifs. ISRE/IRF motifs included motifs for T1ISRE, ISRE, IRF2, IRF1, IRF3, and IRF8. STAT/GAS motifs included motifs for STAT1, STAT3+IL-21, STAT3, STAT4, and STAT5 ( Figure S1B). A variant GAS motif in the promoter of Irf1 (GATTTCCCCGAATG) known to be bound by STAT1 (32) was also included in the search as a GAS motif. The annotatePeaks function in HOMER was used to scan for the presence of motifs within inducible STAT1 peaks. A peak that contained at least one kind of STAT/GAS motif was classified as a GAS peak, and a peak that contained at least one kind of ISRE/IRF motif was classified as an ISRE peak. A peak that contained motifs from both STAT/GAS and ISRE/IRF categories was classified as a BOTH peak, while a peak that did not contain either motif was classified as a No motif peak.

RNA-Seq Analysis
Reads were trimmed, filtered, and deduplicated in the same manner as ChIP-seq data. Processed reads were aligned to mm10 genome using STAR (33), and count tables were generated by the featureCount function in deepTools (34). Counts were normalized using the TMM-normalization method and RPKM values were generated using edgeR (35). Genes below an expression threshold of 2 CPM in all samples were excluded from downstream analysis. IFN-b inducible genes ( Figure 3A) were identified by FDR < 0.05 and fold-change > 2 compared to unstimulated. For fold-change and Enhancement score calculations a pseudocount of 1 RPKM was added. Gene ontology analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (36). The data were visualized using pheatmap and ggplot2 packages (30,37).

Electrophoretic Mobility Shift Assay (EMSA)
After stimulation with IFNs, cells were scraped in cold PBS + 1 mM EDTA, collected by centrifugation, and resuspended in Cytoplasmic Extract Buffer (10 mM HEPES-KOH pH 7.9, 60 mM KCl, 1 mM EDTA, 0.5% NP40, 1 mM DTT, and 1 mM PMSF). After incubating on ice for two minutes, cells were vortexed for 30 seconds, and nuclei were collected by centrifugation and resuspended in Nuclear Extract Buffer (250 mM Tris pH 7.5, 60 mM KCl, 1 mM EDTA, 1 mM DTT, and 1 mM PMSF). The nuclei were snap frozen and stored in -80°C before use.
Nuclear extracts were prepared by disrupting the nuclear membrane with three freeze-thaw cycles of 30 seconds each using 37°C water bath and dry ice, followed by centrifugation at max speed for 15 minutes. The resulting supernatants containing nuclear proteins were then incubated with P 32 -labeled oligonucleotide probes for 15 minutes in binding buffer (10 mM Tris pH 7.5, 50 mM NaCl, 10% glycerol, 1% NP40, 1 mM EDTA, 0.1 mg/mL polydeoxyinosinic-deoxycytidylic acid). The following probes were used for detecting ISGF3 and GAF DNAbinding activity: ISRE probe (GATCCTCGGGAAAGG GAAACCTAAACTGAAGCC) and GAS probe (TACAAC AGCCTGATTTCCCCGAAATGACGC). The reaction mixtures were run on a 5% acrylamide (30:0.8) gel with 5% glycerol and TGE buffer (24.8 mM Tris, 190 mM glycine, 1 mM EDTA) at 200V for 1 hour and 45 mins. The gels were dried and imaged on a Sapphire Biomolecular imager in phosphor mode (Azure Biosystems, Dublin, CA).

Western Blot
Cells were lysed with 1x Laemmli buffer (62.5 mM Tris pH 6.8, 2% SDS, 10% glycerol, 5% b-mercapto-ethanol, and bromophenol blue) and boiled at 95°C for 15 minutes. Denatured proteins were electrophoresed on an acrylamide gel and transferred to nitrocellulose. The following antibodies were used: rabbit anti-IRF1 (sc640), goat anti-actin (sc1615), donkey anti-goat HRP (sc2020) from Santa Cruz Biotechnology, and mouse anti-rabbit HRP (CST 7074) from Cell Signaling T e c h n o l o g y . I m m u n o b l o t s w e r e d e v e l o p e d u s i n g chemiluminescent substrate (Supersignal West Pico Plus, Thermo Fisher) and visualized using a ChemiDoc MP imaging system (Bio-Rad, Hercules, CA).

IFN-b Induced STAT1 Binds to GAS and ISRE Motifs
To characterize the TFs activated by IFN-b in MLE-12 cells, we performed electrophoretic mobility shift assays (EMSAs) with ISRE and GAS probes. Robust protein binding to both ISRE and GAS probes was detected in response to IFN-b (10 U/ml), indicating inducible activity of ISGF3 and GAF, respectively ( Figure 1A). The GAS signal reached a maximum within one hour post stimulation, while the ISRE signal gradually increased and peaked at two hours.
To examine ISGF3 and GAF binding events genome-wide, we performed STAT1 ChIP-seq at various time points after IFN-b stimulation. A total of 723 inducible STAT1 peaks were identified. K-means clustering revealed two clusters of peaks with distinctive kinetics ( Figure 1B). Cluster 1 contained 353 peaks that were induced most strongly at half hour and continued until two hours after stimulation. Cluster 2 contained 370 peaks that were most strongly induced at two and four hours after stimulation.
As STAT1 is a member of both ISGF3 and GAF, we performed de novo motif analysis on the genomic regions in each cluster to infer the identities of STAT1-containing complexes. The peaks in Cluster 1 were highly enriched for a GAS-like motif (p = 10 e-255), while the peaks in Cluster 2 were highly enriched for an ISRE-like motif (p = 10 e-350) ( Figure  1C). To corroborate the result of our motif analysis, we performed STAT1 ChIP-seq in response to IFN-g stimulation, which strongly activates GAF. As expected, peaks in Cluster 1 showed greater IFN-g-inducible STAT1 binding than peaks in Cluster 2 ( Figure 1D). These results confirmed that peaks STAT1 in Cluster 1 represent GAF binding and peaks in Cluster 2 represent ISGF3 binding.
We also performed STAT1 ChIP-seq with and without cycloheximide (CHX) to investigate whether GAF binding is dependent on IFN-b-induced de novo protein synthesis of secondary signals such as IFN-g. We found that the addition of CHX did not reduce the ChIP signal in either Cluster 1 or Cluster 2 ( Figure S1A), indicating that IFN-b signaling directly activates GAF without the need for protein synthesis.
To complement the unbiased clustering approach, we categorized peaks based on known GAS and ISRE binding motifs. We comprehensively scanned the HOMER database (29) for STAT/GAS and ISRE/IRF motifs under the 723 peaks ( Figure S1B). We found that for Cluster 1, 262 out of 353 peaks (74%) had only STAT/GAS motifs, 35 peaks (10%) had both STAT/GAS and ISRE/IRF motifs (BOTH motifs), and 8 peaks (2%) had only ISRE/IRF motifs ( Figure 1E). For Cluster 2, 193 out of 370 peaks (52%) had only ISRE/IRF motifs, 120 peaks (32%) had BOTH motifs, and 18 peaks (5%) had only STAT/ GAS motifs ( Figure 1E). In sum, this motif classification analysis resulted in 280 peaks that had only STAT/GAS motifs (categorized as GAS peaks), 155 peaks that had both STAT/ GAS and ISRE/IRF motifs (BOTH peaks), and 201 peaks that had only ISRE/IRF motifs (ISRE peaks). A re-clustered heatmap using these categories ( Figure 1F) recapitulated the distinct temporal patterns of STAT1 binding seen by K-means clustering. As examples, peaks from BOTH, GAS, and ISRE motif categories were found in the promoters of Isg15, Ripk1, and Usp18 respectively ( Figure 1G). Genome-wide, the peaks with GAS motifs were on average induced more weakly than peaks with BOTH or ISRE motifs ( Figure S1C). Eighty-seven peaks had neither canonical GAS nor ISRE motifs ( Figure 1E), which may represent STAT1 recruited to chromatin through proteinprotein interactions without making direct DNA contacts. These peaks were removed from subsequent analyses.

Low-Dose IFN-b Activates ISGF3 but Not GAF
Having defined that IFN-b activates GAF genome-wide, we wondered if IFN-b dose affected whether GAF is activated. We stimulated cells with with low (1 U/ml) and high (10 U/ml) doses of IFN-b for one hour and performed STAT1 ChIP-seq. We observed that at the low IFN-b dose, ISRE-containing regions showed inducible STAT1 binding (p < 2.2 e-16), yet no significant binding was found at GAS-containing regions (Figures 2A, B) as defined by the motif categorization in Figure 1. For example, the GAS-containing peak at the Ripk1 promoter is bound in a dose-sensitive manner ( Figure 2C), while peaks at the Isg15 and Usp18 promoters that possess BOTH and ISRE motifs, respectively, are induced to a similar magnitude by low and high dose IFN-b. This demonstrated that low-dose IFNb induces STAT1 binding to ISRE but not GAS sites, suggesting that GAF may act as a molecular switch to distinguish low from high doses of IFN-b.

IFN-b-Induced GAF Does Not Induce Expression of Nearby Genes
We next interrogated the function of ISGF3 and GAF binding events by linking the 636 inducible STAT1 ChIP-seq peaks identified in Figure 1 (belonging to BOTH, GAS, and ISRE clusters) to their closest expressed genes ( Figure 3A). This yielded 543 uniquely linked genes, of which 181 had promoter STAT1 peaks (-1000bp to +100bp) and 362 had distal STAT1 peaks. Genes linked to peaks in the BOTH and ISRE clusters showed functional enrichment for defense response to virus and cellular response to IFN-b gene ontology terms, confirming that these are known ISGs. Interestingly, however, the 238 genes linked by proximity to GAS peaks were not enriched for any ontology terms ( Figure S2A).
We then examined the expression of the 543 linked genes by performing RNA-seq over a four-hour time course after stimulation with IFN-b (10 U/ml). Clustering the genes by the motif categories of their linked STAT1 ChIP-seq peaks ( Figure 3B), we found that genes linked to the BOTH and ISRE clusters were induced upon IFN-b stimulation. The magnitude of induction was strongest for genes with promoter STAT1 peaks, while distally linked genes were less likely to be induced. This may be due to the fact that enhancers do not necessarily regulate their closest genes, and the linkage between peak and gene becomes increasingly unreliable with increased genomic distance (38,39).
In striking contrast to the BOTH and ISRE clusters, we found very little induction among genes linked to STAT1 peaks with only GAS motifs, even when the binding events were in their promoters ( Figure 3B, GAS cluster). Only 4% of genes (10 out of 238) in the GAS cluster were induced more than two-fold by IFN-b, compared to 44% and 49% of genes in the BOTH and ISRE clusters, respectively ( Figure 3C). For example, Fastk, Cuta, and Psmd3 all have prominent IFN-b-inducible STAT1 peaks in their promoters that contain GAS motifs ( Figure 3D). However, these genes displayed no IFN-b-inducible change in expression level ( Figure 3E).
Genes linked to GAS peaks were also not inducible by higher doses of IFN-b and IFN-g. Even at 100 U/ml IFN-b and 100 ng/ ml IFN-g, only eight and five out of 238 genes, respectively, were induced greater than two-fold ( Figures S2B, C). One of the genes consistently induced by IFN stimulation was Irf1, which has been implicated as a key factor in amplifying IFN responses (22,40).
We observed that the promoter of Irf1 had a STAT1 peak with a GAS motif ( Figure S3A), and its expression was induced by 10 U/ml IFN-b ( Figures S3B, C). Irf1 may thus be an important exception to the general observation that GAF binding does not induce expression of nearby genes.

In Macrophages, GAF Does Induce Expression of Nearby Genes
Given that IFN-g has a well-described role in macrophage function (20), we wondered if GAF was similarly inactive in a macrophage context. Using publicly available datasets of STAT1, STAT2, and IRF9 ChIP-seq and RNA-seq in bone marrowderived macrophages (BMDMs) (21), we performed an identical analysis to what was done in MLE12 cells. A total of 890 STAT1 peaks were inducible by 90-minute stimulations with either IFNb (200 U/ml) or IFN-g (10 ng/ml). These were categorized into 469 peaks containing GAS motifs only, 38 peaks containing ISRE motifs only, and 335 peaks containing BOTH motifs ( Figure  4A). The low number of ISRE-containing STAT1 peaks was due partly to the relatively weak STAT1 signal and also consistent with prior observations that many ISGF3 binding events in this system are STAT1-independent (21). To overcome this limitation, we used an alternative classification method and categorized BMDM ChIP-seq peaks by patterns of TF binding: locations with STAT1 peaks and no IRF9 peaks were considered GAF binding events, and locations with IRF9 peaks in WT cells and no STAT1 or STAT2 peaks in Irf9 -/cells were considered ISGF3 binding events. Using this approach, we identified 1281 ISGF3 and 677 GAF binding events ( Figures S4A, B). We then linked peaks to their nearest expressed genes to examine IFN-responsive gene expression for both categorization methods. Using the motif-based categorization of STAT1 peaks, we identified 777 uniquely linked genes ( Figure 4B). In contrast to our observation in MLE12 cells, here we found that genes linked to GAS peaks were induced in a similar manner as genes linked to ISRE or BOTH peaks ( Figure 4C), with similar patterns of induction by IFN-b or IFN-g stimulation. We found that 61.8% of genes linked to GAS peaks, 66.7% of genes linked to ISRE peaks, and 70.9% of genes linked to BOTH peaks were induced greater than two-fold ( Figure 4D). The alternative method of classifying peaks by TF binding reproduced this result, with 59.1% of genes linked to GAF binding and 61.6% of genes linked to ISGF3 induced greater than two-fold ( Figure S4C).
These findings indicated that GAF is transcriptionally active in BMDMs but not in MLE12 cells ( Figure 4E), likely due to the presence of collaborating TFs that are expressed in macrophages (41,42). Indeed, promoters of genes linked to GAS peaks were enriched 1.9-fold (p = 10e-7) for motifs of PU.1, a macrophagespecific transcription factor that collaborates with STAT1 in gene expression (42,43).

Combined Activation of GAF and ISGF3 Enhances Expression of ISGs in Respiratory Epithelial Cells
Given the surprising finding in epithelial cells that IFN-b activates substantial GAF binding with no direct transcriptional activity, we hypothesized that GAF might play an indirect role in the regulation of ISGs. To address this hypothesis, we first defined ISGs using our RNA-seq data from MEL12 cells stimulated with high dose IFN-b (10 U/ml). We identified 179 ISGs ( Figure 5A) by stringent statistical thresholds (FDR < 0.05 and fold-change > 2) and used this set of genes for subsequent analysis. To study the effects of GAF on ISG expression, we leveraged the fact that low-dose IFN-b activates ISGF3 but not GAF (Figure 2), while low-dose IFN-g activates GAF but not ISGF3 ( Figure 5B). We stimulated cells with lowdose IFN-b, low-dose IFN-g, a mixture of the two, and high-dose IFN-b, resulting in the TF activation patterns show in Figure 5C. Comparing gene expression in the Mixed condition and the single-stimulus conditions thus allowed for discovery of genes that are co-regulated by ISGF3 and GAF. We performed RNA-seq in response to IFN-b, IFN-g, Mixed, and high dose IFN-b conditions at four hours after stimulation and examined the expression of the 179 previously-defined ISGs. Principal component analysis showed that two replicates of the Mixed condition fell between IFN-b and IFN-g samples and were placed further away from the unstimulated sample compared to the single stimulus conditions ( Figure 5D). This suggested that the Mixed condition resulted in a gene expression profile that is a hybrid of IFN-b and IFN-g profiles and that some ISGs were more strongly induced in the Mixed condition.
Visualizing these ISGs in a hierarchically clustered heatmap provided additional insight into the effects of GAF and ISGF3 ( Figure 5E). We observed that ISGs in the top cluster (C1) were induced by IFN-g alone, and their expression was only mildly enhanced in the Mixed condition. Similarly, genes in the bottom cluster (C3) were induced by IFN-b alone and mildly enhanced in the Mixed condition.
The expression of ISGs in the middle cluster (C2) was of greatest interest. These ISGs were induced weakly in response to either IFN-b or IFN-g alone, but strongly enhanced in the Mixed condition. To quantify the effect of mixing IFN-b and IFN-g, we calculated an Enhancement score for each gene: Enhancement score = RPKM mixed /(RPKM b + RPKM g) (Figure 5E annotation). This objective score confirmed that many genes in the middle cluster are strongly enhanced in the Mixed condition, suggesting that ISGF3 and GAF may co-regulate the expression of some ISGs.
To investigate the properties of highly enhanced genes, we compared ISGs in the top quartile of Enhancement score (n = 45, Enhancement > 1.156) to ISGs in the bottom quartile (n = 45, Enhancement < 0.774)) ( Figure 6A). To confirm that enhanced gene expression occurred over multiple time points, we performed RNA-seq over a six-hour time course. Genes in the Top quartile indeed displayed a pattern of enhancement across multiple time points, and the degree of inducibility in the Mixed condition approached that seen with high dose IFN-b ( Figure  6B). We also observed that genes in the Top category were induced more robustly than genes in the Bottom category for all conditions (Figures 6B, C). This was primarily due to the fact that genes in the Top category had much lower basal expression than genes in the Bottom category ( Figure 6C, p = 2.66 e-15 for difference between basal RPKMs).
We confirmed that genes in the Top category were also substantially more inducible in response to high doses of IFNb (10 and 100U/mL, Figure 6D), suggesting that the Enhancement score in our mixing experiment is predictive of gene expression responses at high doses of IFN-b that activate both GAF and ISGF3.
Prior studies have shown that GC content is anti-correlated with nucleosomal stability, and that genes with high GC content in their promoters are constitutively expressed because they do not require nucleosome remodeling for transcriptional activation (44,45). We calculated GC content in promoters (-300bp to +300bp) of Top vs. Bottom Enhancement score genes and found that promoters of the Bottom category genes showed significantly higher GC content compared to promoters of the Top category genes (p = 0.0001, Figure 6E). We also found that 25 out of 45 (56%) genes in the Top enhancement quartile had STAT1 ChIP-seq peaks in their promoters (-1000bp to +100bp) at any time point, compared to only 10 out of 45 (22%) genes in the Bottom quartile (p = 0.002, Figure 6F). These findings support the idea that genes most strongly enhanced by the combined activation of ISGF3 and GAF are not constitutively expressed but require inducible TF binding at their promoters to establish a chromatin environment that permits maximal expression. Genes in the Top quartile of Enhancement score included the pro-inflammatory cytokine Cxcl10 (Enhancement score = 2.25) and the antiviral effector Ifit1 (Enhancement score = 2.24), both of which demonstrate low basal expression and greatest inducibility in the Mixed condition ( Figure 6G). In contrast, Tor1aip1 was representative of genes in the Bottom quartile, with relatively high expression at the basal state (19 RPKM) and modest inducibility that was similar in IFN-b alone and the Mixed condition. A complete list of Top and Bottom quartile genes is provided in Table S1. TF motif analysis and gene ontology analysis did not reveal any difference between the two categories. Both Top and Bottom genes contained ISRE motifs in their promoters and were enriched for IFN-related ontology terms such as "Viral defense" and "Immune response" (data not shown).

Gene-Specific Mechanisms of GAF-Mediated ISG Enhancement
If GAF binding does not directly induce expression of nearby genes, yet mixing IFN-b and IFN-g enhances expression of ISGs, we inferred that there must be indirect mechanisms by which GAF cooperates with ISGF3. One possible mechanism is that GAF and ISGF3 may work cooperatively at the promoters of enhanced genes ( Figure 7A), where GAF may either recruit general TFs or increase chromatin accessibility to facilitate ISGF3 binding. Consistent with this, genes with promoter STAT1 peaks containing BOTH motifs had higher enhancement scores compared to genes with promoter STAT1 peaks containing only ISRE motifs ( Figure 7B), implying that combinatorial binding of GAF and ISGF3 at a BOTH peak may enhance gene expression. For example, Mx2 had two distinct STAT1 peaks in its promoter, one of which had BOTH motifs, and it had an enhancement score of 1.23 ( Figure 7C). A second possible mechanism is that GAF may bind to distal enhancers of ISGs ( Figure 7D). GAS peaks were distributed further away from annotated genes than BOTH or ISRE peaks ( Figure 7E), and the nearest genes to GAS were not enriched for IFN-related ontology terms ( Figure S2A). This suggested that GAF binds predominantly in intergenic regions functioning as distal enhancers rather than promoters. To explore this hypothesis, we examined promoter-capture Hi-C (Pc-HiC) data from mouse embryonic stem cells (46). Genome-wide analysis of this Pc-HiC data did not identify significant differences between Top and Bottom Enhancement score genes, likely due to the fact that enhancers are highly cell type specific (47). However, we found that the promoter of Cxcl10, a gene with high enhancement score, makes statistically significant Hi-C contacts with three neighboring STAT1 peaks, of which two contain BOTH motifs and one contains an ISRE motif ( Figure 7F).

DISCUSSION
Coordinated activation of both ISGF3 and GAF in response to IFN-b stimulation has been previously observed, but the function of GAF in this context has remained unclear. We found that activation of GAF occurs only at high doses of IFNb, whereas low-dose IFN-b activates only ISGF3 but not GAF. Surprisingly, STAT1 binding to GAS sites was insufficient to induce gene expression in MLE12 cells but not BMDMs. Despite the lack of GAF-mediated transcription, we observed that stimulating MLE12 cells with low doses of IFN-b and IFN-g together, which recapitulates the coordinated activation of ISGF3 and GAF by high doses of IFN-b, enhanced ISG expression compared to either stimulus alone. We posit that GAF enhances ISG expression through a variety of gene-specific mechanisms that enhance the activation potential of ISGF3.
Our finding that IFN-b activates GAF in a dose-sensitive manner carries important implications for immunity and inflammation. Low levels of tonic IFN-b maintain homeostasis in multiple tissues and confer a basal degree of antiviral protection (48)(49)(50)(51). Our study suggests that in tonic IFN-b conditions, the expression of ISGs is restricted to low levels in the absence of GAF. However, in the context of an active infection, IFN-b concentration in the microenvironment increases, and IFN-g is also produced. The increased levels of IFN-b and IFN-g activate GAF, which acts as a molecular switch to enhance expression for a subset of ISGs. These enhanced ISGs include antiviral mediators such as Ifit1, as well as proinflammatory chemokines such as Cxcl10 ( Figure 6G). The absence of GAF in tonic IFN-b conditions may thus restrict the expression of pro-inflammatory genes that would be harmful in homeostatic conditions. We propose a model where ISGs fall into two classes: those that are GAF-independent and weakly induced even at high IFN-b doses, and those that are enhanced when GAF is activated ( Figure S5). Thus the IFN-b-GAF axis controls the ouput of a subset of ISGs in a dose-sensitive manner to limit inflammation at homeostasis but promotes both antiviral and pro-inflammatory gene expression in response to pathogens. A striking conclusion of our study was that STAT1 binding to GAS elements is insufficient to activate gene expression in MLE12 cells. Genes linked to GAF binding events were not induced in response to high-doses of either IFN-b or IFN-g, even when the binding events were in their promoters (Figure 3 and Figure S2). In contrast, using an identical analysis in BMDMs we found that GAF binding is associated with nearby gene expression in macrophages. The likely explanation for this cell type-specificity is the presence of collaborating TFs that are more abundant or basally bound in macrophages, such as IRF1, IRF8, and/or PU.1 (42), thus providing a chromatin environment that better facilitates transcription.
Our findings imply that in certain contexts GAF alone is insufficient for transcriptional activation. This has been suggested by others, and there are two plausible reasons why this may be true. First, there are two isoforms of STAT1: full-length STAT1a and a shorter STAT1b that lacks the S727 residue required for maximum transcriptional activity (52)(53)(54). STAT1b interacts poorly with transcriptional co-activators such as CBP/p300 and thus has reduced transcriptional activity compared to STAT1a (55)(56)(57). It has been suggested that STAT1a preferentially forms ISGF3 complexes while the less active STAT1b forms GAF complexes (52). Second, post-translational modifications are required for STAT1 transactivation. Unphosphorylated STAT1 can participate in transcription as part of the ISGF3 complex (56,58) due to the transactivation activity of STAT2. But as a homodimer, S727 phosphorylation is required for STAT1 transcriptional activity (56). Future ChIP experiments using antibodies against putative collaborating TFs or different STAT1 isoforms and phosphorylation states would shed light on these proposed mechanisms.
Despite its inability to directly activate gene expression in epithelial cells, GAF enhances the expression of ISGs in MLE12 cells. GAF may synergistically collaborate with ISGF3 to induce maximal gene expression, either by colocalizing at promoters or when distally-bound GAF interacts with promoter-bound ISGF3 through chromosomal looping. As GAF alone is transcriptionally inactive, we propose that GAF may modulate the chromatin environment to facilitate ISGF3 binding and transactivation. STAT1 has a well defined epigenetic role downstream of IFN-g in macrophages, priming genes for greater response to other TFs such as NFkB by modulating the chromatin environment (20,(59)(60)(61). A similar mechanism, albeit on a shorter time scale, may underlie GAF's enhancement of ISGF3 in response to high IFN-b doses. Intriguingly, our ChIP-seq time course data show that GAF activity peaks earlier than ISGF3, suggesting that GAF may bind first, altering chromatin states and facilitating subsequent ISGF3 binding and robust ISG induction.
We also found that in MLE12 cells IRF1 was one of only eight IFN-b inducible genes with a GAS-containing STAT1 peak in its promoter. In human hepatocyte cell lines and intestinal organoids, IRF1 is induced by IFN-g and high doses of IFN-b (25 U/ml), and it is required for maximal expression of ISGs including Cxcl10 (22). Our results imply a mechanistic paradigm in epithelial cells where IFN-b activates GAF, which then induces IRF1 to amplify expression of a subset of ISGs. Indeed, IRF1 may be one of the only genes that is a direct target of GAF. In this way, the GAF-IRF1 axis may be analogous to the IRF3-ISGF3 axis, where the secondary TF is a more potent activator of gene expression than the primary TF, whose principal role is to induce the secondary TF (16).
Our study demonstrates that IFN-b activates GAF in a dosedependent manner, which collaborates with ISGF3 to enhance expression of ISGs in MLE12 cells. These results advance our understanding of the regulation of type I IFN signaling while raising a number of questions ripe for further investigation. Whether these mechanisms are seen in other epithelial cell types, including primary cells, is unclear. Further investigation of the role of IRF1 or epigenetic changes induced by GAF binding could include genetic perturbations or measurements of chromatin accessibility and histone modifications at GAF binding sites. Finally, the physiological role of GAF in the IFNb response could be explored using in vivo models of infection in which its coordinated activation is dysregulated.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/geo/, GSE161475.

AUTHOR CONTRIBUTIONS
KK, CW, JB, and QC performed the experiments. KK, MN, CO, and QC analyzed the data. KK and QC wrote the manuscript. AH and QC conceived of the study. AH obtained funding for the study. All authors contributed to the article and approved the submitted version.