Conserved transcriptional connectivity of regulatory T cells in the tumor microenvironment informs n

To explore the impact of Treg cells on the diverse cell states in the TME, we performed single-cell RNA sequencing (scRNA-seq) of sorted CD45− and CD45+ cell populations using the 10X platform (Extended Data Fig. 2a). These populations were isolated from tumor-bearing lungs of KP-DTR mice treated for 48 h with DT or vehicle control 3 months after adenoviral Cre-driven tumor initiation (Fig. 1a). After pre-processing, we clustered cells using PhenoGraph18 and annotated clusters using expression of known markers into major cell types (Extended Data Fig. 2b–d). To ensure our inferences were robust, we focused on the major hematopoietic and non-hematopoietic cell types in the TME that had substantial numbers of cells. The final processed datasets included LECs, VECs, LECs, fibroblasts, lymphoid cells and myeloid cells (macrophages, monocytes, dendritic cells (DCs) and neutrophils; Fig. 2a). Similarly to population-level assessments, scRNA-seq showed that short-term Treg cell depletion had profound effects on transcriptional features of fibroblasts, myeloid and endothelial cells compared to lymphocytes (Extended Data Fig. 2e and Fig. 2a–c). To gain deeper insight into the phenotypic response of accessory cells whose transcriptomes were most affected by Treg removal—endothelial cells, fibroblasts and myeloid cells, we separately clustered and embedded each subtype to ascribe finer-grain identities (Fig. 2d and Extended Data Fig. 3; for annotation strategy see Methods). Furthermore, we used Milo19 to quantify changes in abundance of subpopulations and cell states after Treg cell depletion (Methods). We found several cell states affected by Treg cell depletion, with the most pronounced phenotypic shifts in capillary VECs, mesenchymal stem cells (MSCs), Col14a1 matrix fibroblasts, monocytes and macrophages (Fig. 2d–h and Extended Data Fig. 4a–d). Therefore, Treg cell depletion markedly affected the distribution and abundancies of several cell states and subsets in the TME.

Shared and distinct Treg cell-dependent gene programs

We then sought to characterize genes that respond to Treg cell depletion in these key accessory cell subsets. We used factor analysis to characterize gene expression programs—sets of genes whose expression changes in a coordinated way in a specific set of cells and assessed their differential usage in cell populations from control or DT mice to elucidate the response to Treg cell depletion. Specifically, factor analysis methods are well suited to decompose data into factors, which represent coordinated expression programs across cells and reduce the impact of noise on analysis, which can be dominant at an individual gene level20. We used single-cell hierarchical Poisson factorization (scHPF), designed specifically for scRNA-seq21,22 and applied it to each cell lineage separately to dissect the observed gene expression changes. Each cell and gene present in the expression matrix was assigned a score for each factor, enabling biological interpretation of that factor (see Supplementary Table 2 for factor gene and cell matrices). Factors were robust to random initializations of the model and robust to slight changes in parameters (Methods and Supplementary Fig. 1).

To gain insights into the spatial organization of the identified accessory cell populations, gene programs and their relationship to transcriptional states of tumor cells, we profiled four tissue sections (two control, two Treg cell depleted) using the 10X Visium platform. We used BayesPrism25,26, a Bayesian framework that jointly estimates cell-type fractions and cell-type-specific gene expression using a labeled scRNA-seq reference, to deconvolve each spatial transcriptomics (ST) spot into constituent cell populations. Deconvolution was performed using our scRNA-seq datasets labeled with 26 distinct cell populations selected to optimize granularity, robustness and concordance with underlying histological features in paired H&E-stained sections (Methods, Fig. 4a, Extended Data Fig. 7a–e and Supplementary Table 11). Next, we assessed whether the gene factors that changed upon Treg cell depletion in scRNA-seq were also identified by ST analysis. Consistently, we observed upregulation of endothelial and fibroblast IC and IFN signaling-related gene signatures after Treg cell depletion within spots assigned to the corresponding cell type (Fig. 4b). We also observed increased use of genes associated with the activated VEC factor in capillary aerocyte (aCap) endothelium assigned spots, as well as increased IFN and proliferation related gene signatures in myeloid spots. IC and IFN factors shared many genes across all three analyzed accessory lineages (18 for IC, 103 for IFN), which suggested that similar gene programs were induced across colocalized cell types by common stimuli, indicative of a signaling niche. The spatial behavior of shared genes in these two programs showed localization to two distinct signaling niches in the tissue, with the IC gene program (Cxcl2, Ier3, Fosl1, Il6) localized to the tumor core and the IFN response gene program (Ifit1, Stat1, Isg15, Irf7) localized to the periphery of, or distal to tumor lesions. Inspection of the same H&E-stained sections confirmed dense tumor cell presence with potential hypoxia and neutrophil infiltration at IC foci, and immune cell aggregates at sites with strong IFN response signal (Fig. 4c–e and Supplementary Table 12). Further, ST analysis revealed concordant differential distribution of tumor cells and accessory cell types within these territories with higher frequency of tumor cells, basophils/mast cells, neutrophils and MSCs in IC territories and a high frequency of T cells/type 2 innate lymphoid cells (ILC2s), natural killer (NK) cells, conventional dendritic cells (cDCs), monocytes and alveolar macrophages in IFN territories (Fig. 4f). Taken together, these results point to two primary inflammatory and spatially distinct modes of lung TME response to Treg cell depletion within tumor mass and tumor margin.

