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

doi:10.21203/rs.3.rs-892399/v1

Download PDF

Article

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

https://doi.org/10.21203/rs.3.rs-892399/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 19 Aug, 2022

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

Little is known about the dynamic biological events that underpin therapeutic efficacy in immune checkpoint blockade (ICB) in cancer, due to the inability to frequently sample tumors in patients. Here, we mapped the transcriptional profiles 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 efficacy, but only when IFNβ was targeted, not IFNα. We identified Ly6C+/CD11b+ inflammatory monocytes as the primary source of IFNβ and found that active type-I-IFN signaling in tumor-infiltrating inflammatory 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.

Immunology

Computational Biology

Bioinformatics

Immune checkpoints

PD-1

CTLA4

type I interferon

IFNβ

gene regulatory networks

time-dependent

RNAseq

biomarker

drug discovery

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 specific aspects of the pre-treatment tumor microenvironment, such as PD-L1 expression (Reck et al., 2016), immune cell infiltration ³ or tumor mutation burden ⁴ have been shown to correlate with response, none of these biomarkers are sufficiently robust to guide clinical decisions regarding treatment across cancers ⁵ nor has their characterization yet resulted in approved treatments that improve the efficacy of ICB ⁶. As a consequence, the development of novel combination therapies to improve outcomes is mainly empiric ⁷.

When perturbing complex systems, some important effects may only become apparent over time. For example, in the course of an immune response, some inflammatory mediators need to be switched off, not only to resolve inflammation, but also to mediate a transition from innate to adaptive immunity ⁸. Similar time-dependent mechanisms may underpin an effective anti-tumor immune response ⁵. However, it is not possible to identify these dynamic events in human studies due to the difficulty 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 difficult to identify discrete biological differences between responding and non-responding patient tumors due to inter-individual variation in germline and cancer genetics, tumor microenvironment composition and environmental influences ^11–14.

For these reasons, we optimized murine models with bilateral tumors, derived from syngeneic cancer cell lines ⁶. 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 findings 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–20. Here, we aimed to discover time-critical events in the tumor microenvironment that underlie ICB efficacy. 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.

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 first determined whether there were differences in cellular composition between responders and non-responders at each time point using CIBERSORTx analysis (Figure 1C and Figure S1) ²¹. Although some differences between responders and non-responders were observed, none of these differences were consistent between the two models, except for a significantly higher proportion of NK cells in responders prior to treatment, which we reported previously ¹⁶. 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–24.

To understand gene expression kinetics during treatment with ICB, we clustered genes based on their dynamic expression using TCSeq ²⁵ 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 profile 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) ²⁶. For each gene, the algorithm calculates an importance score reflecting 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 identified 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 ²⁷. 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 confirmed that ISG expression segregated into an early 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 ²⁸, 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) ²⁹. 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 first done in AE17 mesothelioma, which is relatively resistant to ICB ¹⁶. Pre-treatment with poly(I:C) improved the response rate to ICB (0% vs. 26.7% complete responders), which was significantly enhanced by the subsequent blockade of type I IFN (53.3% complete response, p = 0.034), but not type II IFN (Figure 4B). We confirmed these findings in the AB1 mesothelioma model (33.3% vs 53.3% complete responders, Figure S4C). In both models, the beneficial 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 anti-tumor response, confirming 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 first 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β ³⁰, to construct a reference matrix using CIBERSORT ³¹. 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 confirm 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 ³². 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 benefit of blocking IFNAR1 or IFNβ (22.2% and 40% complete response vs 0%), but not IFNγ or IFNα (Figure 4F). We conclude that the beneficial 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 efficacy in AB1-bearing mice, which was dependent on IFNβ (Figure 4G and Figure S5C), and we confirmed these findings 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) first 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 benefit, in contrast to administration 3 days after the first dose of ICB (Figure 4I). These results confirm that temporal restriction of IFNβ activation underlies response to ICB.

Inflammatory 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 identified that particularly monocytes displayed elevated type I IFN signaling in responders (Figure 5B) In particular, a specific monocyte sub-population 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 inflammatory monocytic phenotype ³³. Repeating gene regulatory network inference on this single cell data confirmed 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 confirmed CD11b⁺/Ly6C^hi monocytes attributed the highest expression of Irf1, a key ISG regulator, in tumor samples by flow 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 inflammatory monocytes ³³, this suggests that these transcriptionally dynamic cells are likely infiltrating the tumor and differentiating. To further support this notion, we determined the expression dynamics of the blood-derived monocytic marker CCR2 ³⁴, and found it followed the same kinetics as the Ly6c and on/fast-off gene expression patterns (fig S7E,F).

