Temporally restricted activation of IFNβ signaling determines response to immune checkpoint therapy

Little is known about the dynamic biological events that underpin therapeutic ecacy in immune checkpoint blockade (ICB) in cancer, due to the inability to frequently sample tumors in patients. Here, we mapped the transcriptional proles of 144 responding and non-responding tumors within two mouse models at four time points during ICB. We found that responding tumors displayed on/fast-off kinetics of type-I-interferon (IFN) signaling. Phenocopying of this kinetics using time-dependent sequential dosing of recombinant IFNs and neutralizing anti-bodies markedly improved ICB ecacy, but only when IFNβ was targeted, not IFNα. We identied Ly6C+/CD11b+ inammatory monocytes as the primary source of IFNβ and found that active type-I-IFN signaling in tumor-inltrating inammatory monocytes was associated with T cell expansion in patients treated with ICB. Together, our results suggest that on/fast-off modulation of IFNβ signaling is critical to the therapeutic response to ICB, which can be exploited to drive clinical outcomes towards response. Following we performed differential expression with Sleuth (v0.29.0) We compared responders and non-responders using a model containing time-point and response as covariates using a likelihood ratio test. We aggregated p-values from transcript differential expression to gene-level results with transcript-to-gene mapping relying on the latest Gencode reference M25 (GRCm38.p6) using Lancaster’s method 58 . Genes were deemed differentially at a false-discovery rate of less than 5%, regardless of fold change S1).


Introduction
The response to immune checkpoint blockade (ICB) in cancer is highly variable, with a majority of patients experiencing disease progression. Although the targets of ICB antibodies are known, the downstream therapeutic effector mechanisms are incompletely understood 1,2 . While speci c aspects of the pre-treatment tumor microenvironment, such as PD-L1 expression (Reck et al., 2016), immune cell in ltration 3 or tumor mutation burden 4 have been shown to correlate with response, none of these biomarkers are su ciently robust to guide clinical decisions regarding treatment across cancers 5

nor has
their characterization yet resulted in approved treatments that improve the e cacy of ICB 6 . As a consequence, the development of novel combination therapies to improve outcomes is mainly empiric 7 .
When perturbing complex systems, some important effects may only become apparent over time. For example, in the course of an immune response, some in ammatory mediators need to be switched off, not only to resolve in ammation, but also to mediate a transition from innate to adaptive immunity 8 . Similar time-dependent mechanisms may underpin an effective anti-tumor immune response 5 . However, it is not possible to identify these dynamic events in human studies due to the di culty to sample the same tumor at multiple time points during treatment, usually limiting the number of samples to two; one pre-treatment and one on-treatment sample several weeks after start of ICB 9,10 . In addition, it is exceedingly di cult to identify discrete biological differences between responding and nonresponding patient tumors due to inter-individual variation in germline and cancer genetics, tumor microenvironment composition and environmental in uences [11][12][13][14] .
For these reasons, we optimized murine models with bilateral tumors, derived from syngeneic cancer cell lines 6 . In these models, the response to ICB with antibodies against CTLA4 and PD-L1 is symmetric, leading to either a bilateral response or a failure to respond in both tumors 6,15,16 . This means we can remove one tumor for molecular analysis while tracking the therapeutic fate of the other. Because both responders and non-responders are exposed to identical treatments, this enables the characterization of the time-dependent response to ICB in entire tumors in a highly homogenous background.
We have previously used these bilateral models to identify transcriptomic ICB response signatures in the early tumor microenvironment and extensively validated our ndings in separate patient cohorts of bladder cancer and melanoma patients treated with ICB 6,16 . In addition, our results using these models have been independently validated by other research groups in preclinical models and patient cohorts with various tumor types, underscoring their translational relevance [17][18][19][20] . Here, we aimed to discover time-critical events in the tumor microenvironment that underlie ICB e cacy. We mapped dynamic changes in gene regulatory networks associated with response to ICB and used that information to rationally develop new schedule-dependent combination therapies.