Tumor states associated with response to Treg cell depletion

KP LuAds adopt a range of recurrent transcriptional states with features of differentiated lung ECs, their progenitors or epithelial progenitors from other tissues including the gastrointestinal tract and liver and EMT (epithelial to mesenchymal transition)-associated ones27,28,29,30. We next sought to identify potential associations between tumor states and the identified TME niches, that is, IC-positive, IFN-positive and cold (negative) ones. We first identified tumor cells within our ECs by calling KRAS p.Gly12Asp mutations. Because optimized dissociative TME single-cell analysis protocols are suboptimal for capturing tumor cells, we identified only 239 tumor cells within our mouse scRNA-seq dataset. To enable robust deconvolution of tumor cell states, we substituted tumor cells from our scRNA-seq dataset with those from a published dataset that had better capture of KP LuAd cells (KP-tracer dataset; N = 18,083)28. With this updated reference, we performed an additional spot deconvolution to more accurately capture tumor states in the tissue. TME fractions for other cell types remained relatively unchanged between deconvolutions (Extended Data Fig. 7c).

In spots with tumors, the tumor-state fractions exhibited regional variation in gene expression programs, sometimes within seemingly the same tumor nodule (Fig. 5a). Tumor spots were clustered by their tumor-state fractions and typically showed a dominant tumor state in each spot (Extended Data Fig. 8a) manifested in the expression of corresponding marker genes (Extended Data Fig. 8b), forming continuous spatial patches of similar phenotypes (Fig. 5b and Extended Data Fig. 8c). Spots were grouped into tumor lesion areas, or contiguous patches of tumor cells belonging to the same cluster and quantified across control and Treg cell-depleted conditions. Tumor states were also compared across tumors that had a detectable immune response (>10% of spots in IC or IFN signaling niches) or not in Treg cell-depleted sections. Treg cell depletion resulted in a pronounced increase in tumor lesion areas that corresponded to a high-plasticity cell state, specifically among tumor nodules associated with an immune response (Fig. 5c and Extended Data Fig. 8d,e)29. This was consistent with a significant enrichment of high-plasticity cell-state genes upregulated by tumor cells after Treg cell depletion in scRNA-seq (Extended Data Fig. 8f,g and Supplementary Tables 13 and 14). Therefore, TME response to the removal of Treg cells may elicit a gene program in LuAds that represents an unstable transitional state, which can give rise to other tumor states28,29. While IC and IFN niches were observed in the majority of tumor nodules after Treg cell depletion, there were some nodules, and even areas within individual nodules, that did not (Fig. 4c). In particular, those with a gastric epithelial lineage gene expression program were selectively devoid of IC or IFN responses (Fig. 5c and Extended Data Fig. 8e). We assessed differential gene expression between immune response ‘rich’ and ‘poor’ lesion areas and found increased expression of Gkn2 (gastrokine), Pf4 (platelet factor 4/Cxcl4) and Sox9 among other genes (Fig. 5d,e and Supplementary Table 15). Interestingly, Sox9 expression in lung tumor cells was shown to enable their escape from NK cell-mediated killing in certain cases27, suggesting one potential mechanism of immune evasion. Similarly to a recent analysis of CRISPR-edited tumors31, the observed response to Treg cell depletion was spatially restricted, as even ‘nonresponsive’ areas that were directly adjacent to responsive ones were deprived of immune cell or IC signals (Fig.6a,b). Therefore, regional variation in tumor state appears to define the TME response to Treg cell depletion.

a, H&E staining of representative tumor section characterized by histological and immune response state heterogeneity after Treg cell depletion. Insets at bottom represent a zoomed-in view of gastric (left) and high-plasticity (right) areas. Black arrows highlight neutrophil infiltration in a high-plasticity area. b, Tumor RNA fraction within highlighted high-plasticity and gastric epithelial states (left) and gene expression modules (right) of tumor lesion shown in a. Images are representative of one of two serial sections for each of four samples (DT and Ctrl, two biological replicates each).