To confirm 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 ³⁵. Tumors were analyzed by flow cytometry, one day after intra-tumoral poly(I:C), revealing YFP-positive cells were all CD45⁺, CD11b⁺, MHC-II^-, F4/80^-, CD11c^- with 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 findings, 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 ²⁴. 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 confirmed 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 findings in the murine models of tumour-infiltrating inflammatory monocytes being the predominant source of type I IFN activity and IFNb in particular (Figure 5I, K and Figure S8).

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 first 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 ³⁸. Moreover, inter-patient heterogeneity makes it difficult 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 non-responsive tumor mouse models, or between untreated and treated mice ³⁹, sometimes with suboptimal dosing as not to destroy the biological read-out of the perturbed tumor microenvironment ⁴⁰. We further refined this approach by comparing responders and non-responders to ICB within the same tumor model, sampling entire tumors over time ⁶. Specifically, the bilateral models that we used are fully internally controlled; responders and non-responders are equally exposed to the relevant experimental variables, including potential inflammation 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 ⁴¹. 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 ⁴². We extend these findings by demonstrating in intrinsically responsive tumors that type I IFN only, and more specifically 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 specific for 5 subtypes ³². We cannot exclude that other IFNα subtypes have an equally inhibitory effect on the therapeutic response.

Validation of these preclinical findings in patient samples is difficult, given the aforementioned issues with obtaining repeated samples in patients. However, using our bilateral mouse models, we previously identified a gene expression signature predictive for response ¹⁵, 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–20. In addition, we confirmed the presence of known indicators of clinical efficacy following ICB, such as increased CD8 T cell infiltration ²³. 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 findings, we used single cell RNAseq data from a cohort of breast cancer patients treated with anti-PD-1 ⁴⁵, which confirmed that inflammatory tumor-infiltrating 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 ⁴⁶. And although the validation of the results in clinical samples provides confidence that the therapeutic targeting of IFNb in a time-dependent 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 finding that time-dependent aspects of type I IFN signaling contribute to a powerful anti-tumor immune response resonates with findings 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 ⁴⁷. Yet, paradoxically, blocking IFNAR1 can be beneficial for the control of chronic viral infections, as has been shown for LCMV or coronavirus infections ^32,48–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 findings showing inflammatory 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 ⁵³ Instead, our data suggest a continuous changing landscape of the immune response, where interventions have a time-dependent effect, even to the point that the exact same target must be modified in a diametrically opposite manner; by providing excess recombinant IFNb first, 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 ⁵⁴, these results can be readily translated into the clinic.

Acknowledgments

We thank Michael Small, Ayham Zaitouny, Emma de Jong and Leo Portes dos Santos for discussions and Ian Dick for statistical support. This work was funded by NHMRC grants 1103980 and 1107043. W.L.C. was supported by a scholarship from the ADSA and a WA DoH Cancer and Palliative Care Network Fellowship. R.M.Z. was supported by a Forrest Prospect Fellowship. T.L. was supported by a Feilman Foundation fellowship. W.J.L. was supported by a Simon Lee fellowship, an NHMRC fellowship and a Cancer Council WA fellowship. A.R.R.F. is supported by an NHMRC fellowship, APP1154524. Part of this research was made possible by the Australian Cancer Research Foundation Centre for Advanced Cancer Genomics, a collaborative cancer research grant provided by the Cancer Research Trust and an enabling grant from Cancer Council WA.

Author Contributions

W.L.C, R.M.Z, T.L and W.J.L designed research; R.M.Z., V.S.F., B.W., T.H.C., C.F., C.T., B.G. performed experimental work; W.L.C, R.M.Z, A.B., A.R.R.F., M.J.M., A.K.N., R.A.L., T.L and W.J.L analyzed the data; L.B. contributed new reagents; R.M.Z, W.L.C., T.L. and W.J.L. wrote the manuscript. T.L and W.J.L supervised and managed the overall study.

Declaration of Interests

Patent application pertaining to aspects of this work. W.J.L.: consultancy Douglas Pharmaceuticals, MSD, research funding Douglas Pharmaceuticals, AstraZeneca, ENA therapeutics. A.K.N. advisory boards Boehringer Ingelheim, Bayer, Roche, BMS; research funding from AstraZeneca. M.M.: Advisory boards MSD, BMS, Roche, AstraZeneca.

Data and materials availability

