Lung tumor-infiltrating Treg have divergent transcriptional profiles and function linked to checkpoint blockade response

Abstract

Regulatory T cells (T_reg) are conventionally viewed as suppressors of endogenous and therapy-induced antitumor immunity; however, their role in modulating responses to immune checkpoint blockade (ICB) is unclear. In this study, we integrated single-cell RNA-seq/T cell receptor sequencing (TCRseq) of >73,000 tumor-infiltrating T_reg (TIL-T_reg) from anti-PD-1-treated and treatment-naive non-small cell lung cancers (NSCLC) with single-cell analysis of tumor-associated antigen (TAA)-specific T_reg derived from a murine tumor model. We identified 10 subsets of human TIL-T_reg, most of which have high concordance with murine TIL-T_reg subsets. Only one subset selectively expresses high levels of TNFRSF4 (OX40) and TNFRSF18 (GITR), whose engangement by cognate ligand mediated proliferative programs and NF-κB activation, as well as multiple genes involved in T_reg suppression, including LAG3. Functionally, the OX40^hiGITR^hi subset is the most highly suppressive ex vivo, and its higher representation among total TIL-T_reg correlated with resistance to PD-1 blockade. Unexpectedly, in the murine tumor model, we found that virtually all TIL-T_reg-expressing T cell receptors that are specific for TAA fully develop a distinct T_H1-like signature over a 2-week period after entry into the tumor, down-regulating FoxP3 and up-regulating expression of TBX21 (Tbet), IFNG, and certain proinflammatory granzymes. Transfer learning of a gene score from the murine TAA-specific T_H1-like T_reg subset to the human single-cell dataset revealed a highly analogous subcluster that was enriched in anti-PD-1-responding tumors. These findings demonstrate that TIL-T_reg partition into multiple distinct transcriptionally defined subsets with potentially opposing effects on ICB-induced antitumor immunity and suggest that TAA-specific TIL-T_reg may positively contribute to antitumor responses.

Fig. 1.. Single-cell transcriptomic profiling of T_reg in treatment-naïve and neoadjuvant anti–PD-1–treated NSCLCs.

Coupled scRNA-seq/TCR-seq was performed on T cells isolated from resected tumor (n = 15), adjacent NL (n = 12), TDLN (n = 3), and a resected brain metastasis (n = 1) from patients with NSCLC treated with two doses of neoadjuvant anti–PD-1 as well as resected tumor (n = 10) and paired adjacent NL (n = 8) from treatment-naïve patients with NSCLC. (A) Two-dimensional (2D) UMAP projection of the expression profiles of the 73,882 T_reg that passed QC. T_reg subsets, defined by 10 unique clusters, are annotated and marked by color code. (B) Relative expression [average log₂(fold change)] for top differential genes for each cluster is visualized on a heatmap. Three thousand cells (or all cells in the cluster if cluster size <3000 cells) were randomly sampled from each cluster for visualization. Differential expression tests for T_reg cell subsets were performed with Wilcoxon rank sum test. Genes with >0.25 log₂ fold changes, at least 25% expressed in tested groups, and genes with Bonferroni-corrected P values < 0.05 were regarded as differentially expressed. (C) The expression of canonical T_reg subset marker genes and cell subset selective genes was visualized in red scale using UMAP projection. (D) PCA and canonical correlation of pseudobulk gene expression for individual tumor (yellow, n = 25) and adjacent NL (dark blue, n = 20) samples. Canonical correlation with tissue type = 0.56, P < 0.001. (E) PCA and canonical correlation of pseudobulk gene expression for individual nonresponder (red, n = 9) and responder (blue, n = 6) tumors. Canonical correlation with response status = 0.36, P = 0.50.

Fig. 2.. Activated, OX40^hiGITR^hi T_reg are functionally suppressive and associate with nonresponse to PD-1 blockade.