Conserved Treg cell-dependent features of human and mouse tumor microenvironment

Next, we sought to test whether Treg cell-dependent TME features observed in mice are conserved in human cancer (Fig. 7a) by leveraging variation in Treg cell abundance in 25 primary or local metastatic LuAd samples from 23 individuals, analyzed using scRNA-seq (Supplementary Tables 16 and 17). Despite differences in composition and proportion of accessory cell types in these datasets, we were able to detect all cell populations corresponding to those observed in mice (Fig. 7b and Extended Data Fig. 9a,b). To determine whether the factors induced after Treg cell depletion in mice are present in human LuAd samples with a low abundance of Treg cells, we first determined the frequency of Treg cells among all hematopoietic cells in each sample (Fig. 7c and Extended Data Figs. 9c and 10a,b). Next, we performed scHPF analysis for each of the cell lineages under investigation (Supplementary Table 18) and looked for orthologous genes shared between human and mouse factors to align gene programs between species (Fig. 7d). Then, we assessed the correlation of mean factor usage in single cells to Treg cell frequency across human samples. This identified three factors negatively correlated with Treg cell proportion that corresponded to aspects of the compensatory endothelial response to Treg cell depletion in the KP mouse model (Extended Data Fig. 10c). The latter included factors whose most associated genes were related to activated aCap (CAR4, CD36, IFNGR1, FAS, CX3CL1, TNFRSF11b, EDN1; factors 3 and 5; Fig. 7e), inflammation and hypoxia (VEGFB, PLAUR, SERPINE1, IL6, CXCL1, BCL3, PVR, IRF4, BATF3, TFP12; factors 4 and 5) and angiogenesis factors (factor 3). We used the sum of these three factors as a general Treg cell-responsive endothelial gene program to account for potential sample-specific, cell-type-specific or condition-specific effects that would separate a shared underlying biological program into separate factors during factorization (Extended Data Fig. 10d). Comparing this score to Treg cell proportion, we observed a clear negative correlation—stronger than any factor individually—across tumor samples (Fig. 7f), which suggested conserved Treg cell influence on this gene expression program. To further identify specific components of this shared Treg cell-responsive endothelial gene expression program, we compared the loadings of genes in the factors related to inflammation and hypoxia across species (factors 4 + 5 in human LuAd, factor 15 in KP mouse; Fig. 7g). This identified genes encoding key inflammatory mediators (IL6, CSF3, VCAM1, SELE, PTGS2) and a host of VEGF-induced genes in endothelial cells (RND1, ADAMTS1, ADAMTS4, ADAMTS9, AKAP12) as conserved members of expression programs induced in endothelial cells in Treg cell-poor TMEs across species.

Similar analyses of fibroblasts and myeloid cells also revealed corresponding Treg cell-dependent mouse and human factors. For example, human fibroblast factors 3, 5 and 22 corresponded to IC mouse fibroblast factors 21 and 22, with overlapping genes including IL6, IL1RL1, NFKB1, CCL2 and LIF (Extended Data Fig. 10e,f), while factor 9 (AP1 TF family members, KLF2/4, SOX9, HES1, IRF1) was negatively associated with Treg cell proportion. Additionally, high usage of conserved CSF3R monocyte factor 16 (CSF3R, PROK2, VCAN) in human ‘Treg cell-poor’ LuAd samples was consistent with the hypoxia, angiogenesis and NF-κB signaling related features (VEGFA, HIF1A, CEACAM1, NOTCH1, BCL3, BCL6) of this population in Treg cell-depleted mice (Extended Data Fig. 10g,h and Extended Data Fig. 5d). Notably, several human myeloid factors and corresponding mouse factors showed positive correlation with Treg cell presence, such as an SPP1/FOLR2 macrophage factor, a cell cycle factor and a C1Q+ macrophage factor (C1Q, antigen presentation-related genes), which included genes encoding known negative regulators of innate and adaptive immunity (CFH, CR1L, LAG3, PDCD1LG2, LILRB4, IL18BP; Extended Data Fig. 10i). Interestingly, we observed similarly pronounced downregulation of this gene program upon Treg cell depletion in both lung tumors and bleomycin-induced injury, suggesting that Treg cells within both niches sustain certain immunomodulatory myeloid cell states. Further analysis of correlation between conserved Treg cell-dependent human and mouse factors revealed a set of opposing TME programs (Extended Data Fig. 10j). One factor group in Treg cell-poor or Treg cell-deprived tumors included IL-1β/IL-18 signaling-related genes (IL18RAP, IL1RAP) expressed in angiogenic monocytes and tumor necrosis factor (TNF)/IL-1β-induced genes in fibroblast and endothelial cells involved in monocyte and neutrophil recruitment (CSF3, CXCL1, CXCL2, CXCL8, CCL2). The other, positively associated with Treg cell presence, featured immunomodulatory genes that inhibit IL-1β/IL-18 signaling (TMEM176B, IL18BP; see Supplementary Table 19 for KP/injury/LuAd factors). These findings suggest a conserved role of Treg cells in tuning transcriptional states of principal accessory cell types in the TME.

