**3. Results**

### *3.1. Distinct Induction Profiles of Antiviral and Immunoregulatory Genes in Response to IFN*β*, TNF and IFN*β + *TNF*

First, to validate previous observations, we sought to determine the induction profile by the combination of IFNβ and TNF of a selected panel of immunoregulatory and antiviral genes that had previously been shown to be regulated by IFNβ or TNF alone. A549 cells, in which we previously documented the synergistic action of IFNβ and TNF on the *DUOX2* gene, were stimulated either with IFNβ, TNF, or IFNβ + TNF for various times between 3–24 h and the relative mRNA expression levels were quantified by qRT-PCR. Analysis of the expression of the selected genes revealed distinct profiles of response to IFNβ, TNF, or IFNβ + TNF (Figure 1). *IDO*, *DUOX2*, *CXCL10*, *APOBEC3G, ISG20,* and *IL33* exhibited synergistic induction in response to IFNβ + TNF compared to IFNβ or TNF alone. Expression in response to IFNβ + TNF increased over time, with maximum expression levels observed between 16 and 24 h. While *NOD2* and *IRF1* induction following stimulation with IFNβ + TNF was also significantly higher than upon IFNβ or TNF single cytokine stimulation, they exhibited a steady-state induction profile starting as early as 3 h. *MX1* and *PKR*, two typical IFNβ-inducible ISGs, were found induced by IFNβ + TNF similarly to IFNβ alone. *CCL20* responded to IFNβ + TNF with a kinetic and amplitude similar to TNF, but was not responsive to IFNβ alone. *IL8* expression was not induced by IFNβ, but was increased by TNF starting at 3 h and remained steady until 24 h. In contrast to other genes, *IL8* induction in response to IFNβ + TNF was significantly decreased compared to TNF alone. Overall, these results confirm previous reports that a subset of antiviral and immunoregulatory genes is greatly increased in response to IFNβ + TNF compared to either IFNβ or TNF alone.

**Figure 1.** Expression of a panel of antiviral and immunoregulatory genes in response to IFNβ, TNF, and IFNβ + TNF. A549 cells were stimulated with either TNF (T), IFNβ (I), or costimulated with IFNβ + TNF (I + T) for the indicated times. Quantification of mRNA was performed by qRT-PCR and expressed as fold expression after normalization to the S9 mRNA levels using the ΔΔCt method. Mean +/− SEM, *n* ≥ 3. Statistical comparison of TNF vs. IFNβ + TNF and IFNβ vs. IFNβ + TNF was conducted using two-way ANOVA with Tukey's post-test. *p* < 0.01 (\*\*), *p* < 0.001 (\*\*\*) or *p* < 0.0001 (\*\*\*\*).

### *3.2. Workflow for Genome-Wide Characterization of the Delayed Transcriptional Program Induced by IFN*β+*TNF in the Absence of STAT1*

In a previous study, we provided evidence supporting the existence of a STAT1-independent, but STAT2- and IRF9-dependent, pathway engaged downstream of IFNβ + TNF [14]. Here, using the STAT1-deficient human U3A cell line [15], we aimed to fully characterize the STAT1-independent transcriptional program induced by IFNβ + TNF. The human U3A cell line was derived by mutagenesis from the 2ftGH cells [15], in which a synergistic response to IFNβ + TNF can be observed on a subset of genes (Figure 2A). Two hallmarks of STAT2 and IRF9 activation, i.e., STAT2 Tyr690 phosphorylation and induction of IRF9, were observed in the U3A cells following stimulation with IFNβ + TNF, although to reduced levels compared to the parental 2ftGH cells expressing endogenous STAT1 (Figure 2B). This observation implies that the activation of STAT2 and IRF9 depends to a large extent on the STAT1-dependent canonical pathway, but that a significant response occurs in the absence of STAT1. Therefore, the human U3A cell model is suitable for specifically studying STAT1-independent, but STAT2- and IRF9-dependent, gene expression in response to IFNβ + TNF.