Results
On/fast-off dynamics in IFN signaling are associated with response to ICB To map the dynamic processes underlying the response to ICB, we utilized our bilateral tumor model to remove responsive and non-responsive tumors at 1-hour prior and at 2-, 4-and 6-days following administration of anti-CTLA4/anti-PD-L1 therapy ( Figure 1A, B and Figure S1) and examined the transcriptomes of these tumors using RNA-sequencing. To avoid bias towards one tumor type, we utilized two different tumor models, AB1 mesothelioma and Renca renal cell carcinoma, and explored dynamic trends that were consistently differentially regulated between responders and non-responders in both models.
We rst determined whether there were differences in cellular composition between responders and nonresponders at each time point using CIBERSORTx analysis ( Figure 1C and Figure S1) 21 . Although some differences between responders and non-responders were observed, none of these differences were consistent between the two models, except for a signi cantly higher proportion of NK cells in responders prior to treatment, which we reported previously 16 . We observed an increase in CD8+ T cells after ICB, which was more prominent in responders ( Figure 1D), consistent with CD8+ T cell gene signatures as predictors of response in patient samples reported in the literature [22][23][24] .
To understand gene expression kinetics during treatment with ICB, we clustered genes based on their dynamic expression using TCSeq 25 and analyzed the resulting clusters for enriched biological pathways. We discovered four clusters that showed consistent time-dependent behavior between the two models.
Clusters 1 and 2 contained genes associated with activation of myeloid cells and T cells, including IFNg production. The expression of these genes gradually increased over time in both responders and non-responders, albeit to a greater magnitude in responders (Figure 2A and 2B), which was in agreement with the CIBERSORT results ( Figure 1D). Genes associated with cancer cell signaling (cluster 3) decreased in expression over time, but again in both responders and non-responders ( Figure 2C and S2).
In contrast, cluster 4, demonstrated a kinetic pro le that was strikingly different between responders and non-responders ( Figure 2D and S2). Cluster 4 contained genes associated with IFN signaling, which showed a gradual increase in expression over time in non-responders, while in responders it was initially highly expressed, followed by a rapid decrease in both AB1 and Renca ( Figure 2D).
To understand the transcriptional regulation of these genes, we constructed gene regulatory networks in our ICB responders and non-responders, using the GENIE3 algorithm ( Figure S3A) 26 . For each gene, the algorithm calculates an importance score re ecting the inferred effect of the gene on all other genes. A high importance score denotes a strong effect of a gene (putative regulator) on the dynamics of expression of a downstream gene (target) in the network. We ranked putative regulators by the sum of their outgoing importance scores ( Figure 3) and focused on the top 10 known transcription factors (TF) that formed central hubs in the network. For each of these TFs, we explored downstream genes that were differentially expressed between responders and non-responders and that contained a binding site for the TF within their promoter ( Figure S3B). When we visualized the most connected TFs in these networks, we identi ed multiple gene modules with increased expression early during treatment ( Figure 3A). A module containing IFN related genes exhibiting an on/fast-off dynamic response over time was found in both AB1 and Renca responders ( Figure 3B,C and Figure S3D,E). This suggested that dynamic regulation of the IFN pathway is a determinant of ICB response, with a fast reduction in expression levels over the time course of the experiment.
Next, we investigated candidate upstream regulators of the IFN module. Because functional annotations of this module suggested that both type I and type II IFN contributed to the responder phenotype, we obtained a list of interferon-stimulated genes (ISGs) for both type I and type II IFN from the Molecular Signatures Database hallmark gene sets and analyzed direct edges from TFs to these ISGs weighted by GENIE3 importance scores 27 . This analysis demonstrated that both AB1 and Renca showed similar on/fast-off kinetics of IFN signaling in responders, with on/fast-off changes by day 2 in AB1 and day 4 in Renca ( Figure 3D). In contrast, non-responders displayed a slower and less intense activation of ISGs which remained chronically active over time. We con rmed that ISG expression segregated into an early Page 6/35 and late phase, with the on/fast-off component containing genes stimulated by IFNγ, IFN α/β or a combination of both ( Figure S3F).
For both models, GENIE3 scores indicated that these ISGs were regulated by a common set of TFs, namely IRF1, STAT1, STAT2, IRF7 and IRF9, (Figure 3D), independent of tumor type. When the strength of these regulatory interactions are compared to all other interactions using hive plots 28 , these TF to ISG connections dominate ( Figure 3E, F) in both models, reinforcing the essential role of these IFN-associated TFs to the dynamic response. Taken together, these results show that on/fast-off dynamics in IFN signaling are associated with response to ICB, driven by common transcription factors across different tumor models.
Dynamic on/fast-off targeting of IFNb improves response to ICB Because there is a large overlap in ISGs that are induced by type I or type II IFNs, we were unable to resolve which IFN type was driving the response based on the transcriptomic data alone and therefore interrogated this experimentally. To phenocopy the active IFN signature in vivo, prior to ICB, we pre-treated mice with intra-tumoral injections of poly(I:C), which is known to induce both type I and II IFNs, particularly IFNβ ( Figure S4A,B) 29 . To mimic the subsequent on/fast-off-IFN signature, ICB was followed three days later by functionally blocking either type I or type II IFN signaling, using antibodies against the IFNα/β receptor (IFNAR1), IFNγ, or both ( Figure 4A). These studies were rst done in AE17 mesothelioma, which is relatively resistant to ICB 16 . Pre-treatment with poly(I:C) improved the response rate to ICB (0% vs. 26.7% complete responders), which was signi cantly enhanced by the subsequent blockade of type I IFN (53.3% complete response, p = 0.034), but not type II IFN ( Figure 4B). We con rmed these ndings in the AB1 mesothelioma model (33.3% vs 53.3% complete responders, Figure S4C). In both models, the bene cial effect of blocking type I IFN after administration was negated by blocking type II IFN simultaneously, demonstrating not only that type I IFN was responsible for the observed on/fast-off dynamics of IFN signaling in responders, but that these dynamics indeed played an important mediating role in the therapeutic response. To explore the biological relevance of these IFN dynamics, we assessed the effect of blocking IFNAR1 before rather than after ICB initiation, or concomitantly with poly(I:C) prior to ICB. This treatment completely abrogated both the response to ICB and the priming effect of poly(I:C). However if IFNAR was blocked after treatment with poly(I:C) alone, there was no detriment to the antitumor response, con rming the crucial time-dependent nature of IFN signaling underlying the therapeutic response to ICB ( Figure S4D-G).
To further dissect these on/fast-off type I IFN dynamics, we rst looked for signatures to computationally deduce the IFN subtype present in the responder tumor microenvironment. We used single cell RNAseq data from cells stimulated with a diverse array of cytokines, including IFNβ 30 , to construct a reference matrix using CIBERSORT 31 . Our deconvolution analysis revealed that genes associated with IFNβ signaling, but not other cytokines, followed the on/fast-off IFN pattern in responders in both AB1 and Renca tumor models, suggesting IFNβ was responsible for these observed dynamics ( Figure 4C, D. and Figure S5A,B). To con rm this, we treated AE17 tumor-bearing mice 3 days after administration of ICB with antibodies against either IFNβ, IFNα (subtypes A, 1, 4, 5, 11, and 13) or their shared receptor IFNAR1 32 . Mice treated with the antibody against IFNβ had a similarly increased response rate following ICB as the mice treated with the anti-IFNAR1 antibody (60% and 55.6% complete response vs. 10%), while mice treated with anti-IFNα displayed no increase in response versus controls ( Figure 4E). We repeated these experiments in the Renca model, which exhibited the same bene t of blocking IFNAR1 or IFNβ (22.2% and 40% complete response vs 0%), but not IFNγ or IFNα ( Figure 4F). We conclude that the bene cial effect of on/fast-off kinetics in type I IFN signaling after ICB is entirely dependent on switching off IFNβ.
As intra-tumoral administration of poly(I:C) is not approved for use in clinical practice in combination with ICB, we tested whether similar results could be achieved in the absence of initial IFN induction by poly(I:C). Time staggered blockade of type I IFN again improved ICB e cacy in AB1-bearing mice, which was dependent on IFNβ ( Figure 4G and Figure S5C), and we con rmed these ndings in the AE17 and Renca models ( Figure S5D,E). Although the response could be increased by blocking IFNβ after ICB, priming with poly(I:C) rst to mimic the "on" IFN signature gave optimal results. Having established that blocking IFNAR1 abrogated the effect of poly(I:C) ( Figure S4F), we determined whether IFNa or IFNb activity was driving the "on" signal in respondersusing recombinant cytokines. We found that priming with IFNβ, but not IFNα, increased the response to ICB (p=0.012), and mimicking the on/off IFN signature targeting IFNβ was superior to targeting IFNα (p= 0.027) ( Figure 4H). Notably, the temporal aspect of scheduling the respective treatments was crucial, as treatment with anti-IFNβ concomitantly with ICB did not offer any therapeutic bene t, in contrast to administration 3 days after the rst dose of ICB ( Figure 4I). These results con rm that temporal restriction of IFNβ activation underlies response to ICB.
In ammatory monocytes are the primary source of IFNb in ICB responders To further understand where the on/fast-off IFN signal was derived from, we performed single cell transcriptome sequencing one hour prior to ICB ( Figure 5A). We interrogated responder and non-responder samples using gene set enrichment analysis ( Figure 5B). We identi ed that particularly monocytes displayed elevated type I IFN signaling in responders ( Figure 5B) In particular, a speci c monocyte subpopulation drove the on/fast-off-IFN gene signature found in our bulk RNAseq data ( Figure 5B,C). These monocytes (cluster 1, Figure 5D) displayed high Ly6c expression ( Figure 5E), consistent with an in ammatory monocytic phenotype 33 . Repeating gene regulatory network inference on this single cell data con rmed that IRF1, IRF7, STAT1, IRF9 and STAT2 were major transcription factors driving the response in these cells ( Figure 5F, Figure S6A), supporting our network analysis results in bulk samples.
Additionally, we con rmed CD11b + /Ly6C hi monocytes attributed the highest expression of Irf1, a key ISG regulator, in tumor samples by ow cytometry ( Figure 5G).
We performed RNA velocity analysis and observed a trajectory from cluster 1 to cluster 2 monocytes characterized by a diminished activation of ISGs ( Figure 5H). Individual ISGs velocities along this trajectory are markedly diminished in non-responders compared to responders ( Figure S6B). We therefore examined gene velocities for IFN-related TFs and ISGs in the on/fast-off gene signature. Hierarchical clustering showed separation of responder versus non responder monocytes based on ISG velocities ( Figure S6B). Along this trajectory, velocity analysis showed that monocytes downregulated transcription of ISGs such as Irf1, and this was more pronounced in responders than non-responders ( Figure S6C). Both these results are consistent with the on/fast-off IFN dynamic we observed in bulk RNAseq data ( Figure S6D-F). Additionally, we saw a change in expression from Ly6c hi to Ly6c lo monocytes along this trajectory ( Figure S7A,B) and these changes at single cell level were compatible with Ly6c expression over time observed in the bulk RNAseq data ( Figure S7C,D). As Ly6C1/2 is a marker of blood derived in ammatory monocytes 33 , this suggests that these transcriptionally dynamic cells are likely in ltrating the tumor and differentiating. To further support this notion, we determined the expression dynamics of the blood-derived monocytic marker CCR2 34 , and found it followed the same kinetics as the Ly6c and on/fast-off gene expression patterns ( g S7E,F).
To con rm whether monocytes were indeed the source of IFNβ in the tumor microenvironment, we used B6.129-Ifnb1 tm1Lky /J mice, which co-express IFNβ1 and eYFP from the Ifnb1 locus, bearing AE17 tumors 35 . Tumors were analyzed by ow cytometry, one day after intra-tumoral poly(I:C), revealing YFPpositive cells were all CD45 + , CD11b + , MHC-II -, F4/80 -, CD11cwith the majority expressing Ly6C, consistent with the single cell RNAseq results ( Figure 5I). Together, these results pinpoint CD11b + monocytes as the key IFNβ producing cells in the tumor microenvironment.
On/fast-off IFN kinetics in monocytic cells correlate with response in patients treated with ICB To validate our preclinical ndings, we tested whether activation of IFN signaling also occurred in patients who responded to ICB therapy, using single cell RNAseq data from a cohort of treatment naïve breast cancer patients, subsequently treated with anti-PD-1 and with paired T cell expansion data 24 .
In this dataset, biopsies that were taken prior to treatment with ICB showed global type I/II IFN-responsive genes positively correlated with T cell expansion after treatment. In our analysis, we interrogated expander (responder) and non-expander (non-responder) samples pre-treatment using gene set enrichment analysis which con rmed a globally elevated on/fast-off IFN signature in expanders, which was highest in myeloid cells ( Figure 5J and Figure S8). Within the myeloid cell cluster, we found the Ccr2 + subpopulation in expanders highly expressing the on/fast-off-IFN gene signature, consistent with our ndings in the murine models of tumour-in ltrating in ammatory monocytes being the predominant source of type I IFN activity and IFNb in particular ( Figure 5I, K and Figure S8).