Combinatorial Treg cell depletion therapy

These results highlighted candidate compensatory pathways, whose targeting in combination with current clinical-stage intratumoral Treg cell depletion strategies32,33 could improve therapeutic efficacy. In this regard, the increased expression of VEGF pathway-related genes upon Treg cell deprivation was of particular interest suggesting that heightened VEGF signaling may ‘buffer’ the negative impact of Treg cell depletion on the tumor-supporting TME function and facilitate a rebound in the tumor progression. We tested the above possibility by investigating whether combining short-term Treg cell depletion with VEGF blockade can lead to an improved control of KP tumor progression. We transplanted KP adenocarcinomas into Foxp3GFP-DTR mice and administered them with DT and mouse VEGF-A neutralizing antibody (aVEGF) after tumors became macroscopically detectable (Fig. 8a). While Treg cell depletion and VEGF blockade alone could slow tumor progression, their combination had a markedly more pronounced therapeutic effect (Fig. 8b). Assessment of the overall survival rate, when mice were left untreated after the initial response and killed after rebound (tumor volume reached 1 cm3, maximum allowed by the institutional guidelines) showed that combination therapy improved survival in comparison to either monotherapy or untreated groups (Fig. 8b). While similarly increased numbers and activation level of tumoral T cells were observed in ‘DT + aVEGF’ and ‘DT-only’ in comparison to ‘aVEGF-only’ tumor samples on day 20 of transplantation (Fig. 8c), IFN-γ-producing CD4+ T cells and IFN-γ-producing and TNFα-producing CD8+ T cells were markedly increased in the combination treatment group as were monocyte numbers (Fig. 8d). Moreover, we observed further increases in tumor hypoxia and apoptosis upon combined Treg cell depletion and VEGF blockade in comparison to both monotherapeutic modalities and untreated control groups (Fig. 8e,f). Notably, KP tumors failed to respond to PD-1 blockade, which did not offer additional therapeutic benefits when combined with VEGF blockade in full agreement with a recent study of antiangiogenic, anti-PD-1 and chemotherapy in a KP lung cancer model34. Recent studies revealed high amounts of chemokine receptor CCR8 displayed by Treg cells in human cancers35,36 highlighting their depletion as a therapeutic strategy33,37,38. Thus, we examined the therapeutic potential of short-term VEGF blockade combined with antibody-mediated depletion of CCR8-expressing Treg cells, which represented only a fraction of intratumoral Treg cells in KP tumors (Fig. 8g). While CCR8 antibody treatment alone diminished tumor growth, a markedly more pronounced effect was observed when it was combined with VEGF blockade (Fig. 8h). Notably, this regimen was associated with a mere 15% decrease in overall population of tumor-associated Treg cells in the absence of their noticeable changes in the secondary lymphoid organs (Fig. 8i). Besides VEGF-A, whose neutralization was conducted as a proof-of-concept approach for the discovery of orthogonal combination therapy, we noted additional candidate compensatory pathways enriched in the Treg cell-poor or cell-depleted TME including the CCR2–CCL2 axis, inhibitors of which are currently tested as monotherapies or combination therapies of human cancers. To further test the utility of assessment of early TME responses to Treg cell depletion for identifying combinatory therapeutic modalities, we subjected KP tumor transplanted mice to a similar short-term treatment with CCR8 antibody and a selective CCR2 antagonist RS-504393 (CCR2i). The latter combination provided minimal additional therapeutic benefit in comparison to anti-CCR8 and CCR2i monotherapies contrary to aVEGF/CCR8 combination (Fig. 8j,k). These results suggest that CCR2 blockade and Treg cell depletion may converge on shared or partially overlapping TME states, whereas VEGF blockade offers an orthogonal intervention and highlights potential for discovery of orthogonal cancer therapies through single-cell and spatial analyses of early TME responses to acute perturbation.

Discussion

