Single-cell transcriptome profiling reveals intra-tumoral heterogeneity in human chordomas

Chordoma is a rare and aggressive bone tumor. An accurate investigation of tumor heterogeneity is necessary for the development of effective therapeutic strategies. This study aims to assess the poorly understood tumor heterogeneity of chordomas and identify potential therapeutic targets. Single-cell RNA sequencing was performed to delineate the transcriptomic landscape of chordomas. Six tumor samples of chordomas were obtained, and 33,737 cells passed the quality control test and were analyzed. The main cellular populations identified with specific markers were as follows: chordoma cells (16,052 [47.6%]), fibroblasts (6945 [20.6%]), mononuclear phagocytes (4734 [14.0%]), and T/natural killer (NK) cells (3944 [11.7%]). Downstream analysis of each cell type was performed. Six subclusters of chordomas exhibited properties of an epithelial-like extracellular matrix, stem cells, and immunosuppressive activity. Although few immune checkpoints were detected on cytotoxic immune cells such as T and NK cells, a strong immunosuppressive effect was exerted on the Tregs and M2 macrophages. In addition, the cellular interactions were indicative of enhancement of the TGF-β signaling pathway being the main mechanism for tumor progression, invasion, and immunosuppression. These findings, especially from the analysis of molecular targeted therapy and tumor immune microenvironment, may help in the identification of therapeutic targets in chordomas.


Introduction
Chordoma is a rare, aggressive, and locally invasive bone tumor. It originates from notochordal remnants [1]. Chordomas are treated by bloc excision, which is the standard treatment option. In patients with recurrence or patients who undergo subtotal or piecemeal excision, radiotherapy is recommended [2]. Many patients with chordoma undergo morbid surgeries and eventually die from tumor recurrence; the 5-and 10-year survival rates are 67.6% and 39.9%, respectively [3,4]. Thus, the treatment strategies are limited [5]. Novel and effective treatment methods are needed.
Tumor heterogeneity is depicted by different cellular populations in the complex tumor microenvironment (TME) of chordomas. These include cancer cells, vascular cells, fibroblasts, and immune cells. Intertumoral and intratumoral heterogeneity refers to heterogeneity between patients with the same histological type, and within a single patient, respectively. Tumor heterogeneity is a key challenge in tumor treatment resistance [6]. Clinical studies on molecular targeted therapy (MTT) in chordomas have identified potential therapeutic targets: PDGFR, EGFR, HER2, VEGFR, etc. [7]. However, the therapeutic effect of these MTTs is limited, and the response rate of apatinib and sorafenib is only approximately 3.7%.
Moreover, sorafenib has significant toxicity; 77.8% and 4.8% of patients develop Grades 3 and 4 toxicities, respectively [8]. Immune checkpoint inhibitors are an effective treatment option for cancers such as sarcomas, but their efficacy is limited to mismatch repair-deficient tumors [9]. However, programmed cell death-1 or programmed cell death-ligand 1 (PD-L1) inhibitors have clinical benefits in some patients. Tumor heterogeneity is the main reason for the diverse responses to these drugs [10]. The complex TME explains the high heterogeneity. TMEs comprise cellular populations, such as cancer cells, vascular cells, fibroblasts, and immune cells [11]. The TME of chordomas consists of immune cell populations, including immune checkpoints [12] and vascular endothelial cells [13]. However, it is difficult for bulk gene analyses to provide an accurate reflection of the TME [14]. Therefore, the heterogeneity of chordomas is not well studied, and novel assessment methods are needed.
Single-cell RNA sequencing (scRNA-seq) may elucidate the mechanisms underlying carcinogenesis and the molecular features of cancers. It may provide a clearer description of the TME at the cellular level, allowing more effective therapeutic strategies for tumors [15]. scRNA-seq is used to delineate the complex tumor heterogeneity in several malignancies such as glioblastomas [16], liver cancers [17], lung cancers [18], renal cell carcinoma [19], and head and neck cancers [20]. However, the knowledge on tumor heterogeneity and cellular interactions in chordomas is limited. We used an scRNA-seq sequencing platform (10 × Genomics) to describe the tumor landscape of chordomas. This study aimed to assess tumor heterogeneity and identify potential therapeutic targets in chordomas.