Discussion
Here, we report that on/fast-off activation of IFNb in cancer is required for the therapeutic response to ICB, we dissect the role of IFNβ versus IFNα in the anti-tumor immune response, and we provide a rst example of time-dependent activation and inhibition of a drug target being required to achieve optimal anti-cancer effect.
Clinical studies have shown that an IFN gene signature is associated with treatment with ICB, and that an IFN mediated signature can predict response 12,18,36,37 . However, the dynamic nature of an orchestrated innate and adaptive immune response against cancer cannot be adequately interpreted from a single snapshot 5,8 . Ideally, tumor biopsies would be frequently obtained early during ICB treatment to identify underlying mechanisms of response, but perhaps with the exception of some bone marrow cancers, this is usually not feasible due to anatomical tumor location and the requirement for invasive procedures. In addition, given the presence of intra-patient tumor heterogeneity, repeat biopsies may not provide a consistent representation of the tumor microenvironment 38 . Moreover, inter-patient heterogeneity makes it di cult to identify small yet meaningful biological differences underlying the response to ICB. Preclinical studies using cell line-derived tumors in syngeneic mice can negate some of these issues, including patient and tumor heterogeneity. Typically, comparisons are made between responsive and nonresponsive tumor mouse models, or between untreated and treated mice 39 , sometimes with suboptimal dosing as not to destroy the biological read-out of the perturbed tumor microenvironment 40 . We further re ned this approach by comparing responders and non-responders to ICB within the same tumor model, sampling entire tumors over time 6 . Speci cally, the bilateral models that we used are fully internally controlled; responders and non-responders are equally exposed to the relevant experimental variables, including potential in ammation associated with cancer cell inoculation, tumor size, surgical tumor removal (with sham surgery having no effect on symmetry 6,15 ) and anesthesia, in addition to the inherent genetic and environmental homogeneity of using inbred mice. Yet, mice display a stark dichotomy in response, within the models. This difference in outcome is likely due to the stochastic nature of the induction of an effective anti-cancer immune response, which contains many switches, thresholds and feedforward and feedback loops 41 . Importantly, these models allowed us to identify dynamic immune response-intrinsic changes that determine the therapeutic outcome in a highly homogenous genetic and environmental background.
Our results show that the on/fast-off dynamics of IFNβ signaling are crucial to the response to ICB, which can be therapeutically exploited using antibodies against IFNβ or its receptor IFNAR1, resulting in enhanced tumor clearing. Others have shown that secondary ICB-resistant cancers, which are cancers that initially responded but then relapsed, display chronically active IFN signaling 42,43 . This has resulted in the suggestion to co-treat patients with JAK inhibitors, which block both type I and II IFN signaling 42 . We extend these ndings by demonstrating in intrinsically responsive tumors that type I IFN only, and more speci cally IFNβ only, plays a dual role and that the response rate and depth of response can be improved by therapeutically mimicking these on/fast-off dynamics. We did not observe these effects when we blocked IFNα using an antibody that is speci c for 5 subtypes 32 . We cannot exclude that other IFNα subtypes have an equally inhibitory effect on the therapeutic response.
Validation of these preclinical ndings in patient samples is di cult, given the aforementioned issues with obtaining repeated samples in patients. However, using our bilateral mouse models, we previously identi ed a gene expression signature predictive for response 15 , which was extensively and independently validated by others, by directly comparing the mouse data with multiple patient cohorts across different tumor types and different ICB antibodies [18][19][20] . In addition, we con rmed the presence of known indicators of clinical e cacy following ICB, such as increased CD8 T cell in ltration 23 . Lastly, several other groups have now independently demonstrated that bilateral models identify important and translationally relevant biology in the context of ICB 17,44 .
In order to further validate our preclinical ndings, we used single cell RNAseq data from a cohort of breast cancer patients treated with anti-PD-1 45 , which con rmed that in ammatory tumor-in ltrating monocytes were mainly responsible for the on/fast-off IFN signal in these patients, and that enrichment of these cells was correlated with T cell expansion. Whether the downregulation of IFNb is involved in driving the T cell response, or whether chronic IFNb signaling provides cancer cells with a survival advantage under immune pressure remains to be established. Interestingly, in that context, it has been reported that cancer cells that are unable to downregulate the IFNAR1 receptor during melanoma development respond better to ICB compared to normal melanoma cells 46 . And although the validation of the results in clinical samples provides con dence that the therapeutic targeting of IFNb in a timedependent manner could improve the response to ICB in patients, future clinical trials are needed to assess IFNβ as a biomarker and dynamic drug target to determine the optimal time to start blocking after ICB.
The nding that time-dependent aspects of type I IFN signaling contribute to a powerful anti-tumor immune response resonates with ndings from viral immunology, where type I IFN is required for acute clearance of viral infections such as hepatitis C, and it is even used therapeutically in that context 47 . Yet, paradoxically, blocking IFNAR1 can be bene cial for the control of chronic viral infections, as has been shown for LCMV or coronavirus infections 32,[48][49][50] . We propose that the anti-tumor response following ICB mimics aspects of the acute or chronic anti-viral immune response, resulting in either swift regression or non-response, respectively. This notion is also in line with recent ndings showing in ammatory monocytes as major responders to type I IFN production in viral immunity, as well as our results and those from others in the context of ICB in cancer 51,52 .
The interaction between the immune system and cancer cells is often conceptualized as a cycle, which can be pushed at any level, at any given time, to induce the appropriate momentum 53 Instead, our data suggest a continuous changing landscape of the immune response, where interventions have a timedependent effect, even to the point that the exact same target must be modi ed in a diametrically opposite manner; by providing excess recombinant IFNb rst, followed by blocking antibodies against IFNb later. In oncological treatments in general, drugs, including in combination, are typically administered empirically until they are no longer effective, or toxicities preclude continuation. Our results challenge this approach and could have important implications for drug discovery research, demonstrating that in order to obtain optimal clinical effect, some targets need to be therapeutically modulated in a time-dependent, bidirectional manner. As recombinant IFNb has been FDA approved for multiple sclerosis and antibodies targeting the IFNβ/IFNAR1 pathway have been fully developed in the context of autoimmunity 54 , these results can be readily translated into the clinic. Tumor preparation for RNA sequencing Whole tumors and lymph nodes were surgically resected, the surrounding tissue was removed and immediately submerged in RNAlater (Life Technologies, Australia). Samples were stored at 4°C for 24 hours, after which supernatant was removed and samples transferred to -80°C. Frozen tumors were dissociated in Trizol (Life Technologies, Australia) using a TissueRuptor (QIAgen, Australia). RNA was extracted using chloroform and puri ed on RNeasy MinElute columns (QIAgen, Australia). RNA integrity was con rmed on the Bioanalyzer (Agilent Technologies, USA). Library preparation and sequencing (50 bp, single-end) was performed by Australian Genome Research Facility, using Illumina HiSeq standard protocols.
Alignment and differential expression.
We processed a total of 144 RNA-seq single-end read samples across four time points in two mouse models. After reviewing quality control on all samples using FastQC software, we used Kallisto 56 (v0.43.0) for transcript abundance estimation. Following alignment, we performed differential expression analysis with Sleuth (v0.29.0) 57 We compared responders and non-responders using a model containing time-point and response as covariates using a likelihood ratio test. We aggregated p-values from transcript differential expression to gene-level results with transcript-to-gene mapping relying on the latest Gencode reference M25 (GRCm38.p6) using Lancaster's method 58 . Genes were deemed differentially expressed at a false-discovery rate of less than 5%, regardless of fold change (Data S1).