RNA-Seq data for this analysis are available in the Gene Expression Omnibus under accession number [GEO: GSE153943] for the bulk RNAseq data, and [GEO: GSE153942] for the single cell RNAseq data. Code for the analysis in the main manuscript and the supplementary data is available through GitHub (https://github.com/wlchin/IFNsignalling).

Mice

BALB/cArc, Balb/cAusB, C57BL6/J or Ifnb1^tm1Lky/J mice 8-12 weeks of age were used for all experiments. BALB/cArc or C57BL6/J mice were obtained from the Animal Resource Centre (Murdoch, WA), Balb/cAusB mice were obtained from the Harry Perkins Institute for Medical Research Bioresources Centre South (Murdoch, WA). Ifnb1^tm1Lky/J mice, generated by knock-in of a yellow fluorescent protein (YFP) reporter cassette into the endogenous Ifnb locus ³⁵, were imported from The Jackson Laboratory (Bar Harbour, Maine) and maintained at the Harry Perkins Institute for Medical Research Bioresources Centre South (Murdoch, WA). All mice were housed at the Harry Perkins Institute of Medical Research Bioresources Facility North under specific pathogen free conditions. Mice were fed Rat and Mouse cubes (Specialty Feeds, Glen Forrest, Australia) and had access to water ad libitum. Cages (Techniplast, Italy) were individually ventilated with filtered air, contained aspen chips bedding (Tapvei, Estonia) and were supplemented with tissues, cardboard rolls and wood blocks as environmental enrichment, and were changed every 14 days. Mice were housed at 21-22°C with 12-hour light/dark cycle (06:00 – 18:00). Sentinel mice (n = 3) in the animal facility were screened monthly for a standard panel of bacteria and fungi, ectoparasites, endoparasites, non-pathogenic protozoa and viruses (Cerberus Sciences, Australia). All experiments were conducted in compliance with the institutional guidelines provided by the Harry Perkins Institute for Medical Research animal ethics committee (approval numbers AE047, AE091, AE157, AE201).

Cell culture

Cell lines AB1 and AE17 were obtained from CellBank Australia. Cell line Renca was kindly donated by Dr E. Sotomayor and Dr F. Cheng (University of South Florida, Tampa, FL). Cell lines were maintained in RPMI 1640 supplemented with 20 mM HEPES, 0.05 mM 2-mercaptoethanol, 100 units/ml penicillin/streptomycin (Thermo Fisher), and 10% FCS (Invitrogen, Mulgrave, Australia). Cells were grown to 70-80% before passage and passaged 3-5 times before inoculation. Cells were frequently tested for mycoplasma by PCR and remained negative. Cell lines were validated yearly by flow cytometry for MHC class I molecules H2-Kb (consistent with C57BL/6) and H2-Kd (consistent with BALB/c), and for fibroblast markers E-cad, EpCam and PDGFR⍺ (negative) and by PCR for mesothelin (positive for AB1, negative for Renca).

In vivo treatments

When cell lines were 70-80% confluent, they were harvested and washed 3 times in PBS. 5 x 10⁵ cells in 100 μl were inoculated subcutaneously (s.c.) onto the lower right-hand side (RHS) flank (for single inoculations) or both flanks (for dual tumor inoculations) using a single 26G needle per injection. Mice were randomized when tumors became palpable, approximately 3-5 days after tumor inoculation.

Surgery experiments

A detailed protocol for the surgery experiments has been previously described ¹⁶. For the pre-treatment samples, tumors were resected eight (AB1) or 10 (Renca) days post tumor inoculation, when tumors were ~ 9 mm², and mice were administered ICB 1 hour after surgery. For the post-treatment samples, tumors were resected 2, 4 or 6 days after the first administration of ICB. Mice were dosed with 0.1 mg/kg buprenorphine in 100 μl s.c. (30 min prior) and anesthetized using isoflurane (4% in 100% oxygen at a flow rate of 2 L/min). Whole tumors and the corresponding draining inguinal lymph node on the right-hand side were removed by surgical excision and immediately immersed in RNAlater (Life Technologies, Australia). The wound was closed with staples (Able Scientific, Australia). Mice were placed in a heat box for recovery. The remaining tumor was monitored for response as an indicator of response for the removed tumor. Mice were designated as responders when their tumor completely regressed, and they remained tumor free for up to 4 weeks after treatment. Mice were designated as non-responders if their tumors grew to 100 mm² within 4 weeks after start of treatment, similar to saline-treated controls (Fig. S1). Mice that had a delay in tumor growth or partial regression were designated as intermediate responders and excluded from the analysis. For internal consistency, we only used experiments in which mice displayed a dichotomous response, i.e. in any cage there had to be at least one non-responder amongst responders or vice versa.

In vivo ICB treatment

The anti-PD-L1 hybridoma (clone MIH5) and the anti-CLTA4 hybridoma (clone 9H10) were cultured in IMDM containing 1% of FCS and gentamycin at Bioceros (Utrecht, The Netherlands). Clarified supernatants were used to purify the antibody using affinity chromatography. The antibodies were sterile formulated in PBS. Alternatively, antibodies from the same clones were obtained from BioXcell (New Hampshire, US). Mice received an intraperitoneal (i.p.) dose of 100 μg of anti-CTLA4 and 100 μg anti-PDL1 combined in 100 μl phosphate-buffered solution (PBS). Mice received additional doses of 100 μg anti-PDL1 two and four days later. In previous experiments ⁵⁵, we had not found any difference in effect of control IgG versus PBS, and therefore vehicle controls received PBS alone.

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 purified on RNeasy MinElute columns (QIAgen, Australia). RNA integrity was confirmed 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 ⁵⁶ (v0.43.0) for transcript abundance estimation. Following alignment, we performed differential expression analysis 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 ⁵⁸. 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-responder tumor microenvironment at each time point. The CIBERSORTx ³¹ algorithm was used to estimate the relative proportions of 22 mouse hematopoietic immune cell types based on the transcriptomic profiles of each sample, using the LM22 matrix as a reference. We broadly classified the 22 cell sub-types into 9 major populations by collapsing several related subpopulations as follows: B cells include memory, naïve, and plasma cells; CD4 T cells include CD4 memory resting, memory activated, naïve, follicular helper; macrophages include M0, M1, and M2 phenotypes; NK cells include activated and resting cells; dendritic

cells include activated and resting cells; granulocytes include neutrophils, mast cells resting, mast cells activated and eosinophils. Raw data for all subpopulations is available in Data S2. Before analysis, transcript-level data were library-sized, and gene length normalized to TPM.

We clustered time course RNAseq data using the fuzzy c-means (FCM) clustering algorithm Mfuzz ⁵⁹ in the TCseq package ²⁵. Z-normalised/scaled counts were used in the algorithm and expression profiles 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 ⁶⁰.

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 ²⁶, which achieved the best performance on the DREAM5 network inference challenge ⁶¹. 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 ⁶² 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 ⁶³ (Fig. S2). We identified and retained “direct connections”, defined 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 ⁶⁴. The BEDtools ⁶⁵ “window” function was used to obtain all direct regulatory interactions and GENIE3 importance scores appended to these interactions using the R data.table package (v1.12.2).

Network visualization and analysis of IFN-stimulated genes

We visualized the induced subgraph of the top ten nodes and their first order neighbors in both AB1 and Renca using the R igraph package with the Kamada-Kawai ⁶⁶ 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 confirmed 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 significant 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 ⁶⁸ .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 profiles for ISGs common to both AB1 and Renca networks using the R ComplexHeatmap ⁶⁹ 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 confined 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 ²⁸. 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 axis – non-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’ workflow (version 2 chemistry) as per the manufacturer’s directions. We aimed to capture 9000 cells per sample. Libraries were quantified using the TapeStation D1000 kit (Agilent). Sequencing was performed by Novogene, using NovaSeq S2 flowcell sequencing protocols.

For single cell analysis, we processed FASTQ files 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 ⁷⁰ (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 identification and UMAP embedding and visualization.

To avoid subjective biases in cell identification, we used an automated labelling strategy based on bulk RNAseq references. The R package SingleR ⁷¹ was used in “cluster mode” using species-specific annotation references provided with the package. For annotation of human single cell data, we used the human primary cell atlas reference ⁷². For annotation of murine data, we used the mouse RNAseq dataset from Benayoun et al ⁷³.

Clusters were defined from Seurat’s FindClusters function at default (0.8) resolution. Similarly, labelled clusters were merged. We confirmed that this approach was robust to cluster size by showing that labels were consistent even when cluster size was modified 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 ⁷⁴ 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 non-parametric 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 influenza hemagglutinin ⁷⁵. 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. ⁷⁶ in which the authors identified 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 specific 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 ⁷⁷ 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 ⁷⁸ package, which detects statistically significant coordinated variability at single cell level. Briefly, from the original 38K cells in our AB1 samples, we constructed KNN error models for the 17K surviving default SCDE’s library size filters. 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 ²⁴ was conducted using the escape ⁷⁹ R package (v1.3.1), which provides convenience functions for the GVSA algorithm ⁸⁰. 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 ⁸¹. Filtering, normalization and scaling was performed as outlined in the original paper ²⁴ .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 ⁸². We summarized the result of these analyses using a difference heatmap between responders and non-responders of average binarized TF activity per Seurat cluster. Specifically, 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 quantifies cell transcriptional activity by modelling the time derivative of gene expression states. We used the velocyto ⁸³ package for this analysis. Per-sample loomfiles, containing quantified spliced vs. unspliced transcripts, were combined using the python loompy (v3.0.6) package. After count normalization, filtering and feature selection on these genes, approximately 2.5K genes survived these filtering steps to be used for gene velocity modelling.

We compared “transcriptional momentum” of AB1 monocytes, which we defined as the squared L2 norm for each cell’s embedding vectors with respect to their UMAP ⁸⁴ 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 filtering. 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 ³¹ was used to estimate the relative proportions of 7 cytokine induced T cell signatures based on the transcriptomic profiles of each sample, where the induced T cell gene signature developed by Cano-Gamez et al. ³⁰ was used as a reference (Fig. S5). We broadly classified the 94 samples into 7 major populations by collapsing several related sub-populations by their cytokine treatment: IFNB, Resting, Th17, Th2, Th1, Th0, iTreg to generate the reference file (Data S4). Prior to analysis, gene count data for both AB1 and Renca was normalized to TPM. The data was filtered to retain genes with an TPM value > 0.3 in at least 8 samples (being the smallest experimental group size). CIBERSORT was run on AB1 and Renca separately, with quantile normalization disabled as recommended for RNAseq data.

IFN modulation drug dosing schedules

As we aimed to boost the IFN response using poly(I:C) in the tumor microenvironment prior to administration of ICB to improve response, the timing of administration of ICB antibodies was scheduled late, as to have a low background response rate to ICB. Dosing with drugs commenced on day 12 . Poly(I:C) (HMW, Invivogen) was dosed intratumorally at 50 μg daily for 3 days. Recombinant murine IFNa or IFNb (Biolegend) was dosed intratumorally at 40,000 U daily for 3 days. ICB dosing began 3 days after the final dose of poly(I:C) or recombinant IFN, on day 17. Anti-IFNAR1 (Bioxcell, clone MAR1-5A3, 0.5 mg i.p.), anti-IFNγ (Bioxcell, clone XMG1.2, 0.5 mg i.p.), anti-IFNα (Leinco, clone TIF-3C5, 1 mg i.p., which blocks subtypes A, 1, 4, 5, 11, and 13), anti-IFNβ (Leinco, clone HDß-4A7, 0.6 mg i.p.), or IgG2a isotype (Leinco, clone C1.18.4, 0.6 mg i.p.). Treatment began 3 days after the first day of ICB administration, on day 20, and dosed every 3rd day for a total of 3 doses. For treatment schedules without poly(I:C), treatment with ICB commenced on day 10 for AE17 or day 12 for AB1 and Renca, followed by anti-IFN treatment 3 days after. Treatments were administered by one investigator (RMZ), while tumors were measured at least 3 times weekly using calipers by another researcher (TC) who was blinded for treatment allocation, to guarantee blinded assessment of the primary endpoint.

Flow cytometry of Irf1⁺ cells

For flow cytometric analysis of Irf1 gene expression in different cell subsets, AB1 tumors (n= 6) were harvested 6 days after inoculation and immediately submerged in cold PBS, cut into 1-2 mm pieces with a scalpel blade and dissociated using the GentleMACS system (Miltenyi). Fc block (anti-CD16/CD32, BD) was used for 10 minutes on ice. Cells were stained with Fixable Viability Stain 780 (BD) for 30 minutes at RT, to discriminate live cells. Cells were stained using antibodies for surface markers for 30 minutes at 4^oC (Table S1). To identify Irf1⁺ cells, we used the PrimeFlow Kit (Invitrogen). Briefly, cells were fixed in RNA fixation buffer 1, permeabilized with RNA Permeabilization buffer with RNase inhibitors, then fixed with RNA fixation buffer 2 before using Target Probes against Irf1. The signal was then amplified, followed by addition of fluorescent label probes (Alexa Fluor 647). Data were acquired on a BD Fortessa flow cytometer and analyzed using FlowJo software (TreeStar). Cells were gated on Irf1⁺, followed by CD45⁺ to identify immune infiltrating cells, and CD45^- non-immune cells (e.g. tumor cells). Immune cell populations were analyzed by their expression of CD11b and Ly6C: Ly6C^- Monocytes (CD11b⁺, Ly6C^-); Ly6C⁺ Monocytes (CD11b⁺, Ly6C⁺, also F4/80^-, CD3^-, CD335^-); Other Ly6C⁺ cells (CD11b⁺/^-, Ly6C⁺); and remaining cells (CD11b^-, Ly6C^-). See Supplementary Materials for gating strategies.

Flow cytometry of YFP⁺ cells

For flow cytometric analysis to determine which cells express IFNB, we used B6.129-Ifnb1^tm1Lky/J mice, which co-express IFNβ1 and eYFP from the Ifnb1 locus ³⁵. AE17 tumours (n=3) were treated with poly(I:C) i.t. 21 days after inoculation and harvested 24h later. Tumors were cut into 1-2 mm pieces with a scalpel blade and dissociated using the GentleMACS system (Miltenyi). Fc block (anti-CD16/CD32, BD) was used for 10 minutes on ice. Cells were stained with Fixable Viability Stain 780 (BD) for 30 minutes at RT, to discriminate live cells. Cells were stained using antibodies for surface markers for 30 minutes at 4^oC (Table S1). To detect YFP+ cells, cells were fixed using a cytofix/cytoperm kit (BD), then stained using antibodies against GFP (YFP cross-reactive) in perm buffer overnight. Data were acquired on a BD Fortessa flow cytometer and analyzed using FlowJo software (TreeStar). Cells were gated for single and live cells. Cells were then gated on YFP+, which were also CD45⁺ indicating they are immune infiltrating cells. Immune cell populations were analyzed by their expression of CD11b and Ly6C. YFP+ CD11b⁺ immune cells were MHC-II^-, CD11c^-, F4/80^-, Ly6G^-, CD19^-, CD3^- and NK1.1^-. See Supplementary Materials for gating strategies.

Statistics

The sample size calculation for in vivo mouse experiments was based on prior experiments in which we found that the median survival time on the control treatment (ICB alone) was 35 days ^15,16. Using a proportional hazards model we determined that, if the true hazard ratio (relative risk) of control subjects relative to experimental subjects is 5, we would need to study 10 experimental subjects and 10 control subjects to be able to reject the null hypothesis that the experimental and control survival curves are equal with probability (power) 0.8. The type I error probability associated with this test of this null hypothesis is 0.05. Differences in population frequencies in responders and non-responders using flow cytometry were determined using Mann-Whitney U testing on means. Prism software (GraphPad) was used to analyze tumor growth and to determine statistical significance of differences between groups by applying a Mann-Whitney U test. P-values were adjusted for multiple comparisons using the Benjamini-Hochberg (B-H) method; those < 0.05 were considered significant. The Kaplan-Meier method was used for survival analysis, and p-values were calculated using the log-rank test (Mantel-Cox). For comparison of tumor size when complete response was not achieved, growth curves were analyzed using the TumGrowth package ⁸⁵. For comparison of deconvolution estimations, we used two-way ANOVA with Tukey's multiple comparisons test.

IFNβ RT-PCR on sorted cell populations

To analyze IFNβ gene expression in different cell subsets, AB1 tumors were treated with 50ug of Poly(I:C) (HMW, Invivogen) i.t. (n=5) or untreated (n=5), 5 days after inoculation, then harvested 24 hours later. Tumors were immediately submerged in cold PBS, cut into 1-2 mm pieces with a scalpel blade and dissociated using the GentleMACS system (Miltenyi). Fc block (anti-CD16/CD32, BD) was used for 10 minutes on ice. Cell were stained with UV Zombie live/dead (Biolegend) for 30 minutes at RT, to discriminate live cells. Cells were stained using antibodies for surface markers for 30 minutes at 4ºC. Cells were then sorted into RNAlater (Invitrogen) for the following populations: non-immune cells (CD45-); Ly6C^hi monocytes (CD45⁺ CD11b⁺ Ly6C^hi CD3^- CD335^-), Ly6C^lo/- monocytes (CD45⁺ CD11b⁺ Ly6C^lo/-CD3^-CD335^-), and the remaining immune cell (CD45⁺ CD11b^-). RNA was extracted using the RNAqueous-Micro Kit (Life Technologies). Resulting purified RNA was reverse transcribed using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). Next, we performed RT-PCR using TaqMan Fast Advanced Master Mix (Applied Biosystems) and TaqMan Assay mouse IFNβ1 (Mm00439552_s1, ThermoFisher) or mouse GAPDH (Mm99999915_g1, ThermoFisher) in triplicate for each sample in a MicroAmp optical plate (Applied Biosystems) using QuantStudio 7 Flex Real-Time PCR System (Applied Biosystems). Ifnb1 expression was calculated as dCT of the housekeeping gene GAPDH. ΔCt = Ct ( Ifnb1) – Ct (GAPDH).

1. Hegde, P. S. & Chen, D. S. Top 10 Challenges in Cancer Immunotherapy. Immunity 52, 17–35 (2020).

2. Wei, S. C., Duffy, C. R. & Allison, J. P. Fundamental Mechanisms of Immune Checkpoint Blockade Therapy. Cancer Discov 8, 1069–1086 (2018).

3. Binnewies, M. et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med 24, 541–550 (2018).

4. Yarchoan, M., Hopkins, A. & Jaffee, E. M. Tumor Mutational Burden and Response Rate to PD-1 Inhibition. N Engl J Med 377, 2500–2501 (2017).

5. Lesterhuis, W. J. et al. Dynamic versus static biomarkers in cancer immune checkpoint blockade: unravelling complexity. Nat Rev Drug Discov 16, 264–272 (2017).

6. Zemek, R. M. et al. Bilateral murine tumor models for characterizing the response to immune checkpoint blockade. Nat Protoc (2020) doi:10.1038/s41596-020-0299-3.

7. Lancet Oncology, T. Calling time on the immunotherapy gold rush. The Lancet Oncology 18, 981 (2017).

8. Fullerton, J. N. & Gilroy, D. W. Resolution of inflammation: a new therapeutic frontier. Nat Rev Drug Discov 15, 551–67 (2016).

9. Bi, K. et al. Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma. Cancer Cell (2021) doi:10.1016/j.ccell.2021.02.015.

10. Cloughesy, T. F. et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med 25, 477–486 (2019).

11. Chin, W. L., Zemek, R. M., Lesterhuis, W. J. & Lassmann, T. Functional genomics in cancer immunotherapy: computational approaches for biomarker and drug discovery. Mol. Syst. Des. Eng. 4, 689–700 (2019).

12. Hugo, W. et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 168, 542 (2017).

13. Mariathasan, S. et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–548 (2018).

14. Riaz, N. et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 171, 934-949.e16 (2017).

15. Lesterhuis, W. J. et al. Network analysis of immunotherapy-induced regressing tumours identifies novel synergistic drug combinations. Sci Rep 5, 12298 (2015).

16. Zemek, R. M. et al. Sensitization to immune checkpoint blockade through activation of a STAT1/NK axis in the tumor microenvironment. Sci Transl Med 11, (2019).

17. Chen, I. X. et al. A bilateral tumor model identifies transcriptional programs associated with patient response to immune checkpoint blockade. Proc Natl Acad Sci U S A 117, 23684–23694 (2020).

18. Jerby-Arnon, L. et al. A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. Cell 175, 984-997.e24 (2018).

19. Lee, H.-S. et al. Comprehensive immunoproteogenomic analyses of malignant pleural mesothelioma. JCI Insight 3, (2018).

20. Ock, C.-Y. et al. Genomic landscape associated with potential response to anti-CTLA-4 treatment in cancers. Nature Communications 8, 1050 (2017).

21. Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M. & Alizadeh, A. A. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259 (2018).

22. Chen, P.-L. et al. Analysis of Immune Signatures in Longitudinal Tumor Samples Yields Insight into Biomarkers of Response and Mechanisms of Resistance to Immune Checkpoint Blockade. Cancer Discov 6, 827–837 (2016).

23. Herbst, R. S. et al. Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature 515, 563–567 (2014).

24. Bassez, A. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021).

