Article Text

Original research
Multimodal profiling of chordoma immunity reveals distinct immune contextures
  1. Siddh van Oost1,2,
  2. Debora M Meijer1,2,
  3. Marieke E Ijsselsteijn1,
  4. Jessica P Roelands1,
  5. Brendy E M W van den Akker1,
  6. Ruud van der Breggen1,
  7. Inge H Briaire-de Bruijn1,
  8. Manon van der Ploeg1,
  9. Pauline M Wijers-Koster1,
  10. Samuel B Polak3,
  11. Wilco C Peul3,
  12. Robert J P van der Wal4,
  13. Noel F C C de Miranda1,2 and
  14. Judith V M G Bovee1,2
  1. 1Department of Pathology, Leiden University Medical Center, Leiden, Netherlands
  2. 2Leiden Center for Computational Oncology, Leiden University Medical Center, Leiden, Netherlands
  3. 3University Neurosurgical Center Holland, Leiden University Medical Center, Leiden, Zuid-Holland, Netherlands
  4. 4Department of Orthopaedic Surgery, Leiden University Medical Center, Leiden, Netherlands
  1. Correspondence to Professor Judith V M G Bovee; j.v.m.g.bovee{at}lumc.nl

Abstract

Background Chordomas are rare cancers from the axial skeleton which present a challenging clinical management with limited treatment options due to their anatomical location. In recent years, a few clinical trials demonstrated that chordomas can respond to immunotherapy. However, an in-depth portrayal of chordoma immunity and its association with clinical parameters is still lacking.

Methods We present a comprehensive characterization of immunological features of 76 chordomas through application of a multimodal approach. Transcriptomic profiling of 20 chordomas was performed to inform on the activity of immune-related genes through the immunologic constant of rejection (ICR) signature. Multidimensional immunophenotyping through imaging mass cytometry was applied to provide insights in the different immune contextures of 32 chordomas. T cell infiltration was further evaluated in all 76 patients by means of multispectral immunofluorescence and then associated with clinical parameters through univariate and multivariate Cox proportional hazard models as well as Kaplan-Meier estimates. Moreover, distinct expression patterns of human leukocyte antigen (HLA) class I were assessed by immunohistochemical staining in all 76 patients. Finally, clonal enrichment of the T cell receptor (TCR) was sought through profiling of the variable region of TCRB locus of 24 patients.

Results Chordomas generally presented an immune “hot” microenvironment in comparison to other sarcomas, as indicated by the ICR transcriptional signature. We identified two distinct groups of chordomas based on T cell infiltration which were independent from clinical parameters. The highly infiltrated group was further characterized by high dendritic cell infiltration and the presence of multicellular immune aggregates in tumors, whereas low T cell infiltration was associated with lower overall cell densities of immune and stromal cells. Interestingly, patients with higher T cell infiltration displayed a more pronounced clonal enrichment of the TCR repertoire compared with those with low T cell counts. Furthermore, we observed that the majority of chordomas maintained HLA class I expression.

Conclusion Our findings shed light on the natural immunity against chordomas through the identification of distinct immune contextures. Understanding their immune landscape could guide the development and application of immunotherapies in a tailored manner, ultimately leading to an improved clinical outcome for patients with chordoma.

  • Sarcoma
  • Immunotherapy
  • Tumor Microenvironment
  • T-Lymphocytes
  • Dendritic Cells

Data availability statement

Data are available in a public, open access repository. The RNA sequencing data generated in this study are publicly available on GEO (GSE239531). The imaging mass cytometry data is available on BioImage Archive (S-BIAD830). The RNA sequencing data for the other studied sarcoma subtypes is unpublished and therefore not yet available, with the exception of the osteosarcomas (GSE237033). The used code and processed data will be published on GitLab for reproducibility. All other data are available from the corresponding author on reasonable request.

http://creativecommons.org/licenses/by-nc/4.0/

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.

Statistics from Altmetric.com

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • Chordomas present with conspicuous immune infiltration and can respond to immune checkpoint blockade despite their low mutational burden.

WHAT THIS STUDY ADDS

  • This study has delineated discrete subgroups within chordomas characterized by coordinated antiumor and inflammatory responses, prominently featuring clusters of T cells and antigen presenting cells, particularly dendritic cells. Notably, our findings highlight that the majority of chordomas retain functional antigen processing machinery, rendering them viable candidates for T cell-mediated immunotherapeutic interventions.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • This work not only opens avenues for potential immunotherapeutic strategies tailored to chordomas but also establishes a valuable resource for fellow researchers, offering a unique fusion of transcriptomic insights and spatial immunophenotyping across a uniquely large cohort of patients with chordoma.

Background

Chordomas are cancers displaying notochordal differentiation that arise predominantly in the axial skeleton at the skull base, the craniocervical and mobile spine, the sacrum, and, occasionally, at extra-axial areas.1 General treatment for chordomas consists of surgical resection and/or (neo-)adjuvant proton beam radiation. Since most chordomas are indolent tumors, chemotherapy is considered ineffective for the treatment of this disease.2 The clinical management of patients with chordoma is complicated by the anatomical location of the tumors, their often late detection, as well as treatment-related morbidities.2 3 Aiming at a gross total or en bloc resection involves sacrifice of neurovascular structures with permanent neurological disabilities. As such, disease recurrence is common after treatment and 30–40% of patients with chordoma develop metastases.4 The overall 5-year survival following diagnosis is around 67% and there are currently no curative therapies for patients with chordoma with advanced disease.1 5 This reinforces the need for innovative treatments as local recurrences and distant metastases persist.