RNAseq analysis of dynamic gene expression data
We used a deconvolution approach to deduce the cell subtypes present in the responder and non- We clustered time course RNAseq data using the fuzzy c-means (FCM) clustering algorithm Mfuzz 59 in the TCseq package 25 . Z-normalised/scaled counts were used in the algorithm and expression pro les were grouped clusters (k = 6) based on their dynamic patterns. We used a Pearson correlation score on trend lines to compare trends of each cluster between responders and non-responders, to identify which patterns were unique to responders. Each matching cluster between the two models had overlapping genes extracted (Data S3) and enrichment of per-cluster genes was performed using Enrichr 60 .
To acertain transcription factors involved in these processes, we constructed two networks, one for Renca responders and one for AB1 responders. We used the GENIE3 algorithm 26 , which achieved the best performance on the DREAM5 network inference challenge 61 . To construct each gene regulatory network, we used 36 responder samples across all time points as input to the GENIE3 algorithm. Since GENIE3 requires gene counts as input, we summarized transcript abundances derived from Kallisto as gene counts using the Bioconductor tximport 62 package. We ran GENIE3 (v1.8.0) with default parameters (treeMethod ="RF", k = "sqrt", nTrees = 1000).
We pruned GENIE3 output isolate biologically relevant edges, incorporating DE information and transcription factor binding site predictions from JASPAR 63 (Fig. S2). We identi ed and retained "direct connections", de ned as connections between a transcription factor (TF) to differentially expressed genes only if the TFBS for the TF was situated in a genomic window 400 base pairs upstream or 300 base pairs downstream of a DE gene's transcription start site (TSS). TSS sites for differentially expressed genes were obtained from the UCSC genome browser 64 . The BEDtools 65 "window" function was used to obtain all direct regulatory interactions and GENIE3 importance scores appended to these interactions using the R data.

