Research ArticleCancer

Sensitization to immune checkpoint blockade through activation of a STAT1/NK axis in the tumor microenvironment

See allHide authors and affiliations

Science Translational Medicine  17 Jul 2019:
Vol. 11, Issue 501, eaav7816
DOI: 10.1126/scitranslmed.aav7816

Spot the differences between tumors

Immune checkpoint inhibitors, a form of cancer immunotherapy, have revolutionized cancer treatment in recent years, but unfortunately, not all tumors respond to these drugs. To gain insight into why some tumors might not respond, Zemek et al. compared gene expression patterns and cellular makeup of murine tumors that did or did not respond to immune checkpoint targeting. The authors then identified an immunotherapeutic combination that could be used to convert the microenvironment to a more favorable configuration in mice, suggesting that it may be possible to sensitize patients’ tumors to immunotherapy by a similar approach.


Cancer immunotherapy using antibodies that target immune checkpoints has delivered outstanding results. However, responses only occur in a subset of patients, and it is not fully understood what biological processes determine an effective outcome. This lack of understanding hinders the development of rational combination treatments. We set out to define the pretreatment microenvironment associated with an effective outcome by using the fact that inbred mouse strains bearing monoclonal cancer cell line–derived tumors respond in a dichotomous manner to immune checkpoint blockade (ICB). We compared the cellular composition and gene expression profiles of responsive and nonresponsive tumors from mice before ICB and validated the findings in cohorts of patients with cancer treated with ICB antibodies. We found that responsive tumors were characterized by an inflammatory gene expression signature consistent with up-regulation of signal transducer and activator of transcription 1 (STAT1) and Toll-like receptor 3 (TLR3) signaling and down-regulation of interleukin-10 (IL-10) signaling. In addition, responsive tumors had more infiltrating-activated natural killer (NK) cells, which were necessary for response. Pretreatment of mice with large established tumors using the STAT1-activating cytokine interferon-γ (IFNγ), the TLR3 ligand poly(I:C), and an anti–IL-10 antibody sensitized tumors to ICB by attracting IFNγ-producing NK cells into the tumor, resulting in increased cure rates. Our results identify a pretreatment tumor microenvironment that predicts response to ICB, which can be therapeutically attained. These data suggest a biomarker-driven approach to patient management to establish whether a patient would benefit from treatment with sensitizing therapeutics before ICB.


Immunotherapy has become a remarkable success story in the treatment of cancer (1). In particular, treatment with monoclonal antibodies targeting immune checkpoint molecules such as programmed death receptor 1 (PD-1) and cytotoxic T lymphocyte–associated antigen 4 (CTLA4) improves survival in patients with melanoma, non–small cell lung cancer, and several other cancer types (24). However, a large proportion of patients do not benefit, and some cancer types seem to be less sensitive to immune checkpoint blockade (ICB) than others (5, 6). Some patients have an immediate and complete regression of all their tumors, sometimes within weeks (7), whereas other patients do not experience any therapeutic response whatsoever (6).

Although several correlates of response to ICB have been reported, such as expression of checkpoint ligands (8), mutational load (9), neoantigen expression (10), interferon (IFN) signatures (11), and an inflammatory tumor microenvironment (12), no definitive predictive biomarkers have been identified (6). It is unclear whether it is possible to manipulate a nonresponsive tumor microenvironment toward a responsive phenotype. Consequently, the clinical development of combination therapies is largely empiric and sometimes based on scant preclinical data (13, 14). With more than 2000 currently active clinical trials involving ICB worldwide (15), there is a need to prioritize immunotherapy combinations, preferably based on preclinical data, as to limit patient exposure to futile treatments and potentially severe side effects (5).

Ideally, this would involve therapeutically altering the tumor microenvironment to a favorable phenotype before therapy, allowing subsequent assessment of the tumor state to define a go/no go signal for subsequent checkpoint-targeted treatment. This would not only improve prioritization of immunotherapeutic combinations to test in patients but also allow patient stratification in experimental trials and clinical practice.

Here, we made use of the fact that even in the highly homogeneous setting of inbred mouse strains bearing tumors derived from monoclonal cancer cell lines, there remains a dichotomy in responsiveness to treatment with ICB (6, 1621). Against this highly uniform background, we asked whether we could identify a signature in tumors before treatment that would correlate with response and whether we could use that information to turn nonresponders into responders. Using flow cytometry and bulk and single-cell RNA sequencing (RNA-seq) data from tumors derived from two different mouse cancer models, we mapped the key cellular and molecular networks associated with response, the results of which were corroborated in several datasets of patients with cancer treated with ICB (22, 23). We prioritized upstream regulators for therapeutic targeting with drugs, recombinant proteins, or antibodies and found that we could drive the tumor microenvironment from a nonresponsive into a responsive state.


Tumors derived from clonal cancer cell lines and grown in inbred mouse strains display two distinct gene signatures, predicting sensitivity to ICB

There are many potential reasons why a definitive biomarker for the response to ICB has not emerged, including differences in host genetics, environmental factors, and the diverse genetic and cellular makeup of cancers between patients (22, 24). Even inbred mouse strains bearing transplantable tumors display a dichotomous outcome after immunotherapy (Fig. 1A) (1621). This is unexpected because the genomes of these mice are nominally identical. In addition, the tumors are derived from a clonal cell line, excluding the possibility that a difference in tumor rejection antigen expression caused these disparate responses. Mice are age- and gender-matched, kept under controlled conditions, and receive identical treatment. However, they respond very differently. Differences in outcomes between animals may be related to differences in T cell repertoire, which is not completely encoded in the germ line (25), or to stochastic immunological events (26). Regardless of the cause of this dichotomy, we reasoned that this model would allow us to assess potentially small differences in microenvironmental regulation of therapeutic responses in a controlled background. We inoculated mice bilaterally with either AB1 mesothelioma cells or Renca kidney cancer cells and treated them with anti-CTLA4 and anti–programmed death receptor ligand 1 (PD-L1) ICB antibodies, which resulted in symmetric yet dichotomous responses (Fig. 1, B and C). These models thus allow the biological assessment of a whole tumor at any time point, by surgically removing it, while still being able to infer the therapeutic fate of that tumor if it had been left in situ, by monitoring the fate of the remaining contralateral tumor (16). We used this model to characterize the pretreatment tumor microenvironment. We surgically removed one tumor for analysis using bulk RNA-seq, flow cytometry, or single-cell RNA-seq shortly before administering ICB. The remaining tumor was assessed for therapeutic outcome (Fig. 1D).

Fig. 1 Inbred mouse strains carrying monoclonal tumors display a symmetrical yet disparate response to ICB, associated with a distinctive gene signature.