25. Wu & Gu. TCseq: Time course sequencing data analysis. R package version 1.14.0. (2020).

26. Huynh-Thu, V. A., Irrthum, A., Wehenkel, L. & Geurts, P. Inferring regulatory networks from expression data using tree-based methods. PLoS ONE 5, (2010).

27. Liberzon, A. et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739–1740 (2011).

28. Krzywinski, M., Birol, I., Jones, S. J. & Marra, M. A. Hive plots—rational approach to visualizing networks. Brief Bioinform 13, 627–644 (2012).

29. Salem, M. L., El-Naggar, S. A., Kadima, A., Gillanders, W. E. & Cole, D. J. The adjuvant effects of the toll-like receptor 3 ligand polyinosinic-cytidylic acid poly (I:C) on antigen-specific CD8+ T cell responses are partially dependent on NK cells with the induction of a beneficial cytokine milieu. Vaccine 24, 5119–32 (2006).

30. Cano-Gamez, E. et al. Single-cell transcriptomics identifies an effectorness gradient shaping the response of CD4(+) T cells to cytokines. Nat Commun 11, 1801 (2020).

31. Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457 (2015).

32. Ng, C. T. et al. Blockade of interferon Beta, but not interferon alpha, signaling controls persistent viral infection. Cell Host Microbe 17, 653–61 (2015).

