A stimulus‐contingent positive feedback loop enables IFN‐β dose‐dependent activation of pro‐inflammatory genes

In order to determine whether IFN‐β and IFN‐λ3 induce differential expression of ISGs, we produced replicate RNA‐seq datasets from an extensive time course. The murine MLE‐12 lung epithelial cells were stimulated with either 10 U/ml IFN‐β or 100 ng/ml IFN‐λ3. Concentrations were selected such that the initial activation of the IFN signaling pathway was comparable as measured by ISGF3 activity at 30 min (Fig EV1), thereby sidestepping potential differences in IFNAR and IFNLR receptor expression. RNA was isolated at 11 time points of stimulation (0, 0.5, 1, 2, 3, 4, 6, 8, 12, 24, and 36 h) and subjected to polyA+ mRNA sequencing (RNA‐seq) for global transcriptome analysis. A total of 345 genes were identified as inducible (FDR‐corrected P‐value < 0.05 for IFN‐β dataset; Dataset EV1). Principal component analysis (PCA) of the selected genes yielded two components that accounted for almost 90% of the variance (Fig 1A). While the transcriptomes of IFN‐β and IFN‐λ3 time courses remain similar at early time points (0.5 and 1 h), they diverge on the PCA plot at 2 h and are well‐separated at 4 h. The number of induced genes was similar at the early and late time points, but at 4 h, IFN‐β showed 2.3 times as many induced genes as IFN‐λ3 (FDR‐corrected P‐value < 0.05 and log₂(FC) ≥ 1; Fig 1B). Interestingly, even for genes induced by both stimuli at 4 h, there was a higher magnitude of induction with IFN‐β (Fig 1C). However, by 24 h many genes showed a slightly higher induction with IFN‐λ3.

In order to assess the temporal dynamics at a single gene level, a heatmap was generated with each row representing expression of an individual gene (Fig 1D; Dataset EV1). The selected 345 genes were grouped into two clusters based on induction by only IFN‐β (cluster 1; FDR‐corrected P‐value ≤ 0.05) or by both IFN‐β and IFN‐λ3 (cluster 2; FDR‐corrected P‐value ≤ 0.05). Only a small subset of 27 genes were classified as IFN‐λ3 specific and were not included in further analysis (Table EV3). Cluster 1 had 188 genes that generally showed transient activation with a peak at around 4 h of IFN‐β stimulation. The 157 genes in cluster 2 also showed transient induction in response to IFN‐β but were also induced in response to IFN‐λ3 with gene expression peaking at later time points (12, 24, and 36 h). Gene Ontology (GO) analysis identified potential cluster‐specific biological functions such as “inflammatory and adaptive immune response” as well as “cytokine production” for cluster 1 genes and “regulation of adaptive immune cytotoxicity” as well as “macrophage activation” for cluster 2 genes (Fig 1E).

To inspect the expression time course of specific genes with important cluster‐specific biological functions in more detail, we selected genes based on top quartile loadings of principal component 2 (Fig 1A; Dataset EV1), presence in cluster 1 (Fig 1D) and association with cluster 1‐specific GO functions (i.e., inflammatory and adaptive immune response and cytokine production). Line graphs of these cluster 1 genes confirmed the IFN‐β‐specific induction around the 4‐h time point (Fig 1F). A similar analysis was conducted for the top quartile genes contributing to variance in PC1 that are also found in cluster 2 and contribute to cluster 2‐specific biological functions (i.e., regulation of adaptive immune cytotoxicity and macrophage activation). While many of these genes showed higher expression levels with IFN‐β at 4 h, with IFN‐λ3 they reached almost equivalent expression levels at late time points. Many of these genes (e.g., IRF7, RSAD2, TLR3, DDX58, and IFI35) have been classified as “robust” IFN response genes (Levin et al, 2014).

Since the temporal expression trajectories of cluster 2 genes were distinct when stimulated with either IFN‐β or IFN‐λ3, we hypothesized that IFN‐specific temporal dynamics of the transcription factor ISGF3 were responsible. To examine this question, we measured ISGF3 activity in a time course of MLE‐12 lung epithelial cells incubated with 10 U/ml IFN‐β or 100 ng/ml IFN‐λ3 using an electrophoretic mobility shift assay (EMSA) with an ISRE oligo probe using constitutive transcription factor NFY as control (Fig 2A). We found three phases of ISGF3 activity during the IFN‐β stimulation: an activation phase (0 to 60 min), a potentiated phase (80 min to 4 h), and an inactivation phase (6 to 12 h). IFN‐λ3‐stimulated ISGF3 shared the activation phase, but lacked the potentiated phase, instead showing a trough (40 min to 2 h), and reactivation phase (4 to 12 h). The IFN‐β‐specific potentiated phase is characterized by a more than a twofold increase in activity that peaks at 4 h of stimulation.