Successes in therapeutic targeting of PD-1 and CTLA-4 pathways in T lymphocytes are viewed as clinical evidence supporting the notion of cancer surveillance by cells of the adaptive immune system akin to that of pathogens. On the other hand, a growing realization of the important roles immune cells play in supporting normal tissue function, maintenance and repair suggests an alternative, even if not mutually exclusive view of tumor–immune interactions. Within the latter framework, the TME can be considered as a tissue-supporting multicellular network, which in response to cues emanating from cancerous cells supports their growth. In this regard, cancer represents a special state of a parenchymal cell, whose support by both immune and non-immune cells is guided by common, yet poorly understood principles of tissue organization. Treg cells suppress immune responses directed against self-antigens and foreign antigens to protect tissues from inflammation-associated loss of function1. Besides this indirect tissue-supporting functionality, Treg cells were also implicated in direct responses to injury and other forms of tissue damage through production of tissue repair factors16,39,40,41,42 suggesting that these functions of Treg cells are conserved. Furthermore, Treg cells were shown to support skin and hematopoietic stem cell niches15,43,44. Therefore, it is reasonable to assume that in human solid organ malignancies and in experimental mouse cancers Treg cells likely play similar dual roles supporting tumor growth-promoting accessory cell states.

Here, we showed that Treg cells have a profound impact on states of key accessory cells in a genetic autochthonous mouse model of NSCLC, in an experimental model of lung injury and in human LuAds. Using robust unsupervised data-driven computational analyses, we found that Treg cells support conserved gene programs—factors—across experimental models of lung cancer and injury, suggesting their role in coordinating broad, shared accessory cell functions that extend to various conditions and tissue states. The latter included human immunomodulatory C1Q+ (CFH, CR1L, LAG3, PDCD1LG2, LILRB4, IL18BP) and SPP1/FOLR2 macrophage factors and their mouse counterparts, which were positively associated with the Treg cell presence. A similar macrophage gene program is also reported to be enriched in NSCLC lesions45 and sustained by Treg cells in mouse models of melanoma and breast cancer46.

Our analysis of the distribution of activated cell types and gene expression programs with respect to their localization within and around tumor nodules showed high concordance of characteristic gene programs that were identified by scRNA-seq and ST analyses in situ. Notably, Treg cell depletion induced the IC response program localized to tumor nodule cores, while the IFN response program was most notable in the margins of tumor foci. The display of these programs by multiple cell types present within the same local niche suggests that they are elicited by common signals (‘signal niche’), for example, hypoxia response in the tumor nodule cores and a transient burst of IFN-γ produced by CD8+ T cells and NK cells concentrated in the tumor margins47. We also observed heterogeneity between Treg cell depletion responsive and nonresponsive tumor foci distinguished by the presence or paucity of the IC gene program. Interestingly, tumor nodules that failed to induce the IC gene program in response to Treg cell depletion expressed Sox9 in agreement with a recent study where upregulation of Sox9 in human LuAd conferred resistance to NK cells27.

Among the conserved gene programs negatively associated with Treg cell presence in mouse and human lung cancer, we noted the VEGF signaling pathway. This included increased expression of VEGF signaling-related genes in endothelial cells and increased expression of VEGF-A in myeloid and other cell types. This most likely reinforces the immunosuppressive TME providing support for tumor growth consistent with a recent report of tumor ischemia caused by the transient spike in intratumoral IFN-γ following CD25 antibody photoimmunotherapy-induced Treg cell depletion47. In addition, lung EC-derived VEGF was shown to specify development of CAR4hi endothelial cells and promote vascularization and tissue regeneration following injury48,49. VEGF has also been suggested to exert an immunomodulatory effect on cells of the innate and adaptive immune system50. Considering VEGF targeting being an approved therapy for some human cancers, combined VEGF-A and Treg cell targeting serves as a proof-of-concept for a rational combination therapy instructed by the new knowledge of TME transcriptional connectivity. While near complete loss of the Treg cell pool in KP-DTR mice led to systemic autoimmunity and inflammation, VEGF blockade coupled with CCR8 antibody-mediated depletion of intratumoral Treg cells showed impressive therapeutic efficacy with no adverse effects. The latter owes to the fact that CCR8 expression is selectively enriched in highly activated intratumoral Treg cell subsets in human and mouse malignancies35,36,38. Our observation that a combination of CCR8 antibody-mediated intratumoral Treg cell depletion with CCR2 blockade did not yield additional benefit in comparison to the corresponding monotherapies suggests that the latter either directly or indirectly converge on a shared regulatory node and highlights the utility of preclinical selection of combinatorial therapeutic strategies informed by scRNA-seq and ST analyses of early TME responses.