Deciphering the cellular tumor microenvironment landscape in salivary gland carcinomas using multiplexed imaging mass cytometry
Marcel Mayer, Lisa Nachtsheim, Louis Jansen, Philipp Wolber, Marcel Schmiel, Luc G. T. Morris, Fengshen Kuo, Alan L. Ho, Alexander Quaas, Jens Peter Klussmann, Christoph Arolt
TL;DR
This study uses advanced imaging to map the tumor environment in salivary gland cancers, identifying key cell types and their roles in cancer spread and prognosis.
Contribution
The study introduces a novel spatial analysis of tumor microenvironments in salivary gland carcinomas using imaging mass cytometry and spatial transcriptomics.
Findings
SDC tumors have high levels of Collagen- and matrix-CAFs, linked to poor prognosis.
ACC tumors show strong immune infiltration, suggesting potential for immunotherapy.
mCAF-endothelial interactions are associated with metastasis and poor survival in SDC.
Abstract
To spatially characterize the single-cell tumor microenvironment (TME) of salivary gland carcinomas (SGC) and identify prognostic biomarkers. SGC, including salivary duct carcinomas (SDC), acinic cell, mucoepidermoid, and secretory carcinomas, were analyzed using a 13-marker imaging mass cytometry panel. Multichannel image data from 54 primary cases and nodal metastases were processed to generate single-cell datasets. Cell phenotypes (tumor cells, cancer-associated fibroblasts (CAFs), endothelia, immune cells) were classified using a validated CAF algorithm, followed by spatial analysis and clinicopathological correlation. Clinicopathological results were validated using a previously published independent bulk RNA-sequencing (RNA-seq) cohort of n = 67 SDC cases. Spatial transcriptomics data of three SDC cases was leveraged to understand the molecular mechanisms of spatial interactions.…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5- —German Research Foundation
- —Jean-Uhrmacher Foundation
- —https://doi.org/10.13039/501100011566Marga und Walter Boll-Stiftung
- —Gusyk Program
- —Universitätsklinikum Köln (8977)
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsSalivary Gland Tumors Diagnosis and Treatment · Single-cell and spatial transcriptomics · Cancer Cells and Metastasis
Introduction
Salivary gland carcinomas (SGC) are rare and account for approximately 6% of head and neck malignancies [1]. They form a heterogeneous group of 21 morphologically and prognostically distinct tumor entities [2], some of which have high potential for recurrence and metastasis. Apart from biomarker-guided therapies for salivary duct carcinomas (SDC) and secretory carcinomas (Sec), there is currently no established targeted systemic therapy for SGC.
The tumor microenvironment (TME), which consists of a cellular compartment and the extracellular matrix (ECM), is an essential component of carcinomas and significantly influences their growth and response to chemotherapy [3]. Cancer-associated fibroblasts (CAFs) are activated fibroblasts that form the cellular TME, together with immune cells and tumor vasculature. Previous studies have characterized a variety of CAF subpopulations, each with a distinct marker profile, functionality, and spatial distribution [4–7]. However, because of the requirement for several markers, the precise characterization of CAF subsets remains a challenge when analyzing larger patient cohorts. Various studies have examined the importance of CAFs in carcinogenesis, ECM synthesis, angiogenesis, immune modulation, and chemotherapy resistance [8–10]. As anticipated, CAFs also influence patient outcome [11, 12]. Pro-carcinogenic factors make CAFs promising targets for new therapeutic approaches [13]. Cords et al. recently published two studies that systematically characterized CAF subsets across carcinomas from different organs using an unbiased single-cell RNA sequencing (scRNA-seq) approach and proposed a generalized CAF classification concept. They methodically transferred their classification algorithm to space-agnostic imaging mass cytometry (IMC) and demonstrated that CAF subtypes reside in distinct spatial niches, partly correlate with reduced immune infiltrate, and can determine prognosis in non-small cell lung cancer (NSCLC) [5, 12].
To date, little is known about the occurrence and function of CAFs in SGC. In a pilot study, we demonstrated significant differences in the frequency of COL11A1^+^ CAFs between different SGC entities [14]. More recently, we analyzed a large SGC cohort using deep proteomic profiling and demonstrated that SGC can be subclassified into two large groups based on the ECM, which, among other parameters, differentiates with regard to immune cell infiltration [15]. Given the marked impact of CAF subtypes on patient outcomes and the potential of TME-directed therapies, we aimed to analyze the cellular TME of 54 patients with SGC. This included four different tumor entities without myoepithelial differentiation, as myoepithelial cells may show an overlapping immunophenotype with CAFs when using a restricted marker panel. Using a multiplex IMC panel and a thoroughly validated classification scheme, we analyzed the spatial distribution of CAF subtypes, immune cells, tumor vasculature, and tumor cells, and correlated our findings with clinical data and the ECM landscape to shed light on the cellular TME landscape of SGC.
Methods
Patient cohort
The cohort comprised patients with sufficient formalin-fixed paraffin-embedded (FFPE) material of the primary tumor and/or lymph node metastasis collected at first diagnosis who were treated for primary SGC at the Department of Otorhinolaryngology, Head and Neck Surgery of the University Hospital of Cologne, Germany between January 1 st, 2004, and December 31 st, 2021. Demographics, survival, and histopathological data were extracted from clinical records and histopathological reports with respect to entities and disease stage at the time of diagnosis according to the AJCC TNM staging system (8th edition, 2020) [16]. In cases of missing follow-up data within the clinical records, patients and/or their general practitioners were phoned to follow up on the current tumor status. The study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of the University of Cologne (Approval Code: 13–091). Informed consent was obtained from all participants included in this study. Because this was not a prospective clinical trial, formal power calculations and randomization were not performed.
Antibody panel and conjugation of antibodies with metal isotopes
A panel including 15 markers that distinguish CAFs, immune cells, tumor cells, and endothelia was designed. The following monoclonal antibodies were purchased from Fluidigm®: CKAE1/3-148Nd (dilution 1:100; RRID:AB_3678170), CD31-146Nd (dilution 1:200; RRID:AB_3678176), Vimentin-143Nd (dilution 1:50; RRID:AB_2811069), aSMA-141Pr (dilution 1:200; RRID:AB_2890139), Collagen1-173Yb (dilution: 1:300; RRID:AB_3678205), Ki67-168Er (dilution: 1:50; RRID:AB_2800467), CD34-151Eu (dilution 1:200; RRID:AB_3678086), CD74-166Er (1:50; RRID:AB_3675795), CD45-152Sm (dilution: 1:50; RRID:AB_3678248), CD4-156Gd (dilution: 1:100; RRID:AB_2811051), CD8a-162Dy (dilution 1:200; RRID:AB_2811053), IDO-155Gd (dilution: 1:100; RRID:AB_3678089), and CD73-158Gd (dilution: 1:100; RRID:AB_3678093). The carrier-free CD138 monoclonal antibody (mouse, clone: B-A38; RRID:AB_3676290) was purchased from Origene® and conjugated with 164Dy (dilution 1:100) using the Maxpar X8 Multimetal Labeling Kit (Fluidigm®, 201300), following the manufacturer’s protocol. The intercalator iridium (Fluidigm®, 2011922 B) was included for the identification of nuclei. The conjugated antibodies were titrated and tested for specificity by visual comparisons in SGC samples, as well as in ovarian, breast, lung, and colon carcinoma control samples, resulting in the abovementioned dilutions for optimal signal yield.
Staining
The tissue was freshly sectioned and placed onto tissue microarray (TMA) slides five days before staining. Tissue slides underwent imaging mass cytometry (IMC) staining following the Fluidigm® protocol (PN400322A3). FFPE tumor samples were baked at 65 °C for two hours to remove residual wax. First, deparaffinization was performed. The slides were rehydrated in an ethanol series (100%, 95%, 90%, 80%, and 70%). Antigen retrieval was performed in a water bath at 96 °C using Tris–EDTA buffer (pH 9) for 30 min. Next, the slides were blocked with 3% bovine serum albumin (BSA) in Maxpar phosphate-buffered saline (PBS) for 45 min at room temperature in a hydration chamber. Afterwards, the slides were incubated with the antibody cocktail in the above-mentioned dilutions (0.5% BSA buffer) and stored overnight at 4 °C in a hydration chamber. The following day, the slides were washed four times with 0.2% Triton X-100 in Maxpar PBS and twice with Maxpar PBS. DNA staining was performed by incubating the slides with Intercalator-Ir (Fluidigm®, 201192 B) in Maxpar PBS (dilution 1:400) for 30 min in a hydration chamber. Finally, the slides were washed twice with Maxpar water and air-dried for 30 min.
Imaging mass cytometry
Over a period of four weeks, all cores on all six TMA slides were laser-ablated at 200 Hz with an ablation energy of 3 dB using a Hyperion imaging mass cytometry system connected to a Helios mass cytometer (Fluidigm®). A tuning slide (Fluidigm®) was used prior to each TMA slide. All ROIs were set to 500 μm × 500 μm.
Preprocessing
The IMC raw data were preprocessed and segmented using a combination of Jupyter notebooks (RRID:SCR_018315), CellProfiler pipelines (RRID:SCR_007358), and the Ilastik tool (RRID:SCR_015246) as described previously [12, 17]. Briefly, using published Jupyter notebooks (RRID:SCR_018315), proprietary MCD files were converted into single-channel TIFFs, filtered for “hot pixels,” upscaled, and cropped to facilitate the segmentation step. Then, an Ilastik pixel classifier was trained, which used all marker channels except Ki67. Because CAFs can have extensive cellular processes, cell bodies were reduced to their nuclear core to enable segmentation of individual cells [12]. Next, a CellProfiler pipeline (RRID:SCR_007358) was used to extract single-cell information as a.csv file. All subsequent steps were performed in R using the imcRtools package (RRID:SCR_017132) and published workflow [17]. After reading the data in R, the counts were transformed using an inverse hyperbolic sine function. This was followed by optical and numerical verification of the segmentation quality. T-SNE and UMAP plots revealed no batch effect among non-neoplastic cells that could have been introduced by TMA slides (Supplementary Figure S1A and B). All TMA cores were visualized and those without tumor cells were excluded from subsequent analyses. This led to the exclusion of 58 cores (100,425 cells) and a remainder of 199 cores (408,939 cells). These cores from 54 patients comprised 188 and 11 cores from primaries and metastases, respectively.
Phenotyping
First, the tumor cells were identified using a Gaussian mixture model of cytokeratin expression (Cytokeratin AE1/3 antibody). Next, all ROIs were plotted with tumor cell annotations to ensure morphologically valid classification. All subsequent stepwise clustering of TME cells was performed using a combination of FlowSOM (RRID:SCR_016899) clustering and meta-clustering with ConsensusClusterPlus (RRID:SCR_016954), as described previously [5]. First, all TME cells were clustered into endothelial cells (CD31^+^, Cytokeratin^−^), immune cells (CD45^+^/CD138^+^, Cytokeratin^−^), and CAFs (CD31^−^, Cytokeratin^−^, CD45^−^, CD138^−^). Negative selection of CAFs has been previously carried out in other exploratory studies [18]. Immune cells were then clustered by the expression of CD45, CD4, CD8, and CD138 into plasma cells (CD138^+^), CD8^+^ T cells, CD4^+^ T cells (CD4^high^), and other CD4-expressing immune cells (CD45^+^, CD4^low^). The latter cell population likely corresponds to other CD4-expressing hematopoietic cells, such as macrophages and other myeloid-derived cells, which have been shown to express CD4 to a lesser extent using IMC [5, 12, 19]. A more detailed subclassification was not feasible based on the choice of markers. CD4^+^ T cells were clustered again, using CD45, CD4, CD74, and Ki67, into CD4^+^ T cells, proliferating CD4^+^ T cells (Ki67^high^), and CD4^+^CD74^+^ T cells. This CD4^+^ T cell subset has recently been found to correspond mostly to tumor-associated Tregs and more specifically to Tregs with effector cytokine production [20, 21].
Next, CAFs were subclassified by clustering with the markers Vimentin, aSMA, CD34, IDO, CD73, CD74, Ki67, and Collagen1 into Collagen CAFs, matrix CAFs (mCAF), SMA^+^ CAFs (SMA CAF), developing CAFs (dCAF) and antigen-presenting CAFs (apCAF), largely following a classification scheme recently proposed and validated in different carcinoma types [5, 12]. Cells without specific CAF marker expression (SMA/Collagen1/Ki67/CD74) were split into CAFs with high Vimentin expression and other cells. The latter clustered again and were shown to correspond to other specific CAF and immune cell subsets. The remaining cells that were negative for all markers in the panel were designated as tumor cells because they were located in tumor cell patches. We reasoned that these cells likely express low levels of cytokeratin which could not be detected with the chosen dilution of the anti-cytokeratin antibody after preliminary experiments with a more limited number of tumor samples. Moreover, since the main scope of this work was the TME, this conservative approach minimizes the danger of producing false-positive results.
Spatial analysis
Because we assumed different spatial contexts in metastases and primaries, we included only primaries in the subsequent spatial analyses with imcRtools (RRID:SCR_017132). First, we constructed a cellular interaction graph using k-nearest neighbor detection with k = 20. We then defined the cellular fractions of the 20 nearest neighbors of each cell and clustered each cell based on this information into one of nine clusters, that is, cellular neighborhoods. We also calculated the results for 6 and 12 clusters but found that nine neighborhoods yielded interesting spatial cellular distributions, such as a tumor-stroma interface neighborhood, without splitting the data into clusters with only one cell type. We then carried out tumor patch detection (minimal patch size of 10 cells) to allow for the identification of tumor-infiltrating cells and the calculation of the spatial distances between each TME cell and the tumor border. Finally, we carried out statistical testing of interactions/co-localization by calculating the averaged spatial interaction of each combination of the two cell types and comparing it against a null distribution in each ROI using the testInteractions function from the imcRtools package (RRID:SCR_017132), leveraging a testing strategy proposed by Schapiro et al. [22].
Correlation with SGC ECM modules
We recently proposed three different ECM protein modules, which explain the ECM of SGC and are heavily enriched in CAF-associated terms, basement membrane-associated terms, and those that are linked to peripheral blood components as well as blood-producing and blood-processing organs [15]. We found that these modules describe ECM differences between SGC with and without myoepithelial differentiation. From the single-cell dataset, we calculated the median cell type frequency per primary and correlated these values with the ECM module Eigengenes in 40 samples that were present in both cohorts using the Spearman method. P-values were corrected for multiple testing using the Benjamini–Hochberg procedure.
Clinico-pathological analysis
Statistical analyses were performed using R software (version 4.4.1) in RStudio (version 2024.4.2.764; RRID:SCR_000432). The R packages ggplot2 (RRID:SCR_014601), survfit, survminer (RRID:SCR_021094), survival (RRID:SCR_021137), forestmodel (RRID:SCR_025250), haven, broom (RRID:SCR_026874), lazyeval, tidyverse (RRID:SCR_019186), rstatix (RRID:SCR_021240), and ggpubr (RRID:SCR_021139) were used for analysis and visualization. Categorical data were presented as percentages, whereas continuous data were expressed as means. The Mann–Whitney U test and Kruskal–Wallis test were used for comparisons between two and multiple groups, respectively. The Benjamini–Hochberg procedure was used to correct for multiple tests. The Kaplan–Meier method was used to visualize overall survival (OS), recurrence-free survival (RFS), recurrence-free probability (RFP), and distant control rate (DCR), and the log-rank test was used to test differences between groups. OS was defined as the time interval between the date of diagnosis and date of death. RFS was defined as the time interval between the date of diagnosis and date of recurrence or death. RFP was defined as the time interval between the date of diagnosis and date of recurrence. DCR was defined as the time interval between the date of diagnosis and date of distant recurrence. Cox proportional hazards survival regression was used to determine the influence of different variables on OS, RFS, RFP, and DCR. In the case of multiple significant associations in the univariate Cox regression analysis, the relevant variables as well as other potentially relevant variables were tested in a multivariate Cox regression analysis to identify independent prognostic factors. A p-value < 0.05 was considered statistically significant.
Generation of an mCAF signature from single cell RNA-seq data and its application in a clincally-annotated bulk-RNA seq dataset
Raw RNA-sequencing (RNA-seq) data along with clinical parameters of n = 67 SDC samples previously published by Kohsaka et al. [23] were received from the authors (raw RNA-seq data is publicly available under accession number hum0094 in the Japanese Genotype–Phenotype Archive). Briefly, raw RNA-seq data were aligned against human genome assembly hg19, quality check was performed, and TPM (transcripts per million)-normalized data were generated. Next, we downloaded the results of the differential expression analysis contrasting mCAFs and other CAF types that were generated from an integrated single-cell RNA-sea data set spanning different carcinoma types (breast, colon, head and neck, non-small-cell lung, pancreas) by Cords et al. [5]. We selected the top 20% of overexpressed genes (n_genes_ = 98) for the construction of an mCAF-like signature in the bulk RNA-seq dataset. For validation analysis, also other cutoffs were applied. Even though the underlying data analysis from Cords et al. suggests a generalizability, we cannot exclude a different biology in SGC and therefore term the signature “mCAF-like”. For each sample, the mCAF-like signature was created by multiplying the expression of each signature gene with its respective log2FC value as calculated by Cords et al. [5] Then, all products were summed and averaged. Cases were stratified by median portion of the mCAF-like signature into mCAF_high_ and mCAF_low_ cases.
To validate results from the IMC data analysis, the Kaplan–Meier method was used to visualize RFS and RFP, and one-sided log-rank tests were applied to test differences between mCAF_high_ and mCAF_low_ SDC cases.
Spatial correlation of mCAF signature gene expression and endothelial markers
10 × Visium spatial transcriptomics (ST) data sequenced on the Visium CytAssist Instrument according to the manufacturer’s protocol within a cooperation with Memorial Sloan Kettering Cancer Center, New York City, USA were available for further analysis. Space Ranger outputs from three SDC patients (one primary, one lymph node metastasis, and one lung metastasis sample) were imported into R and analyzed using the Seurat package (RRID:SCR_016341). First, the gene expression data of each sample was normalized and variance-stabilized using the SCTransform algorithm [24]. Then, the samples were merged (total of 14,976 spots, i.e., transcriptomes) and module scores were calculated for mCAFs, pericytes, and endothelia using Seurat’s AddModuleScore function. For this, a list of all mCAF-like signature genes (n_genes_ = 98), the top 20% of marker genes for pericytes as defined by Cords et al. [5] (also by average logFC; n_genes_ = 123), and the canonical endothelial marker genes PECAM1 (CD31) VWF, CDH5 (CD144), KDR (VEGFR2), and MCAM (CD146) served as input for the module scores (log2FC values are not used as input for the AddModuleScore function). Next, we determined the spatial co-localization of mCAFs, pericytes, and endothelia using both vessel-associated module scores and the mCAF-like score (Pearson method).
To characterize areas with high and low interaction between mCAFs and endothelia, for each sample separately, spots with the 20% highest and lowest module score values were defined as mCAF low/high and endothelia high/low. Alternative cutoffs (10%, 30%) were evaluated to ensure the robustness of this approach (Supplementary Figure S6 A and B). The intersecting spots of each of these groups were identified sample-wise and then, the spots were pooled across samples. The following spot categories were created: "mCAF^high^endothelia^high^" (n_spots_ = 964), "mCAF^high^endothelia^low^" (n_spots_ = 416), "mCAF^low^endothelia^high^" (n_spots_ = 291), "mCAF^low^endothelia^low^" (n_spots_ = 1,414), “other” (n_spots_ = 11,891).
Since we were interested to specifically interrogate the metastatic capacity of the interaction between mCAFs and tumor vessels, we downloaded 49 manually curated metastasis-associated signatures from “The Molecular Signatures Database” (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb). Additionally, we systematically searched all “REACTOME”, “HALLMARK”, and GeneOntology gene sets for “VEGF”, “IL6”, “IL-6” and “IL_6”, which resulted in four VEGF- and two IL6-related signatures which were used for further analyses. After excluding lists with a gene coverage of < 40% in the ST data, signature scores were calculated for each gene set and each ST-spot using the AddModuleScore function.
To characterize tumor niches with a high co-localization of mCAFs and endothelia in a broader manner, we then carried out a differential expression analysis contrasting mCAF^high^endothelia^high^ spots with all other ST spots (Seurat’s FindAllMarkers function with default settings) and annotated the results with a gene set enrichment analysis (GSEA; clusterProfiler R package (RRID:SCR_016884), GeneOntology terms).
Results
Clinicopathological data
Overall, 54 patients with SGC were included in this study. The most frequent entity was salivary duct carcinoma (SDC; 43.6%), followed by acinic cell carcinoma (ACC; 24.1%), mucoepidermoid carcinoma (MEC; 22.2%), and secretory carcinoma (Sec; 11.1%). The mean age across the whole cohort was 57.7 (± 17.5) years, and 46.3% of all patients were female. Most patients (57.7%) had T1/2 tumors. Among SDC, 81.1% were androgen receptor (AR)-positive, and 54.5% showed nuclear AR positivity in > 70% of tumor cells (AR_high_), whereas 45.5% showed nuclear AR positivity in ≤ 70% of tumor cells (AR_low_). Furthermore, 36.4% of the SDC were positive for human epidermal growth factor receptor 2 (HER2; score 3 + or HER2-amplified). The primary therapy was surgery in 98.1% of cases, and chemoradiation in 1.9% of cases. Neck dissection was performed in 74.1% of the cases, and 57.4% of the patients received adjuvant (chemo-)radiation therapy. Detailed data for the specific entities are presented in Table 1. Table 1. Clinicopathological data**All n = 54 (%)****SDC n = 23 (%)****ACC n = 13 (%)MEC n = 12 (%)**Sec n = 6 (%)**LocalizationParotid gland50 (92.6)21 (91.3)13 (100.0)11 (91.7)5 (83.3)Submandibular gland2 (3.7)1 (4.3)0 (0.0)0 (0.0)1 (16.7)Small salivary gland2 (3.7)1 (4.3)0 (0.0)1 (8.3)0 (0.0)DemographicsFemale25 (46.3)5 (21.7)8 (61.5)9 (75.0)3 (50.0)Male29 (53.7)18 (78.3)5 (38.5)3 (25.0)3 (50.0)Age57.7 ± 17.566.5 ± 15.454.4 ± 18.348.8 ± 16.849.2 ± 18.8Histopathological parameters T stageT1/230 (57.7)11 (47.8)6 (50.0)9 (75.0)4 (80.0)T3/422 (42.3)12 (52.2)6 (50.0)3 (25.0)1 (20.0)N/A20101N stageN030 (57.7)5 (21.7)9 (75.0)11 (91.7)5 (100.0)N + 22 (42.3)18 (78.3)3 (25.0)1 (8.3)0 (0.0)N/A20101M stageM052 (96.3)21 (91.3)13 (100.0)12 (100.0)6 (100.0)M12 (3.7)2 (8.7)0 (0.0)0 (0.0)0 (0.0)Vascular invasionV044 (81.1)17 (73.9)11 (91.7)11 (100.0)5 (100.0)V16 (9.8)5 (26.1)1 (8.3)0 (0.0)0 (0.0)N/A41111Perineural invasionPn028 (56.0)4 (19.0)9 (75.0)10 (83.3)5 (100.0)Pn122 (34.0)17 (81.0)3 (25.0)2 (16.7)0 (0.0)N/A42101Lymphovascular invasionL037 (74.0)11 (50.0)10 (83.3)11 (100.0)5 (100.0)L113 (26.0)11 (50.0)2 (16.7)0 (0.0)0 (0.0)N/A41111Extracapsular extensionECE-41 (80.4)13 (59.1)11 (91.7)12 (100.0)5 (100.0)ECE + 10 (19.6)9 (40.9)1 (8.3)0 (0.0)0 (0.0)N/A31101AR status (IHC)PositiveN/A18 (81.1)N/AN/AN/ANegativeN/A4 (18.2)N/AN/AN/AN/A1> 70%N/A12 (54.5)N/AN/AN/A≤ 70%N/A10 (45.5)N/AN/AN/AN/A1HER 2 IHC Score0N/A5 (22.7)N/AN/AN/A1 + N/A4 (18.2)N/AN/AN/A2 + N/A5 (22.7)N/AN/AN/A3 + N/A8 (36.4)N/AN/AN/AN/AN/A1N/AN/AN/APrimary therapySurgical resection532213126Chemoradiation11000Neck dissectionYes402111104No121222N/A21000Adjuvant therapyRadiation2114222Chemoradiation107210None222894N/A10100SDC salivary duct carcinoma, ACC acinic cell carcinoma, MEC mucoepidermoid carcinoma, Sec secretory carcinoma, AR androgen receptor, HER2 human epidermal growth factor receptor 2, IHC immunohistochemistry
With a median follow-up of 63.5 months for the entire cohort, SDC patients showed the most unfavorable five-year OS (50.6%), followed by MEC (83.3%), ACC (83.9%), and Sec patients (100.0%). Five-year RFS rates were 30.7% for SDC, 69.2% for ACC, 83.3% for MEC, and 66.7% for Sec. Moreover, five-year RFP was 39.9% for SDC, 69.2% for ACC, 91.7% for MEC, and 66.7% for Sec. Lastly, five-year DCR was 47.0% for SDC, 69.2% for ACC, 100.0% for MEC, and 83.3% for Sec.
A 13-marker IMC panel identifies distinct immune cell and CAF subsets in nonMYO SGC
We used tissue from SGC primaries and metastases from 54 patients with SGC to create a tissue microarray (TMA). Tissue sections from TMA were stained with a cocktail of 13 metal-conjugated antibodies. Ablation of the tissue slices produced single-channel images that were used to generate multichannel images. After segmentation, the cells were categorized into cellular subsets. The cellular frequencies and results from the spatial analyses were then correlated with clinical parameters and patient outcomes (Fig. 1A).Fig. 1. Phenotyping of the cellular SGC TME. A Overview of the study: SGC samples from 54 patients were captured on a TMA, which was incubated with 13 metal-conjugated antibodies and ablated using a Hyperion system. Through image segmentation, a spatially-resolved single cell dataset was obtained and cell phenotyping was carried out in a semi-supervised fashion with established cell type markers. Cell frequencies and cellular neighborhoods were then correlated with clinical data. B tSNE plots of all 408.939 cells colored by cell type and (C) cell category. D Heatmap depicting the mean marker expression by cell type. E Representative images of the tumor architecture by marker expression (top row) and cellular composition after phenotyping (bottom row). 100-micron scale bars
After the thorough exclusion of non-tumor-bearing TMA cores, we analyzed 199 SGC TMA cores from 54 SGC, including 188 primaries and 11 metastases from SDC, ACC, MEC, and Sec. A median cell count of 3,019 (ACC), 2,274 (Sec), 2,104 (MEC), and 2,029 (SDC) per 1 mm^2^ ROI was noted. Using a stepwise Gaussian mixture model of cytokeratin expression and SOM clustering, we measured the expression of 13 different markers in 408,939 cells, (Fig. 1B and C). First, we separated tumor cells (CKAE1/3^+^; n = 266,172) from endothelial cells (CD31^+^, CKAE1/3^−^; n = 1,349), immune cells (CKAE1/3^−^, CD45^+^, or CD138^+^; n = 50,432), and CAFs (CKAE1/3/CD45/CD31/CD138^−^; n = 90,986). Immune cells were then clustered into plasma cells (CD138^+^), CD8^+^ T cells, CD4^+^ T cells, and weakly CD4-expressing cells, which are most likely monocytes or other myeloid-derived cells (other CD4^+^cells; [5, 12, 19]). CD4^+^ T cells were subclassified as CD4^+^CD74^+^ T cells, which most likely correspond to effector Tregs [20, 21], proliferating CD4^+^ T cells (Ki67^+^), and other CD4^+^ T cells. Immune cells solely expressing CD45 were classified as “other immune cells.” To classify CAFs, we adopted an IMC-validated classification scheme proposed by Cords et al. [5, 25], separating mCAFs (SMA^+^Collagen1^+^), Collagen CAFs (SMA-Collagen1^+^), SMA-CAFs (SMA^+^Collagen1^−^), dCAFs (Ki67^+^, Vimentin^+^), and apCAFs (CD74^+^, Vimentin^+^). We also included CD73, IDO, and CD34 to detect rCAFs, ifnCAF/IDO CAFs, and iCAFs, respectively. However, the CD73 and IDO antibodies used did not produce a specific IMC signal. We could not discriminate iCAFs, since all CD34^+^ cells co-expressed CD31, a specific vascular marker, which was used to identify endothelia. Two subsets of cells did not display a specific marker profile: cells without expression of any of the markers used were designated as tumor cells, as they were located in tumor cell patches. Regarding the TME-centered focus of this work, we reasoned that this conservative approach would not impact TME-cell classification and thus minimize the chance of false-positive results. Cells with sole expression of Vimentin were classified as “other CAFs” with Vimentin expression (Fig. 1D). Figure 1E displays representative ROIs of the included tumor entities as multiplex images and the corresponding cellular maps after segmentation and cell phenotyping.
SGC entities are characterized by different TME cell frequencies
When pooling all ROIs, the cell type distribution across SGC tumor types displayed marked differences with a prominent immune cell compartment in ACC (Supplementary Figure S1C and D). No prominent differences in cell category contributions were noted when comparing primaries and metastases or peripheral and central tumor regions (Supplementary Figure S1E and F).
To account for differing numbers of ROIs per patient due to the exclusion of non-tumor-bearing and torn TMA cores, we calculated the median cell frequency per patient and compared the distribution of cells across the different SGC entities. We found a significantly varying abundance of immune cell subsets across tumor types (Fig. 2A), which was largely due to significantly elevated levels of CD4^+^ T cells (p = 0.03), CD8^+^ T cells (p = 0.04), other CD4^+^ cells (p = 0.03), and unclassified immune cells (p = 0.04) in ACC. This resulted in a significantly higher frequency of overall immune cells in ACC (p = 0.04; Fig. 2B and C).Fig. 2. Comparison of median cell frequencies per patient across tumor subtypes (A) Immune cell subsets across SGC entities. B Immune cell subsets in ACC vs. other carcinoma types. C Representative core immune marker expression (top row) and selected immune cell subsets in SDC versus ACC. D Frequencies of CAF subtypes across SGC entities. E Examples of the spectrum of core CAF marker expression (top row) and cell categories (bottom row) in SDC. F Frequency of tumor-infiltrating CAFs dependent on AR expression in SDC. G Frequency of tumor-infiltrating apCAFs in ACC versus other SGC subtypes. H Representative images of CD74^+^ expression in apCAFs in ACC (top row) and corresponding cell categories after segmentation (bottom row). 100-micron scale bars
Although the overall distribution of CAFs did not significantly vary in SGC types (Fig. 2D), we detected the highest mean frequency of mCAFs in SDC, a CAF subset that has been associated with ECM production. However, in direct comparison with other entities, this elevation was not significant (Supplementary Figure S2A), probably because of the broad spectrum of mCAF frequencies within SDC (Fig. 2D and E).
We also analyzed the frequency of tumor-infiltrating CAFs using tumor patch detection (Supplementary Figure S2B). As AR and HER2 are frequently expressed in SDC and are molecular targets for systemic therapy approaches [26], we examined the association between AR/HER2 expression and CAF frequency. We were able to show that the number of tumor-infiltrating CAFs was significantly elevated in AR_high_ SDC compared to AR_low_ tumors (p = 0.049; Fig. 2F). No other significant cellular TME alterations were noted when SDC was subgrouped based on AR (Supplementary Figure S2C-E) and HER2 status (Supplementary Figure S2F-H). In contrast, in direct comparison to other tumor types, ACC displayed elevated levels of tumor-infiltrating apCAFs (p = 0.003; Fig. 2G and H, Supplementary Figure S2I). However, the frequency of these cells was generally very low. No significant differences in tumor-infiltrating immune cells were noted across SGC tumor types (Supplementary Figure S2J).
Spatial interaction analyses reveal a co-localisation of mCAFs with intratumoral vasculature
The spatial composition of the TME is an integral feature of tumor biology and influences the response to immunotherapy and patient outcomes [27]. Therefore, we leveraged the spatial information preserved in the dataset.
After computing a spatial interaction graph (examples in Supplementary Figure S3A), we clustered cells based on their 20 nearest neighbors into 9 cellular neighborhoods (CN; Fig. 3A and B). We found that the tumor cell-rich neighborhoods were largely devoid of CAFs and immune cells (CN1 and CN3). In contrast, CD4^+^ T cells, proliferating CD4^+^ T cells, other CD4^+^ cells, and CD8^+^ T cells formed a distinct lymphocyte-rich cluster (CN9, Fig. 3C). Another immune-related CN mainly consisted of CD4^+^CD74^+^ T-cells, apCAFs, and dCAFs.Fig. 3. Spatial analysis of the TME of SGC and correlation with ECM protein modules. A Exemplary ROIs with cells colored by cellular neighborhood (CN). B Scaled cell abundance per CN. C Representative ROI rich in cells that mainly contribute to CN9 (ACC) and CN8 (MEC). D Frequency of CN across SCG types. E Spatial interaction of cell types colored by the number of ROIs with significant co-localization of two cell types. A red color indicates more ROIs with significant interaction while blue tiles indicate more ROIs with significant avoidance of two cell types. F Correlation of cell types with ECM module eigengenes which were published previously. Symbols within the tiles indicate the test significance (* p < 0.05, ns not significant). 100-micron scale bars
However, as the contributing cell types, the latter CN occurred at very low frequencies (CN2, Supplementary Figure S3B). SMA CAF and mCAFs contributed to two distinct CN with partial overlap (CN 6 and 8, respectively; Supplementary Figure S3C). Interestingly, the matrix-associated CAF-rich (mCAF-rich) CN8 displayed a strong contribution to endothelia (Fig. 3C). No specific marker for vCAFs, another cell type that was found to be associated with endothelia by Cords et al., was available in this panel. However, in contrast to mCAFs, very low expression of SMA and Collagen1 was reported in vCAFs [12], allowing for discrimination between the two cell types.
Although we did not observe significant differences in CN frequencies across SGC types, immune cell-rich CN was distinctively enriched in ACC. In addition, in SDC, we noted higher frequencies of the mCAF-rich CN8, but not of CNs 5, 6, and 7, which have high contributions of Collagen CAFs and SMA CAFs (Fig. 3D, Supplementary Figure S3C and D).
The CN analyses were complemented by a more direct spatial interaction analysis, as proposed by Schapiro et al. [22] (Fig. 3E). This approach compares the spatial interactions of cell type pairs with a null distribution in each ROI. The number of ROIs with significant positive (red) and negative (blue) interactions can then be summarized and graphically displayed. This methodically different analysis largely validated the aforementioned results, indicating potential interactions among immune cells and between mCAFs and endothelia as well as SMA CAFs. Again, mutual exclusion of TME and tumor cells was noted.
In summary, we describe distinct spatial cellular TME neighborhoods in SGC that largely show a mutually exclusive predominance of tumor cells, immune cells, and CAFs. However, mCAFs and Collagen CAFs tend to localize in proximity to the tumor vasculature.
CAFs and particularly mCAFs are associated with a distinct ECM profile
Recently, we dissected the ECM of SGC using an unbiased proteomic approach and discovered three protein clusters (“ECM modules”) that explain ECM differences across SGC tumor types via module Eigengenes (i.e., the resulting vector of the module’s proteomic signature) [15]. Using gene set enrichment analysis, these modules were biologically annotated and found to be enriched for classic CAF-associated (CAF-module), basement membrane-associated (BM-module), and peripheral blood-associated (Hem-module) biological terms. The most significant members of the latter include coagulation-related factors such as kallikrein, kininogen, prothrombin, plasminogen, and angiotensin. Since CAFs are considered the main source of ECM in carcinomas, we leveraged an overlap of 40 patients between that and the present SGC cohort and correlated the module Eigengenes with the cell type frequencies (Fig. 3F), thereby establishing a connection between the ECM and the cellular TME. As anticipated, we observed a strong positive correlation between CAFs with a more classic myofibroblastic phenotype (mCAF, Collagen CAF, SMA CAF) and the CAF protein module. However, only the correlation between mCAFs and overall CAFs was significant. In contrast, all immune cell subsets, except for plasma cells, were negatively correlated with the CAF module. Unexpectedly, very similar associations were found for the Hem-module, including an anti-correlation with endothelia and all immune cell subsets except plasma cells. We previously reasoned that this module consists of blood-related factors that are deposited within the ECM. However, the present data indicate that this does not imply exaggerated vasculature or deposition of cellular components of the peripheral blood. The BM module was mainly expressed in SGC with myoepithelial differentiation (e.g., adenoid cystic carcinomas), which were not analyzed in the present study.
Together, these data clearly link SGC mCAFs, Collagen CAFs, and SMA CAFs to an increase in classical CAF-associated ECM proteins and provide evidence that the ECM components of the CAF and the Hem module participate in an immune-exclusive TME.
Metastasis-associated signatures are enriched in tumor niches with a co-localization of mCAFs and endothelia and the interaction may be mediated by specific gene signatures
A significant positive correlation between mCAF-like and pericyte as well as endothelia signatures found in all three ST SDC data sets additionally supported the previously identified spatial co-localization of mCAFs and intratumoral vasculature (Fig. 5A and B; Supplementary Figure S5A and S5B).
A marked enrichment of metastasis-associated signatures such as “ANASTASSIOU_MULTICANCER_INVASIVENESS_SIGNATURE”, “GILDEA_METASTASIS”, “ROZANOV_MMP14_TARGETS_UP”, and “HEBERT_MATRISOME_TNBC_LUNG_METASTASIS” among others in mCAF^high^endothelia^high^ and a simultaneous depletion of those signatures in mCAF^low^endothelia^low^ ST spots (Fig. 5C) supported our previous clinically-based findings indicating that the spatial interaction of mCAFs and endothelia may be associated with metastasis. Notably, enrichment of metastasis-associated signatures was markedly higher in mCAF^high^endothelia^high^ than in mCAF^high^endothelia^low^ spots. In line with these findings, the invasiveness suppressing signature “RODRIGUES_DCC_TARGETS_UP” [28] was markedly depleted in mCAF^high^endothelia^high^ tumor niches. Alternative cut-offs for spot classification were evaluated and led to similar results, demonstrating the robustness of this approach (Supplementary Figure S6 A and B).
DE analysis followed by GSEA for mCAF^high^endothelia^high^ versus all other spots was performed to identify which factors may mediate the interaction between mCAFs and endothelia and therefore may be associated with metastasis of SDC (Fig. 5D-F and Supplementary Fig. 5 C). Notably, beside several ECM-/fibroblast-/collagen-related signatures, the “metalloendopeptidase activity” signature (enrichment score = 0.59; q-value < 0.001) was among the top-enriched signatures and, accordingly, the pro-metastatic matrix metalloproteinases [29–31] (MMP-28 (avg_log2FC = 2.10; p-adj. < 0.001), MMP-9 (avg_log2FC = 2.01; p-adj. < 0.001), MMP-11 (avg_log2FC = 1.69; p-adj. < 0.001), MMP-14 (avg_log2FC = 1.68; p-adj. < 0.001), and MMP-2 (avg_log2FC = 1.64; p-adj. < 0.001) were identified as strongly upregulated genes. More importantly, the “platelet-derived growth factor binding” (enrichment score = 0.84; q-value < 0.001) signature showed the second highest enrichment score among all signatures. Further, the “insulin-like growth factor binding” (enrichment score = 0.67; q-value = 0.007) signature ranked among the top-enriched signatures within mCAF^high^endothelia^high^ tumor niches. Accordingly, platelet-derived growth factor receptor β (PDGFRβ; avg_log2FC = 2.07; p-adj. < 0.001) and PDGFRα (avg_log2FC = 1.67; p-adj. < 0.001) as well as insulin-like growth factor family member 1 (IGFL1; avg_log2FC = 1.84; p-adj. < 0.001) were strongly upregulated genes within mCAF^high^endothelia^high^ niches (Supplementary Table S1).
A higher frequency of mCAFs is an independent prognostic factor for recurrence and distant metastasis in SDC
Cox regression analysis and log-rank test were performed to identify a potential influence of TME on the prognosis of SDC. The results of the univariate Cox regression model showed that a higher frequency of mCAFs was a prognostic factor for a higher rate of distant metastasis in SDC (p = 0.02; HR (95%CI) = 10.99 (1.37–88.16); Fig. 4A). Regarding RFP, a higher frequency of mCAFs (p = 0.02; HR (95%CI) = 6.23 (1.34–28.93)) and a higher frequency of other CD4^+^ cells (p = 0.04; HR (95%CI) = 3.82 (1.07–13.69)) were prognostic factors for a higher probability of recurrence in the univariate Cox regression model (Fig. 4B). Multivariate Cox regression revealed a higher frequency of mCAFs as an independent prognostic factor for a higher probability of recurrence (p = 0.032; HR (95%CI) = 7.81 (1.19–51.3); Fig. 4C). The five-year DCR was 75.0% for patients with a low frequency of mCAFs and 24.2% for those with a high frequency of mCAFs (p = 0.0049; Fig. 4D). The five-year RFP was 68.2% for patients with a low frequency and 18.2% for those with a high frequency of mCAFs (p = 0.0077; Fig. 4E). A high frequency of CN8 (mCAF-rich) was marginally non-significant as a prognostic factor for a higher probability of recurrence (p = 0.06; HR (95%CI) = 4.54 (0.94–22.01); Supplementary Fig. 4 A) and distant metastasis (p = 0.07; HR (95%CI) = 3.46 (0.91–13.10); Supplementary Fig. 4B) in the Cox regression model. Nevertheless, the five-year DCR was significantly lower in patients with a high frequency (27.7%) than in those with a low frequency of CN8 (68.2%; p = 0.04; Fig. 4F). Further, the five-year RFP was 20.4% for patients with a high frequency and 61.4% for those with a low frequency of CN8 (p = 0.052; Fig. 4G). The association between higher mCAF frequencies and a decreased DCR as well as RFP was consistent across AR subgroups (Supplementary Fig. 4 C and D). Finally, a significantly higher frequency of mCAFs was found in patients with distant metastatic disease than in those without (p = 0.048). Sex did neither show significance as prognostic factor for DCR (p = 0.23; HR (95%CI) = 0.28 (0.03–2.25), nor for RFP (p = 0.36; HR (95%CI) = 0.48 (0.10–2.28). When including survival data, a higher frequency of mCAFs was not a prognostic factor for OS (p = 0.90; HR (95%CI) = 1.10 (0.34–3.63) or RFS (p = 0.16; HR (95%CI) = 2.10 (0.74–5.91) in univariate cox regression. Furthermore, the five-year OS did not differ significantly between patients with a high frequency of mCAFs (45.5%) and those with a low mCAF frequency (58.3%; p = 0.87). There was no difference in the five-year RFS between patients with a high mCAF frequency (18.2%) and those with a low mCAF frequency (43.7%; p = 0.15). Notably, patients with low mCAF frequencies were markedly older (mean age = 69.17 years) than those with high mCAF frequencies (mean age = 63.64 years).Fig. 4. Association of cell types with distant control rate (DCR), recurrence-free probability (RFP), and frequency of distant metastasis in salivary duct carcinoma (SDC) patients. A Univariate cox proportional hazards model for DCR for cell types and T stage, stratified by median proportion or as negative vs. positive in case median equals zero. N stage was excluded due to complete separation, resulting in an unbounded 95% confidence interval for its odds ratio. B Univariate cox proportional hazards model for RFP for cell types and T stage, stratified by median proportion or as negative vs. positive in case median equals zero. N stage was excluded due to complete separation, resulting in an unbounded 95% confidence interval for its odds ratio. C Multivariate cox proportional hazards model for RFP for mCAFs, AR status (AR_high_ > 70% and AR_low_ ≤ 70 of tumor cells positive for AR), HER2 status, T stage, age, and other CD4^+^ cells. N stage was excluded due to complete separation, resulting in an unbounded 95% confidence interval for its odds ratio. D Kaplan–Meier plot with log-rank test for DCR comparing patients stratified as high and low based on the median proportion of mCAFs. E Kaplan–Meier plot with log-rank test for RFP comparing patients stratified as high and low based on the median proportion of mCAFs. F Kaplan–Meier plot with log-rank test for DCR comparing patients stratified as high and low based on the median proportion of CN8. G Kaplan–Meier plot with log-rank test for RFP comparing patients stratified as high and low based on the median proportion of CN8. H Boxplot displaying the proportion of mCAFs samples from patients with vs. without distant metastasis
When testing these results for validity within a previously published independent cohort with RNA-seq data from n = 67 SDC cases [23], we found that patients with low scores of the mCAF-like signature (mCAF_high_) had a significantly lower RFP (median RFP = 10.2 months; two-year survival = 23.5%) than mCAF_low_ patients (median RFP = 21.6 months; two-year-survival = 46.4%; p = 0.0294; Fig. 5G). Accordingly, mCAF_high_ patients had a significantly lower RFS (median RFS = 10.2 months, two-year survival = 23.5%) than mCAF_low_ patients (median RFS = 21.6 months; two-year-survival = 46.4%; p = 0.049; Fig. 5H). RFS was also less favorable in mCAF_high_ compared to mCAF_low_ patients when using mCAF-like signatures with different cut-offs (Supplementary Fig. 5 D-G).Fig. 5. Validation analyses with RNA-seq and Spatial Transcriptomics (ST) data. A Spearman correlation of the mCAF-like module score with the pericytes and the B) endothelial module score, leveraging ST data of three SDC specimens. C Expression of metastasis-, IL-6- and VEGF-associated module scores in mCAF^high^endothelia^high^ ST-spots. IL-6- and VEGF-related signatures were not significantly enriched. D Differential expression analysis contrasting mCAF^high^endothelia^high^ ST-spots with all other ST-spots. The top10 overexpressed genes are highlighted. E Gene set enrichment analysis of DE-results depicted in C) using “Cellular Component” and (F) “Molecular Functions” GO terms. G and H) Kaplan–Meier plots and log-rank tests depicting the prognostic impact of the transcriptomic mCAF-like signature (top20% mCAF marker) after median dichotomization of 67 samples through estimation of the recurrence-free probability (G) and recurrence-free survival (H); LN = lymph node; met. = metastasis
No prognostic tests were performed for entities other than SDC because of the low absolute number of events within those entities.
Discussion
Although studies have described the tumor immune microenvironment (TIME) in a limited number of SGC entities [32, 33], there are no data on the composition of the non-immune cellular TME in SGC. The non-immune cellular TME mainly consists of tumor vasculature and CAFs. CAFs have been shown to directly influence tumor cells but also indirectly interact with them through the ECM and immune cells [34]. Previous research has comprehensively demonstrated the pro-carcinogenic role of CAFs [35] but approaches specifically targeting CAFs have not been successful to date [13]. Therefore, a more precise delineation of CAFs is necessary. Although different CAF subsets have been discovered in various solid tumor entities within the last years [8, 10, 36], a study by Cords et al. recently introduced a classification system for CAFs based on scRNA-seq data and demonstrated its applicability using IMC [5, 12]. In the present study, we are the first to comprehensively decipher the cellular TME and evaluate the influence of cellular TME components on the prognosis of SGC.
We found that the highest mean frequency of mCAFs was observed in salivary duct carcinoma (SDC) samples compared to other entities. Importantly, in SDC, a higher mCAF frequency was an independent prognostic factor for a higher probability of recurrence and distant metastasis, and patients with a higher mCAF frequency had significantly worse DCR and RFP than those with a lower frequency of mCAFs. The unfavorable prognostic influence of high mCAF frequencies was confirmed in an independent SDC cohort leveraging transcriptomic data. In the spatial interaction analysis, the mCAF-rich CN8 showed a strong contribution from endothelia. This is in line with previously published data that showed a positive spatial interaction of mCAFs with “blood,” high endothelial venules, and lymphatic channels in other solid tumors [5, 12]. We also observed a negative prognostic effect of mCAF-rich CN8. However, regarding the stronger prognostic potential of mCAFs, this effect might be driven by mCAFs and diluted by other cells within the cellular neighborhood. These findings suggest that mCAFs are associated with hematogenous metastatic spread in SDC. A prior study has identified an anti-proportional distribution of mCAF and immune cell density in NSCLC and has hypothesized that mCAFs may negatively influence prognosis by shielding tumor cells from immune cells [12]. However, our study is the first to suggest a potential role of mCAFs in hematogenous metastatic spread. More specific, we were not only able to show that mCAFs and the spatial co-localization of mCAFs and endothelia is associated with a significantly worse outcome in SDC, but also that mCAF^high^endothelia^high^ tumor niches show enrichment of metastasis-associated gene signatures using ST data and that enrichment of these signatures is higher in mCAF^high^endothelia^high^ than in mCAF^high^endothelia^low^ niches, thereby supporting this clinical finding with an independent cohort and an additional method on the spatial-biological level. Importantly, we also found factors potentially mediating this interaction. The “platelet-derived growth factor binding” and “insulin-like growth factor binding” signatures were among the top signatures enriched and the respective genes PDGFRβ, PDGFRα, and IGFL1 were among the top-enriched genes in mCAF^high^endothelia^high^ tumor niches. Interestingly, platelet-derived growth factor signaling enriched in stromal cells has previously been shown to mediate brain metastasis in breast cancer [37] and multiple tyrosine kinase inhibitors such as sorafenib and sunitinib inhibiting PDGFRβ and PDGFRα are available and FDA-approved for advanced renal cell carcinoma [38]. Similarly, insulin-like growth factor secreted by CAFs was linked to metastasis in different solid tumor entities [39]. Although several small-molecules and monoclonal antibodies targeting insulin-like growth factor showed promising preclinical results, none of those has reached approval. Overall, our findings deciphering the interaction mediators between mCAFs and endothelia are in line with previous studies. Targeting mCAFs and the pathways which we have found to be associated with a pro-metastatic mCAF^high^endothelia^high^ niche could therefore be a new perspective for the development of targeted therapies for SDC.
SDC is the clinically and biologically most aggressive SGC entity, with a five-year OS and RFS of approximately 50% and 30%, respectively, mainly due to distant recurrence [40–43]. Although systemic therapy options targeting AR and HER2 are available for SDC, the median progression-free survival for drugs targeting both receptors is below nine months [44, 45]. Therefore, new prognostic markers and therapeutic targets are urgently needed in this highly aggressive entity to guide treatment decisions in the curative setting and to improve outcomes in the R/M setting. Given the results of this study, the frequency of mCAFs could serve as a new prognostic marker for SDC. Dual immunohistochemical staining of SMA and Collagen1 could facilitate the quantification of mCAFs in FFPE samples at first diagnosis and thereby the stratification of patients into risk groups. Moreover, a targeted approach leading to a lower frequency of mCAFs in the primary tumor may decrease distant metastatic potential. To date, most clinical trials targeting CAFs have shown negative results [46]. This is most likely because the relevant drugs were targeting CAFs in general rather than specific pro-tumorigenic subpopulations. For example, the FAP-targeting monoclonal antibody sibrotuzumab and the FAP-inhibiting small-molecule inhibitor talabostat did not show any objective response in phase II trials for metastatic colorectal cancer [47, 48]. FAP has been found to be expressed in mCAFs, iCAFs, tCAFs, ifnCAFs, apCAFs, and dCAFs (among others) [5], some of which showed pro-tumorigenic and others anti-tumorigenic features [12]. Identification of tumor-specific pro-tumorigenic subpopulations, such as mCAFs in SDC, and targeting these more specifically may help overcome the challenges faced in previous CAF-targeted therapies and lead to more effective treatment strategies.
Interestingly, we found a significantly higher frequency of tumor-infiltrating CAFs in AR_high_ than AR_low_ SDC cases. Although the impact of CAFs on tumor metabolism and signaling caused by cytokines and chemokines has been proven previously [49], there is currently no evidence showing the influence of CAFs on the AR status of tumor cells. While our data do not allow to draw causal conclusions in terms of the interplay between CAFs and AR expression in tumor cells, SDC patients with AR_high_ tumors could particularly benefit from CAF-targeted therapies or therapies combining androgen blockade and CAF-targeting.
With regard to the TIME, we found significant differences in the abundance of immune cell subsets across entities. This was due to higher levels of CD4^+^ T cells, CD8^+^ T cells, other CD4^+^ cells, and unclassified immune cells, leading to a significantly higher overall immune infiltration in ACC than in the other entities. The propensity for immune infiltration has been described as a characteristic of ACC in H&E diagnostics of these tumors [50]; however, the exact cell types that contribute to this phenomenon have not yet been studied. Moreover, ACC has not been included in studies comprehensively describing the TIME in SGC [32, 33]. Although ACC generally has a favorable prognosis, there is a subset of ACC patients with recurrence, resulting in a five-year RFS of approximately 76% [51] and requiring systemic therapeutic approaches. To date, there are no established therapeutic targets for ACC, and conventional chemotherapy has been shown to have low efficacy in SGC [52]. Of note, programmed death (-ligand) 1 (PD(L)1)checkpoint blockade has not shown promising results in R/M SGC, with an objective response rate (ORR) of approximately 10% across entities [53, 54]. This may be due to the relatively low expression of PD-L1 in SGC tumor cells, with a mean tumor proportion score (TPS) of approximately 4% across entities [55]. According to our results, patients with ACC, compared to other SGC entities, could most likely benefit from new immunotherapeutic approaches that specifically target the TIME.
Furthermore, we found a strong positive correlation between the frequencies of mCAFs, SMA CAFs, and overall CAFs (classic myofibroblastic CAF phenotypes) and the CAF ECM protein module. This finding validates the previously suggested production of CAF ECM proteins by myofibroblastic CAF subpopulations [15]. Additionally, all immune cell subsets, except for plasma cells, showed a negative correlation with the CAF module, which is consistent with the established immune-exclusive features of a stiffened carcinoma-associated ECM [56]. Interestingly, the associations between the Hem module and cellular TME components, including an anti-correlation with endothelia, were very similar to those of the CAF module and its associations. We previously reasoned that the Hem module might consist of blood-related factors deposited within the ECM. However, the present data indicate that this does not imply exaggerated vasculature or deposition of the cellular components of the peripheral blood. The BM module had been shown to be mainly expressed in SGC with myoepithelial differentiation (e.g., adenoid cystic carcinomas) [15] which were not included in the present study. Together, these data clearly link SGC mCAFs, Collagen CAFs, and SMA CAFs to an increase in classical CAF-associated ECM proteins and provide evidence that the CAF and the Hem module are part of an immune-exclusive ECM.
This study has several limitations. First, SGG cases were identified retrospectively, which introduces the possibility of selection bias. Additionally, although relevant markers were included in the panel, we were unable to identify vCAFs, rCAFs, ifnCAFs, and iCAFs because of non-specific IMC signals or the unavailability of appropriate antibodies. Consequently, it is unclear whether these CAFs were not present in SGC or whether the absence of signals was due to the lack of specific metal-conjugated antibodies. Additionally, we did not include SGC types with (partial) myoepithelial differentiation, such as adenoid cystic carcinoma, because we identified a considerable overlap of myoepithelial and CAF marker expression during pretests, which would have resulted in high uncertainty in cell type identification. Lastly, the analyses were based on a limited number of cases and are therefore exploratory, which is a common problem in rare tumors, specifically in SGC showing a high number of biologically distinct entities. The sample size of n = 54 patients reflects the rarity of SGC and the challenge of assembling larger cohorts with available high-quality tumor material and multiplexed imaging data. Due to the limited prior knowledge of effect sizes in cellular tumor microenvironment analyses in SGC, a formal prospective sample size or power calculation was not performed. This limitation may impact the statistical power of some analyses and the generalizability of the findings. Therefore, our results should be considered hypothesis-generating. However, we addressed this limitation by leveraging a second independent clinically-annotated SDC cohort and two other methodologies (ST, bulk RNA-sequencing). Further, it must be mentioned that, due to the limited number of events, the multivariable Cox regression analysis is exploratory and that true independence awaits larger cohorts. Also, we acknowledge that the results of this study are correlative due to no functional experiments performed and therefore causality cannot be inferred from the presented cross-sectional data. Lastly, since only three ST samples were analyzed, these additional analyses and the results derived thereof are exploratory and require validation in future studies with larger collectives.
This is the first study to comprehensively decipher the cellular TME in several of the more frequent SGC entities. To the best of our knowledge, this is also the first study to use multiplexed imaging mass cytometry for spatial single-cell analyses in SGC. We were able to identify a high frequency of mCAFs as a prognostic factor for distant metastasis and further link mCAFs to hematogenous spread by spatial clustering in SDC. Moreover, we were able to validate the prognostic value of a high mCAF frequency in an independent SDC cohort. The frequency of mCAFs may be therefore used to stratify patients with SDC into risk groups. Using ST data, we further validated our results be showing enrichment of metastasis-associated signatures in mCAF^high^endothelia^high^ spots and deciphered potential interaction mediators within these tumor niches. Systemic therapy approaches targeting mCAFs in the TME of SDC could lead to improved prognosis in this highly aggressive tumor entity in the future. Moreover, we were able to demonstrate higher immune infiltration and cellular composition in ACC than in other SGC entities. Therefore, ACC may be amenable to new immunotherapeutic approaches. Overall, our results lay the groundwork for the improvement of diagnostics and therapy for these rare but clinically challenging carcinomas.
Supplementary Information
Supplementary Material 1: Supplementary Figure S1: (A) tSNE and (B) UMAP projection of all analyzed cells colored by TMA. (C) Overall cell type and (D) cell category frequency. (E) Cell category frequency of primaries as well as metastases and (F) central and tumor compartments of the primary (all entities pooled). Supplementary Figure S2: (A) Abundance of CAF subsets in SDC versus other tumor types. (B) Examples of tumor patch detection (tumor patches are colored, other cells depicted in grey). (C-E) Frequency of immune cells, CAFs and tumor-infiltrating CAFs compared between ARhigh and ARlow SDC. (F-I) and between SDC with different HER2 expression levels. (J) Frequency of tumor-infiltrating immune cells in SGC entities. Supplementary Figure S3: (A) Examples of cellular interaction graphs via k nearest neighbor detection with k = 20; colors correspond to cell types. (B-D) exemplary ROIs with co-localization of cells that are characteristic for cellular neighborhoods 2, 6, and 5. 100-micron scale bars. Supplementary Figure S4: Univariate Cox proportional hazards model for (A) DCR and (B) RFP for cluster numbers (Cn) stratified by median proportion or as negative vs. positive in case median equals zero. Kaplan-Meier plots for C) DCR and (D) RFP with log-rank tests comparing patients with high and low mCAF frequencies stratified by AR status (ARhigh >70% and ARlow ≤70 of tumor cells positive for AR). Supplementary Figure S5: A) Spatial expression of the mCAF, pericytes, and endothelia module scores. B) Fraction of spots of each ST-sample within each spot category as defined by expression of the module scores for mCAFs and endothelia. C) Gene set enrichment analysis of DE-results depicted in Figure 5 C using “Biological Processes” GO terms. D-G) Estimation of recurrence-free survival in 67 SDC patients which were median-dichotomized using mCAF signatures derived with different log2FC cutoffs (Cords et al., Nat Com., 2023): (D) top 5 genes, (E) top 5%, (F) top 10% and (G) top 100 genes. Log-rank tests. LN = lymph node; met. = metastasis. Supplementary Figure S6: A) and B) Pathway-enrichment analyses analogous to Figure 5 C (20% cut-off) using 10% (A) and 30% (B) cut-offs for designation of spots as mCAFhigh or endotheliahigh, leading to similar results.Supplementary Material 2.