(A) Representative Renca tumor growth curves showing ICB nonresponders (red), intermediate responders (orange), and responders (blue) to ICB. (B and C) Representative growth curves of mice with bilateral AB1 tumors (B) or Renca tumors (C) treated with ICB (n = 10), color-coded per mouse. (D) Experimental design. (E and F) PCA of responsive (RS) and nonresponsive (NR) AB1 tumors (E) and Renca tumors (F) (n = 12 tumors per group). (G and H) Unsupervised-hierarchical clustering of top differentially expressed genes separating responsive and nonresponsive tumors. For AB1 (G), 10,307 genes were differentially expressed in responders versus nonresponders (top 200 shown) and 127 genes for Renca (H) (all shown; also see data file S1).

Using principal components analysis (PCA) of the bulk RNA-seq data, we found that responsive and nonresponsive tumors clustered separately in both models (Fig. 1, E and F). Unsupervised hierarchical clustering of the top differentially expressed genes also resulted in a clear separation of responsive and nonresponsive tumors (Fig. 1, G and H, and data file S1). Thus, untreated, ostensibly identical tumors in inbred mice displayed an unexpected heterogeneity in gene expression profiles, differentiating animals that were going to be responders from nonresponders even before they were treated with immunotherapy.

ICB responsive tumors are characterized by an inflammatory microenvironment, driven by signal transducer and activator of transcription 1

We aimed to gain insight into the biological relevance of the differentially expressed genes in responsive and nonresponsive tumors. First, we noted a notable difference between the two models in terms of differentially expressed genes between responders and nonresponders, with AB1 having more than 10,000 genes differentially expressed, whereas in Renca only 127 genes were differentially expressed. However, the majority of the 127 genes in Renca overlapped with AB1 (118 genes), resulting in a refined response-associated signature (data file S1).

To further characterize the overarching biological processes associated with this signature, we analyzed the distribution of a well-characterized and curated collection of “hallmark” gene sets using gene set enrichment analysis (GSEA; Fig. 2A and fig. S1) (27). This revealed enrichment of genes associated with a type I and type II IFN response and an inflammatory environment in responsive tumors in both murine models. These gene sets were also enriched in pretreatment tumors from a cohort of patients with urothelial cancer treated with PD-L1 antibody atezolizumab that went on to respond (Fig. 2B) (23), validating the translational relevance of the findings. Pathway analysis of the common differentially expressed genes in the two models identified antigen presentation and T helper 1 (TH1) type immune responses as the most enriched pathways in responsive tumors (Fig. 2C). These data accord with published data in human melanoma (11) and suggest that the activation of inflammatory pathways, particularly the IFN response, renders tumors sensitive to ICB in animal cancer models and patients alike.

Fig. 2 Inflammatory pathways with STAT1 as a key regulator are enriched in ICB responsive tumors in mouse models and patients.

(A) GSEA top hallmark gene sets in responsive versus nonresponsive tumors. TNFα, tumor necrosis factor–α; NF-κB, nuclear factor κB. (B) GSEA of responsive versus nonresponsive tumors from the patient cohort (n = 192) (23). NES, normalized enrichment score; ES, enrichment score. (C) Overrepresented canonical pathways within responsive tumors. (D) WGCNA module enrichment for differentially expressed genes between responsive and nonresponsive tumors. Bars indicate standard deviation with outliers; horizontal dashed line, P = 0.05. (E) Previous knowledge-based graphical reconstruction of the wiring diagram of module 1. (F) GSEA of responsive [complete response (CR)] versus nonresponsive [progressive disease (PD)] tumors from the patient cohort (n = 192) (23), using a STAT1 gene set (54). (G) Representative pSTAT1 immunohistochemistry in nonresponsive and responsive AB1 tumors. Scale bars, 100 μm. (H) Survival curves of mice with pSTAT1 high– and low ICB–treated tumors (n = 10 responsive and 8 nonresponsive; log-rank test, ***P < 0.001).

Next, we used weighted gene correlation network analysis (WGCNA) (28) to map the molecular networks underlying responsiveness to ICB. This identified seven modules of highly coexpressed genes operating within tumors, one of which was significantly (P < 0.05) up-regulated in responsive tumors in both models [module 1 in Fig. 2D and fig. S2, dashed line (P = 0.05)]. This response-associated module was enriched for genes involved in immunity, including natural killer (NK) cell–mediated cytotoxicity, PD-1 signaling, and costimulation (fig. S3). Using NetworkAnalyst (29), we assembled a putative network that identified Stat1 as a hub (Fig. 2E). These findings were corroborated in single-cell analysis of AB1 tumors, which identified greater Stat1 gene expression in both immune and nonimmune cells in a responding tumor compared to a nonresponding tumor (fig. S4). In addition, there was a significant (P = 0.019) enrichment of a signal transducer and activator of transcription 1 (STAT1) signature in responsive tumors from the urothelial cancer/atezolizumab cohort (23), as well as in a second cohort of patients with melanoma treated with nivolumab (P = 0.044; Fig. 2F and fig. S5) (22). Because STAT1 is activated through phosphorylation, which can be assessed by immunohistochemistry, we assessed total and phosphorylated STAT1 abundance in our mouse models and found that responsive tumors had significantly higher percentages of positive cells than nonresponsive tumors (P = 0.001; Fig. 2, G and H, and fig. S6, A and B). Similarly, the absence of a STAT1 gene signature (30) in human pretreatment samples was associated with a decreased progression-free survival in the urothelial cancer/atezolizumab clinical dataset (fig. S6C), which accords with a previous report on the predictive value of tumor STAT1 activation in patients with melanoma treated with anti–PD-1 (30). Together, these data suggest that STAT1 activation is a driver of the ICB response–associated tumor microenvironment and can serve as a potential biomarker to identify patients more likely to respond.

Cellular analyses of resistant and sensitive tumors identify the presence of NK cells as a prerequisite for response to ICB

