Immunostimulatory Cancer-Associated Fibroblast Subpopulations Can Predict Immunotherapy Response in Head and Neck Cancer

This article is featured in Highlights of This Issue, p. 1991

Anti-PD-1 immune checkpoint inhibitors (ICI) are currently the first-line therapy for recurrent/metastatic head and neck squamous cell carcinoma (HNSCC; refs. 1, 2). Yet, overall response rates can be as low as 20%, with increased responses in tumors with elevated PD-L1 expression and tumor-infiltrating T cells (1, 3, 4). While this is partly a T cell–dependent mechanism, there may be additional cellular subpopulations in the tumor microenvironment (TME) responsible for mediating response to ICI (5–8). Recent studies have suggested that cancer-associated fibroblasts (CAF) are associated with this resistance, but their role in T-cell immunomodulation is still unclear in the human TME (9–11).

In human breast cancer, four CAF subtypes (CAF-S1 to S4), were identified on the basis of the expression of six fibroblast markers—fibroblast activation protein (FAP), integrin β1 (CD29), α-smooth muscle actin (αSMA), fibroblast-specific protein-1 (FSP-1), platelet-derived growth factor receptor β (PDGFRβ), and caveolin-1 (CAV1; ref. 12). In their follow-up study, single-cell RNA sequencing (scRNA-seq) further divided CAF-S1 into eight subtypes with the majority of these subtypes linked to immunosuppression and resistance to immunotherapy (11). In contrast, only three molecularly and phenotypically distinct CAF subpopulations were identified in pancreatic cancer, based on spatial location and imputed function, as defined by cytokine and surface marker expression. These include inflammatory CAF (iCAF), myofibroblastic CAF (myCAF), and antigen-presenting CAF (13, 14). In head and neck, three CAF types were previously identified by scRNA-seq corresponding to myCAF and two undefined CAF subtypes (CAF1 and CAF2), but the functionality of these subtypes and their association with immunotherapy response remains unknown (15). In general, these various orthogonal approaches to cancer-specific CAF characterizations have not provided a concordant classification of this important stromal host cell type.

To assess whether CAF-related or other TME subpopulations can regulate clinical responses to nivolumab, we leveraged scRNA-seq from a neoadjuvant trial to longitudinally profile pretreatment and posttreatment human HNSCC and generate a dynamic atlas of the human HNSCC TME. Our novel bioinformatic approach uses the Virtual Inference of Protein-activity by Enriched Regulon (VIPER) algorithm (16, 17) to address limitations imposed by high noise and significant gene dropout rates in most scRNA-seq analysis platforms. Specifically, VIPER leverages knowledge of regulatory networks to allow full quantitative characterization of protein activity by assessing the enrichment of their transcriptional targets in differentially expressed genes. On average, the resulting protein activity profiles outperform antibody-based measurements and dramatically outperform gene expression–based analyses in terms of identifying and characterizing molecularly distinct TME subpopulations (18), thus enabling mechanistic dissection of the HNSCC microenvironment at hitherto unattained resolution. We present the results of these protein activity–based analyses on clinical biospecimens to generate a high-resolution atlas of the human HNSCC immune and stromal microenvironment under ICI pressures.