**Figure 2.** Experimental design used to study the STAT1-independent delayed transcriptional program induced by the combination of IFNβ and TNF. ( **A**) 2ftGH cells were stimulated with either TNF (T), IFNβ (I), or costimulated with IFNβ + TNF (I + T) for 24 h. Quantification of mRNA was performed by qRT-PCR and expressed as fold expression after normalization to the S9 mRNA levels using the ΔΔCt method. Mean +/− SEM, *n* ≥ 5. Statistical comparison was conducted using one-way ANOVA with Tukey's post-test. *p* < 0.05 (\*), *p* < 0.01 (\*\*), *p* < 0.001 (\*\*\*), or *p* < 0.0001 (\*\*\*\*). (**B**) U3A (STAT1-deficient), 2ftGH (parental STAT1-positive) cells and U3A-STAT1 cells (U3A cells stably reconstituted with STAT1) were left untreated or stimulated with IFNβ + TNF for the indicated times. WCE (whole cell extracts) were analyzed by SDS-PAGE followed by immunoblot using anti STAT1-P-Tyr701, total STAT1, STAT2-P-Tyr690, total STAT2, IRF9, or actin antibodies. ( **C**–**E**) U3A cells were transfected with siCTRL, siSTAT2, or siIRF9 before being left untreated (NS) or stimulated with IFNβ + TNF for 24 h. ( **C**) The schematic describes the workflow of sample preparation and analysis. ( **D**) WCE were analyzed by SDS-PAGE followed by immunoblot using anti STAT2, IRF9, and actin antibodies. (**E**) Graph showing the correlation between fold-changes (FC) measured by RNASeq and qRT-PCR for 13 randomly selected genes. Data from siCTRL NS vs. siCTRL IFNβ +TNF, siSTAT2 IFNβ +TNF vs. siCTRL IFNβ + TNF, siIRF9 IFNβ + TNF vs. siCTRL IFNβ + TNF conditions were used.

To profile the genome wide transcriptional program induced by the combination of IFNβ and TNF in the absence of STAT1 and define the role of STAT2 and IRF9, the U3A cells were transfected with Control (Ctrl)<sup>−</sup>, STAT2- or IRF9-RNAi and further left untreated or stimulated with IFNβ (1000 U/mL) + TNF (10 ng/mL) for 24 h (Figure 2C). Efficient silencing was confirmed by immunoblot (Figure 2D). Total RNA was isolated and analyzed by RNA sequencing (n = 3 for each group detailed in Figure 2B) on an Illumina HiSeq2500 platform. To validate the expression profile of genes in RNASeq results, the fold changes (FC) of 13 genes randomly selected were analyzed by qRT-PCR in each experimental groups, i.e., siCTRL non-stimulated (NS) vs. siCTRL IFNβ + TNF, siCTRL IFNβ + TNF vs. siSTAT2 IFNβ + TNF and siCTRL IFNβ + TNF vs. siIRF9 IFNβ + TNF. A positive linear relationship between RNASeq and qRT-PCR FC was observed (Figure 2E).

### *3.3. A Broad Antiviral and Immunoregulatory Transcriptional Signature is Induced by IFN*β + *TNF in the Absence of STAT1*

To identify differentially expressed genes (DEGs) upon IFNβ + TNF stimulation in the absence of STAT1, comparison between non-stimulated (NS) and IFNβ + TNF-stimulated control cells was performed (Figure 3A). In total, 612 transcripts, including protein-coding transcripts, pseudogenes and long non-coding RNA (lncRNA), were significantly different (FC > 1.5, *p* < 0.05, FDR < 0.05) in IFNβ + TNF vs. NS. Among these, 482 DEGs were upregulated and 130 were downregulated (Figure 3B; See Supplemental Table S1 for a complete list of DEGs).

**Figure 3.** Analysis of STAT1-independent IFNβ + TNF-induced DEGs. (**A**) Diagram describing the bioinformatics analysis strategy used to determine STAT1-independent differentially expressed genes (DEGs) and their regulation by STAT2 and IRF9. (**B**) Volcano plots of the fold-change (FC) vs. adjusted *p*-value of IFNβ + TNF (I + T) vs. non-stimulated (NS) siCtrl conditions. (**C**) Volcano plots of the fold-change vs. adjusted *p*-value of siSTAT2 IFNβ + TNF vs. siCTRL IFNβ + TNF (I + T) conditions. (**D**) Volcano plots of the fold-change vs. adjusted *p*-value of siIRF9 IFNβ + TNF vs. siCTRL IFNβ + TNF conditions.