Chordoma samples
Six patients who underwent surgery at Xuanwu Hospital, Capital Medical University, Beijing, China, from July 2019 to April 2020 were included in this study. None of the patients received preoperative radiation therapy, chemotherapy, or other targeted therapy. Using the 2013 World Health Organization classification of bone tumors, all the patients were diagnosed with classic chordomas [21] (Fig. S1A). On immunohistochemistry [22], all chordoma cells had strongly positive TBXT (brachyury) expression (Fig. S1B). The samples were obtained from fresh tumors intraoperatively. Written informed consent was obtained from all the patients. The study was approved by the Ethics Committee of the Xuanwu hospital of Capital Medical University (No.2016033).

sc-RNA sequencing
Single-cell transcriptomic sequencing was performed (Capitalbio Technology Corporation, http:// www. capit albio tech. com). Following the manufacturer's instructions for Single Cell 30 Library and Gel Bead Kit V2 (10 × Genomics), cell suspensions were loaded on a Chromium Single Cell Controller (10 × Genomics, San Francisco, CA) to generate single-cell gel beads in emulsion. After Drop-seq droplet collection, cDNA amplification and sequencing library preparations were carried out as described previously [17], and the libraries were sequenced on Illumina HiSeq X Ten. The libraries from one batch of droplets were sequenced separately for Drop-seq data from chordoma cells.

Cell filtration, clustering, and downstream analysis
Downstream analysis was done using the Seurat package (version 3.2.0) [23] in R software (version 4.0.2), which helped in preprocessing to clustering, dimension reduction, visualization, and differential gene expression. In processing, raw gene expression matrices were imported and filtered using the following standards: (1) cells expressing at least 200 genes and genes expressed by at least 3 cells and (2) cells with a mitochondrial percentage below 25% (Fig. S2A, B). A total of 33,737 high-quality cells were included in the downstream analysis. Furthermore, the Seurat package was used to normalize and integrate expression data and remove the batch effect. The new integrated matrix was used to scale, run the principal component analysis, and visualize the landscape with UMAP.

Cell type annotation
We adopted 19 clusters and identified all cluster markers using the Seurat package. Cell types were annotated to known biological cell types using canonical marker genes. The SingleR package [24] (version 1.3.6), CellMarker website [25], and inferCNV package (version 1.5.0) were used. The TBXT gene, a specific feature of chordomas, was also used to annotate the tumor cells.

Functional enrichment analysis
After the annotation of each cell type, a functional enrichment analysis was performed to identify differentially expressing genes between different clusters for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). This analysis illustrated the biological processes and potential functions of different cells using the cluster-Profiler package (version 3.17.0) [26] and the org.Hs.eg.db 1 3 package (version 3.11.4). The p-value cutoffs of GO and KEGG were both 0.05. The top 10 terms of the results were visualized in a bar or dot plot.

Pseudotime
Trajectory analysis was performed using the monocle package (version 2.17.0) [27]. Tumor cells, mononuclear phagocytes, and fibroblasts were included. The following parameters of each group were analyzed: lowerDetection-Limit = 0.5, min_expr = 0.1, and num_cells_expressed ≥ 10. For visualization, the plot_cell_trajectory function was used to plot the potential trajectory according to pseudotime, Seurat clusters, and the data.

Cellular communications
The CellChat package (version 0.0.1) was used to analyze cell-to-cell interactions [28]. Most ligand-receptor interactions were mainly identified using the KEGG signaling pathway database and recent peer-reviewed experimental studies. The main steps of inference of intercellular communications included: (1) identification of differentially expressed signaling genes, (2) calculation of an ensemble average, and (3) calculation of intercellular communication probability.