Because we observed an enrichment of genes associated with an inflammatory IFN/STAT1-driven environment in responsive tumors and because it has been previously found that “hot” tumors are characterized by increased CD8 T cell infiltration (30), we examined the immune cell infiltrates of responsive and nonresponsive tumors. We used flow cytometry on dissociated tumors and compared the data with a CIBERSORT analysis (31); CIBERSORT is a deconvolution approach for characterizing cell composition of complex tissues that we applied to our RNA-seq data. AB1 tumors were characterized by a predominantly myeloid infiltrate, whereas the infiltrate in Renca tumors was mostly lymphoid (Fig. 3, A and B). Despite the notable differences in cellular infiltrates, these models respond similarly to checkpoint blockade. We did not observe any consistent difference between responders and nonresponders with regard to infiltrating CD8 or CD4 T cells, B cells, macrophages, monocytic cells, or dendritic cells (Fig. 3, A and B, fig. S7, and data file S2). However, there was a greater proportion of NK cells (flow cytometry) and activated NK cells (CIBERSORT) in responding tumors in both models (Fig. 3, C and D). The percentage of overall leukocytes (identified as CD45+ cells by flow cytometry) was the same in responsive and nonresponsive tumors (fig. S8A). The percentage of tumor-infiltrating NK cells (CD335+/CD3) of all tumor-containing cells was significantly increased (P < 0.0001) in responders (fig. S8B). To assess whether tumor NK cell enrichment could be relevant in humans, we interrogated gene expression data from a cohort of patients with urothelial cancer, obtained before treatment with anti–PD-L1 (23). Using CIBERSORT, we found that gene sets specific for activated NK cells significantly correlated with objective radiological response [complete response (CR) and partial response (PR) versus progressive disease (PD); P = 0.034; Fig. 3E and fig. S9], similar to the mouse models. To test whether NK cell infiltration of the pretreatment tumor microenvironment was required for response in our models, we depleted NK cells with a single injection of anti–asialo-GM1 3 days before ICB in both AB1 and Renca models, which resulted in a significantly (P = 0.0185 and P < 0.0001, respectively) diminished response (Fig. 3, F and G). These data show that very different cellular tumor microenvironments between models are still conducive to a response to ICB. Furthermore, pretreatment tumor NK cell infiltration is a common prerequisite for response to ICB.

Fig. 3 NK cells are enriched in ICB-sensitive tumors in mouse models and patients and are required for response.

(A) Flow cytometry of dissociated tumors from responsive (n = 6) and nonresponsive (n = 9) AB1 tumor– or responsive (n = 7) and nonresponsive (n = 11) Renca tumor–bearing mice. (B) CIBERSORT analysis from responsive and nonresponsive tumors (n = 12 per group). (C) NK cells (CD335+) in responsive and nonresponsive tumors shown in (A). (D) NK fraction relative to total leukocyte infiltrate by CIBERSORT (Mann-Whitney U test with B-H correction for multiple comparisons). (E) Activated NK cell fraction (CIBERSORT) in tumors from the patient cohort. PR, partial response; SD, stable disease. n = 298, Mann-Whitney U test; bars indicate standard deviation). (F and G) Survival plots of AB1 tumor– (F) or Renca tumor–bearing mice (G) treated with ICB, with or without anti–asialo-GM1, 3 days before starting ICB. Mice were treated with ICB early (day 5 for AB1 and day 7 for Renca), allowing interrogation of response attenuation (n = 10 mice per group, two independent experiments, log-rank test). *P < 0.05, **P < 0.01, ***P < 0.001.

A rationally developed pretreatment combination therapy sensitizes tumors to ICB

Having identified genes differentially expressed in responsive tumors, we then rationalized that by targeting the regulators of these genes, we could induce nonresponsive tumors to respond to ICB. We performed upstream regulator analysis (URA) (32) to identify higher-level regulators of the inflammatory pathways and networks enriched in responsive tumors. URA identifies transcriptional regulators that, based on previous experimental data, might be expected to account for the genes that are differentially expressed. Despite the marked difference in the number of differentially expressed genes and cellular infiltrates between the models, the regulators associated with response were highly similar, with the most significant (P < 0.0001) positive regulators being IFNγ and STAT1, and interleukin-10 receptor subunit α (IL-10RA) being the top negative regulator [Fig. 4A, dashed line (P = 0.05), and fig. S10]. The patient cohort showed similar results (Fig. 4B and fig. S11). Having identified a regulator signature associated with response, we reasoned that by targeting the regulators predicted to phenocopy this signature before ICB, we could convert nonresponsive to responsive tumors. We focused on therapeutics that are clinically readily available (at least in phase 2 clinical trials) and therefore chose a pretreatment schedule of IFNγ, anti–IL-10, and/or Toll-like receptor 3 agonist poly(I:C), which was also one of the top regulators and is known to activate STAT1 (33). We administered a short course of these treatments, for 3 days only, followed by ICB 3 days later (Fig. 4C). In addition to the models previously described, we also used AE17 mesothelioma and B16 melanoma–bearing C57BL/6 mice, because both are relatively resistant to ICB (34). In all cases, we observed that mice pretreated with the triple combination were sensitized to ICB, with significantly increased response rates, from 0 to 10% with ICB alone to up to 80% for the combination therapy (Renca, P < 0.0001; AB1, P = 0.0362; AE17, P = 0.002; B16, P = 0.0002; Fig. 4, D to G).

Fig. 4 Rational therapeutic modulation of the tumor microenvironment sensitizes tumors to ICB.

(A) URA predicted regulators of the response-associated gene signature in mouse models (red, positive correlation; blue, negative correlation). (B) URA-predicted regulators in responders versus nonresponders of the patient cohort (n = 192) (23). (C) Sensitizing treatment schedule. (D to G) Survival curves of Renca (D), AB1 (E), AE17 (F), and B16 tumor–bearing mice (G) mice treated with the triple combination [poly(I:C) + IFNγ + anti–IL-10], followed by ICB compared to ICB alone (n = 10 per group; n = 15 for AE17). (H) Survival curves of Renca tumor–bearing mice treated with single agent therapy [anti–IL-10, poly(I:C), or IFNγ] versus the triple combination before ICB. (I) Survival curves of Renca-bearing mice treated with double-agent therapy versus the triple combination before ICB (n = 10 mice per group, two independent experiments of five mice per group; log-rank test compared to ICB-only group, *P < 0.05, **P < 0.01, ***P < 0.001). (J) Treatment schedule to test ICB followed by triple-combination therapy. (K) Survival curves of Renca-bearing mice treated with the triple combination followed by ICB, as shown in (C), or with the reverse schedule, as shown in (J) (n = 10 mice per group, two independent experiments of five mice per group; log-rank test. *P < 0.05, **P < 0.01).

To determine whether this was the effect of a single drug or due to the predicted complex treatment combination, mice were pretreated with each selected therapeutic alone, followed by ICB, in which case, we observed no added benefit compared to ICB alone (Fig. 4H). Although pretreatment with a combination of poly(I:C) and either IFNγ or anti–IL-10 resulted in sensitization of some tumors, the triple combination was superior (Fig. 4I). Next, to determine whether the triple combination was sensitizing the tumor to ICB, rather than enforcing the effector response, we changed the scheduling and compared the triple combination followed by ICB (Fig. 4C) to ICB followed by the triple combo (Fig. 4J). When checkpoint blockade was given first, followed by the triple combination, response rates were similar to ICB alone. When the combination was given first, although start of ICB was substantially delayed compared to the controls, the tumors had been sensitized and responded to therapy, in both Renca (Fig. 4K) and AB1 (fig. S12) models. These data show that a rational complex combination of multiple clinically available therapeutics results in marked sensitization to ICB.