A simple model, composed of an ordinary differential equation (ODE), was used to determine whether the temporal dynamics of ISGF3 activity might account for the measured cluster‐specific expression dynamics in response to each stimulus. The ODE model simulates the change in mRNA abundance as a function of mRNA synthesis activated by the ISGF3 transcription factor with an activation constant (K_A), Hill coefficient (n), and an mRNA degradation rate constant (k_deg). The only time‐dependent model variable is the ISGF3 activity, quantified from the EMSA (Fig 2A), which was used as an input. The model simulations were fit to the median mRNA abundances for genes in each cluster at each time point for both stimuli (Fig 2B; Dataset EV1) using a simplex search parameter optimization with 50 random initializations. To ensure robustness, the optimization workflow was conducted five times and the mean and associated standard error of the top 10 parameter sets from the five independent runs is reported. As expected, we were able to identify a parameter set that produced good fits for cluster 2 for both stimuli (RMSD = 0.119 ± 0.0001; Fig 2C). Then, using the time course data of cluster 1, we identified a different parameter set that produced model simulations with a good fit (RMSD = 0.152 ± 0.0002). This suggests that the IFN‐specific temporal dynamics of ISGF3 activity may be sufficient to mediate not only IFN‐specific temporal control of ISG expression (cluster 2) but also activation of an IFN‐β‐specific gene expression program (cluster 1), without the need to invoke other pathways or transcription factors.

The parameter sets for the models for cluster 1 and 2 differed by 1.2× in the activation constant, by 1.5× in the degradation rate constant, and by 2× in the Hill coefficient. Using a sensitivity analysis where one parameter was changed by a multiplier, we found that for the cluster 1 model a fourfold increase or decrease in the activation constant, K_A, diminished the performance of the model only minimally (RMSD = 0.160 and 0.153, respectively; Fig 2D). In contrast, a much smaller range (≥ 0.7× and ≤ 1.4×) for the values of the fitted Hill coefficient, n, was found to be necessary in both parameter sets. The cost function more than doubled (RMSD = 0.342) when decreasing the Hill coefficient in the cluster 1 model to 1.0, and more than tripled (RMSD = 0.380) when increasing it to 2.0 in the cluster 2 model. The mRNA degradation constant showed an intermediate sensitivity: a 0.7‐ and 1.4‐fold change for models 1 and 2, respectively, resulting in a slightly increased RMSD (0.173 and 0.153, respectively). Taken together, this suggests the decoding of ISGF3 temporal dynamics relies on specific values of the Hill coefficient and, to a smaller extent, the mRNA degradation constant.

Given that the cluster 2 model fit poorly to cluster 1 data (RMSD = 0.432), we asked which parameter difference was critically important to distinguish cluster 1 from cluster 2 (Fig 2E). We found that substituting the Hill coefficient alone was sufficient to improve the model performance to more closely simulate the cluster 1 data (RMSD = 0.183). Substituting both the Hill coefficient and degradation rate further improved the fit (RMSD = 0.152). The dose–response relationship of mRNA abundance against ISGF3 activity for the cluster 2 model is close to linear for the experimentally determined range of ISGF3 activity (Fig 2F). For the experimentally measured maximum ISGF3 activity induced by IFN‐λ3 (1.4 A.U., red dotted vertical line), the mRNA abundance is 62% compared with the maximum mRNA abundance seen with the measured IFN‐β‐induced activity (2.5 A.U., blue dotted vertical line). This is in contrast to the cluster 1 model whose nonlinear dose–response results in only 31% mRNA abundance for 1.4 A.U. ISGF3 activity.

The modeling suggested an ultrasensitive dose–response curve for cluster 1 genes versus a linear dose–response curve for cluster 2 genes. Such nonlinear dose–responses may arise from multiple transcription factors binding and synergizing in the activation of the promoter (Johnson et al, 1979; Carey et al, 1990), or from chromatin looping to connect distal enhancers with promoters (Kuhlman et al, 2007; Earnest et al, 2013; Hou et al, 2018). We investigated each possibility by examining STAT1 ChIP‐seq data and identified STAT1 binding motifs within the measured STAT1 peaks. This analysis did not reveal an overrepresentation of multiple motifs in STAT1 peaks associated with cluster 1 genes (Fig 2G). However, many more cluster 1 genes showed larger distances between STAT1 peaks and transcription start sites (TSSs), while for cluster 2 genes STAT1 peaks tended to be promoter‐proximal (Fig 2H). This analysis suggests the possibility that the ultrasensitive dose–response behavior that is driving the type I‐specific expression of cluster 1 genes may be mediated by the need to engage in chromatin looping events that connect distal ISGF3 binding to the TSSs of target genes.