Network visualization and analysis of IFN-stimulated genes
We visualized the induced subgraph of the top ten nodes and their rst order neighbors in both AB1 and Renca using the R igraph package with the Kamada-Kawai 66 layout. We assigned color to genes in the network by their average gene expression across time, normalized by Z-score. From inspection of these graph visualizations, we con rmed in both AB1 and Renca that these subnetworks were enriched for interferon related transcription factors (Irf1, Stat1, Stat2, Irf7, Irf9) and their direct downstream targets (Fig. S3). Pathway analysis on these hubs using Enrichr showed statistically signi cant enrichment of terms relating to both Type I and Type II interferon signaling (Fig. S3), suggesting that dynamic changes in these two pathways were crucial to the checkpoint blockade response.
To identify ISGs involved in Type I and Type II interferon signaling in these direct networks, we used gene sets of interferon alpha/beta signaling and interferon gamma signaling from the Molecular Signature Database (MSigDB) 27,67 converted into their murine equivalents using BioMart 68 .After conversion, we retained 88/97 genes from the alpha dataset and 186/200 genes from the gamma dataset (Fig. S4) We visualized TF-to-ISG GENIE3 scores and expression pro les for ISGs common to both AB1 and Renca networks using the R ComplexHeatmap 69 package. From each direct regulatory network, we extracted the TF-to-ISG matrices, setting any non-existent interactions in the matrix were set to zero prior to visualization. We observed regulation of on/fast-off ISGs con ned to just 5 key interferon-related transcription factors -Irf1, Stat1, Stat2, Irf7 and Irf9, with minimal regulatory impact from other transcription factors in these networks ( Fig. S4 and Fig. 2).