Sensitizing treatment induces the responsive phenotype and is dependent on NK cells

To test whether the sensitizing therapeutic combination induced the response-associated phenotype, we analyzed the tumors by flow cytometry for STAT1 activation and NK cell infiltration after 3 days of pretreatment with IFNγ, poly(I:C), and anti–IL-10 or with vehicle controls (Fig. 5A). We found that the sensitizing pretreatment significantly (AB1, P = 0.079; Renca, P = 0.0002) increased the frequencies of CD335+ cells (Fig. 5B and fig. S13A). In addition, there was an increase in phospho-STAT1 (pSTAT1)–positive and IFNγ-producing leukocytes infiltrating the tumors (Fig. 5, C and D). In addition, the pretreatment increased PD-L1, pSTAT1, and major histocompatibility complex class I (MHC-I) expression on tumor cells (Fig. 5E and fig. S13, B and C). Although the pretreatment-induced increase in IFNγ production was derived from multiple leukocyte subsets in the tumor, this was significant only for the NK cell population (P < 0.0001), not for the other cell subsets (Fig. 5F). This was also the case for phosphorylation of STAT1 (P = 0.027; Fig. 5G). To confirm that the triple combination, rather than a single agent, was driving this response phenotype, we assessed the effect of each agent alone compared to the triple combination. We found that the triple combination was superior over single treatments in terms of increased NK cell infiltration, IFNγ production, and STAT1 phosphorylation (fig. S14). Furthermore, recruitment of NK cells was greatest 1 to 3 days after the final dose of the triple combination, suggesting that it would be optimal to start ICB early after pretreatment (fig. S15). Further phenotypic characterization of these CD335+ NK cells in the tumor revealed that they were conventional NK cells and not tissue-resident CD335+ type 1 or type 3 innate lymphoid cells (ILCs) (fig. S16) (35). High expression of markers CD11b and killer cell lectin-like receptor G1 (KLRG1) confirmed activation and maturity of these conventional NK cells (Fig. 5, H and I). When we depleted NK cells before pretreatment, we completely abolished the sensitizing effect to ICB (Fig. 5J), suggesting that this was mediated through treatment-induced infiltration of circulating NK cells in the tumor microenvironment. To further substantiate the role of STAT1 in the sensitizing effect of the triple combination therapy, we functionally inactivated the STAT1 pathway using blocking antibodies against IFNγ/IFNAR [interferon (IFN) alpha receptor], which are both upstream of STAT1. This significantly diminished tumor NK cell infiltration after triple combination therapy (P = 0.002; fig. S17). Together, these data show that the triple combination therapy attracts conventional NK cells into the tumor microenvironment in a STAT1-dependent manner, resulting in sensitization to subsequent ICB.

Fig. 5 Sensitizing triple therapy activates the STAT1/NK axis in the tumor microenvironment.

(A) Schedule of treatment and harvest for analysis by flow cytometry. (B to E) NK cell fraction (B), STAT1 activation (C), and IFNγ production (D) by CD45+ tumor-infiltrating leukocytes and PD-L1 expression by CD45 nonleukocytes (E) after treatment with IFNγ, poly(I:C), and anti–IL-10 (bars indicate standard deviation; Mann-Whitney U test, **P < 0.01, ***P < 0.001, ****P < 0.0001). (F and G) IFNγ expression (F) and STAT1 phosphorylation (G) in tumor-infiltrating lymphocytes (CD45+ cells). Mann-Whitney U test. (H and I) CD11b (H) and KLRG1 (I) expression by CD335+ NK cells in tumors from Renca-bearing mice after 3-day treatment with PBS (untreated) or IFNγ, poly(I:C), and anti–IL-10 (treated). Bars indicate standard deviation; Mann-Whitney U test, *P < 0.05, **P < 0.01. (J) Survival curves of Renca-bearing mice, receiving sensitization followed by ICB, with or without NK-depleting anti–asialo-GM1 antibody (n = 10 per group; n = 5 PBS controls; log-rank test, ***P < 0.001). (K) Proposal for an adaptive treatment approach.


Together, our data demonstrate that a pretreatment tumor microenvironment dominated by infiltrating activated NK cells and an inflammatory gene expression signature characterized by STAT1 activation are sensitive to ICB and that this profile can be therapeutically induced. These data point to a two-step approach to treating patients with cancer, in which tumor profiling would allow a decision to treat initially with ICB or after pretreatment with sensitizing therapeutics (Fig. 5K).

We have previously argued that the therapeutic response to ICB can be visualized as a critical transition, moving from a pretreatment cancer state toward a normal tissue state (6). Because complex systems that display critical transitions tend to be highly sensitive to initial conditions (36), we set up our mouse models so that in a highly homogeneous background and under controlled conditions, we would still observe dichotomous outcomes, with a proportion of mice displaying a full cure and a proportion no response whatsoever. We presumed that we would find little, if any, difference in pretreatment gene expression between responsive and nonresponsive tumors (6). To our surprise, however, although in the Renca model, only a little over hundred genes were differentially expressed between responders and nonresponders, this was markedly different for AB1, where despite the homogeneous background and clonality of the tumors, several thousand genes were differentially expressed. Our data thus uncover a notable immunological heterogeneity preceding treatment in subcutaneous tumor models derived from clonal cell lines in inbred mouse strains, which are widely used in cancer biology and drug discovery. It is likely that this will affect the results of immunotherapy studies or other treatments such as chemotherapy, which are influenced by immunological mechanisms (37).

In addition, we noted large differences in the cellular makeup of the Renca and AB1 tumors, although they displayed a similar response rate to ICB. Our data indicate not only that ICB-sensitive, inflammatory (hot) tumors are characterized by markedly increased numbers of immune cells such as CD8+ T cells (30) but also that discrete differences in cell populations, such as NK cells or activation state (as measured by STAT1 activation), can be associated with major consequences in gene expression and thus sensitivity to ICB, which can be therapeutically attained.

Limitations of our analysis include the fact that we have focused on transcriptional events and thus will have missed other mechanisms that drive the response. In addition, we have only taken a single snapshot of the pretreatment state, so future studies will need to address the temporal immunological effects associated with the therapeutic response, for example, using bilateral tumor models (6). We acknowledge that some effects appear to be more pronounced in the mouse models compared to the human datasets, such as the correlation between NK numbers and response. This is likely related to the heterogeneity of patients, their cancers, and environmental factors, whereas in our mouse models, those factors are all controlled, allowing the detection of small differences that would otherwise have gone unnoticed. We also cannot exclude the possibility that differences in cancer type and the antibodies tested may have resulted in slight differences in the importance of the identified regulators between the mouse models and human datasets.