Overall landscape of the chordoma tumor samples
The raw sequencing data could be obtained at the Genome Sequence Archive for Human (GSA-Human), and the access ID is HRA000513 (https:// ngdc. cncb. ac. cn/ gsa-human/s/ GokpM dli). Six tumor samples (from four male and two female patients) were included in the scRNA-seq profiles (Fig. 1A). The patients' age range was 17-49 years. The tumor sites included the following: the sacrum (n = 1), mobile spine (n = 3), and skull base (n = 2). Details of the characteristics of six chordomas are shown in Supplemental Table 1. After strict filtering, 33,737 cells were included in the final analysis. After normalization of gene expression and principal component analysis, these cells were divided into 18 clusters using the UMAP method (Fig. 1B). The 18 clusters were equally distributed among the six samples, indicating a smaller batch-corrected effect (Fig. S3A) 1B-D). The marker genes for the seven cellular populations are shown in Fig. 1B, and the characteristic genes are shown in Supplemental Table 2. Downstream analyses of the B cells, neutrophils, and endothelial cells were not performed due to their small numbers. The variation in the cellular populations was investigated between patients ( Fig. 1E) and within patients (Fig. S3) and indicated a potential clinical significance, especially in tumor location. Primary tumor location was a strong prognostic factor and was related to the extent of surgical resection [29].

Transcriptomic tumor heterogeneity of tumor cells in chordoma
A total of 16,052 tumor cells were identified according to the chordoma-specific marker, TBXT. These chordomas cells were further assigned to six clusters with differently expressed genes (DEGs) ( Fig. 2A and B). There was a high degree of within-patient ( Fig. 2A) and between-patient ( Fig.  S4A) tumor heterogeneity; this means that we found six different clusters with different phenotypes in the overall study population. Meanwhile, each patient had a tumor with specific cellular characteristics. These six tumor subclusters varied with tumor location (Fig. S4B). The GO analysis investigated the biological function of DEGs, which mainly influence immune responses, such as leukocyte migration, regulation, chemotaxis, regulation of B-cell activation, and humoral immune response (Fig. 2C). These results suggest that the tumor heterogeneity explained the varied immunologic tumor microenvironment. Next, the expression of major histocompatibility complex class I genes and PD-L1/ PD-L2 expression on chordoma cells (Fig. S4C) were analyzed. However, most tumor cells did not express major histocompatibility complex genes or PD-L1-these findings discourage the use of programmed cell death-1 inhibitors for chordomas.
There were 6950 chordoma cells in subcluster 0; these cells had the largest tumor population (43.3%). They were mainly involved in the extracellular matrix (ECM) organization by GO analysis (Fig. S4D), which was characterized by CTGF, COL2A1, COL6A3, and COL5A2 (Table S3). They were considered matrix tumor cells (mTCs, mTCs-C0-CTGF). Cluster 1 (4408 cells [27.5%]) displayed a high expression of genes related to metabolism of metal ions: MT1X, MT2A, and MT1G. The enrichment analysis revealed metal ion-related pathways (Fig. S4D). Cluster 1 cells were considered epithelial tumor cells (eTCs-C1-IGFBP3). Cluster 4 was involved in nuclear division (Fig. S4D), i.e., active tumor cell proliferation. The biomarkers for cancer stem cells (CSCs) were expressed in this cluster. Biomarkers such as STMN1 [30,31], UBE2C [32], and PTTG1 [33] [34] are associated with tumor cell survival, clonality, and tumorigenicity in several malignancies. This subcluster was considered to consist of stem tumor cells (sTCs, mTCs-C4-STMN1). The pseudotime analysis showed that subcluster 4 was in the early tumor stage. Clusters 2, 3, and 5 did not express CDKN2A, which is absent in several chordomas [35] (Fig. S4E). Next, we conducted a pseudotime analysis of six clusters to investigate their developmental trajectories: clusters 2, 3, and 5 were in a more advanced stage (Fig. 2D, Fig. S4F). A close relationship between CDKN2A loss and tumor immune microenvironment changes was seen in a prior study [36]. Subclusters 2 and 5, which had CDKN2A loss, were involved in immune functions. As mentioned above, chordomas showed great tumor heterogeneity, which indicates that the pathogenesis and phenotype might be varied. Regarding the clinical treatment, different patients with a risk of recurrence responded differently to the same targeted drug for chordomas [37]. Therefore, MTT plays a crucial role in the personalized treatment of chordoma. The MTT target expression in chordomas was analyzed. This included EGFR, HER2, PDGFR, VEGF, VEGFR2, and stem cell factor receptor (KIT) (Fig. 2E). VEGF and VEGFR2 were the most upregulated, indicating their potential as therapeutic targets.

T/NK cell clustering and subtype analysis
Unsupervised clustering was performed for 3944 T/NK cells; they were divided into five subclusters with their unique signature genes (Fig. 3A, B). Cluster 2 was composed of NK cells (374 cells, 9.5%) marked with KRF1, and others were T cells (3570 cells [90.5%]) with CD3E expression (Fig. 3C). High tumor heterogeneity was associated with T/NK cells among patients ( Figure S5a) and tumor location (Fig. S5b) ) cells, respectively. The specific markers for the immune cells identified above were CD4, CD8, and foxp3, respectively (Fig. 3C). A cluster of double-negative T cells was found (Fig.S5C), which was involved in systemic inflammation and tumor damage [38]. Only 160 double-negative T cells were detected, and these were not included in the downstream analysis. The C1-CD8 + T and C2-NK + showed a high level of cytotoxic activity, which was characterized by GZMA, GZMK, GNLY, PRF1, and IFNγ (Fig. 3B, Table S4). Moreover, few immune checkpoints were expressed on these cytotoxic immune cells (Fig. 3D). There was overexpression of CTLA4, TIGIT, and TIM3 on Tregs (Fig. 3D), which are the main immunosuppressive mechanisms and therapeutic targets for chordomas.