Since the IFN‐type‐specific temporal dynamics of ISGF3 appear to control IFN‐specific gene expression levels, we examined the mechanisms for regulating ISGF3 dynamics. We leveraged an ODE model of the JAK–STAT signaling pathway that recapitulates IFN‐α‐induced signaling in hepatocytes (Kok et al, 2020), which contains 41 molecular species connecting ligand‐receptor interactions with nuclear ISGF3 activity (Table EV1). To adapt the model to IFN‐β‐induced ISGF3 signaling in epithelial cells, we produced experimental data from time course measurements of the pathway regulators STAT1, STAT2, IRF9, SOCS1, SOCS3, and USP18 during IFN‐β stimulation using immunoblot, EMSA, and qPCR assays (Fig 3A). We sought to make the least number of parameter changes and thus implemented an iterative algorithm starting with only three systematically selected parameters to produce an optimal fit (Fig 3B; Table EV2). However, key features were not captured, such as the potentiated phase of some of the nuclear active species (ISGF3 and pSTAT1), the late phases of the nuclear total species (STAT1, STAT2, IRF9), as well as the dynamics of most of the mRNA species (total RMSE = 73.26; Fig 3C). To improve the model fits, three additional parameters were systematically selected for optimization. We found that allowing 12 parameters to be estimated was sufficient to capture the potentiated phase of the nuclear active species and the late phase of the nuclear inactive species (total RMSE = 43.72), but the 21‐parameter estimation further improved the overall fit, particularly the STAT1 mRNA data (total RMSE = 35.18). To avoid overfitting, no larger sets of parameters were tested since the error had reached a plateau, showing a similar error when fitting all the parameters concurrently (total RMSE = 33.38; Fig 3B). The parameter changes (Table EV2) suggest that MLE‐12 lung epithelial cells have altered nucleo‐cytoplasmic transport to increase the nuclear residence time of ISGF3 and lower STAT1 and 2 expression levels than in the original model that was parameterized to hepatocytes. With the new MLE‐12 adapted model (i.e., 21‐free‐param model), we could now explore the circuit mechanisms that are key to IFN‐β‐specific ISGF3 dynamics.

Using the signaling model adapted to IFN‐β‐responding MLE‐12 lung epithelial cells, we explored how varying the dose of IFN‐β stimulation would affect the temporal dynamics, and specifically, the potentiated phase of ISGF3 activity. First, we decreased the dose of IFN‐β by two‐, four‐, ten‐, or twenty‐fold (Fig 4A). While 1× IFN‐β (50,000 molecules/cell) induced an initial activation phase up to 1 h followed by a potentiated phase up to 4 h, our simulations predicted a dose‐dependency in the magnitude of the ISGF3 response. We plotted the ISGF3 activity during the initial activation phase at 1 h and during the potentiated phase at 4 h and found that the potentiated phase is reduced more drastically than the initial activation phase (Fig 4B). For example, the model predicted that a 10‐fold lower IFN‐β dose would reduce the 4‐h activity by about 2.3×, whereas the 1‐h activity remained within 1.2× of the original, resulting in close to equal magnitudes at both time points. Experimentally, we tested three doses (10, 1, and 0.1 U/ml) in an ISGF3 time course experiment (Fig 4C). Similar to the model predictions, ISGF3 activity trajectories were diminished at lower doses. The quantified experimental data showed a much more severe reduction in the 4‐h potentiated phase than the 1‐h activation phase such that both time points had close to equal ISGF3 activity at lower doses (Fig 4D). This suggests that the IFN‐β‐specific potentiated phase of ISGF3 is dose‐dependent and requires a higher dose‐range than the initial activation phase, which occurs also at lower doses.

To determine whether the ISGF3 potentiated phase is dependent not only on the dose but also the duration of stimulation, we compared model simulations of the sustained 1× IFN‐β condition and a pulse stimulation where IFN‐β was set to zero after 15 min (Fig 4E). While the activation phase was similar in both conditions, the potentiated phase did not occur with the shorter stimulation time. Instead, ISGF3 activity rapidly decreased to baseline within 6 h. We tested this prediction experimentally in MLE‐12 lung epithelial cells stimulated with 10 U/ml IFN‐β for 15 min (Fig 4F). Similar to the model predictions, the experimental ISGF3 activity did not increase after the activation phase but diminished to basal activity by 6 h. Together, the modeling and experimental results suggest that the IFN‐β‐specific potentiated phase of ISGF3 requires higher doses and longer durations of ligand exposure.