Recent findings have provided evidence for the presence of conspicuous immune cell infiltration in chordomas6 7 and, encouragingly, clinical responses of patients with chordoma to checkpoint blockade immunotherapy targeting cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) and/or the programmed cell death protein 1 (PD-1)/programmed death-ligand 1 (PD-L1) axis have been reported.8–12 This is an intriguing observation as chordomas display a low mutational burden and are, instead, affected by gross chromosomal alterations.13–15 Furthermore, PD-L1 expression is rarely encountered in sarcomas, including chordomas.16–18 Instead, the presence of pre-existing antitumor immune responses perceived by the infiltration of cytotoxic T cells in the tumor tissue or by the presence of tertiary lymphoid structures appear to inform on the sensitivity of soft tissue sarcomas to immune checkpoint blockade.19 20 So far, whether tertiary lymphoid structures are predictive of responsiveness to immune checkpoint blockade has not been investigated in chordomas.

An in-depth immunological characterization of chordomas and their association with clinical parameters is currently lacking. Here, we applied transcriptomic profiling and multidimensional immunophenotyping by imaging mass cytometry and multispectral immunofluorescence to resolve the immune landscape of chordomas. Chordomas presented hallmarks of antitumor immunity and were considered immune “hot” in comparison with other sarcomas. Within our cohort, we identified groups of chordomas with distinct immune contextures, independent from any clinical parameter. Our findings support the rational utilization of immunotherapeutic approaches, tailored to target these distinct microenvironments in chordoma.

Methods

Patient samples

Formalin-fixed paraffin embedded (FFPE) and snap-frozen tissue material was initially collected for 43 and 24 chordomas, respectively, from 32 patients, followed by a second expansion cohort containing FFPE material of an additional 53 patients. Regions of interest (ROIs) representative for the tumor’s immune microenvironment were selected on H&E sections of decalcified FFPE tissue by a bone and soft tissue tumor pathologist (JVMGB). These ROIs were then used to create tissue microarrays (TMAs) using a TMA Master (3DHISTECH). These TMAs were holding two to four cores of 1.6 mm in diameter per tumor, dependent on the source of the material or the morphology of the tumor (two cores: biopsy; three cores: resection/excision, etc; four cores: heterogeneous morphology, or separate areas of dedifferentiation). As controls and for orientation purposes, each TMA contained three cores tonsil and three cores placenta, either decalcified with EDTA or formic acid or non-decalcified.

RNA sequencing data

RNA was isolated from snap-frozen tissue from 20 primary tumors with TRIzol and isopropanol/ethanol. In addition, RNA sequencing data from primary tumors of patients with osteosarcoma (n=7; GSE237033), patients with chondrosarcoma (n=9), patients with undifferentiated soft tissue sarcoma (n=13), and patients with myxofibrosarcoma (n=10) was taken along as a means of comparison. RNA purification was performed with the RNeasy kit which included treatment with DNase. Paired-end 150 base pair (bp) reads were created on a NovaSeq6000 Illumina at GenomeScan (Leiden, The Netherlands). The sequencing data was processed with the RNA sequencing BioWDL pipeline from the SASC with default settings (BioWDL Github, LUMC, The Netherlands). Alignment of the 150 bp reads to hg38 was performed with STAR, after which the gene expression was quantified with HTSeq-count.21 The processed sequencing data was stored as a matrix containing the counts per gene per sample.

Gene expression analysis

Normalization and exploratory analysis of the RNA sequencing data was performed in R (V.4.0.2) by using DESeq2 (V.1.30.1). BiomaRt (R, V.2.46.3) was used to translate the ensemble IDs to known protein coding genes. The most variably expressed genes were determined by calculating the variance per gene across all samples. The immunologic constant of rejection (ICR) gene signature was used to describe the activity state immune microenvironment and to identify distinct immune-related clusters (Th1-related markers: IFNG, STAT1, IL12B, IRF1, TBX21; CD8+ T cells: CD8A, CD8B; immune effector molecules: GZMA, GZMB, GZMH, PRF1, GNLY; chemokine ligands: CXCL9, CXCL10, CCL5; immune suppressive molecules: IDO1, PDCD1, CD274, CTLA-4, FOXP3).22 The mean of the log2 normalized expression of these 20 genes was used to attribute an ICR “score” to each sample for comparison between the sarcoma subtypes, which was visualized with ggplot2 (R, V.3.3.6). For both the heatmaps in figure 1, a Z-score was calculated per gene to visualize the transcriptional profile of the samples by using the Ward variance method with ComplexHeatmap (R, V.2.11.1). Differential gene expression analysis was performed with DESeq2 comparing chordomas to all other sarcomas. Genes were considered differentially expressed when the log2 fold change was above 2 and the adjusted p value was below 0.01. Immune-related genes were selected for visualization in the enhanced volcano plot and additional boxplots. For visualization purposes, the axes of the enhanced volcano plot were manually adjusted as follows: log2 fold change (x-axis) was set between −10 and 10, and the −log10 adjusted p (y-axis) was set to a maximum of 100.

Figure 1

Chordomas display hallmarks of antitumor immunity. (A) Heatmap demonstrating the log2 normalized gene expression Z-score of the top 200 most variably expressed genes in all sarcoma samples (n=59). Samples are annotated according to their histotype. Genes associated with chordoma samples are highlighted on the right. These genes include chordoma-specific genes such as TBXT, CD24 and several keratins. In addition, immune-related genes are highlighted, including human leukocyte antigen (HLA*) and immunoglobulin (IG*) genes. (B) Heatmap presenting the Z-score for the log2 normalized expression of the 20 ICR signature genes in 20 chordoma samples. Based on the expression Z-score, the samples are classified into two groups through unsupervised hierarchical clustering: ICR high (n=17) and ICR low (n=3). (C) Boxplots presenting the immunologic constant of rejection (ICR) score per sarcoma subtype in descending order. (D) Kaplan-Meier curves displaying disease-specific survival for the ICR high and ICR low chordomas. Two-sided log-rank tests were performed for the Kaplan-Meier estimates. ECM, extracellular matrix; USTS, undifferentiated soft tissue sarcoma.

Imaging mass cytometry