Mononuclear phagocyte clustering and subtype analysis
To investigate the heterogeneity among macrophages, 3698 tumor-associated macrophages (TAMs) were clustered into six subgroups (Fig. 4A) and C1 were the most populated, and they were associated with M2 subtype polarization because they exhibited high levels of CD68 + CD163 + CD204 + expression (Fig. 4B). In addition, VEGFA was heavily expressed in these cells in cluster 1 (Table S3). M2-subtype TAMs secrete VEGFA, which promotes tumor growth [39]. In the GO analysis, this cluster was involved in the regulation of angiogenesis (Fig. 4C). Subcluster 1 was composed of inflammatory chemokines such as CCL3, CXCL2, and CXCL3 (Table S5), consistent with the GO analysis findings (Fig. 4C). The release of these proteins by TAMs may promote inflammatory response and tumor development. Clusters 0 and 1 represented an M2-like TAM cluster. Moreover, the genes S100A8, S100A9, and S100A12, which encode calciumbinding proteins, were expressed by the cells in cluster 2 (Table S3) [40]. Activated mononuclear cells can release these proteins, which promote inflammatory responses in vivo. Thus, these data indicate that C3 may have a proinflammatory and anti-tumor role in chordomas. C4 and C5 indicate nuclear division and neutrophil activation, respectively (Fig. 4C). CD47 and TIM3 were expressed on TAMs (Fig. 4D); however, other immune checkpoints were not (data not shown).