With many different immune checkpoint inhibitors currently available, in addition to other antibodies, cytokines, and small-molecule drugs targeting the immune system or cancer cells, there is an abundance of choice for combination therapies. Currently, these combinations are tested in many thousands of patients in the context of clinical trials, not infrequently with little preclinical rationale (5, 38, 39). In particular, when more than two drugs are given together, the number of potential combinations becomes prohibitive for testing in clinical trials. In this context, we showed that effective complex combinations of therapeutics can be identified by systems analysis of gene expression data from responding and nonresponding tumors, and we identified a treatment combination that sensitizes tumors to ICB very effectively in multiple difficult-to-treat preclinical models. Our data will provide a resource for further investigating the biology associated with response to ICB and for identification of combination treatments in the context of immuno-oncology.


Study design

Our primary objective was to define the pretreatment microenvironment associated with an effective outcome by comparing gene expression and flow cytometry data from ICB-responsive and nonresponsive tumors within the same mouse cancer model. Our secondary objective was to test whether we could therapeutically alter the tumor microenvironment toward a responsive phenotype and thus increase the response to ICB. By implanting a tumor subcutaneously on both flanks of the mouse, we were able to remove one entire tumor 1 hour before therapy with ICB and then monitor the remaining tumor for response. The removed tumor was analyzed by flow cytometry, bulk RNA-seq, or single-cell RNA-seq and categorized as a responsive or nonresponsive tumor based on the outcome of the corresponding tumor left in situ. Tumors that had an intermediate response (PR or relapse) were excluded from analysis. We performed these experiments using AB1 mesothelioma and Renca renal cell carcinoma cell lines. To exclude model-specific effects, we focused on genes that were differentially expressed between responsive and nonresponsive tumors in both models.

The sample size for the bulk RNA-seq experiments was estimated using the method developed by Hart et al. (40); for sample sizes of n = 12 and a within group coefficient of variation of 0.3, there is >90% power to detect a 1.7-fold change in gene expression. Differences in population frequencies in responders and nonresponders using flow cytometry and CIBERSORT were determined using Mann-Whitney U test on means.

We then used these transcriptomic data to identify pathways and regulators that could be targeted. We validated the findings using publicly available RNA-seq data from two cohorts of patients with cancer (22, 23). We used in vivo targeting studies with tumor growth as an end point. These studies were performed in BALB/c mice carrying AB1 and Renca tumors and were validated in unrelated subcutaneous AE17 mesothelioma and B16 melanoma models in C57BL/6 mice, to establish robustness across tumor models and mouse strains. Mice were randomized after tumors were established, before tumor removal and therapy. Treatments were administered by one researcher, and tumors were measured by another researcher who was blinded for treatment allocation.

The sample size calculation for in vivo mouse experiments was based on previous experiments, in which we found that the median survival time on the control treatment (ICB alone) was 35 days (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) of 0.8. The type I error probability associated with this test of this null hypothesis is 0.05.


BALB/cArc, BALB/cJAusb, and C57BL6/J mice 8 to 12 weeks of age were used for all experiments. Mice were obtained from the Animal Resource Centre (Murdoch, WA) or the Australian BioResources (Moss Vale, NSW) and housed at the Harry Perkins Institute of Medical Research Bioresources Facility under specific pathogen-free conditions. Mice were fed Rat and Mouse cubes (Specialty Feeds) and had access to water ad libitum. Cages (Tecniplast) were individually ventilated with filtered air, contained aspen chips bedding (TAPVEI), and were supplemented with tissues, cardboard rolls, and wood blocks as environmental enrichment. Cages were changed every 14 days. Mice were housed at 21° to 22°C with 12-hour light/12-hour dark cycle (06:00 to 18:00). Sentinel mice (n = 3) in the animal facility were screened monthly for a standard panel of bacteria and fungi, ectoparasites, endoparasites, nonpathogenic protozoa, and viruses (Cerberus Sciences). 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 and AE091).

Cell culture

Cell lines AB1 and AE17 were obtained from CellBank Australia. Cell line Renca was donated by E. Sotomayor and F. Cheng (University of South Florida, Tampa, FL). Cell line B16 was obtained from American Type Culture Collection (CRL-6475). Cell lines were maintained in RPMI 1640 (Invitrogen) supplemented with 20 mM Hepes, 0.05 mM 2-mercaptoethanol, penicillin (100 U/ml; CSL), gentamicin (50 μg/ml; David Bull Labs), and 10% fetal bovine serum (FBS; Invitrogen). Cells were grown to 70 to 80% confluence before passage and passaged three to five times before inoculation. Cells were frequently tested for mycoplasma by polymerase chain reaction (PCR) and remained negative. Cell lines were validated yearly by flow cytometry for MHC-I molecules H2-Kb (consistent with C57BL/6) and H2-Kd (consistent with BALB/c), and for fibroblast markers E-cadherin, epithelial cell adhesion molecule, and platelet-derived growth factor receptor α (negative) and by PCR for mesothelin (positive for AB1 and negative for Renca).

In vivo treatments

When cell lines were 70 to 80% confluent, they were harvested and washed three times in phosphate-buffered saline (PBS). A total of 5 × 105 cells in 100 μl were inoculated subcutaneously onto the lower right-hand side (RHS) flank (for single inoculations) or both flanks (for dual-tumor inoculations) using a single 26-gauge needle per injection. Mice were randomized when tumors became palpable, about 3 to 5 days after tumor inoculation. Treatments were administered by one investigator (R.M.Z.), whereas tumors were measured at least three times weekly using calipers by a researcher (T.H.C.) who was blinded for treatment allocation, to guarantee blinded assessment of the primary end point.

Surgery experiments

Seven (AB1) or 10 (Renca) days after tumor inoculation, when tumors were ~9 mm2, mice were subcutaneously dosed with buprenorphine (0.1 mg/kg) in 100 μl (30 min before anesthesia) and anesthetized using isoflurane (4% in 100% oxygen at a flow rate of 2 liters/min). Whole tumors and the corresponding draining inguinal lymph node on the RHS were removed by surgical excision and immediately immersed in RNAlater (Life Technologies) for RNA-seq or cold PBS for single-cell RNA-seq or flow cytometry. The wound was closed with staples (Able Scientific). Mice were placed in a heat box for recovery. One hour after surgery, mice were administered ICB. The remaining tumor was monitored 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 at least 4 weeks after treatment. Mice were designated as nonresponders if their tumors grew to 100 mm2 within 4 weeks after the start of treatment, similar to saline-treated controls. 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; in any cage, there had to be at least one nonresponder among responders or vice versa.

In vivo ICB treatment