Conjugation of Bovine Serum Albumin (BSA)-free antibodies and immunodetection with a 40-marker tumor immune microenvironment panel (online supplemental table 1) was performed as previously described.23 In brief, 4 µm sections of the created TMAs were deparaffinized and rehydrated, after which heat-induced antigen retrieval was performed in 1× low pH antigen retrieval solution (pH6, Thermo Fisher Scientific). Subsequently, the first batch of primary antibodies were incubated overnight at 4°C. After washing in phosphate-buffered saline (PBS) supplemented with 1% BSA and 0.05% Tween, conjugated anti-mouse and anti-rabbit antibodies were incubated for 1 hour at room temperature. Next, the sections were incubated with the second batch of primary antibodies for 5 hours at room temperature, followed by an overnight incubation at 4°C with the final batch of primary antibodies. To finish, 1.25 µM of the DNA intercalator iridium (Fluidigm) was incubated for 5 min to stain DNA, after which the sections were washed, air-dried and stored at room temperature. ROIs of 1000×1000 µm were selected with consecutive H&E slides and ablated by using a Hyperion mass cytometry imaging system (Fluidigm) at the Flow Core Facility (LUMC, Leiden) within a month after staining. The imaging data was acquired with CyTOF Software (V.7.0), visually inspected and then exported with MCD Viewer (V.1.0.5). In the end, 127 tumor ROIs were successfully ablated and analyzed.

Supplemental material

Imaging mass cytometry analysis

To correct for variability in signal-to-noise ratios between FFPE sections, the signal intensity of the generated images was normalized at pixel levels by using MATLAB (V.R2021a), after which the images were binarized using semi-automated background removal in Ilastik (V.1.3.3) as described previously.23 Probability masks were generated in Ilastik by using a multicolored image including the DNA channel, CD45 channel, keratin channel and vimentin channel to annotate nucleus, cytoplasm and background. Probability masks were loaded into CellProfiler (V.2.2.0) and used for creating cell segmentation masks. Next, the segmentation masks and binarized images were loaded into ImaCytE24 to generate single cell-FCS files containing relative frequency of positive pixels for each marker, which were then used for phenotyping in Cytosplore (V.2.3.1) by sequentially using hierarchical stochastic neighbor embedding and t-distributed stochastic neighbor embedding. Identified cell phenotypes were loaded into ImaCytE together with the segmentation masks and binarized images, which allows the mapping of the phenotypes back onto the images. Phenotypes per image were extracted to calculate the mean cell density per mm2 for every sample (one image=1 mm2). An overview of the lineage markers used for phenotype identification is listed in online supplemental table 2. Patients were clustered into two groups based on the Z-score normalized mean cell densities per mm2 of all five identified T cell phenotypes by using the Ward variance method with ComplexHeatmap. Ggplot2 was used to visualize the mean cell densities and to compare them between both T cell groups. For statistical testing, a two-sided Student’s t-test was performed followed by a Benjamini-Hochberg false discovery rate (FDR) correction. In addition, an overview heatmap containing all identified phenotypes in the microenvironment of chordomas was visualized with ComplexHeatmap.

Spatial interaction analysis

After phenotype calling, a spatial interaction analysis was performed with ImaCytE for all imaged samples, comparing T cell low with T cell high chordomas. ImaCytE uses a 1000-iteration permutation test which calculates a Z-score for the interaction between adjacent cells. However, permutation testing does not adjust for overabundant cell types, such as tumor cells, and is therefore not able to identify meaningful interactions for tumor cells. Because of the range of the calculated Z-scores, interactions with a Z-score below 10 were considered not significant and therefore removed for visualization only. If an interaction Z-score was higher than 10 in one group and lower than 10 in the other group, this interaction was considered significant for one group. When both Z-scores were higher than 10, this interaction would be considered significant in both groups. A comparative interaction heatmap between all phenotypes excluding tumor cells was visualized with ggplot2.

Multispectral immunofluorescent imaging of T cells