Fig. 3 Clustering and subtype analysis of T and NK cells clustering and subtype analysis, A Clustering for T cells and NK cells, B
Heatmap of the top 10 DEGs of subclusters of T/NK cells, C Expres-sion of known marker genes of T and NK cells. D Expression of the immune checkpoints, NK, natural killer; DEGs, differently expressed genes transcriptomic signatures (Fig. 5B, Table S6). There were 2600 cells in subcluster 0, which had the largest population of fibroblasts (37.4%). The subclusters were characterized by ECM signatures, including collagen molecules (COL1A1, COL3A1, COL4A1), POSTN, and LUM. Furthermore, the GO analysis for this subtype was related to ECM and collagen fibril organization (Fig. 5D). They were considered as matrix cellular-associated fibroblasts (mCAFs, mCAFs-C0-POSTN; Fig. 5B, C). Subcluster 1 consisted of 1185 fibroblasts and was the most populated (17.1%). The significant genes that were upregulated in this subtype were IL1RL1, CCL3, CCL4l2, CCL8, and others (Table S6). The GO analysis confirmed that this subtype was involved in cell chemotaxis and T-cell activation; they were considered inflammatory CAFs (iCAFs, iCAFs-C1-IL1RL1; Fig. 5B, C). The C2 subcluster consisted of 910 fibroblasts and was marked with CRABP2, which was associated with osteogenic differentiation. In the GO analysis, C2 was involved in ossification and cartilage development; therefore, they were considered bone development regulation CAFs (bCAFs, bCAFs-c2-IL1RL1; Fig. 5B, C). Subcluster 3 was made of 796 fibroblasts. The marked genes were CST1 and ANGPTL4. There was a relationship between the GO analysis of this C3 and regulation of vasculature development. The C3 was then regarded as vascular CAFs (vCAFs, vCAFs-c3-MCAM; Fig. 5B, C). However, C5 and C3 had similar functions because the marker genes of C5 were characterized by microvasculature signature genes (MCAM, MYH11, GJA4, and RGS5). Subcluster 5 fibroblasts expressed mainly epithelium-specific marker genes such as KRT19 and KRT8; they were considered epithelial-to-mesenchymal transition (EMT)-like CAFs (eCAFs, eCAFs-c5-KRT19; Fig. 5B, C).

The TGFβ signaling pathway was enriched by the interplay among the CD4 + T cells, fibroblasts, and macrophages
We first analyzed the overall cellular interactions based on ligand-receptor pairs. The most active pathway in the TME of chordomas was the TGFβ signaling pathway (Fig. 6A,  Fig. S6). TGFB1-TGFBR1/TGFBR2 and TGFB1-ACVR1/ TGFBR1 contributed most to this communication network (Fig. 6D), which was consistent with the results of previous studies in which TGFB1 was a potent TGFβ ligand in cancers [41]. Moreover, TGFβ signaling played a key role in tumor progression by its different effects on multiple cell types, including malignant and non-cancerous cells within the TME [42]. This showed that TGFβ was mainly produced and released by fibroblasts and macrophages in chordomas (Fig. 6B). However, this TGFβ signaling pathway network was very active among macrophages, fibroblasts, tumor cells, and CD4 + T cells (Fig. 6C). Although few Tregs were detected, they expressed some immune checkpoints that exerted strong immunosuppressive activity (Fig. 3D). The TGFβ signaling pathway was one main reason for transforming CD4 + T cells into Tregs in the TME of chordomas. Cytotoxic immune cells, such as CD8 + T cells and NK cells, did not express immune checkpoints, resulting in no effect on TGFβ (Fig. 3D) [43,44]. The chordoma cells exhibited properties of the ECM, epithelial cells, and stem cells (Fig. S4D). The enhanced TGFβ signaling was associated with malignant biological processes of ECM, tumor stem cell properties, and EMT [45][46][47]. Therefore, the malignant cells were mainly influenced by the TGFβ pathway in chordomas.

