Conserved transcriptional connectivity of regulatory T cells in the tumor microenvironment informs n
Issuing time:2023-05-05 11:29
Abstract
While regulatory T (Treg) cells are traditionally viewed as professional suppressors of antigen presenting cells and effector T cells in both autoimmunity and cancer, recent findings of distinct Treg cell functions in tissue maintenance suggest that their regulatory purview extends to a wider range of cells and is broader than previously assumed. To elucidate tumoral Treg cell ‘connectivity’ to diverse tumor-supporting accessory cell types, we explored immediate early changes in their single-cell transcriptomes upon punctual Treg cell depletion in experimental lung cancer and injury-induced inflammation. Before any notable T cell activation and inflammation, fibroblasts, endothelial and myeloid cells exhibited pronounced changes in their gene expression in both cancer and injury settings. Factor analysis revealed shared Treg cell-dependent gene programs, foremost, prominent upregulation of VEGF and CCR2 signaling-related genes upon Treg cell deprivation in either setting, as well as in Treg cell-poor versus Treg cell-rich human lung adenocarcinomas. Accordingly, punctual Treg cell depletion combined with short-term VEGF blockade showed markedly improved control of PD-1 blockade-resistant lung adenocarcinoma progression in mice compared to the corresponding monotherapies, highlighting a promising factor-based querying approach to elucidating new rational combination treatments of solid organ cancers.
Main
Diverse stromal cell types found within the tumor microenvironment (TME) can support cancer initiation and progression by acting as accessory cells, yet their relationships and interdependencies remain poorly understood. Cells of the innate and adaptive immune system, when mobilized by immunotherapeutic agents, have been implicated in limiting cancer progression, yet some of the very same cell types can support tumor growth either directly or indirectly by facilitating tumor-promoting functions of other accessory cell types. Treg cells, expressing the transcription factor Foxp3, are highly enriched in human solid organ cancers and their experimental animal models, and at sites of inflammation and injury, where they exert both their essential immunosuppressive function and distinct tissue repair-promoting modalities1,2,3,4. Depletion of Treg cells results in restraint of tumor growth in numerous experimental cancer models5,6,7,8,9. Nevertheless, some tumors eventually progress after an initial response to Treg depletion5. The latter can be due to waning functionality of effector T cells due to negative regulation by co-receptors, foremost PD-1, expected to occur primarily in PD-1 blockade-responsive tumors expressing PD-L1. An alternative, yet not mutually exclusive, explanation, is that Treg cell depletion induces compensatory modulation of key accessory cell types in the TME, which may affect predominantly PD-1 nonresponsive cancers. Thus, early changes in diverse cellular components of the TME upon short-term Treg cell depletion may directly and indirectly impact its overall effect on tumor growth. Thus, we sought to elucidate the interplay between Treg cells and other cellular components of the TME by investigating early changes in their features upon Treg depletion in experimental cancer settings. Specifically, we wished to use a genetically engineered mouse model that is characterized by natural evolution of the TME, pronounced Treg cell presence, resistance to PD-1 blockade and close resemblance to human disease. Therefore, we used KrasLSL-G12D/WTTrp53fl/fl mice harboring a Foxp3GFP-DTR allele (KP-DTR), in which intratracheal infection with a Cre-expressing replication-deficient adenovirus induces lung adenocarcinoma (LuAd) formation. These mice offer a well-established model of non-small cell lung cancer (NSCLC) in humans, a disease where only some respond to PD-1/PD-L1 blockade-based therapies7,10,11,12. Our studies revealed that Treg cells profoundly affect the transcriptional programs of key accessory cells including endothelial cells, fibroblasts, monocytes and macrophages in the TME. Moreover, these Treg cell dependencies of the transcriptional states of accessory cells are largely conserved in human lung cancer.
Results
Early responses of tumor microenvironment cells to Treg cell depletion
To enable temporally controlled Treg cell depletion in KP adenocarcinomas, we generated KrasLSL-G12D/WTTrp53fl/flFoxp3GFP-DTR mice, in which all Treg cells express the diphtheria toxin receptor (DTR)13. We reasoned that since Treg cells are typically found in the tumor margins, early compensatory responses of key accessory cell types—tumor-associated fibroblasts, vascular endothelial cells (VECs) and lymphatic endothelial cells (LECs), and macrophages (Mac)—to Treg cell depletion likely precede effects on the tumor growth. Because the expansion of activated self-reactive T cells, observed 72–96 h after DT-mediated Treg cell depletion in Foxp3GFP-DTR mice, induces pronounced inflammatory responses13, we sought to minimize these confounding factors by analyzing early transcriptional responses of KP tumor cells, lung epithelial cells (ECs), VECs, LECs, macrophages and T cells 48 h following DT administration to tumor-bearing KP-DTR mice (Fig. 1a,b and Extended Data Fig. 1a,e). As expected, pronounced local and systemic T cell activation and inflammation, typically elicited by an extended Treg cell depletion regimen, were not observed (Fig. 1c and Supplementary Fig. 1c), and tumor volume was unaffected at this early time point (Fig. 1d) even though neutrophils were moderately increased (Extended Data Fig. 1c,k). Highly efficient tumoral Treg cell depletion in situ was confirmed by immunofluorescence (IF) microscopy of DT-treated as compared to control (Ctrl) mice, in which Treg cells were found mainly at the boundaries of tumor foci (Fig. 1e). Bulk RNA-sequencing (RNA-seq) analyses of cell subsets purified by fluorescence-activated cell sorting (FACS) from the lungs of DT-treated KP adenocarcinoma-bearing KP-DTR mice showed pronounced changes in gene expression in LECs, macrophages and fibroblasts, while T cells, which are considered the main targets of Treg cell suppression, changed the least (Fig. 1f and Supplementary Table 1). Among accessory cells, the most pronounced transcriptional responses were observed in fibroblasts, endothelial cells and CD11c+ myeloid cells, highlighting Treg cell ‘connectivity’ to these cell types in tumor-bearing lungs (Extended Data Fig. 1f–h and Supplementary Table 1). Importantly, DT-induced Treg cell ablation in tumor-free control KP-DTR mice resulted in minor, if any, changes in gene expression in all lung cell populations analyzed with the exception of VECs (Extended Data Fig. 1j,k). This was consistent with the predominantly intravascular localization of Treg cells in the lung of unchallenged mice in contrast to their heavy presence in the cancerous lung parenchyma (Extended Data Fig. 1f,g)14. These results suggest that the observed transcriptional changes in accessory cells in cancerous lungs are not due to a systemic response to Treg cell depletion. Next, we investigated whether shared groups of genes underwent modulation in different accessory cell types and observed correlated gene expression changes in endothelial cells and fibroblasts (Extended Data Fig. 1f, g). These included programs related to endothelial-to-mesenchymal transition (EndMT)-related genes (Id2, Itgav and Cxcl12), which were previously shown to be modulated by Treg cells in the hair follicles15, and inflammation-related genes (Il6, Ccl5, Acacb, Ccl22, Arg1 and Tnfrsf18), whose expression is affected by Treg cells in adipose tissue in the context of metabolic inflammation and muscle injury16,17 (Extended Data Fig. 1g). Cell-type-specific gene expression changes confirmed the shared gene expression changes were not due to sample impurities (Extended Data Fig. 1h). Considering the early transcriptional response of several TME cell types to Treg depletion, we assessed whether Treg cells were found in the proximity of these ‘first responders’ using IF analysis of tumor-bearing lungs. Indeed, GFP-DTR+ Treg cells were found in markedly closer proximity to Lyve-1+ LECs, GP38+ fibroblasts and F4/80+ macrophages within and near tumor nodules than in areas further away from tumor nodules in the same tumor-bearing lung (Fig. 1g,h). Collectively, we have shown that Treg cells are highly connected in the KP TME.
a, Schematic of the experimental design. b,c, Quantification of Treg (CD4+Foxp3+) one-tailed unpaired t-test P = 12.87, d.f. = 7 ****P < 0.0001 and Tcon (TCRβ+CD4+ and TCRβ+CD8+) cell populations; left, one-tailed t-test P = 0.3799, d.f. = 7, not significant (NS) P = 0.3576; right, one-tailed t-test P= 0.1925, d.f. = 7, NS P = 0.4264, in tumor-bearing lungs 48 h after diphtheria toxin (DT) or PBS (Ctrl) administration. d, Quantification of lung weight in tumor-free and tumor-bearing mice 48 h after DT-induced Treg cell depletion. One-way analysis of variance (ANOVA) followed by Sidak’s multiple-comparisons test. Tumor-free PBS versus tumor-free DT, P = 0.004037, d.f. = 10 NS P > 0.9999; tumor PBS versus tumor DT, P = 0.7450, d.f. = 10, NS P = 0.9787. e, Representative IF staining of Foxp3+ cells in tumor-bearing lungs of Ctrl and DT-treated mice. f, Numbers of upregulated (red) or downregulated (blue) DEGs (P < 0.05) 48 h after DT or PBS administration identified by bulk RNA-seq analysis of the indicated cell subsets. Fib, fibroblasts; Neu, neutrophils; Mac, macrophages; CD4 and CD8, effector CD4+ and CD8+ T cells. g, Representative IF staining of the indicated cell types. h, Quantification of distances between Treg cells and the indicated cell types. One-way ANOVA, alpha = 0.05, followed by Tukey’s multiple-comparison test Treg-Fib tumor-free zone versus tumor nodule, q = 8.041, d.f. = 2544 ****P < 0.0001. Treg-LEC tumor-free zone versus tumor nodule q = 10.08, d.f. = 2544, ****P < 0.0001, Treg-Mac versus tumor-free zone versus tumor nodule q = 17.79, d.f. = 2544, ****P < 0.0001. At least 200 cells were counted in each comparison. Three independent sections per mouse were analyzed. Three and four mice were used in each group in two independent experiments. Data are presented as the mean ± s.e.m. (b–d) (b and c) N = Ctrl-5, DT-4, (d) N = 3 tumor-free PBS, 3 tumor-free DT, 4 tumor PBS, 4 tumor DT. Data are presented as the mean ± s.e.m.
Single-cell analysis of tumoral Treg cell ‘connectivity’
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.
a,b, t-distributed stochastic neighbor embedding (t-SNE) plots (27,000 cells) representing cell populations from major cell lineages isolated from 48 h DT-treated or PBS-treated (Ctrl) tumor-bearing lungs (three mice per group) colored by cell type (a) and condition (b). c, A density plot showing the distribution of cells between experimental conditions. d,e, t-SNE plots (2,815 cells) representing distribution of the VEC populations colored by subtype (d) and condition (e). f, A density plot of endothelial cells showing the distribution of cells between experimental conditions. g, Graph of neighborhoods of endothelial cells computed using MiloR and t-SNE embedding. Each dot represents a neighborhood and is color coded by the false discovery rate (FDR)-corrected P value (alpha = 1) quantifying the significance of enrichment of DT cells compared to control in each neighborhood. The size of the dot represents the number of cells in the neighborhood. h, Swarm plot depicting the log fold change of differential cell-type abundance in DT-treated versus control samples in each neighborhood across different endothelial cell types. Each dot represents a neighborhood and is color coded by the FDR-corrected P value (alpha = 1) quantifying the significance of enrichment of DT cells compared to control in each neighborhood. A neighborhood is classified as a cell type if it comprises at least 80% of cells in the neighborhood, otherwise it is called ‘mixed’. i, Heat map showing average factor cell score in each cluster for each experimental condition in the VEC population. The scores were row normalized between 0 and 1. Each row represents a factor, and each column represents a cluster in a specific experimental condition. The clusters are grouped based on their phenotype. j, Gene expression heat maps showing the top 200 genes that correlated the most with the imputed activated VEC factor indicated (Methods). Each column represents a cell; cells are ordered based on their factor score (in ascending order from left to right) indicated by the green bar. Select genes of interest are noted on the right. b, e, h and i; Ctrl, PBS, gray; DT, red.
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).
We reasoned that gene programs most affected by Treg cell presence would have differential factor cell scores between the control and DT conditions. To evaluate this systematically, we computed the average cell score of every factor in each cluster for each condition (Fig. 2i) and identified those that have higher averages in DT compared to control. In the endothelial lineage, we identified four major gene programs that were robust to random initializations (Supplementary Fig. 1), were biologically relevant and had significantly differential cell scores (Mann–Whitney test; Methods) following Treg cell depletion compared to control in at least one of the endothelial cell subtypes (Fig. 2i). We then visualized expression of the genes with the highest factor loadings in the relevant cell subtype (Fig. 2j). We observed several notable patterns, including the inflammatory or activated capillary VEC factor (factor 3), a highly Treg cell-dependent factor characterized by cytokine/chemokine-, Notch and nuclear factor-κB (NF-κB) signaling-, and co-stimulation pathway-related gene expression (Fig. 2j; see Supplementary Table 3 for endothelial factors of interest). Other highlighted factors enriched following Treg cell depletion in the endothelial cell population included genes related to the NF-κB signaling pathway (Nfkbia, Rel, Hbegf), inflammation/hypoxia (Klf6, Serpine1, Plaur; factor 15) and vascularization (Vegfa, Thbd and Slco2a1; factor 8), and genes linked to transforming growth factor-beta-induced EndMT (Emp3, Timp1 and Tgm2; factor 17). Besides cancer, the latter process is induced in aberrant tissue remodeling and fibrosis23,24. These observations indicate that Treg cells impact specific features of certain endothelial cell subsets in the TME.
Notably, the observed transcriptomic perturbations were not unique to endothelial cells. The Treg cell depletion-induced gene programs related to interferon (IFN) response, inflammatory cytokines (ICs) and chemokines, STAT3 and interleukin (IL)-6 signaling appeared to be shared across accessory cell populations. The three most differentially expressed gene (DEG) programs observed in fibroblasts following Treg cell depletion included an inflammatory secretory phenotype (Ccl2, Hif1a, Rel, Cxcl1; factor 22), IFN response (Irf7, Ifit3, Isg15; factor 9) and ECM-related genes (Fbn1, Fn1, Lamc2, Notch2; factor 14; Extended Data Fig. 5a,b and Supplementary Table 4). On the other hand, several factors in monocytes (factors 2, 5, 7, 13, 17, 21 and 22) and macrophages (factors 15, 17 and 23) including IFN and hypoxia response emerged as differentially abundant (Extended Data Fig. 5c,d and Supplementary Table 5; for all significant factors across cell subsets, see Supplementary Table 6). These results suggested that Treg cell communication with various cells in the TME imparted both shared and distinct transcriptional features across and within specific cell populations in either a direct or indirect manner.
Treg cell dependency of accessory cell states in lung injury
To test whether the Treg cell ‘connectivity’ to key accessory cell types observed in lung cancer represents a generalizable facet of tissue organization, we examined perturbations of their transcriptional states upon identical short-term Treg cell depletion in a setting of bleomycin-induced fibrotic lung inflammation using scRNA-seq analysis (Fig. 3a,b and Extended Data Fig. 6a–d). Not only were all cell populations detected in tumor-bearing lungs also present in inflamed lungs, Treg depletion in this setting also generated similar transcriptional responses (Fig. 3a,b and Extended Data Fig. 6c,d). Independent analysis of the gene programs in the inflamed lung using scHPF (see Supplementary Table 7 for factor matrices) identified Treg cell depletion-associated endothelial factors (Fig. 3c and Supplementary Table 8). We correlated gene scores associated with each factor from lung tumors to lung injury to identify similarities. We found that the activated VEC factor in the lung injury (factor 15) correlated strongly (Pearson correlation > 0.70) with its counterpart in the tumor setting (factor 3), indicating that the same set of genes responded to the loss of Treg cells in both challenges. In fact, 72 of the top 200 genes associated with factor 3 specific to the tumor endothelial cell inflammatory capillary subset were shared with the top 200 genes associated with factor 15 specific to the same subset of cells in the injury model (Fig. 3d and Supplementary Table 9). Other endothelial cell factors, namely inflammation/hypoxia (factor 13), NF-κB signaling and EndMT (factor 12), that were observed in the inflamed lung upon short-term Treg cell depletion also correlated positively, even if weakly, with related tumor factors 15 and 8, respectively (Extended Data Fig. 6e). Consistently, factor analyses of other lineages revealed overlapping differential gene programs between tumor and injury models, including Treg cell depletion-induced gene programs in Arg1+ macrophages (Fig. 3e,f and Supplementary Table 10) and IC signatures in Col14a1 matrix fibroblasts. These findings suggested that Treg cell-dependent transcriptional programs are not limited to the TME and can be shared across pathological conditions.
a, t-SNE plots (24,592 cells) representing cell populations isolated from the lungs of mice administered with diphtheria toxin (DT) or PBS (Ctrl) for 48 h. Lung injury-induced inflammation was induced in both groups of mice upon bleomycin treatment 21 d before DT/PBS administration. The data represent analysis of three mice per group colored by cell type (left) and condition (middle), and a density of the distribution of cells between conditions (right). b, t-SNE embedding of endothelial cells isolated from Ctrl and DT after bleomycin administration color coded by cell type (left) or experimental condition (middle), and density plots of the distribution of endothelial cells between conditions (right). c, Heat map showing average factor cell score in each cell type for each experimental condition for endothelial cell subsets. The scores were row normalized between 0 and 1. Each row represents a factor, and each column represents an endothelial cell subset in a specific experimental condition. Factors of interest are highlighted by a red box. d, Heat map showing the 72 shared genes specific to activated VEC factor in both lung challenge models (Methods and Supplementary Table 9). Each column represents a cell; cells are ordered based on their factor score (in ascending order from left to right), indicated by the green bar. e, Heat map showing factor cell score across experimental conditions averaged over each myeloid cluster in each experimental condition for bleomycin-administered cells. The rows are factors and columns are clusters for each experimental condition. The clusters are grouped based on the cell type they are associated with. The heat map shows row-normalized scores from 0 to 1. The left color bar shows the average factor cell score. f, Heat maps showing the 54 shared genes between mouse lung tumor and injury-induced inflammation in the Arg1+ macrophage factor (tumor factor 23 corresponding to injury-induced inflammation factor 0; Supplementary Table 10). Each column is a cell; cells are ordered based on their factor score (in ascending order from left to right) indicated by the green bar. The treatment condition for each cell is indicated by gray for PBS and red for DT bars. Select genes of interest are shown.
Spatial distribution of Treg cell-dependent tumor microenvironment gene programs
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.
a, Tumor region identification in KP LuAd sections using Visium ST. The fraction of tumor cell RNA in each Visium spot (top right) was determined by BayesPrism deconvolution, binarized (bottom right; Methods), and compared to histological H&E images (left). b, Factor scores and Bonferroni-adjusted two-sided t-test P values differentially expressed factors between control and Treg cell-depleted conditions in ST. c,d, Representative tissue sections from control (left) or Treg cell-depleted (right) conditions. Tumor regions are outlined, and spots are colored by factor score. Scores represent IC (c; 18 genes) or IFN (d; 103 genes) gene programs shared across all lineages (Br, bronchi; A/V, artery/vein; LV, lymphatic vessel; IAs, immune cell aggregates). e, ST analysis revealed distinct signaling niches. Spots were assigned to niches based on thresholding a gamma distribution fitted to IC or IFN signaling module scores across all spots (Methods). f, Enrichment of cell-type RNA fractions in signaling niches. Adjusted empirical P value corresponds to the probability of obtaining the mean observed RNA fraction for that cell type (Methods). Fractions with adjusted P > 0.01 are not shown. In a and c–e, images are representative of, and analysis performed on (b and f), one of two serial sections for each of four samples (DT and Ctrl, two biological replicates each).
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, ST analysis of tumor states. BayesPrism deconvolution using additional labeled tumor cells from Yang et al.28 was performed to assign tumor-state-specific RNA fractions. Correspondence of regions with highlighted differential tumor states (middle) to H&E section is shown (right). Dashed lines denote regions with the indicated dominant tumor states (red, high plasticity; yellow, EMT; black, lung progenitor-like). b, Spots labeled by tumor-state cluster. In a and b, images are representative of, and analysis performed on (c and d), one of two serial sections for each of four samples (DT and Ctrl, two biological replicates each). c, Quantification of tumor lesion area types across Treg cell depletion and control conditions (left) or between tumors with or without detectable immune response in Treg cell-depleted condition (right; N = 85 lesion areas). d, Differential gene expression (two-sided Wilcoxon test Benjamini–Hochberg adjusted) of tumor spots in lesions with and without immune response to Treg cell depletion. e, Log-normalized expression of Sox9 and Pf4 (Cxcl4) in a representative tumor-bearing lung section after Treg cell depletion. Inset at top left indicates immune response status of tumor lesion areas.
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.
a, Schematic of the experimental design. b, t-SNE plot of all cells (82,991 total cells) from 25 primary human LuAd or local metastases labeled by lineage. c, t-SNE of T/NK cell lineage colored by unique molecular identifier (UMI) counts of Treg cell marker genes (maximum of two). d, Jaccard similarity between genes associated with mouse and human factors in tumor endothelial cells. Factors of interest with high correlation are highlighted by a green box. e, Conservation of activated VEC signature genes. Normalized gene loading (fraction of gene score across all factors) of genes within the mouse activated VEC signature across all human endothelial factors. Upper and lower notches of the box plot correspond to the 75th and 25th quartiles, respectively, and the middle notch corresponds to the median. Whiskers extend to the farthest data point no more than 1.5 times the interquartile range from the hinge, with outliers beyond that displayed as individual points. Select genes with high loadings of factors 3 and 5 are highlighted (N = 45 genes). f, Mean log2 sum of inflammation/angiogenesis associated human endothelial factor (3, 4 and 5) cell loadings plotted against log2 Treg cell proportion in each human sample. Spearman correlation estimate (R) and P value are listed. Trend line represents a linear model fit between the two and shading indicating the 95% confidence interval (N = 19 human samples). g, Normalized gene scores (fraction of gene scores across all factors) in orthologous genes between mouse and human inflammation/hypoxia factors. Genes significantly attributed to both human factors and mouse factors are highlighted as conserved. VEGF-induced genes in endothelial cells were derived from the CytoSig database.
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.
a, Schematic of the experimental design; s.c., subcutaneous. b, Tumor growth dynamics upon the indicated therapeutic interventions. The data represent mean values of tumor volume measurements (left). Adjusted P values for day 20 measurements: PBS-IgG versus DT-IgG P < 0.0001; PBS-IgG versus PBS-αVEGF P = 0.0004; PBS-IgG versus DT-αVEGF P < 0.0001; DT-IgG versus PBS-αVEGF P = 0.0328; DT-IgG versus DT-αVEGF P = 0.0109; PBS-αVEGF versus DT-αVEGF; P = 0.0005. Representative image of tumor volumes at day 20 (center). Kaplan–Meyer survival curves followed by log rank (Mantel–Cox) of KP tumor-bearing mice (right). The ‘survival’ time reflects the end point of the experiment when tumor volume in individual mice reached 1 cm3; adjusted P values: PBS-IgG versus DT-IgG P = 0.0012; PBS-IgG versus PBS-αVEGF P > 0.05 (NS), PBS-IgG versus DT-αVEGF P = 0.0078; DT-IgG versus PBS-αVEGF P > 0.05 (NS); DT-IgG versus DT-αVEGF P = 0.05; PBS-αVEGF versus DT-αVEGF P = 0.0186. c,d, Quantification of the indicated immune cell subsets and frequencies of activated (CD44hi CD62lo), proliferating (Ki67+) and IFN-γ-producing TCRβ+ CD4+ and TCRβ+ CD8+ cells in tumor samples shown in Fig. 8b in the indicated experimental groups of mice analyzed on day 20. e, Representative HIF1α and TUNEL staining of KP tumor sections. f, Quantification of HIF1α expression and apoptosis (TUNEL staining) in KP tumor sections; staining areas and signal intensity normalized by the total area and mean background intensity, respectively. 3–5 tumors from each experimental group were analyzed. (PBS-IgG N = 5; DT-IgG N = 4; PBS αVegf N = 3; DT αVegf N = 3) with four sections per individual tumor sample. Data represent the mean ± s.e.m. g, Proportion of intratumoral Treg cells on day 20 after KP tumor transplantation. Data represent the mean ± s.e.m. of one of two independent experiments; N = 8. h, Tumor growth dynamics upon the indicated therapeutic interventions. Gray arrows indicate days of neutralizing antibody administration. The data represent mean values of tumor volume measurements (left). Adjusted P values for day 20 measurements: IgG versus αCCR8 P < 0.0001; IgG versus αVEGF P < 0.0001; IgG versus αCCR8-αVEGF P < 0.0001; αCCR8 versus αVEGF P = 0.0434; αCCR8 versus αCCR8-αVEGF P = 0.0044; αVEGF versus αCCR8-αVEGF P < 0.0001. i, Quantification of proportion and absolute numbers of intratumoral and splenic Treg cells following treatment (left) and the corresponding Treg cell numbers in spleens in the treated animals (right). Data in h and i represent the mean ± s.e.m. of one of two independent experiments, IgG N = 10, CCR8 N = 10, αVegf N = 8, CCR8-αVegf N = 8. j, Tumor growth dynamics upon the indicated therapeutic interventions (left). Gray and black arrows indicate timing of neutralizing antibody and CCR2 inhibitor (CCR2i) administration, respectively. The data represent the mean ± s.e.m. values of tumor volume measurements. Adjusted P values of day 20 measurements: IgG versus αCCR8 P = 0.0009; IgG versus αVEGF P < 0.0001; IgG versus CCR2i P < 0.0001; IgG versus αCCR8-αVEGF P < 0.0001; IgG versus αCCR8-CCR2i P < 0.0001; αCCR8 versus αVEGF P = 0.9982; αCCR8 versus CCR2i P = 0.6138; αCCR8 versus αCCR8-αVEGF P < 0.0001; αCCR8 versus αCCR8-CCR2i P = 0.0041; αVEGF versus CCR2i P = 0.9551; αVEGF versus αCCR8-αVEGF P = 0.0003; αVEGF versus αCCR8-CCR2i P = 0.0363; CCR2i versus αCCR8-αVEGF P = 0.0018; CCR2i versus αCCR8-CCR2i P = 0.2271; αCCR8-αVEGF versus αCCR8-CCR2i P = 0.4530. Plots include data from two independent experiments combined with nine animals in each group in experiment 1 (IgG N = 9, αCCR8 N = 9, αVEGF N = 9, CCR2i N = 9, αCCR8 N = αVEGF-9, αCCR8-CCR2i N = 9) and 4–6 animals per group in experiment 2 (IgG N = 4; CCR2i N = 6; CCR8-CCR2i N = 6). k, Kaplan–Meyer survival curves followed by Log-rank (Mantel–Cox) of KP tumor-bearing mice. The ‘survival’ time reflects the end point of the experiment when tumor volume in individual mice reached 1 cm3. Adjusted P values: IgG versus αCCR8 ***P < 0.0001; IgG versus αVEGF ***P < 0.0001; IgG versus CCR2i ***P < 0.0001; IgG versus αCCR8-αVEGF ***P < 0.0001; IgG versus αCCR8-CCR2i ***P < 0.0001; αCCR8 versus αVEGF P = 0.5687 (NS); αCCR8 versus CCR2i P = 0.7411 (NS); αCCR8 versus αCCR8-αVEGF ***P = 0.0002; αCCR8 versus αCCR8-CCR2i P = 0.0342; αVEGF versus CCR2i P = 0.8054 (NS); αVEGF versus αCCR8-αVEGF ***P = 0.0006; αVEGF versus αCCR8-CCR2i P = 0.0666 (NS); CCR2i versus αCCR8-αVEGF ***P = 0.0003; CCR2i versus αCCR8-CCR2i P = 0.6749 (NS); αCCR8-αVEGF versus αCCR8-CCR2i *P = 0.0489. Plots include data from two independent experiments combined with 5–11 animals in each group in experiment 1 (IgG N = 9, αCCR8 N = 9; αVEGF N = 9; CCR2i N = 11; αCCR8-αVEGF N = 9; αCCR8-CCR2i N = 5) and 4–10 animals per group in experiment 2 (IgG N = 7; CCR2i N = 10; CCR8-CCR2i N = 4). In b–d, h and i, plots are representative of one of two experiments with 8–10 mice per group each, at day 20 after transplantation. Number of mice per group in b and c: PBS-IgG N = 10; DT-IgG N = 10; PBS-αVEGF N = 9; DT-αVEGF N = 9; number of mice per group in h and i: IgG N = 10; αCCR8 N = 10; αVEGF N = 8; αCCR8-αVEGF N = 8.
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.
Our study highlights a generalizable approach where perturbation of a given cell population in an engineered genetic cancer model enabled computational learning of its ‘connectivity’ and influence on the TME and other diseased tissue states, which could then be compared to the human clinical settings. A surfeit of secreted and cell surface molecules has been implicated in Treg cell-mediated immunosuppressive and tissue-supporting functions. However, none of these individual modalities can predominantly account for the bulk of these functionalities. Combinatorial targeting of these putative mediators will enable elucidation of the molecular mechanisms of the observed Treg cell dependencies in the TME. Our results suggest that Treg cells serve as an essential component of a complex network of accessory cells of both hematopoietic and non-hematopoietic origin. Shared perturbations in their transcriptional states observed across the three different settings imply that the identified interdependencies of Treg cells and other components of tissue-supporting cellular networks are conserved and can be exploited to develop new strategies for rational therapies of cancer and other diseases.