(A) Functional analysis of OX40^loGITR^lo and OX40^hiGITR^hi TIL-T_reg–mediated suppression of conventional CD4⁺T cell (T_conv) proliferation. TIL-T_reg from patients MD017–0157, MD017–0092, MD017–0115, MD017–0124, and MD017–0116 were combined to ensure sufficient T_reg numbers for this experiment (n = 1). (B) Volcano plot showing differential expression between the Activated (1)/OX40^hiGITR^hi cluster (right) versus all other T_reg (left) when analyzing all tumor samples (n = 25). Each dot represents one gene. False discovery rate (FDR) < 0.05 is considered significant (blue/red dots). LAG3, TNFRSF18, and TNFRSF4 represent the top three most differentially expressed genes in the Activated (1)/OX40^hiGITR^hi T_reg (red dots). (C) Overlay of the Activated T_reg score on the TIL-T_reg UMAP for each response/treatment group. Red indicates higher expression; blue indicates lower expression. (D) The frequency of T_reg with a high T_reg activation score is shown for tumor (left) and adjacent NL (right) from nonresponders (NR; red), responders (R; blue), and treatment-naive patients (Tx-naive; green). Comparisons were performed at the individual patient level using Wilcoxon rank test. (E) Clonotype sharing pattern across T_reg subsets. The frequencies of T_reg TCR clones that were detected in at least two T_reg clusters were calculated and are shown on a heatmap. (F) Diffusion plot with RNA velocity for the Activated (1)/ OX40^hiGITR^hi and Activated (3) clusters (among which most clonotype sharing was observed). (G) Boxplots showing the relative proportion of Activated (1)/OX40^hi GITR^hi and Activated (3) clusters by R (blue), NR (red), and treatment-naïve status (green). (H) A violin plot shows TNFSF4 (OX40L) expression by neoantigen-specific (red) and flu-specific (blue) CD8 TILs from the same neoadjuvant-treated patients. (I) Violin plot comparing TNFSF4 (OX40L) expression by neoantigen-specific CD8⁺ TIL between nonresponding (red)and responding (blue) tumors. (J) Gene set enrichment analysis to evaluate differing biological functions of ex vivo OX40L, 41BBL, and GITRL agonism in sorted human TIL-T_reg. TIL-T_reg from patients MD043–011 and MD01–019 were combined to ensure sufficient T_reg numbers for this experiment (n = 1). NF-κB, nuclear factor κB.

Fig. 3.. Defining the scRNA-transcriptome associated with T_reg tumor antigen reactivity using a transgenic TCR mouse model.