33. Ma, R.-Y. et al. Monocyte-derived macrophages promote breast cancer bone metastasis outgrowth. J Exp Med 217, (2020).

34. Serbina, N. V. & Pamer, E. G. Monocyte emigration from bone marrow during bacterial infection requires signals mediated by chemokine receptor CCR2. Nat Immunol 7, 311–317 (2006).

35. Scheu, S., Dresing, P. & Locksley, R. M. Visualization of IFNbeta production by plasmacytoid versus conventional dendritic cells under specific stimulation conditions in vivo. Proc Natl Acad Sci U S A 105, 20416–20421 (2008).

36. Ayers, M. et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest 127, 2930–2940 (2017).

37. Keenan, T. E., Burke, K. P. & Van Allen, E. M. Genomic correlates of response to immune checkpoint blockade. Nat Med 25, 389–402 (2019).

38. Jamal-Hanjani, M., Quezada, S. A., Larkin, J. & Swanton, C. Translational implications of tumor heterogeneity. Clin Cancer Res 21, 1258–1266 (2015).

39. Mosely, S. I. S. et al. Rational Selection of Syngeneic Preclinical Tumor Models for Immunotherapeutic Drug Discovery. Cancer Immunol Res 5, 29–41 (2017).

40. Wei, S. C. et al. Distinct Cellular Mechanisms Underlie Anti-CTLA-4 and Anti-PD-1 Checkpoint Blockade. Cell 170, 1120-1133.e17 (2017).