The anti–PD-L1 hybridoma (clone MIH5) and the anti-CTLA4 hybridoma (clone 9H10) were cultured in IMDM (Iscove’s Modified Dulbecco’s Medium) containing 1% of FBS and gentamycin. Clarified supernatants were used to purify the antibody using affinity chromatography. The antibody was sterile formulated in PBS. Mice received an intraperitoneal dose of 100 μg of anti-CTLA4 and 100 μg of anti–PD-L1 combined in 100 μl of PBS. Mice received additional doses of 100 μg of anti–PD-L1 2 and 4 days later. In previous experiments (41), we had not found any difference in effect of control immunoglobulin G (IgG) versus PBS, and therefore, vehicle controls received PBS alone. We varied the time of treatment initiation after tumor inoculation as to obtain a suitable background response rate to ICB, high for experiments in which responses were to be negated and low in experiments in which the response rate was to be improved.

Tumor preparation for RNA-seq

RNA-seq, flow cytometry, and single-cell RNA-seq were performed on tumors from separate experiments. For RNA-seq, whole tumors and lymph nodes were surgically resected, and the surrounding tissue was removed and immediately submerged in RNAlater (Life Technologies). Samples were stored at 4°C for 24 hours, after which supernatant was removed and samples were transferred to −80°C. Frozen tumors were dissociated in TRIzol (Life Technologies) using a TissueRuptor (QIAGEN). RNA was extracted using chloroform and purified on RNeasy MinElute columns (QIAGEN). RNA integrity was confirmed on the Bioanalyzer (Agilent Technologies). Library preparation and sequencing (50–base pair single-end) were performed by Australian Genome Research Facility, using Illumina HiSeq standard protocols.

Single-cell RNA-seq

For single-cell RNA-seq, tumors were harvested and submerged in cold PBS, cut into 1- to 2-mm pieces with a scalpel blade, and dissociated using the gentleMACS system (Miltenyi Biotec). Single-cell suspensions were stored at −80°C in 10% dimethyl sulfoxide until they were processed for single-cell profiling. Cryo-stored cells were rapidly thawed and diluted in PBS and pelleted. Pellets were resuspended in PBS and passed through a 40-μm filter to remove cell clumps. About 5000 cells per sample were then loaded onto a 10x Genomics Chromium Controller to generate Chromium Single Cell 3′ Libraries. Sequencing was carried out by the Australian Genome Research Foundation in Melbourne, Australia.

Primary analysis was carried out using the 10x Genomics Cell Ranger software suite. Raw sequencing data were demultiplexed (Cell Ranger v2.1.1 and bcl2fastq v2.20.0.422), reads were aligned to the reference (GRCm38, Ensembl 84 build), and gene expression was quantified using the Cell Ranger count command. Using default parameters, some droplets containing real cells were discarded. We therefore used the expected cells 6000 flag to include these cells. Responder and nonresponder samples were combined using Cell Ranger aggr function. We annotated all cells in our samples using the SingleR software package (42).

Flow cytometry

For flow cytometry, tumors were harvested and immediately submerged in cold PBS, cut into 1- to 2-mm pieces with a scalpel blade, and dissociated using the gentleMACS system. Fc block [anti-CD16/CD32, Becton Dickinson (BD)] was used for 10 min on ice. Zombie UV live/dead (BioLegend) was used to discriminate live cells. Cells were permeabilized and fixed using FoxP3 Fix/Perm buffer kit (BioLegend) before addition of antibodies for 20 min at room temperature (RT). See table S1 for antibody details. Cells were kept in stabilizing fixative until acquisition. Data were acquired on a BD Fortessa flow cytometer and analyzed using FlowJo software (TreeStar).

Cell populations were defined as follows: monocytes (CD11b+MHCII+/−Ly6C+F4/80CD11c), macrophages (CD11b+ MHCII+F4/80+CD11c+/−Ly6CLy6G), immature myeloid cells (CD11b+MHCIIF4/80+Ly6CLy6G), neutrophils/granulocytes (CD11b+Ly6G+MHC Ly6CintF4/80), CD8 T cells (CD3+CD8+), CD4 helper T cells (CD3+CD4+FoxP3), regulatory T cell (CD3+CD4+FoxP3+), NK cells (CD335+CD3), and B cells (CD19+CD3). See fig. S18 for gating strategies.

For detection of IFNγ, cells were incubated with brefeldin A at 37°C for 4 hours. Fixable Viability Stain 620 (BD) was used to identify dead cells. Cells were surface-stained with antibodies for 20 min at RT. After permeabilization with Cytofix/Cytoperm (BD), cells were stained with anti-IFNγ (Thermo Fisher Scientific) in perm/wash [PBS/0.1% saponin (Sigma-Aldrich)] for 30 min on ice.

For detection of pSTAT1, dissociated tumors were stained with surface antibodies, fixed with 1.5% formaldehyde for 10 min, and then permeabilized with ice-cold methanol at 4°C for 10 min. Cells were stained with a rabbit–anti-pSTAT1 antibody (Tyr701)–PE (Cell Signaling Technology) for 30 min at RT.

For detection of intracellular NK/ILC1 markers, cells were permeabilized using FoxP3 buffer kit (Thermo Fisher Scientific). See table S1 for antibodies used.

ICB-sensitizing drug dosing schedules

Dosing with pretreatment drugs commenced on day 15 for AB1 and Renca, day 10 for AE17, or day 8 for B16. IFNγ (Shenandoah Biotechnology) was dosed subcutaneously into tumor area at 50,000 U daily for 3 days. Poly(I:C) [high molecular weight (HMW), Invivogen] was dosed subcutaneously into tumor area at 50 μg daily for 3 days. The anti–IL-10 hybridoma (clone JES5.2A5) was cultured in IMDM containing 1% of FBS and gentamicin. Clarified supernatants were used to purify the antibody using affinity chromatography. The antibody was sterile-formulated in PBS. Anti–IL-10 was dosed intraperitoneally at 0.5 mg per mouse daily for 3 days. ICB was dosed 2 days after the pretreatment drugs: day 20 for AB1 and Renca, day 15 for AE17, or day 13 for B16; for the ICB-only group, dosing began at the same time as the pretreatment drug dosing commenced in the other arms.

NK cell depletion

To investigate the impact of NK cell depletion, ICB antibodies were administered early (day 5 for AB1 or day 7 for Renca) to give a high-background response rate. Anti–asialo-GM1 (Wako Chemicals) was dosed at 20 μg in 50 μl of saline and injected intravenously 3 days before ICB administration (day 2 for AB1 or day 4 for Renca). NK cell depletion was verified by flow cytometry analysis of peripheral blood using antibodies for CD45 and CD335; we also used antibodies for CD3, CD4, CD8, ICOS, and Ki67 to detect CD8 and CD4 T cells and their activation state.

Bulk RNA-seq data analysis