Coupled scRNA-seq/TCR-seq was performed on T_reg isolated from GP-expressing MC38 (MC38-GP, n = 5 mice per experiment) or parental MC38 (MC38WT, n = 10 mice per experiment) tumors, tumor-draining inguinal lymph node, and spleens on day 14 of tumor growth. (A) 2D UMAP projection of the expression profiles of the 30,325 T_reg that passed QC. T_reg subsets, defined by five unique clusters, are annotated and marked by color code. (B) Relative expression of three to five differential genes for each cluster is visualized on a heatmap. (C) Comparison of the frequency of SMARTA FoxP3⁺ T_reg of total T_reg defined by expression of CD45.1, FoxP3, and the SMARTA TCR between MC38WT (blue) and MC38-GP (red) tumor-bearing mice in each tissue compartment. Means with SEM error bars are shown for each mouse (n = 5). P values were obtained using Mann-Whitney test. (D) Homology to human T_reg clusters is shown by the average expression of genes scores built on the top 20 differentially expressed genes for each murine cluster queried in human clusters. (E) Quantification of cluster designation for all SMARTA clones from WT and GP-expressing MC38. (F) Cluster localization for SMARTA clones is shown with SMARTA TCR (red) overlaid on full 2D UMAP (gray) of T_reg isolated from MC38-GP (left) or MC38WT (right). (G) Averaged log₂-transformed and library size–normalized FoxP3 expression values are shown. Visualized in red scale using full UMAP projection. (H) TBX21 (Tbet) expression is shown for TR-T_reg (blue) and all other T_reg (red) from the GP-expressing MC38 tumor after 14 days in vivo. (I) Five thousand TR-T_reg and non–TR-T_reg were sorted and cocultured ex vivo with 1000 dendritic cells loaded with the LCMV GP peptide for 48 hours. Supernatants were harvested and assayed for IFN-γ using the Meso Scale Discovery platform. Data are shown as IFN-γ production above background, defined as T_reg cocultured with dendritic cells without peptide. Mann-Whitney test was used to compare average IFN-γ production between peptide-loaded and peptide-unloaded conditions (P = 0.33; n = 2). (J) Violin plot showing the TR-T_reg score of adoptively transferred SMARTA TCR^pos T_reg before adoptive transfer (pre-AT) and at 1, 7, and 14 days of tumor residence. (K) Violin plot showing FoxP3 expression by adoptively transferred SMARTA TCR⁺ T_reg before adoptive transfer (pre-AT) and at 1, 7, and 14 days of tumor residence. P values obtained by Student’s t test. The P value of the trend (P_trend) was calculated for the TR-T_reg score and normalized Foxp3 expression over time by testing the linear regression slope with null and alternative hypotheses: H0: β1 = 0 (the slope is equal to zero), HA: β1 ≠ 0 (the slope is not equal to zero). N = 5 to 10 mice per time point. (L) Volcano plots showing differentially expressed genes of SMARTA TCR⁺ T_reg in day 7 versus day 1. X axis shows log2(fold change), and y axis shows −log₁₀(FDR). Differentially expressed genes higher at day 7 are shown in red, and differentially expressed genes higher in 1 day are shown in blue. (M) TSDR methylation profile of sorted SMARTA TCR⁺CD25^hiRFP⁺ T_reg from female SMARTA;RIP-GP mice compared with control values (low, medium, and high). Residue numbers denote individual CpG motifs in reference to the transcription initiation site of FoxP3. Data are an average of n = 2 experimental replicates with cells from n = 3 mice pooled per experiment. Bar colors represent 0 to 100% methylation.

Fig. 4.. Murine TR-T_reg gene signature defines an orthologous subset among human NSCLC TIL-T_reg that is enriched in anti–PD-1 R.

(A) Red scale overlay of the TR-T_reg gene score, with red indicating higher expression and blue indicating low expression. Red dotted line represents the UMAP region with highest expression. Refined clustering was performed on the T_H1-like/cytotoxic subset, and 2D UMAP projection of six unique subclusters (SC0-SC5) is visualized by UMAP and marked by color code. (B) Relative expression of the top 10 most differential genes for each subcluster is visualized on a heatmap. (C) The expression of biologically relevant genes from the TR-T_reg score is visualized in red scale. (D) Boxplots showing differences in expression of FoxP3 (P = 6.8 × 10⁻⁹), IL2RA (CD25; P = 1.9 × 10⁻¹), and TBX21 (Tbet, P = 0.00024) by TIL-T_reg in SC0 versus all other T_reg. Each dot represents an individual tumor sample per patient (n = 14). Comparisons were performed at the individual patient level using paired t test. Patient-averaged log₂-transformed and library size–normalized expression values are shown. (E) Cell density plots of the T_H1-like/cytotoxic T_reg subclusters stratified by response/treatment status. The TR-T_reg score^hi population is indicated with a dotted line. Increased density is represented by red scale, and decreased density is represented by green/blue. (F) Boxplots showing the frequency of SC0 in the tumors (left) and adjacent NLs (right) of NR (red) and responders (blue). Comparisons were performed at the individual patient level using Wilcoxon rank test. (G) Diffusion map with RNA velocity between SC0, Activated (1)/OX40^hiGITR^hi, and LN-homing clusters for NR (left) and R (right). (H) Cross-cluster sharing of T_reg TCR clonotypes detected in the SC0 subcluster with non-SC0 T_reg and T_conv clusters (as shown in fig. S1C). Red indicates a higher frequency of the clone in the relevant subcluster. T_cm, central memory T cells; T_em, effector memory T cells.