41. Germain, R. N. The art of the probable: system control in the adaptive immune system. Science 293, 240–245 (2001).

42. Benci, J. L. et al. Tumor Interferon Signaling Regulates a Multigenic Resistance Program to Immune Checkpoint Blockade. Cell 167, 1540-1554 e12 (2016).

43. Jacquelot, N. et al. Sustained Type I interferon signaling as a mechanism of resistance to PD-1 blockade. Cell Res 29, 846–861 (2019).

44. Qu, Y. et al. Baseline Frequency of Inflammatory Cxcl9-Expressing Tumor-Associated Macrophages Predicts Response to Avelumab Treatment. Cell Rep 32, 107873 (2020).

45. Bassez, A. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021).

46. Katlinskaya, Y. V. et al. Suppression of Type I Interferon Signaling Overcomes Oncogene-Induced Senescence and Mediates Melanoma Development and Progression. Cell Rep 15, 171–180 (2016).

47. Park, S. H. & Rehermann, B. Immune responses to HCV and other hepatitis viruses. Immunity 40, 13–24 (2014).

48. Channappanavar, R. et al. Dysregulated type I interferon and inflammatory monocyte-macrophage responses cause lethal pneumonia in SARS-CoV-infected mice. Cell Host Microbe 19, 181–193 (2016).

