Imaging examinations and further diagnostic work-up of MFH
Here, we report a rare case of MFH-B in the lower femoral segment of the left thigh in a 44-year-old male patient. The patient complained of pain and discomfort in the distal left thigh and went to a local hospital for symptomatic conservative treatment. The result of treatment was not good in the local hospital. The patient suffered from severe pain and found a hard mass at the distal end of his left thigh for one month. Therefore, he came to our hospital for further diagnosis and treatment. X-ray showed irregular thawing bone destruction at the lower end of the left femur, which invaded the bone marrow cavity. Multiple strip like high-density dead bone shadows were seen inside (Fig. 1A). CT and MRI showed mass of soft tissue at the lower end of the left femur and its surrounding (Supplementary Fig. S1B). HE and Immunohistochemical staining showed SMA (+), CD68 (+), Ki67 (+), CD34 (+), Desmin (-), EMA (-), S-100 (-), Bcl-2(-)and the pathological result was malignant fibrous histiocytoma (Supplementary Fig. S1A) . In order to reduce the burden of tumor and relieve the symptoms of leg pain, tumor resection was performed to remove the tumor. The patient was initially diagnosed with MFH at first resection of the tumor.
Two months later, the results of the patient's MR imaging showed that patchy abnormal signal was seen in the left side of the L5 spinous process and the posterior edge of the vertebral body, which was slightly more advanced than before, and the enhancement was slightly increased. The abnormal signal at the left posterior border of sacrum was similar to the anterior range and slightly reduced in signal. The possibility of metastasis was considered in the combination with medical history (Fig. 1E). After the patient developed local recurrence and metastasis, he was treated with PD-1 tirelizumab in his second round of chemotherapy.
Single-cell RNA-seq identified MFH associated cellular constitution in tumor tissues and adjacent muscle tissues
To discover the alteration of gene expression in MFH, we performed scRNA-seq of a patient referred to the SIR RUN RUN HOSPITAL NANJING MEDICAL UNIWERSITY. Thus, we dissociated tumor tissues into a single-cell suspension and performed scRNA-seq analysis. Adjacent muscle tissues served as controls. The harvested cells from different groups were sequenced on a 10 × Genomics platform (Fig. 1F), we obtained a total of 18433 single cell transcriptomes (10532 adjacent muscle tissues; 7901 tumor tissues) from two groups of organizations. We conducted preliminary quality control and evaluation on the sequencing results, removed reads with low sequencing quality, compared reads with the reference genome using cellranger, annotated reads as specific genes, and then corrected unique molecular indexes (UMI) and counted them (Supplementary Fig. S2, S3E).
Unbiased clustering of the cells identified 19 clusters in parallel, based on t-distributed stochastic neighbor embedding (t-SNE) and uniform manifold approximation and projection (UMAP) analyses according to their gene profiles and canonical markers (Supplementary Fig. S3F). Our initial goal was to visualize and ultimately define the various cell subsets in the dataset. Subsequently, we used t-SNE visualization of the cells to reveal 10 major clusters, including CD4 + T cell (2 cell clusters), CD8 + T cell (4 cell clusters), NKT cell (1 cell cluster), Proliferative T cell (1 cell cluster), Macrophage (3 cell clusters), Osteoclast (1 cell cluster), Fibroblast (2 cell clusters), Proliferative Fibroblast (3 cell clusters), B cell (1 cell cluster), Mast cell (1 cell cluster) (Fig. 1G, H). In particular, we identified the marker genes for each clusters as follows [8, 11, 12]: (1) the CD4 + T cell highly expressing CD4, but low expressing CD8; (2) the CD8 + T cell highly expressing CD8, but low expressing CD4; (3) the NKT cell highly expressing NK cell markers and T cell markers GNLY, GZMB; (4) the Proliferative T cell highly expressing T cell markers and cell proliferating markers TUBA1B, and MKI67; (5) the Macrophage with high expression of the markers C1QC, C1QA; (6) the Osteoclast specifically express the markers CTSK and MMP9; (7) the Fibroblast expressing COL1A1 and FN1; (8) the Proliferative Fibroblast with high expression of proliferating markers and fibroblast markers; (9) the B cell specifically expressing IGHM, CD79A; and (10) the Mast cell highly expressing TPSAB1 and TPSB2 (Fig. 1H). The results were consistent with the similarity analysis of each cell population (Supplementary Fig. S3A). Subsequently, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis with up-expressed genes identified specific processes relevant to them (Supplementary Fig. S4), which were coincident with their cell type.
By comparing adjacent muscle tissues, we obtained the relative abundance of tumor cells and infiltrating immune cells in tumor tissues, which first reveal the landscape of tumors and invasive immune cells in MFH. These results provide deeper insights into the mechanisms of tumor clearance by immune cells in the surrounding microenvironment and contribute to improved targeting and immunotherapy for MFH [13, 14]. Next, we attempted to discern the cellular differences between tumor tissues and adjacent muscle tissues. We noticed that almost all types of cell populations were present in the tumor tissues and adjacent muscle tissues except for part of the Proliferative Fibroblast, some of them were almost only recognized in adjacent muscle tissues, and some of them were predominantly recognized in the tumor tissues (Fig. 1G). The proportion of the cellular clusters varied significantly among the tissues, suggesting the intertumoral heterogeneity as well as the consistency among the tissues (Supplementary Fig. S3B).
GO and KEGG analyzed the differentially expressed genes (DEGs) in tumor tissues and adjacent muscle tissues, and we found that DEGs were enriched in different biological processes, indicating that MFH is related to the regulation of immune system process, cell activation, blood vessel development, cell migration and cell motility (Fig. 1B). We found numerous distinctions in the composition ratio and structure of the cell types, as well as in the expression of genes within the clusters, suggesting that the biological features of tumor tissues differed from those of adjacent muscle tissues (Supplementary Fig. S3C). GO and KEGG analyzed the DEGs of non-immune cells in tumor tissues and adjacent muscle tissues, and we revealed that the DEGs were enriched in Pathway in cancer (Fig. 1C). The scRNA-seq data quantified the expression of genes associated with disease developmental pathways, including TGF-β, MAPK, NF-κB, and JAK-STAT, as well as signaling pathways associated with epithelial-to-mesenchymal transition (EMT), in various MFH cell populations (Supplementary Fig. S5, Fig. 1D). In our tissue sample, we observed upregulation of ligands and target genes of these pathways in Osteoclast, but few changes were found in immune cells. We also detected downregulated expression of MAPK pathway genes in several Proliferative Fibroblast populations. Some genes related to EMT, NF-κB, and JAK-STAT pathways were upregulated in Osteoclast and Fibroblast within MFH lesions. EMT programs have been indicated to play an important role in cancer invasion, metastasis and acquiring drug resistance. The analysis of gene signatures associated with EMT programming exhibited that the procedure consisted of EMT markers was significantly highly expressed in the Osteoclast and Fibroblast , suggesting that most of the Osteoclast and Fibroblast in this sample were undergoing an active EMT process. Interestingly, immune cells showed few EMT-related genes except Macrophage, which expressed abundant EMT-related genes (Supplementary Fig. S5, Fig. 1D).
In consistent with the high program of EMT, the MFH sample also showed the significant high levels of invasion and metastasis as well as the high angiogenesis genes signature score (Supplementary Fig. S3D), indicating that the MFH in this patient might have an increased capability for high-grade metastasis, which highly correlate with poor prognosis . Indeed, this patient just presented with metastasis based on clinical examination about two months after the surgery. MFH exhibits the highest number of tumor-infiltrating lymphocytes (TILs), suggesting that MFH could benefit from immune checkpoint blockade (ICB). However, not all MFH patients respond to neoadjuvant ICB. An outstanding question is therefore to identify which underlying mechanisms and associated markers determine therapeutic response. In fact, TILs represent a heterogeneous population of cells with respect to cell type composition, gene expression and functional properties. So far, TIL scores and tumor PD-L1 expression have been proposed to predict clinical outcome, but their ability to act as predictors for MFH remains unclear [16, 17].
Gene expression heterogeneity of T cell subsets was identified in the MFH
T cells are the key elements of cancer immunotherapy. However, their high heterogeneity regarding their cell-type compositions, gene expression patterns and functional properties significantly influence the outcomes of the T cell-based immunotherapy. Interestingly, we found that T cell clusters were present at high levels in the immune cells (Supplementary Fig. S3E), and the presence of infiltrating T cell in tumor tissues was confirmed using immunohistochemistry (IHC) with CD3E, CD4, CD8 antibodies, which was consistent with the results of the scRNA-seq data (Supplementary Fig. S6A), thus demonstrated that the T cell-based immunotherapy might be efficient in the MFH patients (Fig. 1G). Notably, we observed that the overall number of T cell in the tumor tissues was much lower than that in the adjacent muscle tissues , indicating that T cell infiltration was inefficient (Supplementary Fig. S6D). There are many differences in gene expression within T cell clusters, suggesting that T cell biology in tumor tissues differs from that in adjacent muscle tissues (Supplementary Fig. S6B).
To reveal the intrinsic structure and potential functional subtypes of the overall T cell populations, we performed unsupervised clustering of all T cells via t-SNE and UMAP algorithm, subsequently identified 12 distinct subclusters of T cell. These 12 stable clusters include 8 clusters for CD8 + and 4 clusters for CD4 + cells. Each cluster has its unique signature genes (Fig. 2A, 3A and Supplementary Fig. S8F).
Cells of the first CD8 + cluster, C1_Exhausted CD8 + T, expressed high levels of exhaustion markers LAG3, PDCD1, and HAVCR2. Cells of the second CD8 + cluster, C2_Naïve CD8 + T, specifically expressed ‘‘naïve’’ marker genes such as LEF1, TCF7 and CCR7. The third cluster, C3_Proliferative CD8 + T, characterized by specific expression of MKI67, TOP2A, and PCNA. The fourth cluster, C4_Effector CD8 + T, was identified by the high expression of the CX3CR1, FCGR3A, and FGFBP2, which commonly associated with T cell possessing effector functions. The fifth cluster, C5_Cytotoxic CD8 + T, was characterized with relatively high expression of the cytotoxic molecules IL7R, TNF and CD69. The sixth cluster, C6_Proliferative APOE + CD8 + T, shared the proliferative marker with cluster3 and highly expressing APOE. C8_APOE + CD8 + T, highly expressed APOE [8, 11, 18]. The remaining cells were primarily T helper cells, characterized by the high expression of GZMA and CCL5, and fallen into the seventh cluster C7_CD8 + helper T (Fig. 2E and Supplementary Fig. S7A, B).
We counted the cell number and proportion of each cellular subcluster (Fig. 2C, D and Supplementary Fig. S8B, C). The percentage of C1_Exhausted CD8 + T cell within CD8 + T cell isolated from tumor tissues and adjacent muscle tissues was much higher than other cell types, revealing potential enrichment of exhausted CD8 + T cell in the TME (Fig. 2C, D). The percentage of C1_Exhausted CD8 + T cell was increased significantly in tumor tissues, which was consistent with previous findings. GO enrichment analysis showed that C1_ Exhausted CD8 + T showed a state of loss of function, persisted in the tumor tissues but replied poorly to the tumor cells [16, 17]. These exhausted CD8 + T cells expressed high levels of PDCD1 in TME and can be rescued from unresponsive and depleted state by ICBs (Fig. 2H). At the same time, we observed that the overall number of CD8 + T cells in the tumor tissues was lower than that in the adjacent muscle tissues , indicating that the infiltration efficiency of CD8 + T cells was low and might be associated with poor prognosis (Fig. 2B).
Next, we generated computationally imputed pseudotime trajectories using Partition based graph abstraction (PAGA)  and RNA velocity analysis to infer the course of maturation of CD8 + T cells in MFH . Most cells from each cluster were gathered based on the similar gene expression, and variant subsets formed into a relative process in pseudotime. We observed three distinct CD8 + T cell trajectories that began with the C2_Naïve CD8 + T, which we considered as the root of the trajectories (Fig. 2F and Supplementary Fig. S7C). One of these trajectories formed C3_Proliferative CD8 + T and ended with C1_Exhausted CD8 + T. Consequently, those exhausted CD8 + T cells were highly enriched at the late period of the pseudotime, demonstrating that the CD8 + T cell state transformed from activation to exhaustion, Monocle 3 algorithm confirmed these trajectories (Fig. 2G). We analyzed gene expression patterns involved in CD8 + T cell-state transitions. The genes regarding “positive regulation of immune system process” decreased significantly along the pseudo-time axis, while the genes related to “cellular response to cytokine stimulus” increased significantly. Genes referred to “positive regulation of cell adhesion” increased initially and then decreased along the pseudo-time axis (Fig. 2I).
We identified 8 groups of DEGs along the trajectory of C1_Exhausted CD8 + T cells. Firstly, Naïve T cell markers CCR7, LEF1 and TCF7 were reduced following the trajectory . The subsequent groups increased at the end of the trajectory and were characterized by effector TNFSF9, cytotoxicity GZMK, GZMA, NKG7, and early markers of general exhaustion PDCD1, CTLA4, and TIGIT. In the last two groups, the early-activating genes TGFB1, GZMM and TNF increased midway through the trajectory, but decreased thereafter. In each gene set, we authenticated several genes that are previously unidentified as T cell markers (for example, COTL1, CCDC141 as differentiated exhausted CD8 + T markers) (Fig. 2J).
Cells from the first CD4 + cluster, C1_Cytotoxic CD4 + T, overexpressed cytotoxic molecules IL7R, TNF and CD69. Cells in the second CD4 + cluster, C2_Exhausted CD4 + Treg, with high expression of FOXP3, IL2RA, the canonical exhausted gene signatures such as PDCD1, CTLA4 were also significantly over expressed. The third cluster, C3_Effector CD4 + helper T, characterized by specific expression of NKG7, GZMK, representing effector T cells, and also highly expressed CCL5 and GZMA, describing helper T cells. The fourth cluster, C4_Naïve CD4 + T [8, 11], specifically expressed “naïve” marker genes such as LEF1, TCF7 and CCR7 (Fig. 3B and Supplementary Fig. S8A, S8E).
Importantly, we also found that C2_Exhausted CD4 + Treg and C3_Effector CD4 + helper T positively expressed the inhibitory receptors and ligands (IRs) including TIGIT, CTLA4, PDCD1 and LAG3, IRs were associated with the exhaustion program of dysfunctional tumor-infiltrating lymphocytes (TILs), suggesting that these cells became exhausted after the initial activation in the MFH. Treg showed relatively high levels of immune inhibitory molecules TIGIT, CTLA4, PDCD1 and TNFRSF18, which may contribute to Treg-mediated suppression of anti-tumor immune responses in the MFH. Recently, the anti-TIGIT therapeutics have drawn great attention in treating colorectal cancer, breast cancer, and melanoma through modulating the activities of CD8 + T, Treg, and NK cells. We also noticed that TIGIT was widely expressed in CD8 + T cells corresponding to gradual loss of responsiveness (Supplementary Fig. S6C). In the new era of immunotherapies, ICIs including antibodies against PD-1 (nivolumab) and CTLA4 (ipilimumab) are widely used for cancer treatment. ICIs act by blocking the inhibitory receptors of immune system on T cells (PD-1 and CTLA4), and thereby activate tumor-specific T cells to destroy tumor cells. Recently, Nivolumab in combination with ipilimumab were reported to have survival benefits in patients suffering from hepatic melanoma of unknown primary origin. These messages illustrate that TIGIT, CTLA4 and PDCD1 blockade could be an effective therapy for MFH [16, 17].
According to the tips of the patient's pelvic plain scan on August 20th, 2021, abnormal imaging signal was found around the left iliac vessels and left inguinal region, considering lymph node enlargement. The soft tissue thickening accompanied by abnormal enhancement signal was detected around left L5 lamina. Therefore, the possibility of distant metastasis was considered (Fig. 1E). In mid-August, due to the outbreak in Nanjing and personal economic reasons, the patient underwent a second round of chemotherapy in Nanjing Gaochun people's Hospital by using epirubicin hydrochloride 100mg, ifosfamide 2.0g, tirelizumab (PD-1 antibody) 200mg and anlotinib 10mg. On September 18th, 2021 patient reexamined blood routine, the examination of male tumor marker showed that alpha fetoprotein, carcinoembryonic antigen, carbohydrate antigen 19 − 9 and carbohydrate antigen 12 − 5 returned to normal.
Next, we analyzed C2_Exhausted CD4 + Treg and obtained a list of 54 exhaustion-specific genes by comparing exhausted and non-exhausted CD4 + T cells. Multiple known exhaustion markers, such as HAVCR2, PDCD1, ENTPD1, CTLA4, TIGIT, TNFRSF9 and CD27, were selected. The 54-gene list also contained several little-described genes , such as MYO7A, TOX and CXCL13, as well as novel exhaustion makers such as LAYN, PHLDA1, and SNAP47 (Fig. 3H). Based on the TCGA data, high levels of LAYN associated with poor prognosis (Fig. 3I, J). LAYN, encoding layilin, was recently reported to be highly expressed in Treg isolated from hepatocellular carcinoma. Besides, LAYN is linked to the suppressive function of tumor Treg and exhausted CD8 T cells. Thus, our data not only confirmed previously identified genes associated with exhausted CD4 + T cells and Treg, but also discovered additional markers for these cell types .
The percentage of C1_Cytotoxic CD4 + T cells within CD4 + T cells isolated from tumor tissues and adjacent muscle tissues was much higher than other cell types, displaying potential enrichment of C1_Cytotoxic CD4 + T cells in the tumor microenvironment (Supplementary Fig. S8B, C). GO enrichment analysis showed that C2_Exhausted CD4 + Treg was enriched with functions involving “lymphocyte activation” and “response to cytokine” (Fig. 3G) .
After defining the CD8 + T and CD4 + T cells in our dataset, we recognized the DEGs in tumor tissues and adjacent muscle tissues of MFH. We identified the differential gene set for these CD8 + T and CD4 + T cells that allowed a more in-depth analysis of regulatory pathways (Supplementary Fig. S6B). KEGG enrichment profile analysis demonstrated that the majority of the upregulated genes of the CD8 + T and CD4 + T cells focused on IL-17 signaling pathway , Osteoclast differentiation and TNF signaling pathway (Supplementary Fig. S7D, S8G).
Similarly, we generated computationally imputed pseudotime trajectories using PAGA and RNA velocity analysis [19, 20] to infer the course of maturation of CD4 + T cells in MFH. We observed three distinct CD4 + T cell trajectories that began with the C4_Naïve CD4 + T, which was considered as the root of the trajectories. One of these trajectories developed C1_Cytotoxic CD4 + T and ended with C2_Exhausted CD4 + Treg. Accordingly, those C2_Exhausted CD4 + Treg were highly enriched at the late period of the pseudotime, demonstrating that the CD4 + T cell state transformed from activation to exhaustion (Fig. 3C, D). Monocle 3 algorithm confirmed these trajectories (Fig. 3E), and the profiling of marker genes confirmed their functional annotation (Supplementary Fig. S8D).
We analyzed gene expression patterns involved in CD4 + T cell-state transitions . The genes associated with “cytokine-mediated signaling pathway” and “leukocyte differentiation” decreased significantly along the pseudo-time axis, while the genes related to “protein localization to membrane” and “establishment of protein localization to organelle” increased markedly. Genes referred to “positive regulation of metabolic process” and “positive regulation of biological process” were increased at first and then decreased along the pseudo-time axis (Fig. 3F).
Two Distinct States of Tumor-Enriched Macrophages
We detected a total of 3972 macrophages that formed 2 clusters (Fig. 4A). Although genes upregulated in C2-Macrophage-THBS1 were enriched for signatures of myeloid-derived suppressor cell (MDSC), those in C1-Macrophage-C1QA simultaneously resembled the signatures of tumor-associated macrophage (TAM) and M1 and M2 macrophage (Fig. 4E). The co-existence of M1 and M2 signatures indicated that TAMs were more complex than the classical M1/M2 model, which was consistent with the previous study. A diffusion map of their global transcriptomes showed that C2-Macrophage-THBS1 (MDSC-like macrophage) and C1-Macrophage-C1QA (TAM-like macrophage) formed a continuum but with distinct expression features (Fig. 4G). Specifically, MDSC-like cells highly expressed the S100A family genes FCN1 and VCAN, whereas they possessed low levels of HLA-related genes . In contrast, TAM-like cells expressed a set of genes (APOE, C1QA, C1QB and TREM2), which is found previously in TAMs of lung cancer. Furthermore, two additional genes, SLC40A1, encoding ferroproteins, and GPNMB, encoding type I transmembrane glycoprotein, showed high levels in TAM-like cells. The transcription factors of these two clusters were diverse, that is TAM-like cells preferentially expressed MITF, RUNX2 and MAF and MDSC-like cells harbored high levels of NR4A1, RXRA and TCF25 (Fig. 4E).
The number of TAM and MDSC isolated from tumor tissues were higher than that of adjacent muscle tissues, demonstrating that TAM and MDSC were preferentially recruited in tumor microenvironment (Fig. 4B-D) . Indeed, macrophages were increased in tumor tissues compared with the adjacent muscle tissues (Fig. 4J). The slight increase of macrophages in tumor tissue illuminates that macrophage immunotherapy may be effective in patients with MFH . Enrichment analysis of up-regulated gene subsets showed that the function of TAM was mainly “regulation of cell migration” and “regulation of cell motility”. And the function of MDSC was primarily “inflammatory response” and “regulation of immune system process” (Fig. 4F).
After defining the macrophage in our dataset, we identified the DEGs in tumor tissues and adjacent muscle tissues of MFH. We authenticated the differential gene set for these macrophages, allowing for a more in-depth analysis of regulatory pathways (Fig. 4K). We focused on the two-fold upregulated or two-fold downregulated genes from macrophage in tumor tissues compared with the adjacent muscle tissues in MFH. Heatmap for the profile of CTSK, MMP9, CKB, CCL18 and COL6A2 demonstrated that these genes were upregulated in macrophage of tumor tissues compared with adjacent muscle tissues (Fig. 4I). GO enrichment analysis demonstrated that the majority of the down-regulated genes in the macrophage subsets were related to “inflammatory response” and “regulation of response to external stimulus” (Fig. 4H).
The t-SNE analysis is instrumental in revealing the heterogeneity of distinct macrophage clusters. However, it is possible that the clusters share common differentiation trajectories. Most macrophages were arranged into a major trajectory with two minor bifurcations by pseudotime ordering. Macrophages from different samples widely distributed in the pseudotime space. The macrophages in adjacent muscle tissues occupied the lower part, while macrophages in the tumor tissues were located in the higher position, indicating that the cells in lower part may be the origin of differentiated macrophages. MDSCs were mainly distributed in the left lower branche, while TAM basically occupied the right upper part, which also suggests that the cells in lower part may be the starting point of differentiation of macrophages (Fig. 4L). We speculate that macrophages have the tendency to transform MDSC phenotype to TAM phenotype .
Two Distinct States of Tumor-Enriched Osteoclasts
Osteoclast plays a vital role in osteolysis and tumor growth in tumor tissues. Based on t-SNE algorithm (Fig. 5A), two individual subclusters of osteoclast were identified with distinct levels of myeloid markers such as CD74 and/or mature osteoclastic markers such as CTSK and ACP5 . The subcluster described as C1_progenitor_OC with high levels of myeloid markers CD74 and CD27 and low levels of OC markers CTSK and ACP5. C2_mature_OC expressed high levels of CTSK and ACP5, and low levels of CD74 (Fig. 5G).
Both progenitor osteoclasts and mature osteoclasts isolated from tumor tissues were higher than those from adjacent muscle tissues, suggesting that progenitor osteoclasts and mature osteoclasts may be enriched in tumor microenvironment (Fig. 5B, C). Enrichment analysis showed that the functions of C1_progenitor_OC were mainly “type I interferon signaling pathway” and “response to cytokine”. The functions of C2_mature _OC were mainly “bone resorption” and “regulation of bone resorption” (Fig. 5E).
Osteoclasts are specialized cells derived from the monocyte/macrophage hematopoietic lineage. They develop and adhere to bone matrix and then secrete acid and lytic enzymes to degrade the bone matrix in a specialized extracellular compartment. Increased bone resorption is the result of osteoclast formation induced by tumor cells and osteoclast formation facilitates bone resorption. Bone is a heterogeneous environment that is benefit for the growth of tumor cells. Among the different cell types present in the bone, osteoclasts are crucial players in the so called “vicious cycle”. This phenomenon is triggered by tumor cells and eventually leads to both tumor proliferation as well as bone deregulation, which promotes the development of bone metastasis .
Subsequently, we performed the trajectory analysis of the osteoclasts to infer the osteoclast maturation process in MFH. Thus, those C2_mature_OC were highly enriched at the late period of the pseudotime, demonstrating that the osteoclast state transformed from progenitor to mature (Fig. 5D). Then we analyzed the trajectory of macrophages and osteoclasts. MDSCs mainly occupied the left upper branch, while TAMs were primarily located in the right lower part. Both clusters of osteoclasts were concentrated at the end of the lower branch, suggesting that macrophages had the tendency to differentiate into osteoclasts (Fig. 6B) [27, 28]. RNA velocity analysis, Monocle 3 and PAGA algorithm confirmed these trajectories (Fig. 6A, C, D).
We analyzed the gene expression patterns involved in osteoclasts and macrophages. Genes related to “regulation of programmed cell death” decreased observably along the quasi-time axis. Genes associated with “mitotic cell cycle process” increased dramatically along the pseudo-time axis. Genes related to “myeloid leukocyte mediated immunity” increased initially and decreased subsequently along the quasi-time axis. (Fig. 6E).
After defining the osteoclasts in our dataset, we recognized the DEGs in tumor tissues and adjacent muscle tissues. We identified the differential gene set for these osteoclasts that allowed a more in-depth analysis of regulatory pathways (Fig. 5F, H). Gene Set Enrichment Analysis (GSEA) demonstrated that osteoclasts in MFH negatively regulated humoral immune response (Fig. 5I). KEGG enrichment analysis showed that most of the up-regulated genes of osteoclasts in MFH were related to “PD-L1 expression and PD-1 checkpoint pathway in cancer” and “T cell receptor signaling pathway” (Fig. 5J).
Gene Expression Heterogeneity in Fibroblast Subsets was Identified in the MFH
As the important cell component in the disease lesion, Fibroblasts within adjacent muscle tissues and tumor tissues were compared and various DEGs were identified (Fig. 7I). We found abundant DEGs in tumor and adjacent normal tissues, so we inferred that cancerous fibroblasts exist in fibroblasts . Therefore, fibroblasts were identified by CopyKAT [12, 29] and divided into tumor-associated fibroblasts and normal fibroblasts (Fig. 7G). To determine the intrinsic structure and potential functional subtypes of the entire fibroblast population, we performed unsupervised clustering in these two types of cells to examine their heterogeneity (Fig. 7A and Supplementary Fig. S9D). According to the DEGs, normal fibroblasts were divided into 7 clusters (Supplementary Fig. S9E): C1_normal_Fibroblast, C2_normal_Fibroblast, C3_normal_Fibroblast, C4_normal_Fibroblast, C5_normal_Fibroblast, C6_normal_Fibroblast, C7_normal_Fibroblast. Tumor-associated fibroblasts are divided into 8 clusters (Fig. 7E, F): C1_malignant_Fibroblast, C2_malignant_Fibroblast, C3_malignant_Fibroblast, C4_malignant_Fibroblast, C5_malignant_Fibroblast, C6_malignant_Fibroblast, C7_malignant_Fibroblast and C8_malignant_Fibroblast [12, 15, 30, 31]. Besides, we found that normal fibroblasts are mainly distributed in adjacent muscle tissues, while malignant fibroblasts are primarily distributed in tumor tissues (Fig. 7B-D and Supplementary Fig. S9A-C).
Enrichment analysis of up-regulated gene subsets showed that the functions of C1_malignant_Fibroblast were principally “blood vessel development” and “cell migration”. The functions of C6_malignant_Fibroblast were largely “immune system process”, “leukocyte activation” and “T cell activation” (Fig. 7H). In addition, GO enrichment analyses revealed that within the malignant cells of MFH, fibroblasts were enriched for genes associated with “cell activation” and “inflammatory response” (Fig. 7J).
Cell Communication Networks in MFH
We used CellphoneDB to predict receptor-ligand interactions. First, we calculated the interactions in the cell types from tumor tissues and adjacent muscle tissues separately. We observed that cells from tumor tissues had more potential for interaction than those from adjacent normal tissues, especially in malignant fibroblasts, osteoclasts, macrophages and several kinds of T cells . Interestingly, we found a cellular communication network between tumor cells and immune cells through immune checkpoint ligand-receptor interactions (Fig. 8A-C).
The pro-apoptotic interaction between C2_Naïve CD8 + T cells and mature osteoclasts, osteoclast progenitors or C8_malignant_Fibroblast (TNFSF12_TNFRSF12A, TNFRSF1A_GRN and TNFRSF1B_GRN) was increased. Mature osteoclasts, osteoclast progenitors or C8_malignant_Fibroblast express NECTIN2 and NECTIN3, thereby transmitting inhibitory signals to TIGIT on C2_Naïve CD8 + T cells, respectively. Some new communications have been observed between mature osteoclasts, osteoclast progenitors or C8_malignant_Fibroblast and CD8 + T cells in the malignant tissues Indicating that T cells were recruited in tumor tissues (CCL4L2_VSIR and CCL5_CCR1). C2_Naïve CD8 + T cells showed inhibitory interaction with mature osteoclasts or osteoclast progenitors (HLA-DPA1_GAL and HLA-FAM3c), or with C8_malignant_Fibroblast (HLA-DPB1_NRG1). The pro-apoptotic interactions of C1_Exhausted CD8 + T cells with mature osteoclasts, osteoclast progenitors or C8_malignant_Fibroblast were increased (TNFSF12_TNFRSF12A, TNFRSF1B_GRN and TNFRSF1A_GRN), but T cell-homing communications weakened (CCL4_CCR5, CCL5_CCR5 and CXCR3_CCL20). Mature osteoclasts, osteoclast progenitors and C6_malignant_Fibroblast express NECTIN2 and FAM3C, which transmitted inhibitory signals to TIGIT and PDCD1 on C1_Exhausted CD8 + T cells, respectively. The co-stimulatory interaction between C2_Exhausted CD4 + Treg and mature osteoclasts or osteoclast progenitors (TNFRSF1B_GRN, TNFSF12_TNFRSF12A and MIF_TNFRSF14) increased. NECTIN2, CD80 and CD86, expressed on mature osteoclasts, osteoclast progenitors C5_malignant_Fibroblast or C6_malignant_Fibroblast, transferred the suppressive signals to TIGIT and CTLA4 on C2_Exhausted CD4 + Treg, respectively. It is noteworthy that C2_Exhausted CD4 + Treg possessed relatively high levels of adhesion molecules, including CD2 and ICAM3. The corresponding receptors, including CD58 and aLb2 complex, are widely expressed by mature osteoclasts, osteoclast progenitors, C5_malignant_Fibroblast and C6_malignant_Fibroblast, which can enhance the adhesion and growth of MFH. Some new interactions have been observed between MDSCs and mature osteoclasts or osteoclast progenitors in tumor tissues, suggesting that MDSCs were recruited in tumor tissues (CCL3L1_CCR1, CXCL2_DPP4 and CCL3L1_DPP4). Furthermore, angiogenic signals increased between MDSCs and mature osteoclasts, osteoclast progenitors or C8_malignant_Fibroblast (NRP2_VEGFA, NRP1_VEGFA and NRP1_VEGFB). In addition, angiogenic signals (NRP1_VEGFB, IGF1_IGF1R, NRP1_VEGFA and NRP2_VEGFA) and costimulatory effect (TNFRSF1A_GRN and TNFRSF1B_GRN) increased between TAMs and mature osteoclasts, osteoclast progenitors or C8 malignant fibroblasts (Fig. 8D-G) [7, 15, 32].