Hive plot generation
To overcome visual biases from traditional network layouts, we used a hive plot visualization 28 . For each network, we used a 4-axis hive plot in the HiveR package, allowing us to partition edges to visualize interferon-related signaling in our networks. The following axes were used, color coded in the following way (Fig. 2D, 2E): Red axis -interferon related TFs (IRF1, STAT2, STAT1, IRF7 and IRF9); Green axisnon-interferon related important TFs, comprising the union set of top TFs from both networks in Fig. S3; Purple -Downstream gene targets in the "on/fast-off ISG set"; Blue -Other genes in the direct networks not in 1, 2, or 3. For axis 3, the fast-ISG set was derived from k-means (k = 2) clustering on time course expression data from the AB1 responder data set. This node topology allows us to more easily visualize the "quadrant" of graph edges from IFN-related TFs to on/fast-off-ISGs, demonstrating that this quadrant contained edges with high value GENIE3 scores (above 0.9 quantile) denoting important dynamic regulatory links in the network.
Single cell sample pre-processing AB1 and renca tumors were surgically removed 1 hour prior to ICB administration and immediately submerged in cold PBS, cut into 1-2 mm pieces with a scalpel blade and dissociated using the GentleMACS system (Miltenyi). Cell suspensions were frozen in RPMI medium containing 50% FCS and 10% DMSO. Cryopreserved single cell suspensions were rapidly thawed in a 37°C water bath and prepared for single cell library construction as previously described (Denisenko et al., 2020). Libraries were constructed using the 10X Chromium 3' work ow (version 2 chemistry) as per the manufacturer's directions. We aimed to capture 9000 cells per sample. Libraries were quanti ed using the TapeStation D1000 kit (Agilent). Sequencing was performed by Novogene, using NovaSeq S2 owcell sequencing protocols.
For single cell analysis, we processed FASTQ les from 6 AB1 and 6 Renca samples using cellranger v3.0 (10X genomics). For each sample, we performed demultiplexing and read alignment using the cellranger count function, using cellranger's pre-supplied mm10 reference with an expect-cells parameter of 6000.