Biospecimens were harvested from a window of opportunity trial of patients with locally advanced HNSCC (oral cavity, oropharynx, larynx, hypopharynx) who were candidates for primary surgical intervention with curative intent (NCT03238365). This study was conducted in accordance with the Declaration of Helsinki guidelines and approved by the Vanderbilt Institutional Review Board (IRB #171883, IRB #030062). Informed written consent was obtained from all patients. All enrolled patients were treated with 1 month of 240 mg nivolumab every 2 weeks for two doses prior to definitive surgery (N = 50). Half of the patients received tadalafil, and no differences were noted in response rates between the two cohorts (19). Consented patients were required to have fresh prenivolumab biopsy as well as preimaging and postimaging. Metaclinical annotation using pathologic criteria was used to delineate paired subject specimens as responders versus nonresponders. For both pretreatment and posttreatment timepoints, fresh specimens were collected for frozen fixation, paraffin-embedded fixation, and processed for both bulk and single-cell transcriptomic sequencing. Because of insufficient RNA or data quality, only 32 of the 50 patients were used for bulk RNA sequencing (bulkRNA-seq) and scRNA-seq analyses.

Fresh HNSCC tumor specimens were collected in DMEM supplemented with streptomycin (200 μg/mL), penicillin (200 U/mL), and amphotericin B (250 μg/mL). Tumor specimens were minced to 2–4 mm sized pieces in separate 6-cm dishes and digested to single-cell suspension using the Miltenyi Biotec human tumor dissociation kit (Miltenyi Biotec, #130-095-929) on the Miltenyi gentleMACS Octo dissociator (Miltenyi Biotec, #130-096-427) following manufacturer's instructions. Dissociated cells were aliquoted for single-cell sequencing, flow cytometry analysis, or CAF culture.

Samples were processed to generate single-cell gene expression profiles (scRNA-seq) using the 10× Chromium 3′ Library and Gel Bead Kit (10× Genomics), following the manufacturer's user guide. After GelBead in-Emulsion reverse transcription reaction, 12–15 cycles of PCR amplification were performed to obtain cDNAs used for RNA-seq library generation. Libraries were prepared following the manufacturer's user guide and sequenced on the Illumina NovaSeq 6000 Sequencing System. Gene expression data were processed with “cellranger count” on the prebuilt human reference set of 30,727 genes to generate counts matrices. Cell Ranger performed default filtering for quality control, and produced for each sample a barcodes.tsv, genes.tsv, and matrix.mts file containing counts of transcripts for each sample, such that the expression of each gene is in terms of the number of unique molecular identifiers (UMI) tagged to cDNA molecules corresponding to that gene. These data were loaded into the R version 3.6.1 programming environment, where the publicly available Seurat package v3.0 was used to further quality-control filter cells to those with fewer than 10% mitochondrial RNA content, more than 1,500 unique UMI counts, and fewer than 15,000 unique UMI counts.

Gene expression UMI count matrices for each sample were normalized and scaled in R using the Seurat SCTransform command to perform a regularized negative binomial regression based on the 3,000 most variable genes. Scaled data from each patient were batch corrected by Seurat using the functions FindIntegrationAnchors and IntegrateData, with default parameters. The resulting dataset included 22,906 high-quality cells (mean UMI count 4802) across 4 patients, including both pretreatment and posttreatment timepoints for each patient (Patient1: 5,857 pretreatment, 7,360 posttreatment; Patient2: 4,938 pretreatment, 1,550 posttreatment; Patient3: 487 pretreatment, 1,741 posttreatment; Patient4: 401 pretreatment, 572 posttreatment). The batch-corrected dataset was projected into its first 50 principal components using the RunPCA function in Seurat, and further reduced into a two-dimensional visualization space using the RunUMAP function with method umap-learn and Pearson correlation as the distance metric between cells. The data were clustered by the Louvain algorithm with silhouette score resolution-optimization selecting the resolution with maximum bootstrapped silhouette score in the range of resolution from 0.01 to 1.0 incremented by 0.01 (18). This resulted in an initial coarse clustering on gene expression. Within each cluster metaCells were computed for downstream regulatory network inference by summing SCTransform-corrected template counts for the 10 nearest neighbors of each cell by Pearson correlation distance.

For each single cell, inference of cell type was performed using the SingleR package and the preloaded Blueprint-ENCODE reference, which includes normalized gene expression values for 259 bulkRNA-seq samples generated by Blueprint and ENCODE from 43 distinct cell types representing pure populations of stroma and immune cells (Martens and Stunnenberg, 2013; the ENCODE Consortium, 2012). The SingleR algorithm computes correlation between each individual cell and each of the 259 reference samples, and then assigns both a label of the cell type with highest average correlation to the individual cell and a P value computed by Wilcoxon signed-rank test of correlation to that cell type compared with all other cell types. Highest-frequency SingleR labels within each cluster among labels with P < 0.05 are projected into the Gene Expression Uniform Manifold Approximation and Projection (UMAP) space, such that localization of SingleR labels is highly concordant with the unsupervised clustering. Unsupervised Clusters determined by the resolution-optimized Louvain algorithm are therefore labeled as a particular cell type based on the most represented SingleR cell type label within that cluster.

Differential gene expression analysis between single-cell clusters throughout the article is computed by the MAST hurdle model, as implemented in the Seurat FindAllMarkers command, with a log-fold change threshold of 0.5 and minimum fractional expression threshold of 0.25, indicating that the resulting gene markers for each cluster are restricted to those with log fold change greater than 0 and nonzero expression in at least 25% of the cells in the cluster.

From the integrated dataset, metaCells were computed within each gene expression–inferred cluster by summing SCTransform-corrected template counts for the 10 nearest neighbors of each cell by Pearson correlation distance. A total of 200 metaCells per cluster were sampled to compute a regulatory network from each cluster. All regulatory networks were reverse engineered by the ARACNe algorithm. ARACNe was run with 100 bootstrap iterations using 1,785 transcription factors (genes annotated in gene ontology molecular function database as GO:0003700, “transcription factor activity,” or as GO:0003677, “DNA binding” and GO:0030528, “transcription regulator activity,” or as GO:0003677 and GO:0045449, “regulation of transcription”), 668 transcriptional cofactors (a manually curated list, not overlapping with the transcription factor list, built upon genes annotated as GO:0003712, “transcription cofactor activity”, or GO:0030528 or GO:0045449), 3,455 signaling pathway-related genes (annotated in GO biological process database as GO:0007165, “signal transduction” and in GO cellular component database as GO:0005622, “intracellular” or GO:0005886, “plasma membrane”), and 3,620 surface markers (annotated as GO:0005886 or as GO:0009986, “cell surface”). ARACNe is only run on these gene sets so as to limit protein activity inference to proteins with biologically meaningful downstream regulatory targets, and we do not apply ARACNe to infer regulatory networks for proteins with no known signaling or transcriptional activity for which protein activity may be difficult to biologically interpret. Parameters were set to zero DPI (data processing inequality) tolerance and MI (mutual information) P-value threshold of 10⁻⁸, computed by permuting the original dataset as a null model. Each gene list used to run ARACNe is available on github (18).

Protein activity was inferred for all cells by running the metaVIPER algorithm, using all cluster-specific ARACNe networks, on the SCTransform-scaled and Anchor-Integrated gene expression signature of single cells from each patient. Because the SCTransform Anchor-Integrated scaled gene expression signature is already normalized as an internal signature comparing all cells with the mean expression in the dataset, VIPER normalization parameter was set to “none.” The resulting VIPER matrix included 1,239 proteins with activity successfully inferred from ARACNe gene regulatory networks. VIPER-inferred protein activity matrix was loaded into a Seurat Object with CreateSeuratObject, then projected into its first 50 principal components using the RunPCA function in Seurat, and further reduced into a two-dimensional visualization space using the RunUMAP function with method umap-learn and Pearson correlation as the distance metric between cells. Clustering was performed by resolution-optimized Louvain algorithm, as for gene expression, and SingleR-inferred cell type labels were carried over to identify cluster-by-cluster cell type labels. Differential protein activity between clusters identified by resolution-optimized Louvain was computed using bootstrapped t test, run with 100 bootstraps, and top proteins for each cluster were ranked by P value. The entire pipeline is implemented as in Obradovic and colleagues (18). Cluster cell counts were normalized to a fraction of the total sample separately for each patient and separately for pretreatment and posttreatment samples, with differences in pretreatment versus posttreatment frequency distribution.

Fibroblast clusters including 5,414 cells from overall VIPER clustering of all cells were further isolated and subclustered, with differential protein activity and frequency pretreatment versus posttreatment compared as in the analysis of initial clustering for all cells. A proteomic gene set for each head and neck CAF (HNCAF) cluster was defined on the basis of proteins differentially upregulated in each cluster. In the dataset of clinical trial patients profiled by bulkRNA-seq that had been annotated with subsequent response (n = 9) or nonresponse (n = 19) to αPD-1 immunotherapy (19), we applied VIPER transformation using the single-cell ARACNe networks on z-score scaled log₁₀(TPM) counts from pretreatment bulkRNA-seq data, and computed a differential protein activity signature ranking proteins by most upregulated in responders to most downregulated in responders. Enrichment of each HNCAF cluster marker set in the VIPER-transformed signature of responders versus nonresponders from bulkRNA-seq was computed by gene set enrichment analysis (GSEA), with normalized enrichment score and P value determined by 1,000 random permutations of gene labels. Insufficient number of human papillomavirus–positive (HPV+) samples prevented CAF enrichment analysis in HPV+ versus human papillomavirus–negative (HPV−) groups.

Clinical association of HNCAF signatures with outcome was further tested in The Cancer Genome Atlas (TCGA) head and neck cancer cohort processed by ARACNe and VIPER as above. Sample-by-sample normalized enrichment scores were computed ranking VIPER-inferred protein activity in each patient sample from highest to lowest activity and then applying GSEA. Normalized enrichment scores for HNCAF cluster signatures were binarized to less than zero (low) or greater than zero (high), and Kaplan–Meier curve showing association with overall survival time was plotted along with the log-rank P value, such that HNCAF-0 enrichment is associated with improved overall survival (P = 0.0095, median survival time = 602 vs. 1,671 days) and HNCAF-1 enrichment is associated with worse overall survival (P = 0.0041, median survival time = 1,718 vs. 773 days). We further plotted the sample-by-sample enrichment of these HNCAF populations among different TCGA tumor types with high stromal involvement [head and neck squamous cell carcinoma (HNSC), pancreatic adenocarcinoma (PAAD), sarcoma (SARC), uterine carcinosarcoma (UCS), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), liver hepatocellular carcinoma (LIHC)] and plotted the distribution of these enrichment scores by tumor type to assess relative tumor-type specificity of the identified HNCAF signatures.

NanoString GeoMX digital spatial profiling (DSP) was further applied, profiling IO360 immune gene panel expression among three regions of interest (ROI) from 1 patient and four ROIs from another. Anti-CD8, anti-αSMA, anti-PanCK, and DAPI stains were used for morphology identification and ROIs were selected on the basis of high abundance of tumor (PanCK), cytolytic T cells (CD8), and fibroblasts (αSMA). ROIs were split into PanCK-positive and PanCK-negative components, with gene expression evaluated separately in each. To further assess spatial colocalization of HNCAF subtypes with functionally exhausted T cells, we applied segment-by-segment gene set enrichment of HNCAF-0 and HNCAF-1 markers as well as enrichment of a published T-cell exhaustion signature (20), and correlate normalized enrichment scores for these populations between spatial segments.

To assess phenotypic concordance between prior fibroblast categorizations, including CAF-S1/S2/S3/S4 subtypes described in the setting of breast cancer (11, 12) and iCAF/myCAF subtypes described in the setting of pancreatic cancer (13), we have performed pairwise gene set enrichment of fibroblast phenotype marker gene sets among our HNCAF clusters identified by scRNA-seq. For CAF-S1/S2/S3/S4 phenotype matching, we sorted S1/S2/S3/S4 cells by FACS using the gating strategy described by Costa and colleagues (12) and subsequently performed bulkRNA-seq of each sorted population, applied VIPER using single cell–derived ARACNe networks, and computed differential protein activity of each population against the mean to define population-specific signatures, with genes ranked from most differentially upregulated protein activity to most differentially downregulated protein activity in CAF-S1/S2/S3/S4. We then performed pairwise GSEA of HNCAF cluster marker gene sets (by VIPER protein activity) among CAF-S1/S2/S3/S4 gene signatures. We highlight the findings that CAF-S1 gene signature was significantly enriched for the gene sets of HNCAF-0 (NES = 7.43, P = 1.1 × 10⁻¹³), HNCAF-1 (NES = 6.54, P = 6 × 10⁻¹¹), and HNCAF-3 (NES = 6.24, P = 4.4 × 10⁻¹⁰), CAF-S2 gene signature was not significantly enriched for any HNCAF population, CAF-S3 signature was significantly enriched for HNCAF-4 gene set (NES = 3.09, P = 2 × 10⁻³), and CAF-S4 signature was significantly enriched for HNCAF-2 gene set (NES = 6.7, P = 2.2 × 10⁻¹¹). Published CAF-S1 subtype marker gene sets (11) were directly tested by GSEA for enrichment in each single cell, with resulting enrichment scores plotted by HNCAF cell cluster. Published iCAF/myCAF VIPER-inferred marker gene sets (13) were directly tested by GSEA for enrichment in each single cell, with resulting enrichment scores plotted by HNCAF cell cluster, such that cells in HNCAF-1 are enriched for iCAF signature and cells in HNCAF-2 are enriched for myCAF signature. This phenotypic classification scheme highlights the distinction between our HNCAF categorization observed from scRNA-seq and prior CAF classification paradigms.

Receptor–ligand interactions were inferred between coarse-grained cell types using 2,557 high-quality receptor–ligand interactions reported in the RIKEN FANTOM5 database (21). This list of receptor–ligand pairs was filtered to identify pairs where the ligand was significantly upregulated, at the gene expression level, in at least one subpopulation, across patients, and the corresponding receptor was significantly activated in another subtype, based on VIPER protein activity analysis, as proposed in (18). We further filtered these interactions to those detected in CAFs and plotted the number of unique receptor–ligand interaction pairs inferred between fibroblasts and each other detected subpopulation.

Fresh HNSCC tumor specimens were processed to single-cell suspension as described above. For HNCAF-0/3, tumor single-cell suspension was cultured in DMEM supplemented with 10% FBS, streptomycin (100 µg/mL), and penicillin (100 U/mL) for 2–3 weeks at 37°C until fibroblasts grew out. For HNCAF-1, tumor single-cell suspension was cultured in pericyte medium (ScienCell #1201) supplemented with 2% FBS, streptomycin (100 µg/mL), and penicillin (100 U/mL) for 2–3 weeks at 37°C until fibroblasts grew out. To verify CAF identity, RNA was isolated from CAF lysates using TRIzol (Invitrogen #10296010) and sent for bulkRNA-seq. GSEAs for the HNCAF subtype protein activity signatures were then performed on the bulk sequencing data, along with inference of cell type proportions by CIBERSORTx. Fibroblasts were passaged when cultures reached approximately 80% confluence and all experiments were performed with CAF under 10 passages.

CD3⁺ T lymphocytes were isolated from the peripheral blood of healthy human donors. Peripheral blood mononuclear cells (PBMC) were isolated using Ficoll-Paque Plus, following manufacturer's instructions. CD3⁺ T cells were isolated from PBMCs using magnetic bead sort with the MojoSort Human CD3 T-Cell Isolation Kit (BioLegend #480022) according to manufacturer's instructions. For isolation of CD3⁺ tumor-infiltrating lymphocytes (TIL), fresh HNSCC tumor specimens were processed to single-cell suspension as described above. CD3⁺ TILs were isolated from the tumor single-cell suspension using magnetic bead sort with the MojoSort Human CD3 T-Cell Isolation Kit.

A total of 25,000 primary CAF were plated in DMEM supplemented with 10% FBS in 96-well plates. After CAF had attached to the plate, 50,000 CD3⁺ T cells were added to the CAF with or without CD3/CD28 activation beads (Gibco # 11131D) and cocultured at 37°C for 5–7 days with or without 20 ng/mL TGFβ. Media was renewed on days 3 and 5. Cocultures with TILs were only cultured for 3 days to preserve TIL viability. Following incubation, T cells were harvested and stained with Live/Dead Aqua (1:1,600) for 15 minutes in PBS. Cells were then washed and stained for 15 minutes with an antibody cocktail containing anti-CD4-APC/Fire 810 (1:1,000, SK3), anti-CD8-BB515 (1:200, RPA-T8), anti-PD-1-BV421 (1:100, EH12.2H7), anti-TIM-3-BV786 (1:100, F38-2E2), anti-NKG2A-PE (1:200, S19004C), and anti-CD103 (1:400, Ber-ACT8). Cells were then washed, fixed, permeabilized, and stained with an intracellular antibody cocktail containing anti-Perforin-PerCP/Cy5.5 (1:100, B-D48) and anti-Granzyme B-Alexa Fluor 700 (1:100, QA16A02). Cells were subsequently analyzed by flow cytometry using the Cytek Aurora.

A total of 100,000 primary CAF were plated in DMEM supplemented with 10% FBS in the lower chamber of the transwell (0.4 µm pore size, Corning Polycarbonate Membrane Transwells #3401). A total of 200,000 CD3⁺ T cells were plated in DMEM supplemented with 10% FBS in the upper chamber of the transwell. Cells were incubated at 37°C for 7 days. Media was renewed on days 3 and 5. Following incubation, T cells were stained and analyzed by flow cytometry using the Cytek Aurora as described above.

The level of IFNγ in cell culture supernatants was measured using an ELISA MAX Deluxe kit (BioLegend #430104) following manufacturer's instructions. Supernatants were collected from CAF T-cell cocultures as described above.

All quantitative and statistical analyses were performed using the R computational environment and packages described above with the exception of CAF composition and coculture experiments. Statistical analyses of these assays were performed using Prism 9 software (GraphPad). Differential gene expression was assessed at the single-cell level by the MAST single-cell statistical framework as implemented in Seurat v3 (22), and differential VIPER activity was assessed by t test, each with Benjamini–Hochberg multiple-testing correction. Comparisons of cell frequencies were performed by nonparametric Wilcoxon rank-sum test, and survival analyses were performed by log-rank test. In all cases, statistical significance was defined as an adjusted P value less than 0.05. Details of all statistical tests used can be found in the corresponding figure legends.

All sequencing data are deposited and publicly available in a Mendeley data repository (doi:10.17632/yk8wj7xgdg.1) as well as the Gene Expression Omnibus database (GSE195832).

Longitudinal scRNA-seq of patient tumors, pretreatment and posttreatment with nivolumab, and gene expression clustering with Seurat revealed 12 broadly distinct cellular populations, consistently expressed across all the tumors sequenced (Supplementary Fig. S1A and S1B). To achieve higher resolution of cellular subpopulation characterization, we used scRNA-seq profiles from each cluster to infer subpopulation-specific gene regulatory networks with the ARACNe algorithm (23), followed by protein activity analysis using VIPER. Protein activity–based reclustering identified two additional, previously undetected clusters for a total of 14 distinct cellular populations which were also consistently expressed across all 4 patients (Fig. 1A). To visualize key differences between these cellular populations, we generated a heatmap showing the activity of the five proteins most differentially activated in each cluster (Fig. 1B). The full list of differentially active proteins (P < 0.05, FDR corrected) and differentially expressed genes (P < 0.05, FDR corrected) are included in the Supplementary Data (Supplementary Table S1). We first assessed the ability of these VIPER-generated populations to accurately respond to expected treatment-induced changes. As expected, both gene expression and protein activity analyses revealed increased T cells and IFNγ protein activity following nivolumab treatment [Fig. 1C (clusters 1 and 5); Supplementary Fig. S1C and S1D]. When we interrogated changes in the abundance of other cell populations, VIPER clustering revealed heterogeneity among fibroblast cells not discoverable from gene expression clustering alone, with two clusters (clusters 4 and 9) presenting highly statistically significant posttreatment cellular fraction increase (Fig. 1C). Cell lineage inference, using SingleR (24), identified both clusters as fibroblasts, suggesting that PD-1–targeted immunotherapy in head and neck cancer was associated with CAF upregulation (Fig. 1D). Accordingly, imputed receptor-ligand interactions between cell types (18) suggested strong interplay between CAF and CD8 T cells (Supplementary Fig. S1E). Furthermore, two additional clusters (clusters 6 and 7), also characterized as fibroblasts by SingleR (Fig. 1D), did not show significant fractional representation differences following immunotherapy (Fig. 1C), thus suggesting the existence of functionally distinct CAF subpopulations within the HNSCC TME.

To further evaluate functional differences between the distinct CAF subpopulations in the HNSCC TME, we performed protein activity–based subclustering focusing only on fibroblast cells using ARACNe and VIPER. The analysis identified five molecularly distinct CAF clusters, preliminarily termed HNCAF-0 to HNCAF-4 (Fig. 2A). Importantly, gene expression–based subclustering of fibroblasts only identified two distinct CAF populations corresponding to the two fibroblast populations seen in Supplementary Fig. S1A (clusters 3 and 6). As expected, these two clusters phenotypically match the two CAF populations previously identified in HNSCC such that cluster 3 corresponds to CAF1 and cluster 6 corresponds to CAF2 (Supplementary Fig. S1F; ref. 15). Among the five HNCAF populations identified by protein activity–based clustering, cell fractional representation increased for HNCAF-0 and HNCAF-3, decreased for HNCAF-1 and HNCAF-2, and was unaffected for HNCAF-4 (Fig. 2B). The top 10 most differentially active proteins, presented as a ranked list of differentially active transcription factors and signaling molecules, in each of the five clusters highlight their potential functional properties (Fig. 2C) to help define the molecular biology of each HNCAF phenotype (Supplementary Table S2). To assess the associations of each HNCAF subpopulation with clinical response to αPD-1 immunotherapy, we used the HNCAF molecular signatures to analyze the bulkRNA-seq profiles from the 28-patient parental cohort annotated with clinical response to nivolumab. For this purpose, we first used VIPER to generate protein activity profiles from each bulk profile, using fibroblast-specific regulatory networks generated at the single-cell level, and then evaluated the enrichment of the most differentially active proteins in each HNCAF subpopulation (marker protein sets) among proteins differentially active in responders versus nonresponders. The analysis revealed statistically significant enrichment of HNCAF-0 and HNCAF-3 marker genes in pretreatment samples of patients who subsequently responded to αPD-1 immunotherapy (Fig. 2D). This result indicates that the HNCAF-0 and HNCAF-3 populations, which also expand following nivolumab treatment, may be predictive of clinical response in human HNSCC patients. In contrast, HNCAF-1, HNCAF-2, and HNCAF-4 cells did not expand following therapy and their markers were not significantly enriched in responders versus nonresponders.

Because of scant literature on CAF in human HNSCC, we next quantified the extent of CAF infiltration from surgical HNSCC specimens using flow cytometry. CAF abundance—as defined by CD45⁻ EpCAM⁻ CD31⁻—ranged between 12% and 58% of the total live cells (Fig. 3A). Having confirmed significant abundance of CAF in human HNSCC, we proceeded to assess the presence of novel HNCAF subpopulations predicted by the VIPER analysis. Distinct CAF subpopulations termed CAF-S1, CAF-S2, CAF-S3, and CAF-S4 have been previously identified in breast cancer based on the expression of CD29 and FAP by flow cytometry (12). Kieffer and colleagues showed that CAF-S1 can be found in HNSCC but the presence of other CAF-S populations in HNSCC remains unknown (11). Hence this protein-based framework was initially used to assess how the VIPER-imputed HNCAF align with breast cancer CAF scheme. Following the same gating strategy employed by Costa and colleagues (Supplementary Fig. S2), we confirmed the presence of all four breast cancer CAF-S populations in HNSCC (Fig. 3B; ref. 12). Interestingly, CAF-S1 and CAF-S2 were most abundant, while CAF-S3 and CAF-S4 abundance was quite minimal (Fig. 3C). We then sorted CAF-S1 to S4 cells from HNSCC tumors (Supplementary Fig. S2) and performed bulkRNA-seq of each subpopulation to assign these sorted cells to the VIPER-generated HNCAF populations. Pairwise GSEA of the HNCAF protein activity signatures in the bulk transcriptome indicated that the gene sets representative of HNCAF-0, HNCAF-1, and HNCAF-3 were all enriched in the same breast subtype (CAF-S1), while HNCAF-2 and HNCAF-4 were both enriched in CAF-S4 with HNCAF-4 also enriched in CAF-S3 (Supplementary Fig. S3A). CAF-S2—primarily defined as double negative for FAP and CD29—did not significantly align with any HNCAF. These analyses importantly showed that VIPER-clustered HNCAF provide much greater resolution of functionally distinct CAF phenotypes compared with the flow-based CAF-S1/S2/S3/S4 framework. Specifically, the CAF-S1 subtype matched three distinct HNCAF subtypes, which have opposing association with clinical responses to immunotherapy. In addition, the HNCAF subtypes did not clearly correlate with any of the CAF-S1 subclusters later identified by gene expression in breast cancer apart from ecm-myCAF (Supplementary Fig. S3B). HNCAF-0, 1, 2, and 3 were all significantly enriched for the ecm-myCAF signature, with HNCAF-0 and HNCAF-3 being more enriched than HNCAF-1 and HNCAF-2. However, ecm-myCAF are associated with immunosuppression and resistance to immunotherapy leading us to believe that HNCAF-0 and HNCAF-3 differ from ecm-myCAF in terms of functionality (11).

We also tested for concordance of our classification schema with previously defined gene set markers of iCAF and myCAF, as first described in pancreatic cancer (13). Cell-by-cell enrichment of iCAF and myCAF gene signatures revealed an enrichment for the iCAF signature in HNCAF-1 cells and for myCAF in HNCAF-2 cells (Supplementary Fig. S3C and S3D). Correlations between our HNCAF populations and the breast or pancreatic CAF phenotypes is summarized in the classification scheme shown in Fig. 3D. While presenting some similarity to CAF-S1 cells, we conclude that HNCAF-0 and HNCAF-3 represent novel, molecularly distinct fibroblast subpopulations, potentially unique to head and neck cancer and predictive of patient outcome (Fig. 2D), which do not completely match the iCAF/myCAF classification. Furthermore, in contrast to the HNCAF-0 and HNCAF-3 subtypes, iCAF/myCAF and CAF-S1/S2/S3/S4 are not significantly enriched in responder cohorts, suggesting that these classification schemes do not accurately depict CAF function in HNSCC (Supplementary Fig. S4). Cumulatively, these data underscore the greater resolution of HNCAF from our VIPER analysis compared with previous efforts, and more critically, the HNCAF populations identified through VIPER allow prognostic correlations of CAF cells in HNSCC.

To evaluate the prognostic relevance of the CAF populations identified in a setting free from external immunotherapeutic pressures, we quantified the enrichment of HNCAF protein activity signatures in TCGA HNSCC cohort (Fig. 4; Supplementary Fig. S5). GSEA (25), on a patient-by-patient basis, revealed significant enrichment of the HNCAF-0 signature in patients with better overall survival in TCGA (Fig. 4A), suggesting that these cells may not only be important regulators of immunotherapy response but may also play a key role in mounting clinically relevant, endogenous immune responses against HNSCC. In contrast, the HNCAF-1 protein activity signature was associated with early worse overall survival when there were sufficient number of patients (Fig. 4B). Prognosis in HNSCC has been associated with higher TIL level in TCGA cohort (26), and these results were consistent with the differential immunomodulatory functional roles of distinct VIPER-derived HNCAF cells.

Prompted by these intriguing clinical findings (Fig. 2D), we studied the potential interactions of HNCAF-0 and HNCAF-3 cells with other TME subpopulations. Interactome analysis showed that HNCAF have more receptor–ligand interactions with CD8 T cells than with any other cell subtype in the TME (see Materials and Methods; Supplementary Fig. S1E). Therefore, we next interrogated the relationship between HNCAF-0 and human CD8 T cells in situ and in vitro. DSP (NanoString) of immune-related transcripts and protein markers was performed on HNSCC tissue from patients prior to nivolumab treatment. We first analyzed global CAF infiltration pattern in HNSCC tissue using αSMA, αCD8, and αcytokeratin antibodies. These multiplexed immunofluorescent images exhibited colocalization of CAF with CD8⁺ T cells in the stromal compartment (Fig. 5A, white arrow). Reliable validated antibodies for each of the VIPER-derived HNCAF subpopulations are currently unavailable, preventing histologic analysis of each HNCAF population at this time.

To test whether HNCAF-0 cells can directly affect the biology of the CD8 T cells, we performed in vitro coculture assays with HNCAF harvested from surgical resection mixed with either naïve T cells or tumor-infiltrating T cells. Primary fibroblasts enriched for HNCAF-0/3 were isolated from human HNSCC samples and their transcriptional identity was verified by RNA-seq and protein activity analysis (Supplementary Fig. S6). Because of the phenotypic similarity between HNCAF-0 and HNCAF-3, we were unable to enrich solely for HNCAF-0 and proceeded with a heterogeneous population enriched for both HNCAF-0 and HNCAF-3. When cocultured with T cells isolated from PBMCs of healthy donors, HNCAF-0/3 cells reduced the PD-1⁺TIM-3⁺ exhaustion phenotype among exogenously activated CD8 T cells and increased the CD103⁺NKG2A⁺ tissue-resident memory (Trm) phenotype, as well as their cytolytic function, based on Perforin and granzyme B expression (Fig. 5B). It is important to note that reduced exhaustion and increased Trm phenotypes caused by HNCAF-0/3 was not due to PD-1:PD-L1 signaling (Supplementary Fig. S7A and S7B). Transwell coculture assays revealed that HNCAF-0/3–mediated T-cell activation increase and induction of Trm phenotypes—but not T-cell exhaustion phenotype mitigation—depends on cell-to-cell contact (Fig. 5C). In addition, coculture of HNCAF-0/3 cells with tumor-infiltrating T cells isolated directly from human HNSCC specimens resulted in an equivalent increase in Trm phenotype and cytotoxicity among CD8 T cells, but HNCAF-0/3 cells could not rescue the exhaustion phenotype of terminally exhausted tumor-infiltrating T cells (Fig. 5D). However, the duration of coculture for these experiments was shorter than with PBMCs due to issues with tumor-infiltrating T-cell viability and downregulation of PD-1⁺TIM-3⁺ exhaustion phenotypes likely requires longer interaction times. Regardless, HNCAF-0/3 cells strongly promoted production of the activation markers, Perforin, granzyme B, and IFNγ, in tumor-infiltrating T cells (Fig. 5D and E), suggesting that HNCAF-0/3 may prevent terminal exhaustion in early activated T cells as well as reinvigorate already exhausted T cells in the TME. Notably, we found that HNCAF-0/3 completely rescued TGFβ-mediated PD-1/TIM-3 induction in culture, without inhibiting total TGFβ signaling, as defined by CD103 induction (Fig. 5F).

To evaluate CAF influences on T-cell exhaustion in situ without a validated antibody for each HNCAF population, we leveraged the DSP data to evaluate colocalization of HNCAF-0 and HNCAF-1 protein activity signatures in regions enriched for the T-cell functional exhaustion signature. Indeed, the HNCAF-1 signature enrichment significantly correlated with increased T-cell exhaustion signature enrichment (r = 0.94, P = 0.0014). In sharp contrast, the HNCAF-0 signature was not significantly associated with the T-cell exhaustion signature in the TME region of interest (Fig. 5G and H). Given the association of HNCAF-1 cells with an immunosuppressive environment, we aimed to evaluate the direct impacts of isolated HNCAF-1 cells on T-cell phenotypes in coculture, as performed for HNCAF-0/3 cells. However, despite repeated experiments, T cells cocultured with HNCAF-1 rapidly died, leaving an insufficient number of viable cells for further analyses (Fig. 5I; Supplementary Fig. S7C–S7E). T-cell death induction was not observed when T cells were cultured in isolation or with HNCAF-0/3 cells, suggesting HNCAF-1 mediates accelerated T-cell apoptosis ex vivo.

To test for the potential generalizability of these HNCAF subpopulations, we next performed enrichment of HNCAF protein activity signatures across TCGA, by tumor type, focusing on tumors with high stromal cell content. Enrichment analyses revealed that HNCAF-0 enrichment is relatively specific to HNSCC while HNCAF-1 enrichment is more broadly observed (Supplementary Fig. S8A and S8B). Intriguingly, HNCAF-1 enrichment is highest in pancreatic adenocarcinoma, which is known to be unresponsive to PD-1–based immunotherapy (Supplementary Fig. S8B). HNCAF-1 phenotypically matches the previously defined iCAF population from pancreatic cancer (Fig. 3D).

To externally validate our HNCAF-0/3's potential for clinical response prediction, we tested for enrichment of protein activity signatures in another cohort of patients with HNSCC treated with αPD-1 immunotherapy, pembrolizumab (27). The analysis revealed statistically significant enrichment of HNCAF-0 and HNCAF-3 marker genes in pretreatment samples of patients who subsequently responded to pembrolizumab, validating the association of HNCAF-0/3 with immunotherapy response in an independent cohort (Fig. 6A). Although not as enriched as HNCAF-0 and HNCAF-3, analysis of this cohort also revealed a significant enrichment of HNCAF-2 in responders pretreatment (Fig. 6A). HNCAF-3, which showed the strongest overall enrichment in responders compared with nonresponders (Fig. 6A), also showed the highest predictive area under receiver operating characteristic (AUROC), with an AUROC of 0.8 (Fig. 6B), an improvement as compared with reported AUROC of PD-L1 tumor proportion score (AUROC = 0.66) and PD-L1 combined positive score (AUROC = 0.61; ref. 28).

In this study, we used protein activity profiles, as measured by the VIPER algorithm analysis of a longitudinal single-cell transcriptomics HNSCC dataset, to identify five molecularly distinct CAF subtypes. We took advantage of the longitudinally harvested biospecimens from a neoadjuvant clinical trial to show that two subtypes, HNCAF-0 and HNCAF-3, are predictive of favorable clinical responses to PD-1 checkpoint blockade therapy. Moreover, we discovered HNCAF-0/3 cells have an immunostimulatory effect on CD8 T cells while HNCAF-1 cells are associated with immunosuppression. From a functional perspective, we have shown that HNCAF-0/3 fibroblasts prevent induction of an exhausted T-cell phenotype and are associated with increased CD8 T-cell cytotoxicity. HNCAF-1 fibroblasts, in contrast, correlate with increased T-cell exhaustion in situ and T-cell apoptosis in vitro, suggesting contrasting roles for these CAF subtypes. Functional significance of HNCAF-2 and HNCAF-4 subtypes has yet to be defined.

Protein activity–based clustering with VIPER has previously been shown to outperform gene expression–based clustering and we have corroborated this with the identification of five HNCAF subtypes compared to two with gene expression–based clustering (18). Although we are not the first to identify CAF in HNSCC by single-cell sequencing, we were able to achieve greater resolution with VIPER allowing for a more granular picture of the CAF composition in HNSCC (11, 15). We were also able to integrate the major subclasses of CAF identified in breast cancer and pancreatic cancer into our HNSCC CAF data. We showed through GSEA that CAF-S1, identified in breast cancer and in our HNSCC samples as CD29⁺FAP⁺ fibroblasts, correspond to three of the HNCAF groups we identified, HNCAF-0, HNCAF-1, and HNCAF-3 (12). Kieffer and colleagues further dissected CAF-S1 in breast cancer and identified eight subgroups but only focused on the five most abundant groups (11). We performed GSEA of these five subgroups but were only able to identify one with significant enrichment in our HNCAF subpopulations: ecm-myCAF. We did find that the ecm-myCAF signature was significantly enriched on four of our five HNCAF groups. The discrepancies between our subclustering of CAF-S1 and Kieffer and colleagues may be due to tumor-specific differences in the CAF heterogeneity between breast and HNSCC as well as the use of a protein activity–based algorithm in contrast to gene expression. In addition, we have found a CAF population (HNCAF-0) that is relatively specific to HNSCC amid the other CAF types seen across different cancers. At this time, we hypothesize that tumor intrinsic factors may be influencing the differentiation of tumor-specific CAF subpopulations (HNCAF-0) while nonspecific CAF subpopulations like HNCAF-1 could be derived from mesenchymal stem cells (29–33).

Despite the incomplete biological understanding of the HNCAF, our work introduces a novel clinically actionable biomarker for HNSCC. ICIs have revolutionized the field of cancer immunotherapy with mAbs targeted against CTLA-4, PD-1, and PD-L1 being recently approved for use as frontline therapies for HNSCC and other cancer types (1, 3, 34). The factors driving resistance to ICI remain largely unknown, making it difficult to select those who will respond and who will not. Accordingly, there remains an unmet need for reliable biomarkers predictive of response to guide patient selection and optimization of ICI treatment. In recent studies, CAF have been implicated to influence resistance to checkpoint immunotherapy. A preclinical model of pancreatic ductal adenocarcinoma showed that depletion of CAF expressing high levels of FAP improves response to αPD-L1 blockade (9). Similarly, scRNA-seq revealed a LRRC15⁺ CAF population associated with worse response to αPD-L1 immunotherapy in a clinical trial for bladder cancer (10). Furthermore, distinct CAF populations identified in breast cancer were also shown to be associated with poor αPD-1 immunotherapy response in melanoma and lung cancer (11). All these studies have implicated CAF primarily as contributors to resistance; however, the precise nature of molecularly distinct CAF subtypes and their role in mediating the effect of immunotherapy remains poorly investigated. With our ability to provide a higher resolution CAF repertoire, we show that the presence of two unique HNSCC-specific CAF subtypes is predictive of clinical response to immunotherapy in HNSCC. In particular, our functional findings suggest that HNCAF-0/3 fibroblasts are active participants in the immune response elicited by PD-1–directed immunotherapy through T-cell modulation. For HNSCC, we are currently evaluating whether these CAF subtypes behave differently in virally associated HPV+ HNSCC versus HPV− tumors.

Previous studies have typically shown CAF as promoters of immunosuppression. CAF have been shown to prevent T-cell infiltration and to kill T cells in an antigen-dependent manner, via PD-L2 and FasL (9, 35). CAF have also been shown to suppress T cells through upregulation of PD-L1 and PD-L2 and through recruitment of regulatory T cells (12, 36). While confirming the immunosuppressive role of some HNCAF subtypes, our work has also established a novel proinflammatory role for distinct CAF subtypes, which act as a promoter of T-cell activation and cytotoxicity. In light of this immunostimulatory function, we have termed the HNCAF-0/3 phenotype as T cell–stimulating CAF (tsCAF). On the basis of previous studies identifying TGFβ1 signaling through SMAD3 as a regulator of PD-1 expression, we hypothesize that tsCAF may repress SMAD3 to transcriptionally inhibit PD-1/TIM-3 expression (37). Although our data strongly suggest that HNCAF-0/3 are immunostimulatory, the inability to sort these cells to obtain pure populations is a major limitation of this study as the CAF cells used in our functional studies do not encompass a pure population at this point. Proteomic-based approach to select and sort CD29⁺FAP⁺ CAF into tsCAF and HNCAF-1 is an active area of investigation. Recently, tumor restrictive CD105⁻ CAF have been demonstrated in murine models of PDAC which are mediated by the adaptive immune system (38). CD105⁻ CAF highly overlaps with the myCAF gene signature, which has also previously been demonstrated to be tumor constraining (39). Because the tsCAF we describe here are distinct from the myCAF population, both molecularly and by surface marker expression, we believe they represent a distinct CAF population from CD105⁻ CAF.

Interestingly, we found that coculture of HNCAF-0/3 with CD8 T cells induced a Trm phenotype that coexpressed NKG2A. NKG2A is an inhibitory receptor that we and others found to be highly enriched in tumor-infiltrating Trm⁺ CD8 T cells in HNSCC (40, 41). NKG2A ligation with its ligand HLA-E reduces cytotoxicity and effector function and is therefore a novel immunotherapy target (42). Clinical trials combining NKG2A blockade (monalizumab) and other checkpoint inhibitors have shown promising results in HNSCC (40). While it is not clear why NKG2A is upregulated on tumor-infiltrating CD8 T cells, our findings suggest HNCAF characterization may also inform clinical development of NKG2A inhibition along with other ICIs. We found that upregulation of NKG2A required contact between the HNCAF-0/3 and activated CD8 T cells, which suggests that induction is mediated by either a surface ligand on HNCAF-0/3 or a component in the extracellular matrix produced by the CAF. To our knowledge, NKG2A expression on CD8 T cells has never been associated with fibroblasts and our finding here provides a clear link between intratumoral NKG2A expression on CD8 T cells and the tumor stroma. Future studies will need to be performed on HNSCC specimens to determine whether NKG2A expression on CD8 T cells is localized to stromal regions or associated with increased HNCAF-0/3 cells.

Plasticity of CAF subtypes have also been well demonstrated (14, 43), and our CAF atlas provides an excellent framework to develop strategies to force CAF differentiation toward the proinflammatory tsCAF phenotype rather than the immunosuppressive HNCAF-1 phenotype, to be combined with immunotherapy. Our characterization of each CAF subpopulation characterization through their master regulatory network lends itself toward this strategy. Further investigation of the signals that induce tsCAF formation and activation, possibly by targeting master regulators of the two subtypes, either genetically (44, 45) or pharmacologically, via the OncoTreat algorithm (46) is currently underway. Critically, this study highlights a much greater molecular heterogeneity of CAF subtypes than previously appreciated and demonstrates the critical need to functionally characterize their pleiotropic effects in terms of cancer progression, outcome, and response to immunotherapy and other cancer treatments.

A. Obradovic reports grants from NIH F30 CA260765-01 during the conduct of the study; in addition, A. Obradovic has a patent for (TH Docket No. 222230-8210) “Unique Cancer Associated Fibroblast Subsets Predict Response to Immunotherapy” pending. D. Graves reports grants from DoD USAMRAA and NIH, as well as other support from Barry Baker Biorepository Fund and James Rowen Fund during the conduct of the study; in addition, D. Graves has a patent for 63/212,927 pending. M. Korrer reports grants from DoD USAMRAA and NIH, as well as other support from Barry Baker Biorepository Fund and James Rowen Fund during the conduct of the study; in addition, M. Korrer has a patent for 63/212,927 pending. A. Naveed reports grants from DoD USAMRAA and NIH Grant, as well as other support from Barry Baker Biorepository Fund and James Rowen Fund during the conduct of the study. J. Curry reports other support from Rakuten Medical outside the submitted work. P. Hurley reports grants from NIH, American Cancer Society, and Jim Rowen Fund during the conduct of the study, as well as personal fees from Horizon Discovery outside the submitted work. X.S. Liu is a cofounder, board member, SAB member, and consultant of GV20 Oncotherapy and its subsidiaries. R. Uppaluri reports personal fees and other support from Merck Inc outside the submitted work. C.G. Drake reports a patent for LAG-3 for Cancer Immunotherapy licensed from Johns Hopkins University to BMS, and is currently employed by Janssen Research and Development. A. Califano reports other support from DarwinHealth Inc. outside the submitted work; in addition, A. Califano has a patent for VIPER algorithm (US 20170076035 A1)) issued, licensed, and with royalties paid from DarwinHealth. Y.J. Kim reports grants from NCI and NIDCR during the conduct of the study, as well as personal fees from Aduro, AstraZeneca, Sanofi, Takeda, and Mersanna outside the submitted work. In addition, Y.J. Kim has a patent for CAF signature pending, and is currently an employee with Regeneron Pharmaceuticals. No disclosures were reported by the other authors.

A. Obradovic: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Graves: Conceptualization, data curation, formal analysis, validation, investigation, methodology, writing–original draft, writing–review and editing. M. Korrer: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Y. Wang: Data curation, software, formal analysis, methodology. S. Roy: Data curation, validation, investigation, visualization, methodology. A. Naveed: Data curation, software, formal analysis, validation, investigation, methodology. Y. Xu: Resources, data curation, software, formal analysis, funding acquisition, validation, writing–review and editing. A. Luginbuhl: Resources, data curation, funding acquisition, writing–review and editing. J. Curry: Resources, data curation, funding acquisition, investigation, writing–review and editing. M. Gibson: Resources, data curation, formal analysis, supervision, investigation, methodology, writing–review and editing. K. Idrees: Resources, data curation, formal analysis, supervision, investigation, project administration. P. Hurley: Data curation, formal analysis, supervision, investigation, methodology, writing–review and editing. P. Jiang: Data curation. X.S. Liu: Conceptualization, resources, data curation, supervision, validation, investigation, visualization, methodology, project administration, writing–review and editing. R. Uppaluri: Conceptualization, data curation, software, formal analysis, supervision, validation, investigation, methodology, project administration, writing–review and editing. C.G. Drake: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. A. Califano: Conceptualization, resources, data curation, software, formal analysis, supervision, validation, investigation, methodology, project administration, writing–review and editing. Y.J. Kim: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by the Barry Baker Biorepository Fund (Y.J. Kim), James Rowen Fund (Y.J. Kim, P. Hurley), NIH R01 CA178613 (Y.J. Kim), R01 DE027749 (Y.J. Kim), Dept. of Defense CDMRP Breakthrough Award (Y.J. Kim), the American Cancer Society 131356-RSG-17-160-01-CSM (P. Hurley), and NIH RO1CA211695-01A1 (P. Hurley). A. Califano was supported by an NIH/NCI Outstanding Investor Award (R35 CA197745), and two NIH Instrumentation grants (S10 OD012351 and S10 OD021764). Fresh HNSCC tumor tissue was harvested and processed by Drs. James Netterville, Kyle Mannion, Sara Rohde, Robert Sinard, and Brandee Brown.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.