The above‐described stimulation conditions allowed us to assess the functional consequence of the ISGF3 potentiated phase for gene expression. We isolated RNA from MLE‐12 lung epithelial cells incubated with either low IFN‐β (1 U/ml) or pulse IFN‐β (15 min of 10 U/ml) for various durations and assessed the expression of the previously identified 345 IFN‐β inducible genes by RNA‐seq. Principal component analysis revealed that the response profile for neither condition converged with the response profile for sustained IFN‐β after 1 h, but clustered near the IFN‐λ3‐induced profile until 4 h (Fig 4G, dashed circles). After 4 h, the pulse and low IFN‐β‐induced transcriptomes converge toward the unstimulated condition. When examining the previously identified IFN‐β‐specific gene expression program in cluster 1, the lower dose and shorter duration of IFN‐β was not sufficient to induce these genes (Fig 4H). Line graphs of cluster 1‐specific inflammatory & adaptive immune response and cytokine production genes (e.g., Cd274, Il6, TNFSF10, CCRL2, Tlr2, Myd88, and IRF5) confirmed that only the sustained and high dose of IFN‐β was able to induce them (Fig 4I). This demonstrates that the ISGF3 potentiated phase induced with high and sustained IFN‐β stimulation is needed to activate expression of specific genes involved in the regulation of inflammatory and adaptive immune responses, and cytokine production.

Since the ISGF3 potentiated phase has biological importance in regulating a subset of inflammatory genes, we investigated the mechanisms controlling this phase of the ISGF3 response. Given the timing of the response and the dependency on dose and duration, we hypothesized that a stimulus‐contingent positive feedback loop was regulating the ISGF3‐potentiated phase. To test this, the protein synthesis inhibitor cycloheximide (CHX) was incubated with MLE‐12 lung epithelial cells stimulated with IFN‐β to block induced feedback loops (Fig 5A). Time course measurements of ISGF3 activity revealed that while the initial phase of ISGF3 activity was not affected by the CHX treatment, the potentiated phase, which starts after the first hour, was blocked. The level of ISGF3 activity induced during the activation phase was sustained throughout the remainder of the time course, suggesting that a positive feedback loop is required to amplify the ISGF3 response. To identify potential gene candidates for the positive feedback amplifier, we measured STAT1, STAT2, and IRF9 mRNA abundances, known positive regulators of IFN signaling, during an IFN‐β stimulated time course using RT–qPCR (Fig 5B). By 60 min, IRF9 mRNA had already increased about eightfold, while STAT1 and STAT2 mRNA had doubled at that time. Measuring the nuclear STAT1, STAT2, and IRF9 protein, we also observed that IRF9 had the largest increase (sevenfold) during the potentiated phase (100 min) with respect to the activation phase (30 min) compared with STAT1 (40% increase) and STAT2 (30% increase; Fig 5C). This suggests that IRF9 induction, and to a lesser extent STAT1 and STAT2 induction, acts as a critical positive autoregulation circuit for ISGF3 and is needed for the stimulus‐contingent positive feedback loop.

Using the MLE‐12 IFN signaling model, we tested how IRF9 as well as STAT2 positive autoregulation during IFN‐β stimulation affects ISGF3 activity (Fig 5D). Model simulations predicted that when IRF9 induction is blocked the potentiated phase of the ISGF3 response is reduced, while the activation and late phase are similar to the control. Blocking the STAT2 positive autoregulation reduced the potentiated phase but also the late phase. Eliminating STAT1 autoregulation in the model simulation also eliminated the potentiated and late phases (Fig EV2) but this was not further explored experimentally as STAT1 also participates in other complexes (e.g., GAF) making downstream gene expression effects difficult to interpret. In order to test the model predictions, the ISGF3 binding site within the IRF9 and STAT2 promoters in the MLE‐12 lung epithelial cells was mutated using CRISPR‐Cas9 gene editing (Fig EV3A). This resulted in mutant variants in 79% of the IRF9 promoter region amplicons and 91% of the STAT2 promoter region amplicons (Fig EV3B). This significantly diminished IRF9 and STAT2 inducibility by IFN‐β, as assessed by RT–qPCR (Fig 5E). While the IRF9 promoter mutant and lentiviral control cells had similar basal levels of IRF9, at 1 h of stimulation the abundance of IRF9 was reduced in the mutant cells compared with the control (P‐value < 0.01). A similar reduction of STAT2 mRNA expression levels was detected in the STAT2 promoter mutant by 1 h of stimulation (P‐value < 0.05), suggesting diminished induction of the positive autoregulation in the IRF9 and STAT2 promoter mutants.