49. Teijaro, J. R. et al. Persistent LCMV infection is controlled by blockade of type I interferon signaling. Science 340, 207–211 (2013).

50. Wilson, E. B. et al. Blockade of chronic type I interferon signaling to control persistent LCMV infection. Science 340, 202–207 (2013).

51. Griffiths, J. I. et al. Circulating immune cell phenotype dynamics reflect the strength of tumor-immune cell interactions in patients during immunotherapy. Proc. Natl. Acad. Sci. U.S.A. 117, 16072–16082 (2020).

52. Lee, J. S. et al. Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19. Sci Immunol 5, (2020).

53. Chen, D. S. & Mellman, I. Oncology meets immunology: the cancer-immunity cycle. Immunity 39, 1–10 (2013).

54. Ciechomska, M. & Skalska, U. Targeting interferons as a strategy for systemic sclerosis treatment. Immunol Lett 195, 45–54 (2018).

55. Lesterhuis, W. J. et al. Synergistic effect of CTLA-4 blockade and cancer chemotherapy in the induction of anti-tumor immunity. PloS one 8, e61895 (2013).

56. Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34, 525–527 (2016).

57. Pimentel, H., Bray, N. L., Puente, S., Melsted, P. & Pachter, L. Differential analysis of RNA-seq incorporating quantification uncertainty. Nature Methods 14, 687–690 (2017).