To identify the Biological Processes (BP) and Molecular Functions (MF) induced by IFNβ + TNF independently of STAT1, we further analyzed the upregulated DEGs. The top 40 upregulated DEGs are shown in Figure 4A. We subjected the complete list of upregulated DEGs (Supplemental Table S1) through Gene Ontology (GO) enrichment analysis. The top enriched GO BP (*p* < <sup>10</sup>−10) and MF, are depicted in Figure 4B (See Supplemental Table S2 for a complete list of enriched GO). The majority of the top enriched BPs were related to cytokine production and function (response to cytokine, cytokine-mediated signaling pathway, cytokine production, and regulation of cytokine

production), immunoregulation (immune response, immune system process, innate immune response, and regulation of immune system process) and host defense response (defense response, response to other organism, 2'-5'-oligoadenylate synthetase activity and dsRNA binding). Fourteen MF categories were found enriched in IFNβ + TNF. The top enriched MFs were related to cytokine and chemokine functions (cytokine activity, cytokine receptor binding, chemokine activity, and Interleukin 1-receptor binding). Other enriched MF included peptidase related functions (endopeptidase inhibitor activity, peptidase regulator activity, and serine-type endopeptidase activity).

**Figure 4.** Functional characterization of STAT1-independent IFNβ + TNF-induced DEGs. (**A**) Top forty IFNβ + TNF- upregulated DEGs. (**B**) Gene ontology (GO) enrichment analysis of the differentially upregulated genes in IFNβ + TNF vs. non-stimulated siCtrl conditions based on the Biological Processes and Molecular Function categories. Top enriched terms are shown and the full list is available in Supplemental Table S2. (**C**) Modular transcription analysis of upregulated DEGs. Eighteen enriched modules are shown. The full list of enriched modules is available in Supplemental Table S3. (**D**) U3A and STAT1-rescued U3A-STAT1 cells were stimulated with IFNβ (I) or IFNβ + TNF (I + T) for 30 h before infection with VSV at a MOI of 5 for 12 h. The release of infectious viral particles was quantified by plaque forming unit (pfu) assay. The left graphs show dot-plots of all stimulations. Statistical comparisons were performed on the "before and after" plots (displayed on the right of dot-plot graphs) using ratio paired *t*-tests.

To gain deeper insight into the relevance of the STAT1-independent IFNβ + TNF-induced transcripts towards immunological and host defense responses, we conducted a modular transcription analysis of upregulated DEGs against 606 immune-related functional modules. These modules were previously defined from co-clustered gene sets built via an unbiased data-driven approach as detailed in the material and methods section [30,31]. STAT1-independent IFNβ + TNF-induced DEGs showed significant enrichment in 37 modules (See Supplemental Table S3 for a complete list of enriched modules). Six of these modules were associated with virus sensing/Interferon antiviral response, including LI.M75 (antiviral IFN signature), LI.M68 (RIG-I-like receptor signaling), LI.M127 (type I interferon response), LI.M111.0 and LI.M111.1 (IRF2 target network), and LI.M150 (innate antiviral response) (Figure 4C). Additionally, six modules were associated with immunoregulatory functions, including LI.M29 (proinflammatory cytokines and chemokines), LI.M27.0 and LI.M27.1 (chemokine cluster I and II), LI.M38 (chemokines and receptors), LI.M115 (cytokines receptors cluster), and LI.M37.0 (immune activation - generic cluster) (Figure 4C). Module analysis also showed enriched AP-1 transcription factor-related network modules, LI.M20 and LI.M0, and cell cycle and growth arrest LI.M31 module.

GO enrichment and modular transcription analyses were suggestive of an antiviral response being induced by IFNβ + TNF. To demonstrate the physiological relevance of STAT1-independent gene expression, we evaluated the capacity of IFNβ + TNF to restrict virus replication in the absence of STAT1. Previous study has shown that infection by Vesicular Stomatitis Virus (VSV) is sensitive to IFNβ, but not to TNF, in 2ftGH-derived cells [35]. Thus, U3A cells were stimulated with IFNβ alone or IFNβ + TNF and further infected with (VSV). While VSV replicated similarly in untreated U3A cells and in cells treated with IFNβ, treatment with IFNβ + TNF significantly limited VSV replication (Figure 4D). As a control of the e fficiency of IFNβ alone treatment, we observed a significant antiviral response when STAT1 expression was stably restored to endogenous levels in U3A cells (U3A-STAT1) (Figures 2B and 4D). Collectively, these data demonstrate that the combination of IFNβ + TNF induces an e fficient antiviral and immunoregulatory transcriptional response in the absence of STAT1.