Clustering, visualization and cell annotation
We used the Seurat 70 (version 3.14) R package to combine samples for downstream analysis. Gene counts were normalized against both sequencing depth and also against the percentage of mitochondrial DNA in each cell using negative binomial regression. The resulting Pearson residuals from these processing steps were used for downstream PCA, cluster identi cation and UMAP embedding and visualization.
To avoid subjective biases in cell identi cation, we used an automated labelling strategy based on bulk RNAseq references. The R package SingleR 71 was used in "cluster mode" using species-speci c annotation references provided with the package. For annotation of human single cell data, we used the human primary cell atlas reference 72 . For annotation of murine data, we used the mouse RNAseq dataset from Benayoun et al 73 .
Clusters were de ned from Seurat's FindClusters function at default (0.8) resolution. Similarly, labelled clusters were merged. We con rmed that this approach was robust to cluster size by showing that labels were consistent even when cluster size was modi ed by changing resolution parameter in the FindClusters function. We performed annotation diagnostics by checking cell cluster identities in our AB1 and Renca samples against the ImmGen 74 reference. We found both references to be in agreement .

Single cell differential expression analysis and conserved marker analysis
We performed differential expression using the FindMarkers function in Seurat. Genes were deemed differentially expressed at an absolute log-fold change of 0.5 and a q-value of below 0.05 using the nonparametric Wilcoxon test.. In Renca, we observed far fewer differentially expressed genes between responders and non-responders at Day 0 across all cell types, consistent with our bulk RNAseq data (Data S1).

Label transfer and interspecies integration of single cell data samples
We used Seurat's "label transfer" functions, allowing cell identities from a reference sample to be projected onto our target dataset to identify tumor cells and to compare monocytes from murine and human single cell data. To construct a reference sample of mesothelioma cells, we sequenced samples from AB1 mesothelioma tumors in which the tumor cells were tagged with in uenza hemagglutinin 75 .
The cellranger index was rebuilt to incorporate the tag sequence. After alignment, any cell containing the HA-tag was labeled as a tumor cell and used as a reference for label transfer. Clusters containing more than 10% tumor cells were deemed to be putative tumor clusters.
To label Renca tumor cells, we used a reference dataset of mouse kidney single cells. 76 in which the authors identi ed gene markers mapping to anatomical elements of the mouse renal tubular system, since these Renca tumors recapitulate renal tubule elements. We selected clusters which expressed the highest average expression of markers speci c proximal and distal tubules (Lrp2, Slc27a2). Cells in these reference clusters had their identities projected onto the Renca dataset.

Copy number variation analysis in tumor cells
Tumor cells usually display evidence for somatic large-scale chromosomal copy number alterations, such as gains or deletions of entire chromosomes or large segments of chromosomes. We used the inferCNV R package 77 to check tumor cell identity after label transfer. In both AB1 and Renca, tumor clusters labeled by Seurat's projection strategy were in good agreement with clusters of cells deemed to be tumor clusters based on the existence of copy number changes inferred by inferCNV.
Gene set enrichment analysis on single cell data We used the SCDE/pagoda 78 package, which detects statistically signi cant coordinated variability at single cell level. Brie y, from the original 38K cells in our AB1 samples, we constructed KNN error models for the 17K surviving default SCDE's library size lters. For gene sets, we tested GO terms for interferon production and interferon response extracted from the org.Mm.eg.db Bioconductor package and also our "custom" gene set of fast-ISGs derived from bulk expression data. We visualized enrichment scores using python's Seaborn scatter plot with a color scale mapped to enrichment score intensity.
Single cell analysis on the validation dataset of human breast cancer patients 24 was conducted using the escape 79 R package (v1.3.1), which provides convenience functions for the GVSA algorithm 80 .
Random subsampling (without replacement) was performed to create input batches of 20,000 cells to the GVSA algorithm and the results were pooled for visualization with the R ggplot2 package 81 . Filtering, normalization and scaling was performed as outlined in the original paper 24 .Cell clusters were labelled using the labels assigned by the authors.