A T cell panel comprising five markers for five channels was generated for imaging with a Vectra (Akoya Biosciences) as previously described25 for validation purposes. This panel included FoxP3 (236A/E7, 1:1000, Invitrogen, 14-4777-82), CD8α (D8A8Y, 1:100, Cell Signaling Technology, #85336), CD3ε (EP449E, 1:50, Abcam, ab52959), pan-cytokeratin (AE-1/AE-3 & C11, both 1:50, Novus & Santa Cruz, NBP2-33200AF647 & SC-8018) and Ki-67 (8D5, 1:200, Cell Signaling Technology, #9449) and 4′,6-diamidino-2-phenylindole (DAPI) as a nuclear staining. In short, slides were deparaffinized and incubated in 0.3% H2O2 in methanol for 20 min to block endogenous peroxidase. After antigen retrieval in citrate (pH 6.0), slides were incubated for 30 min with SuperBlock Blocking Buffer (Thermo Fisher Scientific). Next, slides were incubated with primary antibody (FoxP3) for 1 hour at room temperature and subsequently with BrightVision DPVO-HRP (Immunologic) for 30 min at room temperature. Opal 520 reagent (Akoya Bioscience) was diluted in Opal amplification diluent and then added to the slides for 1 hour at room temperature, after which the slides were heated in citrate buffer (pH 6.0) for 15 min. The second round of primary antibodies (CD8α and Ki-67) were incubated overnight at room temperature and then stained with their respective fluorescent secondary antibodies (CF555, 1:400, Biotum and Alexa 680, 1:200, Invitrogen) for 1 hour at room temperature. Finally, the slides were stained with directly conjugated antibodies (CD3ε-Alexa 594 and pan-cytokeratins-Alexa 647) for 5 hours at room temperature and then mounted. All TMAs were imaged at 20×, after which the images were processed with inForm (V.2.4), which comprises background removal for each antibody by normalization of the spectra, and analyzed with QuPath (V.0.3.1). Object classification was performed by using the following classes: CD3+FOXP3+=regulatory T cells, CD3+CD8+=CD8+ T cells, CD3+=CD4+ T cells, CD3+Ki-67+=proliferating T cells, CD3=rest. Counts per image (one image=1 mm × 1.3 mm) were normalized per patient to counts per mm2 tissue. Patients were divided into a T cell high and low cluster based on the unsupervized clustering of the Z-scores of all T cell phenotypes in R.

Survival analysis

For all patients with chordoma, Kaplan-Meier curves for the disease-specific survival and recurrence-free survival were generated with the survminer and survival packages (R, V.0.4.9 and V.3.1–12, respectively). The disease-specific survival was considered the survival of a patient until death due to the disease, whereas the recurrence-free survival was considered the time until a patient was diagnosed with a recurrence, thereby excluding patients who did not receive surgery. Since the cohort included many censored cases, the log-rank p value was used to estimate the significance of the associations. Known prognostic factors such as the anatomical tumor location excluding the only extra-axial patient, age and type of treatment, surgical margin excluding patients who did not receive surgery and disease recurrence were investigated for their prognostic value, as well as the newly identified T cell infiltration groups. Univariate and multivariate analyses were performed in R according to the Cox proportional hazard model. Only clinical parameters that were found significant in the univariate analysis were taken along in the multivariate analysis.

Immunohistochemical detection of human leukocyte antigen class I and β-2m

Immunohistochemistry was performed for the human leukocyte antigen (HLA) class I heavy chain (HCA2, 1:1000, Nordic Bio, MUB2036P; HC10, 1:3200, Nordic Bio, MUB2037P) and β−2m (EPR21752-214, 1:4000, Abcam, ab218230) on all generated TMAs as described previously.26 Briefly, antigen retrieval was performed in citrate (pH 6.0), after which the TMA sections were preincubated with PBS-5% non-fat dry milk for 30 min at room temperature. Primary antibody dilutions were incubated overnight at 4°C, after which the sections were incubated with BrightVision DPVO-HRP (Immunologic) for 30 min at room temperature. Finally, TMA sections were further developed with liquid 3,3’-diaminobenzidine (DAB)+chromogen (Dako) for 10 min and subsequently counterstained with hematoxylin. Samples were classified as HLA class I positive, weak, or defective depending on the expression pattern of HC10 and HCA2 in the tissue and the contrast between the expression of the tumor cells and the immune cells. Tumors were classified as defective if one or more of the HLA class I antibodies were negative in the tumor cells. The distribution of the expression patterns among the identified T cell groups was visualized in R using ggplot2.

T cell receptor β-chain sequencing

DNA was isolated from snap-frozen tissue from 27 chordomas from 23 patients for sequencing of the variable region of the TCRB locus. In brief, DNA libraries were generated with the AmpliSeq Library Kit Plus and the OncoMine TCR-β-SR Assay following the manufacturer’s instructions. Sample libraries were measured with the ION Library TaqMan Quantitation Kit and then pooled to a final concentration of 25 pM. Next, the pooled libraries were templated on ION 540 Chips with the Ion Chef System. The 80 bp reads were created on an Ion GeneStudio S5 Prime Sequencing Technology at GenomeScan (Leiden). The number of clones and the Shannon diversity, a measure for clonal enrichment, were calculated in R by using productive counts only, meaning that a clone needed at least 10 plus and 10 minus counts. Because recurrences showed similar levels of T cell receptor (TCR) diversity as their respective primary tumor, we used the TCR repertoire of recurrences for three patients for which there was no frozen material of primary tumors available. Finally, the data was visualized in R using ggplot2.

Results

Patient characteristics

In this study, two cohorts were established to investigate the immune microenvironment of chordomas. The first cohort included primary tumors of 32 patients and was explored using RNA sequencing (n=20) and imaging mass cytometry (n=32). The median age was 59 years (range 1–80 years) and the median follow-up time was 66 months (range 0–418 months). In total, the primary tumors in the cohort comprised 8 skull base chordomas (25%), 7 spinal chordomas (22%) and 17 sacrococcygeal chordomas (53%). The second cohort included primary tumors of an additional 53 patients. The median age of this cohort was 60 years (range 22–80) and the median follow-up time was 41 months (range 0–182). In total, this cohort comprised 17 skull base (32%), 15 spinal chordomas (28%), 20 sacrococcygeal chordomas (38%) and 1 extra-axial chordoma of the ribs (2%). The second cohort also comprised patients solely treated with radiotherapy, for which the treatment-naïve biopsies were used. All profiled primary tumor samples were annotated for their anatomical location as well as whether they underwent neoadjuvant radiotherapy. Further patient characteristics of both cohorts are listed in table 1.

Table 1

Patient characteristics of chordoma cohorts. This table contains an overview of each cohort separately as well as both cohorts combined

Chordomas present transcriptional hallmarks of anticancer immunity

The transcriptomic profiles of the 20 chordomas were compared with those of osteosarcomas (n=7), chondrosarcomas (n=9), undifferentiated soft tissue sarcomas (n=13), and myxofibrosarcomas (n=10). All chordoma samples clustered separately from the remaining sarcomas following unsupervized clustering based on the 200 genes with most variable expression among all samples (figure 1A). Sample clustering was largely driven by features related to the line of differentiation of the different sarcoma samples, with known chordoma-related genes such as TBXT, CD24, and several keratins enriched in the chordomas27–29 (figure 1A, online supplemental table 3). Moreover, genes associated with the extracellular matrix and adhesion-related genes were enriched in chordomas and chondrosarcomas as compared with the other sarcomas (figure 1A, online supplemental table 3). Interestingly, several immune-related genes, including HLA class I, HLA class II, PLA2G2D and immunoglobulin genes, were among the most variably expressed genes and were enriched in specific chordomas. Furthermore, differential gene expression analysis revealed enrichment for other immune-related genes in chordomas in comparison with the other sarcomas (online supplemental figure 1A,B, online supplemental table 4).

Supplemental material

Supplemental material

To substantiate that immune-related processes are enriched in chordomas, we evaluated the expression of genes included in the ICR gene signature to compare between the distinct sarcoma subtypes.22 This 20 gene signature reflects the activation of type 1 T (Th1) cell signaling, expression of CXCR3/CCR5 chemokine ligands, cytotoxicity and counter-activation of immunoregulatory mechanisms and was previously shown to be an independent prognostic factor for metastasis-free survival in soft tissue sarcomas.30 Most chordomas are characterized by a high expression of ICR genes (ICR high) (figure 1B), particularly when comparing the ICR score (average of the 20 ICR genes) with samples from other sarcoma types (figure 1C). Furthermore, the ICR high immune subtype was associated with improved disease-specific survival in this cohort, although only three cases were classified as ICR low (log-rank p=0.0055; figure 1D).

T cell infiltration in chordomas correlates with ICR score and separates patients into distinct groups

To confirm and expand our initial observations, we performed imaging mass cytometry to characterize the immune cell contexture of 32 primary tumors. Twenty-nine cell populations were identified including five T cell phenotypes comprising CD4+ T cells, CD8+ T cells, Ki-67+ T cells, regulatory T cells and (unspecified) T cells including low numbers of CD57+ and gamma delta T cells (online supplemental figure 2). Since the ICR transcriptional signature is largely reminiscent of T cell-driven immune responses, we investigated the association between T cell infiltration and the ICR score. Indeed, T cell infiltration was more abundant in ICR high chordomas as compared with ICR low samples (online supplemental figure 3A). Of note, T cells were detected in all chordomas although the levels of infiltration were highly variable, both within and between lesions, ranging from 6.5 to 284.6 cells/mm2 (median=51 cells/mm2) between tumors (figure 2A,B). Although the vast majority of T cells were confined to the stromal compartment surrounding the tumor lobules, most tumors also show infiltration of a low number of T cells within the cancer cell compartment. Other commonly observed immune cell populations included various subsets of myeloid cells, such as monocytes, macrophages, dendritic cells, and granulocytes (online supplemental figure 2).

Figure 2

Chordomas present distinct groups based on T cell infiltration. (A–B) Heterogeneity in T cell infiltration in chordomas. The images show T cell-related markers (CD4=red, CD8=blue), a proliferation marker (Ki-67=yellow) and a tumor marker (keratin=white) as detected by imaging mass cytometry. (A) chordoma with low T cell infiltration (ICR low). (B) Chordoma with high T cell infiltration and a magnified image of the chordoma with high T cell infiltration (ICR high). (C) Heatmap showing the mean cell density Z-score of all T cell phenotypes in all chordoma samples (n=32), hierarchically clustered into a T cell high (n=18) and T cell low (n=14) group. Samples are annotated with their respective ICR score from the transcriptional analysis. In addition, the samples are annotated for their anatomical location, the follow-up status of the patients as well as whether the sample had received neoadjuvant radiotherapy. (D) Boxplots presenting the mean cell densities per mm2 per classified T cell group (T cell high: n=18; T cell low: n=14) for the identified T cell phenotypes. (E) Kaplan-Meier curves interrogating the association between T cell infiltration and disease-specific survival. Two-sided log-rank tests were performed for the Kaplan-Meier estimates. ICR, immunologic constant of rejection; Tregs, regulatory T cells.

Next, chordoma samples were clustered according to the extent of infiltration of the different T cell populations leading to the definition of two major groups (high vs low T cell infiltration; figure 2C,D). Of note, infiltration by CD4+ and CD8+ T cells was highly correlated across the cohort (online supplemental figure 3B). All three chordoma samples previously classified as ICR low clustered in the T cell low group but samples that were previously classified as ICR high were equally distributed between high and low T cell groups (figure 2C), demonstrating the additional value of image-based evaluation of immune cell infiltration in solid tumors. Although it appeared that patients with higher T cell infiltration have an improved disease-specific survival, the observed association was not statistically significant (log-rank p=0.11, figure 2E). There was no association between T cell infiltration and recurrence-free or metastasis-free survival (online supplemental figure 4A,B). Moreover, unsupervized clustering of sequential tumors of six patients revealed that the immune contexture of chordomas can change on disease recurrence or metastasis (online supplemental figure 4C).

Antigen-presenting myeloid cells are abundant in T cell high chordomas and co-localize with T cells

Further exploration of the imaging mass cytometry data revealed that high T cell infiltration was accompanied by the presence of dendritic cells (T cell high vs T cell low, p<0.001, Benjamini-Hochberg corrected FDR=0.0052, figure 3A). In addition, we found that T cell high chordomas were characterized by increased numbers of immune cell aggregates mainly composed of T cells and myeloid cells (p<0.001, FDR=0.0052, figure 3A). Other phenotypes that were significantly associated with T cell infiltration include CD11c+HLA-DR+ macrophages (p=0.014, FDR=0.048), CD163+HLA-DR+ macrophages (p=0.0079, FDR=0.035), granulocytes (p=0.0082, FDR=0.035), podoplanin+ stromal cells (p<0.001, FDR=0.0052) and vessels (p=0.0088, FDR=0.035; online supplemental figure 5A,B). No specific phenotypes were associated with low T cell infiltration. Instead, this group exhibited a significantly lower density of immune cells (p<0.001) and stromal cells (p<0.01), resulting in a higher relative proportion of tumor cells (online supplemental figure 5C). HLA class II expression (HLA-DR) was only detected on immune cell populations and never on tumor cells.

Figure 3

Pro-inflammatory interactions shape the immune microenvironment of T cell high chordomas. (A) Boxplots comparing cell densities between the T cell high (n=18) and T cell low groups (n=14). Dendritic cells (CD11c+/−HLA-DR+CD68) and immune cell aggregates (CD45+ and multiple differentiation markers, eg, CD3+CD68+) are significantly higher in the T cell high group as compared with the T cell low group. Two-sided t-test was performed in R and the p value was adjusted with the Benjamini-Hochberg false discovery rate. *=p<0.05; **=p<0.01; ***p<0.001. (B) An image generated with imaging mass cytometry displaying co-localization of dendritic cells (HLA-DR=red) and T cells (CD3=blue). Tumor cells are only shown in the overview image (top-left image, keratin=white). Other myeloid and T cell-related markers are presented in separate images around the overview image in green (CD4, CD8, CD11c, CD14 and CD163), to show co-localization of different markers. (C) Heatmap displaying the spatial interaction analysis results comparing the microenvironment of T cell high with T cell low chordomas. Only significant interactions are visualized and colored for the subgroup(s) in which they are enriched. Permutation Z-scores below 10 were considered not significant and thus removed for visualization. Interactions of interest are indicated by black arrows. PD-L1, programmed death-ligand 1; PDPN, podoplanin.

A comparative neighborhood analysis on the imaging mass cytometry data was conducted to reveal distinct interaction patterns between both T cell groups independent of their overall abundance (figure 3C). Naturally, most interactions involving T cells were enriched in the T cell high group, where most T cell phenotypes interacted with each other, with the exception of regulatory T cells which were only enriched in the neighborhood of CD39+ stromal cells (figure 3C) The pro-inflammatory nature of the microenvironment of T cell high chordomas was further reinforced by the observed interactions between dendritic cells and CD8+ and Ki-67+ T cells as well as CD38+ plasma B cells and CD4+ and unspecified T cells (figure 3C). Most notably, the co-localization of Ki-67+ T cells with CD4+ T cells and dendritic cells suggests interaction between these antitumor immune cells in a coordinated manner. In contrast, the most dominant interactions in the microenvironment of T cell low chordomas involved macrophage and stromal cell subsets including CD163+HLA-DR+ macrophages, stromal cells, vessels and CD39+ vessels (figure 3C). These observations suggest that the presence of T cells in T cell high chordomas accompanies a coordinated inflammatory response that also involves antigen-presenting myeloid cells such as dendritic cells, whereas the T cell low chordomas appear to be dominated with myeloid and stromal cell interactions.

The extent of T cell infiltration in chordomas is independent of clinical parameters

To further explore the relation between the microenvironment of chordomas and the clinical behavior of the disease, we performed multispectral immunofluorescence on the complete patient cohort (n=76). A T cell-orientated antibody panel was designed to evaluate the expression of CD3, CD8, FoxP3, Ki-67 and pan-cytokeratin. Following cell phenotyping and automated counting, patients were clustered into T cell high and low groups (figure 4A,B). Also in this cohort, a positive correlation was observed between CD4+ and CD8+ T cell infiltration (online supplemental figure 6A). No significant association was found between T cell infiltration and disease-specific survival (log-rank p=0.17, figure 4C). Further, we have investigated a potential effect of specific T cell phenotypes but not association with survival was found (online supplemental table 5). Of note, a higher proportion of KI-67+ cells was found in the CD8+ T cell population when compared with CD4+ T cells (p=0.006, online supplemental figure 6B).

Figure 4

T cell infiltration is independent from clinical parameters. (A) Heatmap displaying the mean cell density Z-score per T cell phenotype as identified with multispectral immunofluorescence in the entire cohort (n=76). The samples are clustered in an unsupervized manner and thereby marked as T cell high (n=34) or T cell low (n=42). The samples are annotated for their respective anatomical location and whether they had received radiotherapy prior. (B) Images displaying different T cell infiltration patterns as generated by multispectral immunofluorescence. The original and magnified images show T cell-related markers (CD3=red, CD8=blue, FOXP3=cyan) and a tumor marker (keratin=yellow) in two chordomas, each annotated for their respective T cell group (T cell high=red rectangles; T cell low=blue rectangles). (C) Kaplan-Meier curves of the disease-specific survival comparing the T cell high and T cell low groups of the entire cohort. Two-sided log-rank tests were performed for the Kaplan-Meier estimates. Tregs, regulatory T cells.

Clinical features, such as anatomical location of the tumor or neoadjuvant radiotherapy, were not associated with T cell infiltration, suggesting that the presence of a coordinated immune response is not dependent on neoadjuvant radiation or the tumor’s location (figure 4A). Further associations between clinical parameters and prognosis were explored through Cox proportional hazard models for the disease-specific, recurrence-free and metastasis-free survival, which are all presented in online supplemental table 6. In short, previously established clinical parameters such as anatomical location of the primary tumor or the surgical resection margin were confirmed as prognostic factors.

HLA class I expression is maintained in chordomas, but strong clonal expansion of T cells is rare

Strong immune selective pressure by T cells can lead to immune evasion by tumor cells, most notoriously through loss of HLA class I expression. Since high expression of HLA class I-related genes on RNA level does not preclude the loss of HLA class I expression due to, for instance, mutations in the B2M gene, we investigated the prevalence of HLA class I defects by immunodetection (figure 5A). In the complete cohort, HLA class I and β−2m membranous expression was observed in the majority of chordomas (73%). However, 22% of chordomas displayed defective HLA class I expression, which was accompanied by either loss or low expression of β−2m. The remaining 5% of the cases showed very weak expression of HLA class I and/or β−2m, which was identified by stronger stromal expression of HLA class I/β−2m but weak expression in the tumor cells. Thus, up to 27% of chordoma samples presented some form of impaired or altered antigen presentation. Interestingly, patients in the T cell low group more often had defective or weak HLA class I expression in comparison to T cell high patients (figure 5B). Positivity for HLA class I was relatively more frequent in T cell high patients (figure 5B). HLA class I expression was also examined specifically in the context of CD8+ T cell infiltration but no association between HLA class I phenotypes and CD8+ T cell infiltration was found (online supplemental figure 6C). Intriguingly, we observed intratumoral γδ T cell infiltration in one T cell high chordoma lacking HLA class I expression (online supplemental figure 7A,B). γδ T cells have recently been proposed as effectors of immunotherapy in HLA class I-defective cancers.31 These results indicate that HLA class I expression and, therefore, antigen presentation is maintained in the majority of chordomas.

Figure 5

Antigen presentation and T cell clonal enrichment are present in chordomas. (A) Human leukocyte antigen (HLA) class I (HCA2 and HC10) and β−2m expression patterns as validated by immunohistochemistry presenting the distinct expression patterns. From top to bottom, these images display positive expression of HLA class I and thus also β−2m, weak expression of HLA class I and β−2m and defective expression of HLA class I accompanied by weak expression of β−2m. (B) Stacked bar plots displaying the percentual distribution of the observed HLA class I (both HCA2 and HC10) expression patterns in the chordoma samples in both T cell groups (T cell high: n=34; T cell low: n=40). (C, D) Boxplots displaying the Shannon diversity (C) and the number of unique clones (D) of the profiled T cell receptors (TCR) from the immune microenvironment of chordomas (n=23). Both the Shannon diversity and the number of unique clones are compared between the T cell high (n=10) and the T cell low (n=13) groups. Two-sided t-tests were performed to compare the T cell high and the T cell low group. *=p<0.05; ns=not significant (p>0.05). (E) Boxplots presenting the Shannon diversity in relation to the observed HLA class I expression patterns in the subset that underwent TCR profiling (n=22).

Presentation of tumor antigens can lead to the specific enrichment of T cell clones with antitumor specificity in the tumor tissues. To investigate the clonality of T cell responses in chordomas, we sequenced the variable region of the TCRB locus in 27 chordoma samples from 23 patients. Intriguingly, T cell high chordomas presented a more focused TCR repertoire than T cell low chordomas, which was demonstrated by the lower number of unique clones (p=0.068) and Shannon diversity (p=0.047) in T cell high chordomas compared with T cell low chordomas (figure 5C,D). Additionally, tumors displaying a Shannon diversity below 2.5 were only contained within the T cell high group, suggesting the presence of strong clonal expansion in two T cell high chordomas (figure 5C). To substantiate the evidence for enriched TCR clonality, the Shannon diversity was compared with the HLA class I expression pattern of these 23 tumors, which revealed that patients with positive HLA class I expression, presented a more enriched TCR repertoire than patients with either defective or weak HLA class I expression (figure 5E).

Discussion

Through a multimodal approach comprising RNA sequencing, imaging mass cytometry, multispectral immunofluorescence and TCR profiling, we studied the immune microenvironment of in total 77 patients with chordoma. We found that chordomas present many hallmarks of antitumor immunity. We observed the majority of chordomas to have an “immune hot” microenvironment based on their immune-related transcriptomic profile, especially compared with other bone sarcomas (eg, chondrosarcoma and osteosarcoma). Spatial analysis of immune cell infiltration further defined groups of patients with chordoma based on the prevalence of T cell phenotypes. Increased T cell infiltration was associated with an apparently coordinated antitumor immune response as indicated by the presence of dendritic cells and other immune cells in immune cell aggregates. In contrast, chordomas with lower T cell infiltration were characterized by a microenvironment dominated by interactions comprising stromal cells and myeloid cells. Neither group was associated with clinical outcome nor with any other clinical parameter.

Recently, several studies have contributed to the improved understanding of the chordoma microenvironment, based on approaches including immunohistochemistry, multispectral immunofluorescence or single cell-RNA sequencing.6 7 29 32 33 These studies identified myeloid cells and T cells as the most predominant immune cell populations in chordomas.7 29 32 33 Similar to our findings, the majority of these immune cells were confined to the stroma, although both cell types could be found in the tumor nodules.6 7 33 Moreover, other lymphoid cells such as natural killer cells (or innate lymphoid cells) and B cells were considered rare in chordomas.7 29 33 Using imaging mass cytometry, we expanded on the previously identified immune cell populations and were able to establish immunologically relevant interactions within the chordoma microenvironment.

Studies focused on methylation and RNA sequencing have identified groups of patients with either an inflamed or non-inflamed phenotype and associated these groups to clinical outcome, all with varying implications.34–36 For instance, Zuccato et al found through methylation profiling that a group of patients presenting an inflamed phenotype had a worse prognosis than patients without these inflammatory features,35 whereas Huo et al observed that patients from the non-inflamed group had a worse clinical outcome.34 Furthermore, Baluszcek et al described two similar groups based on methylation and transcriptomic profiling, without finding an association between these groups and survival.36 All these studies, however, either only include skull base chordomas or skull base chordomas combined with spinal chordomas. The variable composition of the cohorts based on institutional referral patterns complicates the identification of common prognostic markers across studies. Here we included tumors from all major anatomical locations as well as different treatment strategies, providing an extensive characterization of a uniquely large series of chordomas as representative of the general patient population as possible. Presently, the cause of the apparent disconnect between the immune microenvironment and clinical outcomes in chordoma remains uncertain. It is conceivable that the extensive and intricate surgical interventions required for chordoma treatment exert a more significant influence on clinical outcomes than the natural immune response. Future investigations might delve into this by examining the prognostic significance of the immune microenvironment in patients exclusively undergoing radiotherapy. Alternatively, it is plausible that an inadequate cytotoxic T cell response stems from a lack of available CD8+ T cell epitopes due to the low mutation burden of chordomas, offering another potential explanation.

Nevertheless, most studies, including ours, support the separation of distinct groups of chordomas based on the abundance of pro-inflammatory immune cells in their microenvironment. Little is known, however, regarding potential drivers of these different immune contextures. In other tumor types (eg, melanoma, carcinomas) spontaneous immune responses can often be explained by the mutation burden of a cancer where tumors with high mutation burden more often experience spontaneous antitumor immune responses as a consequence of high neoantigen abundance. In contrast, similar to most other sarcomas, the number of somatic, coding mutations in chordomas is generally low and stable between patients making it unlikely that distinct immune responses can be explained by the overall mutation burden of the tumors.13–15 37 Although we only found evidence for strong clonal expansion in two T cell high chordomas, the TCR repertoire was less diverse in T cell high chordoma than T cell low chordomas, which would support that antigen recognition is a potential driver of immunity in chordomas. Naturally, this raises the question whether there is some specificity for the known chordoma-related transcription factor brachyury (TBXT) in these tumors, as it has increasingly been demonstrated that it can be recognized by T cells as well as that it can be targeted by means of vaccination in chordomas.9 38 Another possible underlying factor to these distinct groups could be the extent of immunosuppressive mechanisms in the chordoma microenvironment. For instance, we identified abundant CD39 expression on vessels and stromal cells, in the vicinity of immune cells in both T cell high and T cell low chordomas, suggesting that the adenosine pathway might be a potential mediator of immune suppression in chordomas. Although the immunosuppressive functionalities of this pathway have mainly been described in carcinomas,39 40 expression of CD39 was recently found to associate with increased inflammation in primary tumors as well as with increased T cell densities and PD-1 expression in recurrent undifferentiated soft tissue sarcomas.41 Further research should determine the extent to which the adenosine pathway can modulate immune responses in this disease.

Going forward, it will be crucial to understand how the chordoma microenvironment impacts the outcomes of immunotherapy, in particular immune checkpoint blockade therapy. Intuitively, and as observed in other tumor types, responses should be more likely to occur in patients where spontaneous antitumor immune responses are observed. Encouragingly, our data also supports the notion that most patients with chordoma may be eligible to immunotherapeutic approaches as: (1) despite differences in the extent of immune infiltration between tumors, the majority of chordomas display enhanced hallmarks of antitumor immunity as compared with other sarcoma types and (2) we determined that the majority of chordoma cases retain HLA class I expression and, therefore, antigen presentation. Therefore, other approaches including antigen vaccination or T cell-based immunotherapies are, in theory, applicable for patients with chordoma and may be particularly important in chordomas that do not present with a conspicuous antitumor immune response. We speculate that providing immunotherapy in the neoadjuvant setting42 to patients with chordoma, either with or without radiotherapy, might have a strong impact on the clinical management of these tumors as it might improve surgical outcomes through tumor shrinkage. This could contribute to less extensive and complicated surgical procedures which could either reduce the risk of surgery-related morbidities (eg, by saving critical nerve structures in the axial skeleton) or possibly facilitate complete surgical resection.

In conclusion, through a multimodal approach we provided evidence that chordomas present multiple hallmarks of antitumor immunity. On a transcriptional level, chordomas appear more enriched for inflammatory signaling than other genetically complex sarcomas. Within our chordoma cohort, we identified distinct immune contextures, either defined by a coordinated, possibly antigen-driven, immune response or a microenvironment dominated by stromal and myeloid cell interactions. Furthermore, antigen presentation was maintained in the majority of chordomas. Our study provides a rationale for using immunotherapeutic strategies to aid in the clinical management of chordomas.

Data availability statement

Data are available in a public, open access repository. The RNA sequencing data generated in this study are publicly available on GEO (GSE239531). The imaging mass cytometry data is available on BioImage Archive (S-BIAD830). The RNA sequencing data for the other studied sarcoma subtypes is unpublished and therefore not yet available, with the exception of the osteosarcomas (GSE237033). The used code and processed data will be published on GitLab for reproducibility. All other data are available from the corresponding author on reasonable request.

Ethics statements

Patient consent for publication

Ethics approval

Patient consent: A waiver of consent was obtained from the medical ethical evaluation committee (Medisch-Ethische Toetsingscommissie Leiden Den Haag Delft; protocol number: B17.039 and B21.009). Ethical approval: All specimens in this study were pseudoanonymised and handled according to the ethical guidelines described in ‘Code for Proper Secondary Use of Human Tissue in The Netherlands’ of the Dutch Federation of Medical Scientific Societies. As described above, a waiver of consent was given by the ethical evaluation committee (Medisch-Ethische Toetsingscommissie Leiden Den Haag Delft; protocol number: B17.039 and B21.009).

Acknowledgments

The authors would like to thank the Flow Core Facility and the Information Technology & Digital Innovation Department of the Leiden University Medical Center for their service.

References

Supplementary materials

Footnotes

  • NFCCdM and JVMGB contributed equally.

  • Contributors SvO contributed to the conceptualization, methodology, investigation, data curation, data analysis, visualization, writing—original draft and writing—review editing. DMM contributed to the conceptualization, investigation, data curation, data analysis and writing—review editing. MEI contributed to the experiments, data analysis and writing—review editing. JPR contributed to the data analysis and writing—review editing. BEMWdA contributed to the investigation and writing—review editing. RvdB contributed to the investigation and writing—review editing. IHB-dB contributed to the investigation and writing—review editing. MvdP contributed to the investigation and writing—review editing. PMW-K contributed to the investigation and writing—review editing. SBP contributed to the resources and writing—review editing. WCP contributed to the resources and writing—review editing. RJPvdW contributed to the resources and writing—review editing. NFCCdM contributed to the conceptualization, methodology, investigation, supervision, writing—original draft and writing—review editing. JVMGB is the guarantor and contributed to the conceptualization, methodology, investigation, supervision, writing—original draft and writing—review editing.

  • Funding This work was financially supported by the intramural Leiden Center for Computational Oncology strategic fund. NFCCdM is funded by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (grant agreement no. 852832) and by the VIDI ZonMW (project number: 09150172110092). JVMGB is the recipient of a Cancer Research Institute/Chordoma Foundation CLIP Grant (CRI5106).

  • Competing interests No, there are no competing interests.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.