### *3.4. Clustering of STAT1-Independent IFN*β + *TNF Induced DEGs According to their Regulation by STAT2 and IRF9*

Having shown that IFNβ + TNF induce a broad antiviral and immunoregulatory transcriptional response independently of STAT1, we next sought to gain insight into the role of STAT2 and IRF9. To do so, we compared transcripts levels in siSTAT2\_IFNβ + TNF vs. siCTRL\_IFNβ + TNF and siIRF9\_IFNβ + TNF vs. siCTRL\_IFNβ + TNF conditions (Figure 3A and Supplemental Table S1). Volcano plots revealed that a fraction of IFNβ + TNF-induced DEGs were significantly (FC > 1.5, *p* < 0.05, FDR < 0.05) downregulated or upregulated upon silencing of STAT2 (Figure 3C) or IRF9 (Figure 3D). Nine distinct theoretical categories of DEGs can be defined based on their potential individual behavior across siSTAT2 and siIRF9 groups (Categories A–I, Figure 5A): a gene can either be downregulated upon STAT2 and IRF9 silencing, indicative of a positive regulation by STAT2 and IRF9 (Category A); conversely, a gene negatively regulated by STAT2 and IRF9 would exhibit upregulation upon STAT2 and IRF9 silencing (Category B); genes that do not exhibit significant di fferential expression in siSTAT2 and siIRF9 groups would be classified as STAT2 and IRF9 independent (Category C); IRF9-independent genes could exhibit positive (Category D) or negative (Category E) regulation by STAT2; conversely, STAT2-independent genes might be positively (Category F) or negatively (Category G) regulated by IRF9; lastly, STAT2 and IRF9 could have opposite e ffects on a specific gene regulation (Category H and I). Based on a priori clustering of RNASeq data we found that in the absence of STAT1 IFNβ + TNF-induced DEGs clustered into only seven of the nine possible categories (Figure 5B). The top 15 upregulated DEGs of each category are shown in Figure 5C. The complete list of genes in each category is available in Supplemental Table S1. Amongst the 482 DEGs, 163 genes exhibited either inhibition or upregulation following silencing of STAT2 and/or IRF9 (Categories B–G). A large majority of upregulated DEGs, i.e., 319 out of the 482 DEGs, were not significantly a ffected by either STAT2 or IRF9 silencing, and were therefore classified as STAT2/IRF9-independent (Figure 5B). No genes were found in Category H and only one gene was found in Category I.

**Figure 5.** Clustering of IFNβ + TNF-induced DEGs according to their regulation by STAT2 and IRF9. (**A**) Theoretical categories in which IFNβ + TNF-induced DEGs can be segregated based on their potential individual regulation by STAT2 and IRF9. (**B**) Hierarchical clustering of the categories of DEG responses according to their regulation by STAT2 and IRF9. Euclidean distance metric is used for the construction of distance matrix and the categories are used as a priori input into clustering algorithm as detailed in Materials and Methods. (**C**) Top 15 induced DEGs (Log2FC, siCTRL NS vs. siCTRL IFNβ + TNF) in each category identified in (**B**) are displayed. Note that category D contains only twelve genes, so all genes are shown. The full list of genes is available in Supplemental Table S1. (**D**) Diagram showing enriched transcription modules in each gene category.

*Cells* **2019**, *8*, 919