SCENIC network analysis
We performed network inference to analyze our single cell data using SCENIC 82 . We summarized the result of these analyses using a difference heatmap between responders and non-responders of average binarized TF activity per Seurat cluster. Speci cally, the regulon binarization scores were averaged across clusters, separated by response and the difference between responder and non-responder averages were visualized using the R pheatmap (v1.0.12). To visualize these results, we created a heatmap showing the difference in percentage of TF activation per cluster for the IRF related TF (Irf1, Stat1, Stat2, Irf7 and Irf9) (Fig. 3F). The full heat map of differential transcription factor activation in AB1 monocytes is displayed in Fig. S6.

Velocity analysis
RNA velocity quanti es cell transcriptional activity by modelling the time derivative of gene expression states. We used the velocyto 83 package for this analysis. Per-sample loom les, containing quanti ed spliced vs. unspliced transcripts, were combined using the python loompy (v3.0.6) package. After count normalization, ltering and feature selection on these genes, approximately 2.5K genes survived these ltering steps to be used for gene velocity modelling.
We compared "transcriptional momentum" of AB1 monocytes, which we de ned as the squared L2 norm for each cell's embedding vectors with respect to their UMAP 84 coordinates. We compared KDE distributions for momentum in various AB1 clusters, separated by response (Fig. S6). To compare Renca and AB1 monocytes, we repeated our velocity analysis on using an embedding of all 12 samples from both AB1 and Renca single cell data. On this common embedding, momentum calculations show that AB1 responders have a higher transcriptional velocity than Renca monocytes in responders (Fig. S6).
As an additional check that differences in transcriptional momentum were, in part, due to differences in interferon signaling, we examined gene velocities for IFN-related TFs and ISGs in the on/fast-off component gene set which survived the above-mentioned ltering. We extracted the normalized velocities of 42 of these genes which survived data preprocessing. Hierarchical clustering showed separation of responder versus non responder monocytes based on ISG velocities (Fig. S6).

Cytokine stimulation estimation
We aimed to determine whether there was evidence of an IFNb-induced signature, irrespective of cell type involved. We used the only available dataset which compared known cell differentiation cytokines under co-stimulatory conditions including IFNb stimulation. We used a deconvolution approach to deduce the IFN subtype present in the responder tumor microenvironment to share weighting of genes conserved across the different stimulatory conditions. The CIBERSORT algorithm 31 was used to estimate the relative proportions of 7 cytokine induced T cell signatures based on the transcriptomic pro les of each sample, where the induced T cell gene signature developed by Cano-Gamez et al. 30 was used as a reference (Fig. S5)   Responders and non-responders have similar reactions to ICB, but responders shut off genes related to IFNβ signaling. TCseq analysis was used to cluster genes with similar expression over time, identifying clusters shared by both AB1 and Renca. Gene expression over time and pathway analysis of overlapping genes between AB1 and Renca was per-formed for cluster 1 (A), cluster 2 (B), cluster 3 (C) and cluster 4 (D). for AB1 (E) and Renca (F) with TF-to-ISG edges situated in the right upper quadrant. The top 10% of (high valued) edges by GENIE3 score are highlighted in red.

Figure 4
Targeting IFNβ in a directionally opposite, time-dependent manner improves the re-sponse to ICB. (A) Treatment strategy. (B) Survival curves of AE17-bearing mice treated with poly(I:C), ICB followed by antibodies blocking type I and/or II IFN (n = 15 per group). (C and D) Deconvolution analysis on bulk RNA seq data from AB1 (C) and Renca (D) using an IFNβ-stimulated T cell signature between responders (red) and non-responders (blue). Bars represent standard deviation. *p ≤ 0.05, **p ≤ 0.01, ****p < 0.0001 from two-way ANOVA with Tukey's multiple comparisons test. (E) Survival curves of AE17 bearing mice treated with poly(I:C), ICB followed by antibodies against the type I IFNs, IFNα or IFNβ compared to blocking their receptor IFNAR (n = 5 per group). (F), Survival curves of Renca bearing mice treated with poly(I:C), ICB followed by antibodies against the type I and/or II IFN (n = 5 or 10 per group). (G) Survival curves of AB1 bearing mice treated with ICB followed by antibodies against type I and/or II IFN (n = 5 to 15 per group).