To examine the functional consequences of diminished autoregulation, the IRF9 and STAT2 promoter mutant cells were stimulated with IFN‐β, and ISGF3 activity at various time points measured using EMSA (Fig 5F). Similar to the model simulation, the potentiated phase was found to be diminished in both the IRF9 and STAT2 promoter mutants compared with the lentiviral control cells. To determine gene expression consequences, the mRNA of inflammatory instigator genes (e.g., IL6 and TNFSF10) as well as inflammatory regulator genes (e.g., CD274, TLR2, CCRL2, GBP5, CGAS, and MYD88) from cluster 1 during IFN‐β stimulation was measured with a time course using qRT–PCR (Fig 5G). For the IRF9 promoter mutant, the amount of mRNA abundance was decreased in at least one time point for TNFSF10, CD274, CCRL2, GBP5, and CGAS compared with the lentiviral control (P‐value < 0.05). For many of the genes including TLR2, MYD88, and GBP5, the peak expression level was diminished and delayed to 8 h. The STAT2 promoter mutant also had a reduction in mRNA abundance for CD274, GBP5, and CGAS (P‐value < 0.05), with shifts in peak expression for many of the genes compared with the lentiviral control. Taken together, this suggests a role for IRF9 and STAT2 positive autoregulatory feedback loops in controlling inflammatory gene expression responses during IFN‐β stimulation.