To functionally interpret these clusters, we applied the modular transcription analysis to each of the categories to assess for the specific enrichment of the functional modules found associated with IFNβ + TNF-upregulated DEGs (Figure 5D). First, most modules, except LI.M31 (cell cycle and growth arrest), LI.M38 (chemokines and receptors), LI.M37.0 (immune activation - generic cluster), and LI.M53 (inflammasome receptors and signaling), were enriched in the category of DEGs positively regulated by STAT2 and IRF9 (Category A). Conversely, the cluster negatively regulated by STAT2 and IRF9 (Category B) exclusively contains enriched LI.M38 (chemokines and receptors), LI.M37.0 (immune activation - generic cluster) and LI.M115 (cytokines receptors cluster). The cluster of IRF9-independent genes that are negatively regulated by STAT2 (Category E) only exhibited enrichment in the virus sensing/IRF2 target network LI.M111.0 module, while the IRF9-independent/STAT2-positively regulated cluster (Category D) encompasses antiviral and immunoregulatory functions. The STAT2-independent but IRF9 positively regulated transcripts (Category F) were mainly enriched in modules associated with the IFN antiviral response, including LI.M75 (antiviral IFN signature), LI.M68 (RIG-I-like receptor signaling), LI.M127 (type I interferon response), and LI.M150 (innate antiviral response). Finally, the STAT2-independent but IRF9 negatively regulated cluster (Category G) was mostly enriched in modules associated with immunoregulatory functions, including LI.M29 (proinflammatory cytokines and chemokines), LI.M27.0 and LI.M27.1 (chemokine cluster I and II), LI.M78 (myeloid cell cytokines), but also with cell cycle and growth arrest (LI.M31) and inflammasome receptors and signaling (LI.M53). Of note, all modules were found enriched in the cluster of genes induced independently of STAT2 and IRF9 (Category C), pointing to a broad function of this pathway(s) in the regulation of the antiviral and immunoregulatory program elicited by IFNβ + TNF. Altogether, these observations reveal that STAT2 and IRF9 are involved in the regulation of a subset of the genes induced in response to the co-stimulation by IFNβ and TNF in the absence of STAT1. Importantly, our results also reveal that STAT2 and IRF9 act in part in a concerted manner, but also independently in distinct non-canonical pathways, to regulate specific subsets of the IFNβ + TNF-induced DEGs.

### *3.5. Di*ff*erential Regulation of CXCL10 in Response to IFN*β *and IFN*β+*TNF*

The canonical ISGF3 complex mediates ISGs transcriptional regulation through binding to ISRE consensus sequences [7]. Identification of DEGs upregulated by IFNβ+TNF in a STAT1-independent, but STAT2 and IRF9-dependent, manner raised the question of the ISRE site usage compared to the ISGF3 pathway. To address this question, we studied the regulation of the *CXCL10* gene promoter that was found induced by STAT2 and IRF9 in the absence of STAT1 in response to IFNβ and TNF (Supplemental Table S1 and Figure 5C), but is also inducible by IFNβ alone in the presence of STAT1 (Figures 1 and 2A). The *CXCL10* promoter contains three ISRE sites. We used full-length wild-type (972bp containing the three ISRE sites), truncated (376bp containing only the ISRE (3) site) or mutated (972bp containing a mutated ISRE (3) site) *CXCL10* promoter luciferase (*CXCL10*prom-Luc) reporter constructs (Figure 6A) to determine the ISRE consensus site(s) requirement. U3A and STAT1-rescued U3A-STAT1 cells were transfected with the *CXCL10*prom-Luc constructs and further stimulated with IFNβ or IFNβ+TNF to monitor the canonical ISGF3-dependent and the STAT1-independent responses. In STAT1-deficient U3A cells, IFNβ was unable to activate the *CXCL10*prom reflecting the dependence on the ISFG3 pathway. In contrast, induction of *CXCL10*prom was significantly induced when cells were stimulated with IFNβ +TNF. This induction in the absence of STAT1 was not altered by the deletion of ISRE (1) and (2) sites, but was significantly impaired when the ISRE (3) site was mutated (Figure 6B). STAT1 expression rescue led to a higher induction of *CXCL10*prom by IFNβ + TNF. Additionally, IFNβ-mediated induction of *CXCL10*prom was restored, albeit to a much lower extent than IFNβ +TNF. The activation of the *CXCL10*prom in the presence of STAT1 involved both the distal ISRE (1) and/or ISRE (2) sites and the proximal ISRE (3) site (Figure 6B). Altogether, this shows that the STAT1-independent, but STAT2/IRF9-dependent, pathway mediates gene expression through a restricted ISRE site usage compared to the ISGF3-dependent regulation.

**Figure 6.** Analysis of the *CXCL10* promoter regulation in response to IFNβ + TNF vs. IFNβ. (**A**) Schematic representation of the *CXCL10* promoter (CXCL10prom) luciferase constructs used in this study indicating the main transcription factors consensus binding sites. (**B**) U3A and U3A-STAT1 cells were transfected with the indicated CXCL10prom-luciferase reporter constructs and either left untreated or stimulated with IFNβ or IFNβ + TNF. Relative luciferase activities were measured at 16 h post-stimulation and expressed as fold over the corresponding unstimulated condition. Mean +/− SEM, n = 6. Statistical analyses were performed using an unpaired t-test comparing each promoter to the CXCL10prom-972bp construct. *p* < 0.01 (\*\*), *p* < 0.001 (\*\*\*), and *p* < 0.0001 (\*\*\*\*).