Discussion
In this study, we first delineated the transcriptomic landscape of chordomas by scRNA-seq. The cellular populations of the TME in chordomas included tumor cells, fibroblasts, and immune cells; this varied within and between tumors, indicating high tumor heterogeneity. Moreover, tumor heterogeneity varied with tumor location, suggesting novel mechanisms for different clinical outcomes based on the chordoma site. Molecular therapy targets and immunotherapy biomarkers were also analyzed, and potentially valuable therapeutic targets were shown. The cellular interactions enhanced the TGFβ signaling pathway: this was the main mechanism for tumor progression and invasion and immunosuppression in chordomas. This study is the first to investigate tumor heterogeneity at the single-cell level, and it provides new insights into mechanisms of tumor development and possible therapeutic strategies for chordomas. Chordoma cells accounted for approximately half of the total cellular populations and were made up of six subclusters. These subclusters conferred malignant cell properties such as ECM, epithelial-like, stem cell, and immunosuppression. There were 4408 (27.5% of total tumor cells) chordoma cells in subcluster 1, which were related to epithelial cells. Previous studies revealed a dual epithelial-mesenchymal differentiation of chordomas [48], and these subclusters were likely to resist chemotherapy and radiotherapy. Biomarkers (CD133 + , CD15 + , WNT5, ABTG2, and MYCBP) of CSCs are found in chordomas [49], but these markers were absent in the tumors in the present study. Instead, the biomarkers of CSCs, such as STMN1, UB2C, and PTTG1 in hepatocellular carcinoma, breast cancer, and prostate tumors were detected [30][31][32][33][34].
The pathological mechanisms of these markers should be investigated in future studies. In addition, VEGF and VEGFR-2 were the most upregulated in chordomas. These results are consistent with the clinical evidence that chordomas respond objectively to sorafenib and apatinib [7]. Our results showed a complex immune TME of chordomas, consisting of macrophages, T cells, and NK cells. The CD8 + T cells and NK cells showed a high cytotoxic activity, characterized by GZMA, GZMK, GNLY, PRF1, and IFNγ. Two cytotoxic immune cells were devoid of immune checkpoints. However, a strong immunosuppressive activity was observed, which was mainly exerted by Tregs and macrophages. CTLA4, TIGIT, and TIM3 were highly expressed by Tregs in our study. Tregs have been detected in human chordoma samples and associated with poor clinical outcomes due to their strong immunosuppressive activity [50,51]. TAMs mainly exhibited the M2-subtype, which is involved in angiogenesis and immunosuppression. As CD47 and TIM3 are highly expressed by TAMs, immunotherapy may suppress the pro-tumor development role of TAMs in chordomas.
Cellular interactions showed that the TGFβ signaling pathway played a key role in tumor development, tumor invasion, immunosuppression, and CSCs. According to previous reports, fibroblasts and macrophages are the main sources of TGFβ in chordomas [52][53][54]. Chordoma cells are strongly influenced by TGFβ in several malignant biological factors such as ECM, CSC properties, and epithelial-like characteristics [55,56]. Furthermore, in the present study, the TGFβ signaling pathway also negatively affected the immune cells, which were CD4 + T cells. The CD4 + T cells were the largest   immune cell population in chordomas. CD4 + T cells may be transformed into Tregs or Th0 cells under the influence of certain cytokines or chemokines [57]. CD4 + T cells were affected by the TGFβ signaling pathway and developed into Tregs in this study. Together, these results indicate that the TGFβ signaling pathway is a fundamental mechanism for tumor development and immunosuppression in chordomas. This study had some limitations. First, the results and conclusions were based on sequencing data and were not validated by experiments and clinical trials. However, the primary study objective was to provide a landscape of chordomas, and our results showed this tumor's heterogeneity. Second, the TGFβ signaling pathway was found to play a critical role in chordoma progression by promoting tumor invasion, immunosuppression, and CSCs. There are no existing reports on TGFβ signaling pathway in chordomas. Therefore, the results of the TGFβ signaling pathway in chordomas may not be conclusive and need further research.
Our findings provide a large transcriptomic landscape and details of the single-cell resolution of chordomas. This study is an established resource for elucidating chordoma diversity. However, details of the mechanisms and efficacy of therapeutic options still need to be further explored in experimental and clinical studies.