MLE‐12 lung epithelial cells were obtained from ATCC (#CRL‐2110) and cultured in DMEM with 4.5 g/l glucose, l‐glutamine, and sodium pyruvate (Corning #10‐013‐CV) supplemented with 10% fetal bovine serum, 100 IU/ml penicillin, 100 μg/ml streptomycin (Corning #30‐002‐CI), and 2 mM l‐glutamine (Corning #25‐005‐CI) at 37°C at 5% CO₂. For stimulations, MLE‐12 cells were incubated with 1 U/ml or 10 U/ml recombinant mouse IFN‐β (PBL Assay Science #12401‐1) and 100 ng/ml recombinant mouse IL‐28B/IFN‐λ3 (R&D Systems #1789‐ML‐025/CF). For select experiments, cells were also incubated with 10 μg/ml cycloheximide (Sigma‐Aldrich #C7698) or a vehicle (ethanol) control.

Protein from extracts was measured (BioRad DC™ Protein Assay Kit II #5000112), and equal amounts of protein were loaded into a 4–15% gradient polyacrylamide gel (BioRad TGX™ Precast Midi Protein Gel #5671085) and proteins separated based on molecular weight using electrophoresis at 110 V. Proteins were transferred onto a PVDF membrane (Fisher Scientific Immobilon®‐P IPVH00010) at 100 V for 1 h. The membrane was incubated for 1 h in 5% bovine serum albumin (BSA; Millipore Sigma #A3059) followed by overnight incubation at 4°C in a primary antibody solution. The following antibodies were used in this study: phospho‐STAT1 (pY701.4A) mouse monoclonal IgG2aĸ targeting to phosphorylated Tyr701 (Santa Cruz Biotechnology #sc‐136229), STAT1 p84/p91 (E‐23) rabbit polyclonal IgG with epitope matching near the C terminus (Santa Cruz Biotechnology #sc‐346), STAT1 p84/p91 rabbit polyclonal IgG antibody (Cell Signaling #9172), phospho‐STAT2 (Tyr689) rabbit polyclonal IgG antibody with epitope matching at phosphorylated Tyr689 (EMD Millipore #07–224), STAT2 (D9J7L) rabbit monoclonal IgG antibody with epitope corresponding near Leu706 (Cell Signaling #72604), IRF9, clone 6F1‐H5, mouse monoclonal IgG_2aĸ antibody with epitope matching the N terminus (EMD Millipore #MABS1920), nuclear matrix protein p84 [EPR5662(2)] rabbit monoclonal antibody mapping with aa350‐450 (abcam #ab131268), and αTubulin (B‐7) mouse monoclonal IgG2aĸ antibody raised against epitope with aa149‐448 (Santa Cruz Biotechnology #sc‐5286). The membrane was washed in TBS‐T (0.5% Tween‐20) and incubated in a secondary antibody solution for 1 h at room temperature. The following horseradish peroxidase (HRP)‐conjugated secondary antibodies were used for this study: mouse anti‐rabbit IgG (L27A9) monoclonal antibody (Cell Signaling #5127), goat anti‐rabbit IgG polyclonal antibody (Cell Signaling #7074S), horse anti‐mouse IgG polyclonal antibody (Cell Signaling #7076), mouse anti‐rabbit IgG monoclonal antibody (Santa Cruz Biotechnology #sc2357), and bovine anti‐mouse IgG polyclonal antibody (Santa Cruz Biotechnology #sc2380). The protein bands were resolved by chemiluminescence (Thermo Fisher Scientific SuperSignal West Pico and Femto Chemiluminescent Substrate #34080 and 34095).

The protocol was used as described previously (Hoffmann et al, 2002; Werner et al, 2005). In brief, nuclear fractions were collected by extraction using a Tris (250 mM Tris, 60 mM KCl, 1 mM EDTA) buffer solution. Equal amounts of nuclear protein were incubated for 30 min in a binding buffer (10 mM Tris‐Cl, 50 mM NaCl, 10% glycerol, 1% NP‐40, 1 mM EDTA, 0.1 mg/ml Poly(deoxyinosinic‐deoxycytidylic) acid sodium salt) with 0.01 pmol P‐32 radiolabeled oligo containing the ISRE consensus sequence (ISG15 loci; 5′‐GATCCTCGGGAAAGGGAAACCTAAACTGAAGCC‐3′; 5′‐GGCTTCCAGTTTAGGTTTCCCTTTCCCGAGGATC‐3′) or a NFY binding sequence (5′‐GATTTTTTCCTGATTGGTTAAA‐3′; 5′‐ACTTTTAACCAATCAGGAAAAA‐3′). Samples were loaded onto a 5% polyacrylamide gel (5% glycerol and TGE buffer [24.8 mM Tris base, 190 mM glycine, 1 mM EDTA]). Bands were resolved using electrophoresis at 200 V. Gels were dried and imaged using an Amersham Typhoon 5 laser scanner (GE Healthcare).

RNA was collected from cells lysed in TRIzol Reagent (Fisher Scientific #15–596‐018) and isolated using the Direct‐zol RNA Extraction Kit (Zymo Research #R2052). Equal amounts of RNA were reverse transcribed into cDNA using the Iscript™ Reverse Transcription Supermix (BioRad #1708841). Equal volumes of diluted (1:8) cDNA were loaded for qPCR analysis along with the SsoAdvanced Universal SYBR Green Supermix (BioRad #1725272) and 0.5 μM forward and reverse primers for target genes (listed below). Amplification curve thresholds were determined and selected for each primer set.

mSTAT1 (5′‐GGCCTCTCATTGTCACCGAA‐3′; 5′‐TGAATGTGATGGCCCCTTCC‐3′);

mSTAT2 (5′‐GTAGAAACCCAGCCCTGCAT‐3′; 5′‐CTTGTTGCCCTTTCCTGCAC‐3′);

mIRF9 (5′‐TCTTTGTTCAGCGCCTTTGC‐3′; 5′‐CTGCTCCATCTGCACTGTGA‐3′);

mSOCS1 (5′‐CAACGGAACTGCTTCTTCGC‐3′; 5′‐AGCTCGAAAAGGCAGTCGAA‐3′);

mSOCS3 (5′‐CTTTTCTTTGCCACCCACGG‐3′; 5′‐CCGTTGACAGTCTTCCGACA‐3′);

mUSP18 (5′‐CTCACATGTTTGTTGGGTCACC‐3′; 5′‐TGAAATGCAGCAGACAAGGG‐3′);

mTNFSF10 (5′‐GATGAAGCAGCTGCAGGACAAT‐3′; 5′‐CTGCAAGCAGGGTCTGTTCAAG‐3′);

mCD274 (5′‐TCGCCTGCAGATAGTTCCCAAA‐3′; 5′‐GTAAACGCCCGTAGCAAGTGAC‐3′);

mCCRL2 (5′‐GAGCAAGGACAGCCTCCGAT‐3′; 5′‐CCACTGTTGTCCAGGTAGTCGT‐3′);

mGBP5 (5′‐TGCTGACATGAGCTTCTTCCCA‐3′; 5′‐TCATCGCTACCTTGCTTCAGCT‐3′);

mCGAS (5′‐TGGGCACAAAAGTGAGGACCAA‐3′; 5′‐CGCCAGGTCTCTCCTTGAAAACT‐3′);

mMYD88 (5′‐TCTCCAGGTGTCCAACAGAAGC‐3′; 5′‐TGCAAGGGTTGGTATAGTCGCA‐3′);

mTLR2 (5′‐GGGGCTTCACTTCTCTGCTT‐3′; 5′‐AGCATCCTCTGAGATTTGACG‐3′);

mIL6 (5′‐ACCAGAGGAAATTTTCAATAGGC‐3′; 5′‐TGATGCACTTGCAGAAAACA‐3′);

mGAPDH (5′‐AGCTTGTCATCAACGGGAAG‐3′; 5′‐TTTGATGTTAGTGGGGTCTCG‐3′).

RNA was collected from cells lysed in TRIzol Reagent (Fisher Scientific #15–596‐018) and isolated using the Direct‐zol RNA Extraction Kit (Zymo Research #R2052). RNA libraries were generated using the KAPA Stranded mRNA‐Sequencing Kit with KAPA mRNA Capture Beads (Roche #07962207001). Samples were enriched for poly(A) + mRNA strands using mRNA Capture Beads. The captured mRNA was synthesized into complimentary DNA (cDNA). cDNA libraries were ligated using Illumina TruSeq single index adapters (Roche #KK8700), amplified, and quantified using a Qubit dsDNA Broad Range Assay Kit (Life Technologies #Q32850) and a Qubit 2.0 fluorometer. Sequencing was conducted on an Illumina HiSeq 2500 with single‐end 50 bp reads through the UCLA Broad Stem Cell Research Center. Sequencing reads were aligned and mapped onto the mm10 genome. The EdgeR package (Bioconductor 3.7) was used for differential gene expression analysis of raw read counts, which were normalized using trimmed mean of the mean (TMM) and filtered by the number of mapped reads and gene length (RPKM) of greater than 2. Normalized counts were applied to a negative binomial generalized linear model, which was used to calculate biological coefficients of variation. Pairwise dispersions between the normalized counts were calculated using the Cox‐Reid Likelihood profile method. The differential expression testing was initiated using the likelihood ratio test (LRT), a form of ANOVA testing. Using the TREAT Method and Benjamini–Hochberg multiple test correction, significantly induced genes were selected based on differences in RPKM values at each time point compared with values in unstimulated conditions using a false discovery rate (FDR) of less than 0.05 based on replicate data and fold change greater than 2. PCA plots were generated using the prcomp function in R and plotted with ggplot2 (Gómez‐Rubio, 2017). Correlation plots were plotted using ggplot2. Heatmaps were generated using the pheatmap package.

The protocol was used as described previously (Kishimoto et al, 2021). In brief, cells were fixed with 2 mM disuccinimidyl glutarate for 30 min followed by a 10‐min incubation in 1% formaldehyde. Cells are lysed (Buffer 1: 50 mM HEPES‐KOH, 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); Buffer 2: 10 mM Tris–HCl, 200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, and 1X protease inhibitor cocktail; Buffer 3: 10 mM Tris–HCl, 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 to shear the chromatin. The DNA‐bound complexes were incubated with STAT1 p84/p91 rabbit polyclonal IgG antibody (Cell Signaling #9172) overnight at 4°C. STAT1 was immunoprecipitated by incubating for 5 h at 4°C with Protein‐G DynaBeads (Invitrogen #10004D). DNA bound to STAT1 complexes was extracted using 1% SDS and 0.6 mg/ml of Proteinase K (New England Biolabs #P81075) at 60°C overnight, followed by purification with Agencourt AMPure XP solid‐phase reversible immobilization (SPRI) beads (Beckman Coulter #A63881). DNA was quantified using a Qubit dsDNA High Sensitivity Assay Kit (Life Technologies #Q32851) and a Qubit 2.0 fluorometer. Libraries were prepared using NEBNext Ultra II DNA Library Prep Kit NEBNext Multiplex Oligos (New England Biolabs) according to the manufacturer instructions and sequenced with a length of 50 bp on an Illumina HiSeq 2500 at the UCLA Broad Stem Cell Research Center. Sequencing data were aligned to mm10 genome as described previously (Kishimoto et al, 2021). MACS2 version 2.1.0 was used to identify peaks for each sample using pooled input samples as control with FDR < 0.01 and extension size of average fragment length (Zhang et al, 2008). STAT1 motifs were identified by searching the HOMER database for position‐weighted matrix files of STAT/GAS and ISRE/IRF motifs (Heinz et al, 2010). 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. A variant GAS motif in the promoter of Irf1 (GATTTCCCCGAATG) known to be bound by STAT1 was also included in the search as a GAS motif (Decker et al, 1991). Peak‐to‐gene distances were calculated with bedtools, using the TSS of the longest isoform (Quinlan & Hall, 2010).

The maximum RPKM value and gene names for the selected 188 genes in cluster 1 were analyzed using the Panther Classification System Database DOI: https://doi.org/10.5281/zenodo.4495804 Released 2021‐02‐01 (Mi et al, 2019). The genes were mapped to the mm10 genome and analyzed using the statistical overrepresentation test, which determines overrepresentation of gene ontology annotations in cluster 1 compared with all expressed genes using the Fisher's exact test and FDR correction. A list of GO biological process terms and IDs as well as associated genes was generated. Each GO term was categorized into 10 groups based on function. To determine whether a function was enriched in a gene cluster the fraction of GO terms for each biological function compared with the total number of terms for the cluster was calculated and plotted as a pie graph. Similar analysis was conducted on the 157 genes in cluster 2.

The cloning strategy for the lentiCRISPRv2 plasmid was used (Sanjana et al, 2014; Shalem et al, 2014). Briefly, the following oligos were annealed and cloned into a BsmBI (Thermo Fisher Scientific #FD0454)‐digested lentiCRISPRv2 vector: STAT2 promoter mutant (5′‐CACCGAGGGAAAGGAAACTGAAACC‐3′; 5′‐AAACGGTTTCAGTTTCCTTTCCCTC‐3′); IRF9 promoter mutant (5′‐CACCGACTCAGACCACGTGGTTTCT‐3′; 5′‐AAACAGAAACCACGTGGTCTGAGTC‐3′). Plasmids were transformed into One Shot™ Stbl3™ Chemically Competent bacteria (Thermo Fisher Scientific #C737303). Single colonies were selected and inoculated in LB Broth. Plasmids were isolated from expanded bacterial colonies using the ZymoPURE™ II Plasmid Kit (Zymo Research #D4203). Ten μg of isolated plasmid was transfected into HEK293T cells using Lipofectamine® 2000 Transfection Reagent (Invitrogen #11668–019). Two days following transfection, viral particles were isolated from the HEK293T‐conditioned media. MLE‐12 lung epithelial cells were infected with the purified viral particles for 24 h, followed by 2 μg/ml puromycin antibiotic selection.

Sequence mutations in the IRF9 and STAT2 promoter mutant cells were confirmed using Sanger sequencing. In brief, 50 ng of genomic DNA (gDNA) was amplified using PCR with 0.5 μM of forward and reverse primers: STAT2 promoter mutant (5′‐CCTATCCTAAGGGACGCAGG‐3′; 5′‐GCCCAACTAAAGTCTTAGCC‐3′); IRF9 promoter mutant (5′‐CAAGGTGCTACTGCTGACTG‐3′; 5′‐TTTCTGATCTCTCCTCCCC‐3′). PCR amplicons were resolved on 2% agarose gels and DNA extracted using the QIAquick Gel Extraction Kit (Qiagen #28706). Ten ng amplicons with 5 μM primers were submitted for sanger sequencing to Laragen, Inc. To quantify the number of sequence variants in the IRF9 and STAT2 promoter mutant cells, PCR amplicons (500 ng) were submitted for Amplicon‐EZ sequencing to Genewiz®.

An ordinary differential equation (ODE) mathematical model was built to simulate RNA abundance levels as a function of transcription factor binding and RNA degradation. Briefly, the median RPKM values for all genes within each cluster from the RNA sequencing data were calculated for each time point. The temporal gene expression trajectory of the median values was generated using a modified Akima interpolant fit (MATLAB). The model was composed of four parameters (e.g., k_act, activation rate; k_deg, mRNA degradation rate, K_a, association constant, and n, Hill coefficient). The model was parameterized by model fitting to the gene expression trajectories normalized to max peak. For optimization, the simplex search method (e.g., fminsearch) was implemented using an objective function of the sum of RMSD values calculated at each time point. Fifty multiple random seeds for the initial parameter values were used to optimize over a larger parameter space. The parameter optimization was run 5 times, and the top 10 parameter sets with the lowest RMSD values were averaged and used as the optimized model parameter values. Heatmaps and line graphs were plotted using MATLAB.

We used the previously published model of interferon signaling (Kok et al, 2020), with minor modifications. We replaced the equation for the active receptor complex formation by instead of to be able to simulate the model in presence of cycloheximide which blocks protein synthesis, leading to SOCS protein decay, ultimately converging to 0. To better fit our data, we moved the delay term prior to mRNA synthesis instead prior to protein synthesis.

There are 41 molecular species in the mathematical model (Table EV1). The model parameters (Table EV2) were fit to the data, optimizing parameter value as well as which parameter to change iteratively 3 by 3 in order to adjust the minimal number of parameters allowing for a good fit. Parameters were selected at each iteration by the algorithm based on optimizing the best fit. For comparison, a particle swarm optimization algorithm (740 particles) was used to fit all parameters concurrently.

The modeling and parameter optimization were done in MATLAB R2018a using the particleswarm function implemented in MATLAB's optimization toolbox.