58. Yi, L., Pimentel, H., Bray, N. L. & Pachter, L. Gene-level differential analysis at transcript-level resolution. Genome Biol. 19, 53 (2018).

59. Futschik, M. E. & Carlisle, B. Noise-robust soft clustering of gene expression time-course data. J Bioinform Comput Biol 3, 965–988 (2005).

60. Kuleshov, M. V. et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res 44, W90–W97 (2016).

61. Marbach, D. et al. Wisdom of crowds for robust gene network inference. Nat. Methods 9, 796–804 (2012).

62. Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res 4, 1521 (2015).

63. Khan, A. et al. JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework. Nucleic Acids Res 46, D260–D266 (2018).

64. Karolchik, D. et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res. 32, D493-496 (2004).

65. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).

66. Kamada, T. & Kawai, S. An algorithm for drawing general undirected graphs. Information Processing Letters 31, 7–15 (1989).

67. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550 (2005).

68. Smedley, D. et al. The BioMart community portal: an innovative alternative to large, centralized data repositories. Nucleic Acids Res 43, W589-98 (2015).

69. Gu, Z., Eils, R. & Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849 (2016).

70. Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).

71. Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20, 163–172 (2019).

72. Mabbott, N. A., Baillie, J. K., Brown, H., Freeman, T. C. & Hume, D. A. An expression atlas of human primary cells: inference of gene function from coexpression networks. BMC Genomics 14, 632 (2013).

73. Benayoun, B. A. et al. Remodeling of epigenome and transcriptome landscapes with aging in mice reveals widespread induction of inflammatory responses. Genome Res. 29, 697–709 (2019).

74. Heng, T. S. P., Painter, M. W., & Immunological Genome Project Consortium. The Immunological Genome Project: networks of gene expression in immune cells. Nat. Immunol. 9, 1091–1094 (2008).

75. Marzo, A. L., Lake, R. A., Robinson, B. W. & Scott, B. T-cell receptor transgenic analysis of tumor-specific CD8 and CD4 responses in the eradication of solid tumors. Cancer Res. 59, 1071–1079 (1999).

76. Park, J. et al. Single-cell transcriptomics of the mouse kidney reveals potential cellular targets of kidney disease. Science 360, 758–763 (2018).

77. Tickle, T., Tirosh, I., Georgescu, C., Brown, M. & Haas, B. InferCNV: Inferring copy number alterations from tumor single cell RNA-Seq data. GitHub repository (2019).

78. Fan, J. et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nat Methods 13, 241–4 (2016).

79. Borcherding, N. & Andrews, J. escape: Easy single cell analysis platform for enrichment. (Bioconductor version: Release (3.13), 2021). doi:10.18129/B9.bioc.escape.

80. Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics 14, 7 (2013).

81. Wickham, H. ggplot2: Elegant Graphics for Data Analysis. (Springer-Verlag New York, 2016).

82. Aibar, S. et al. SCENIC: Single-cell regulatory network inference and clustering. Nat Methods 14, 1083–1086 (2017).

83. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

84. McInes, L., Healy, J., Saul, N. & Großberger, L. UMAP: Uniform Manifold Approximation and Projection. Journal of Open Source Software 3, 861 (2018).

85. Enot, D. P., Vacchelli, E., Jacquelot, N., Zitvogel, L. & Kroemer, G. TumGrowth: An open-access web tool for the statistical analysis of tumor growth curves. Oncoimmunology 7, e1462431 (2018).

Yes there is potential Competing Interest. Patent application pertaining to aspects of this work. W.J.L.: consultancy Douglas Pharma-ceuticals, MSD, research funding Douglas Pharmaceuticals, AstraZeneca, ENA therapeu-tics. A.K.N. advisory boards Boehringer Ingelheim, Bayer, Roche, BMS; research funding from AstraZeneca. M.M.: Advisory boards MSD, BMS, Roche, AstraZeneca

Data4StimulatedTCellReference.xlsx
Data 4
Data2CIBERSORTabsolute.xlsx
Data 2
Data1DEgenesbulk.xlsx
Data 1
Data3TCseqclusters.xlsx
Data 3
TableS1.docx
supplementaryfiles.docx

Download PDF

Journal Publication

published 19 Aug, 2022

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

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

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results

Discussion

Declarations

Star Methods

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1