Read libraries were quality-assessed using FastQC (v0.11.3) (43) and mapped to the mouse genome (mm11) at both the transcript and gene levels using HISAT2 (v2.0.4) (44). Gene-level quantitation (counts) of aligned reads was performed using summarizeOverlaps (45), and transcript discovery and quantification were performed using StringTie (v1.3.0) (46) and Ballgown (47).

PCA was performed to visualize the major contributions of variation within the data using the prcomp() function within R (v3.3.3). Gene count data were transformed for PCA using the variance-stabilizing transformation (48). The top 1000 most variable genes across samples (those with the greatest median absolute deviation) were selected for PCA. Differentially expressed genes were identified between immunotherapy nonresponders and responders using DESeq2 (49). P values were adjusted for multiple comparisons using the Benjamini-Hochberg (B-H) method; P < 0.05 was considered significant (50). A full list of differentially expressed (DE) genes is in data file S1.

Differentially expressed genes were uploaded to InnateDB ( along with associated gene expression data. A list of pathways mapping to the uploaded genes was returned, and pathway analysis was undertaken to determine which pathways were significantly overrepresented (P < 0.05) in the up- and down-regulated gene datasets. InnateDB simultaneously tests for overrepresentation of DE genes in more than 3000 pathways, from which we looked at the KEGG (Kyoto Encyclopedia of Genes and Genomes; and Reactome ( databases. The B-H false discovery rate (FDR) correction was applied to correct for multiple testing.

The WGCNA algorithm was used to construct a signed network across all samples and identify clusters (modules) of genes with highly correlated patterns of gene expression. To prepare the data for WGCNA, we (i) applied a variance stabilizing transformation to normalize the data, (ii) filtered out genes expressed at a low level (only those with at least 10 counts per sample were retained), (iii) removed genes without an official Mouse Genome Informatics (MGI) symbol, and (iv) removed genes with low variability by applying the varianceBasedfilter() function within the DCGL package (significance threshold set to 0.01). The resulting set of 6026 genes was used as input for network construction. Network modules of coexpressed genes identified by WGCNA were tested for enrichment of differentially expressed genes between immunotherapy nonresponders and responders by plotting the –log10 P values derived from the DESeq2 analysis, on a module-by-module basis. Protein-protein interactions of the differentially expressed genes were created using InnateDB network analysis tool (unfiltered) and network created using NetworkAnalyst (29), with STRING interactome, confidence cutoff of 900 (51). Rendering of module interactions was performed with Cytoscape software. Network modules of interest and differentially expressed genes (filtered to include those with an absolute fold change in expression of >2) were also analyzed within Ingenuity Systems (32) to identify associated upstream transcriptional regulators, using right-tailed Fisher’s exact tests and default settings for other options. P < 0.05 was considered significant. An activation z score was also calculated for each upstream transcriptional regulator by comparing their known effect on downstream targets with observed changes in gene expression. Those with activation z scores of ≥2 or ≤−2 were considered “activated” or “inhibited,” respectively.

The CIBERSORT (31) algorithm was used to estimate the relative proportions of 25 mouse hematopoietic immune cell types based on the transcriptomic profiles of each sample, where the 511 mouse gene signature developed by Chen et al. (52) was used as a reference. We broadly classified the 25 cell types into 12 major populations by collapsing several related subpopulations as follows: B cells include memory, naïve, and plasma cells; CD8+ T cells include memory, naïve, and activated cells; CD4 T cells include memory, naïve, follicular, TH1, TH2, and TH17 cells; macrophages include M0, M1, and M2 phenotypes; NK cells include activated and resting cells; and dendritic cells include activated and immature cells. Before analysis, transcript-level data were library-sized, and gene length was normalized using the Ballgown package (53), resulting in fragments per kilobase of transcript per million mapped reads (FPKM). Individual transcripts were collapsed to gene-level data based on the mean FPKM value using the aggregate() function. Last, the data were filtered to retain genes with an FPKM value of >0.3 in at least eight samples (the smallest experimental group size).

The Broad Institute GSEA software (54) was used to run analyses on normalized gene expression data with prefiltering for low counts. We used the hallmarks gene set database, which includes 50 MSigDB hallmarks gene sets (27) and a STAT1 signature. The STAT1 signature was derived from Care et al. (55) and defined on the basis of variance-ranked Spearman correlations of gene expression across 11 diffuse large B cell lymphoma (DLBCL) datasets, selecting strongly correlated genes (cutoff median correlation score of >0.6). A total of 1000 permutations were performed. All other default parameters were used. Gene sets enriched at a nominal P < 0.05 and FDR < 0.25 were considered significant.

Survival analysis for patients with STAT1 signature enrichment

To define the subset of patients with STAT1 pathway activation for survival analysis, we closely followed the Classification Algorithm Based on a Biological Signature (56). We used the previously defined STAT1 gene set from Care et al. (55) and leveraged RNA-seq data from the IMvigor210 trial (23), available as an R data package (IMvigor210CoreBiologies). The dataset contains gene count data of 348 patients with urothelial cancer treated with anti–PD-L1 antibody atezolizumab. Patients were classified by best response: 25 with CR, 43 with PR, 63 with stable disease (SD), and 167 with PD. Fifty patients did not have evaluable disease. From these RNA-seq data, two prototype vectors were defined on the basis of mean expression values of genes in this STAT1 gene set; a “responder” prototype from radiological complete responders and partial responders and a “progressor” prototype based on progressor samples. For each patient, expression profiles of the STAT1 gene set were compared to the two prototype vectors. Similarity was calculated on the basis of Pearson correlation coefficient, and a decision score was formulated on the basis of the ratio of correlations. This decision score was used to classify all patients in the IMvigor210 dataset regardless of response, with a decision score of >1 denoting higher activation of the STAT1 pathway. Survival analysis was performed using the Log-rank test, dividing patients based on the decision score, with a decision score of >1 denoting higher activation of the STAT1 pathway; a Kaplan-Meier survival curve was constructed.

Immunohistochemistry for STAT1 and pSTAT1

For immunohistochemical staining of STAT1 and pSTAT1, slides of 4-μm thickness were cut from formalin-fixed, paraffin-embedded tissue blocks. Slides were deparaffinized in two changes of xylene, followed by rehydration in changes of 100, 96, 70, and 40% ethanol and distilled water. For STAT1 staining, antigen retrieval using 10 mM citrate buffer (pH 6.0) was performed for 10 min at 121°C using a pressure cooker (Cuisinart). For pSTAT1 staining, slides were sub-boiled in 1 M EDTA buffer (pH 8.0) for 10 min in a microwave. Slides were rinsed in tris-buffered saline with Tween 20 (TBST), and the endogenous peroxidase was blocked using 3% hydrogen peroxidase (Sigma-Aldrich) in distilled water. Sections were washed again in TBST and blocked with goat serum (VECTASTAIN) diluted in PBS as per the manufacturer’s instructions. Primary antibody against STAT1 or pSTAT1 (dilution of 1:800 and 1:200, respectively; Cell Signaling Technology) was incubated for 60 min at RT. Sections were washed, and secondary antibody, goat anti–rabbit IgG-horseradish peroxidase (Santa Cruz Biotechnology), was incubated for 30 min at RT, before washing again. Betazoid DAB chromogen (Biocare Medical) was prepared as per the manufacturer’s instructions. Chromogen was applied, sections were monitored for the development of a brown color, and the reaction was stopped with H2O. Sections were counterstained with hematoxylin and rehydrated before mounting coverslips with Pertex. An experienced pathologist (C.L.) scored the sections by assigning an estimated proportion of positive tumor cells (positivity defined as moderate or strong nuclear staining for pSTAT1 and moderate or strong nuclear and cytoplasmic positivity for STAT1), while blinded to treatment outcome.


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 B-H method; P < 0.05 was considered significant. The Kaplan-Meier method was used for survival analysis, and P values were calculated using the log-rank test (Mantel-Cox).


Fig. S1. Individual GSEA graphs of the top 8 hallmark gene sets in AB1 and Renca.

Fig. S2. WGCNA showing expression of modules up-regulated in responders in AB1 and Renca separately.

Fig. S3. Reactome and KEGG pathway enrichment of response-associated module (module 1) in overlap of AB1 and Renca.

Fig. S4. Single-cell analysis of responsive/nonresponsive AB1 tumors.

Fig. S5. STAT1 expression and correlation with response in human patient cohorts.

Fig. S6. STAT1 immunohistochemistry on AB1 and Renca tumors.

Fig. S7. CIBERSORT of individual samples of responders and nonresponders in AB1 and Renca before treatment with anti-CTLA4/anti–PD-L1.

Fig. S8. Analysis of flow cytometry cell populations as a percentage of all events.

Fig. S9. CIBERSORT stacked graph of overall cell populations in atezolizumab patient cohort.

Fig. S10. URA of differentially expressed genes in AB1, Renca, and module 1 (combined AB1 and Renca) by P value and z score.

Fig. S11. Negative regulators of response of the atezolizumab cohort identified by URA.

Fig. S12. Sequencing of treatments and antitumor efficacy of combination immunotherapy.

Fig. S13. Flow cytometry analysis after treatment with IFNγ + poly(I:C) + anti–IL-10.

Fig. S14. Flow cytometry analysis of single-agent– compared to triple-agent–treated tumors.

Fig. S15. Analysis of tumors 1, 3, or 6 days after triple combination by flow cytometry.

Fig. S16. Conventional CD335+ NK cells in tumors.

Fig. S17. The IFN-STAT1 pathway and NK cell infiltration.

Fig. S18. Flow cytometry gating strategy.

Table S1. Flow cytometry antibodies.

Data file S1. Lists of differentially expressed genes between responders and nonresponders.

Data file S2. Cell proportion data for Fig. 3 (A and B).


Acknowledgments: We acknowledge the facilities and the scientific and technical assistance of the Australian Microscopy and Microanalysis Research Facility at the Centre for Microscopy, Characterisation and Analysis, The University of Western Australia, a facility funded by the University, State and Commonwealth Governments. Cell line Renca was donated by E. Sotomayor and F. Cheng (University of South Florida, Tampa, FL). We thank C. To, C. Tilsed, C. Rinaldi, P. Deng, and D. Hope for technical assistance and I. Dick for statistical advice. Funding: This work was funded by project grant 1103980 from the National Health and Medical Research Council (NHMRC). W.J.L. is funded by an NHMRC RD Wright Fellowship and a Fellowship from the Cancer Council of Western Australia. W.L.C. is funded by a Western Australian Cancer and Palliative Care Network fellowship. M.A.D.-E. is funded by an NHMRC Principal Research Fellowship. A.Z. is funded by discovery project DP180100718 by the Australian Research Council. T.L. is supported by a fellowship from the Feilman Foundation. The National Centre for Asbestos Related Diseases is supported by NHMRC Centre of Research Excellence grant 1107043. Part of this research was made possible using equipment funded by the Australian Cancer Research Foundation to establish the ACRF Centre for Advanced Cancer Genomics and with the support of a collaborative cancer research grant provided by the Cancer Research Trust. Author contributions: R.M.Z. performed and interpreted all mouse experiments, flow cytometry experiments, GSEA, and ingenuity pathway analyses and co-wrote the paper. E.D.J. performed and interpreted the bulk RNA-seq analyses, URA, and WGCNA. W.L.C. performed and interpreted the computational analyses of the human RNA-seq data. I.S.S. and M.A.D.-E. developed, designed, and performed the IFNγ and NK/ILC flow cytometry panels, interpreted the data, and aided in manuscript preparation. V.S.F., C.F., and T.H.C. provided technical assistance with surgery and in vivo treatment experiments. S.J.D. and C.L. designed, performed, and analyzed the pSTAT1 immunohistochemistry. L.B. provided critical reagents and aided in manuscript preparation. A.R.R.F. and D.O.M. designed and performed the single-cell RNA-seq experiments. R.A.L., M.J.M., and A.K.N. were involved in experimental design and interpretation of all in vivo studies, the flow cytometry experiments and the human dataset analysis, and critical revision of the manuscript. T.L. and W.L.C. performed and interpreted all analyses pertaining to single-cell RNA-seq data. T.L. supervised and interpreted the computational analysis of the human RNA-seq dataset and critically revised the manuscript. A.B. was involved in experimental design of the murine RNA-seq experiments, supervised, and, together with T.L., A.Z., and M.S., interpreted the computational analyses of bulk RNA-seq data and critically revised the manuscript. W.J.L. designed the study, developed the models, interpreted all data, and supervised the project. R.M.Z and W.J.L. wrote the manuscript with input from all authors. All authors gave final approval of the version to be published. Competing interests: Patent application pertaining to aspects of this work: R.M.Z., W.J.L., R.A.L., and A.B. (PCT/AU2019/050259, “Method for immunotherapy drug treatment”). A.K.N., M.J.M., A.B., R.A.L., and W.J.L.: Research funding and consultancy for Douglas Pharmaceuticals. W.J.L.: Research funding from AstraZeneca. A.K.N. advisory boards for Boehringer Ingelheim, Bayer Pharmaceuticals, Roche Pharmaceutics, and Bristol-Meyers Squib; research funding from AstraZeneca. M.M.: Advisory boards for Merck Sharp & Dohme, Bristol-Myers Squibb, Roche, and AstraZeneca. Data and materials availability: The raw RNA-seq data are available from the Gene Expression Omnibus data repository (accession no. GSE117358). All data associated with this study are present in the paper or the Supplementary Materials.

Stay Connected to Science Translational Medicine

Navigate This Article