Fibroblasts identified in the tumor microenvironment

We acquired gene expression profiles of three MESO samples (ADU-S100_10um, ADU-S100_50um, and control) at the single cell level from the Gene Expression Omnibus (GEO) database. A total of 14,347 single-nucleus RNA sequencing data were obtained after stringent filtration and control for subsequent analysis. The detailed procedure of the study was visualized in Fig. 1A and Additional file 1: Fig. S1. To identify the key stromal cells that may have regulatory functions in tumor progression, we utilized dimension reduction analysis and discovered 15 seurat clusters, which were categorized into seven cell types (B cell, endothelial cell, epithelial cell, fibroblast, mast cell, myeloid cell, and NK/T cell) in accordance with cell markers, in the tumor microenvironment (Fig. 1B, C). The expression conditions of cell type marker genes, which assisted in the classification of cell types, were shown in Fig. 1E. What’s more, the average number and cell proportion of the seven cell types in the three samples were depicted in Fig. 1D, illustrating the highest number of the epithelial cell, followed by the NK/T cell. In the differential expression analysis of the 15 seurat clusters, we found that some clusters that were recognized as the same cell types showed a similar differential expression landscape, exhibiting the rationality of the cell type identification (Fig. 1F). We further explored the differential expression conditions of the seven cell types, visualizing the expression of four classic cell markers (CD3D, LYZ, CD79A, and VWF) and listing the top five differentially expressed genes in each cell type, respectively, in Additional file 1: Fig. S2A, B. As CD3D, LYZ, CD79A, and VWF were highly expressed in NK/T cells, myeloid cells, endothelial cells and B cells, respectively, their expression was low in the other cells, illustrating the accuracy of the cell type identification. Subsequently, to figure out the relationship and interactions between fibroblasts and other cells, a cell communication network and ligand-receptor interactions among seven cell types were established. As shown in Additional file 1: Fig. S2C, D, fibroblast lay at the heart of the communication network and possessed the most frequent connections with the other six cell types, emphasizing its momentous role in regulating the biological behavior of MESO. Moreover, cell cycle score and cell cycle distribution based on the tSNE showed fibroblasts were uniformly distributed across the G1, S, and G2M phases of the cell cycle (Additional file 1: Fig. S2E, F), partly reflecting the heterogeneity among fibroblasts.

Fig. 1

The gene expression landscape of cells in MESO microenvironment. A The workflow of bioinformation analysis in the study. B The tSNE plot of single cell sequencing data from three samples (ADU-S100_10um, ADU-S100_50um, and control). C Identification of 15 cell clusters and seven cell types (B cell, endothelial cell, epithelial cell, fibroblast, mast cell, myeloid cell, and NK/T cell) in the tSNE analysis with cell markers. D The histogram of average number and proportion of each cell type in the three samples. E Respective cell markers of each cell type. F Gene differential expression pattern of the 15 cell clusters

Among the three MESO samples, 1969 fibroblasts were identified (Fig. 2A), from which we derived six fibroblast subtypes through dimension reduction analysis and clustering (Fig. 2B). Figure 2C displayed the expression level of four cell markers in each subtype. In addition to the archetypal fibroblast markers, DCN and COL1A2, the expression of APOE was generally high. We also listed the top four or five differentially expressed genes in each subtype in the form of a heatmap (Fig. 2D). The cell communication network and ligand-receptor interactions among six subtypes showed the central position of KRT19_FIB and MMP2_FIB among all the subtypes and the tight link between EFGR and both CALM2 and ANXA1 (Fig. 2E, F). In order to further explore the functions of each subtype, we performed GSVA based on variation scores of gene sets in each fibroblast, utilizing the Hallmark pathway database, revealing that all the subtypes were prevalently enriched in hallmark KRAS signaling dn, and both KRT19_FIB and MMP2_FIB had relatively high enrichment of immune-related hallmark gene sets (hallmark interferon gamma response and hallmark interferon alpha response) and proliferation-related hallmark gene sets (hallmark MYC targets V1 and hallmark MYC targets V2) [25], which implied that these two subtypes might have higher proliferative potential, compared to other subtypes (Fig. 2G). Therefore, among all six subtypes, KRT19_FIB and MMP2_FIB deserved more attention.

Fig. 2

Identification and features’ capture of fibroblast subtypes in MESO environment. A The tSNE plot of fibroblasts which were identified in the MESO environment in the last step. B Recognition of seven fibroblast clusters and six fibroblast subtypes (IGF_FIB, KRT19_FIB, LARP6_FIB, MMP2_FIB, MYL9_FIB and PTPRC_FIB) in the tSNE analysis. C The expression levels of four classic cell markers (DCN, COL1A2, PTPRC and APOE) and proportion of the three samples in each fibroblast subtype. D The top four or five differentially expressed genes in each fibroblast subtype. E The cell communication network of fibroblast subtypes. F The prediction of the interactions between each fibroblast subtype based on the ligand-receptor interactions. G The heatmap of pathway enrichment conditions of the six fibroblast subtypes based on GSEA analysis

The differentiation states of fibroblast subtypes and differentiation-related fibroblast differential genes (FDGs)

Monocle 2 and tSNE were utilized in the process. With 4 tSNE cell clusters recognized (Additional file 1: Fig. S3A), the differentiation trajectory of the fibroblast was visualized. In Fig. 3A, which illustrated the chronological order of the trajectory, the gradation of color represented the time cells took to differentiate. A total of three differentiation states were captured, of which state 1 was the earliest stage and was then separated into states 2 and 3 (Fig. 3B). The distribution of the six subtypes in the trajectory was presented in Fig. 3C, D. While IGF2_FIB, LARP6_FIB, and MYL9_FIB were mainly located in state 3, which implied relatively mature differentiation, KRT19_FIB and PTPRC_FIB were mostly distributed in state 1, a hypodifferentiated stage. The distribution of MMP2_FIB was dispersive, suggesting multiple differentiation states. Consistent with the distribution of KRT19_FIB in state 1, the expression of the top five differentially expressed genes of KRT19_FIB, including C19orf33, ITLN1, KRT19, SLPI, and UPK3B, was extremely concentrated in state 1, especially at the origin of the trajectory (Additional file 1: Fig. S3B). What’s more, the expression of these five genes in the origin was primarily derived from KRT19_FIB (Additional file 1: Fig. S3C), and their expression conditions in the three MESO samples were exhibited in Additional file 1: Fig. S3D. Therefore, the differentiation degree of fibroblasts in KRT19_FIB might be low.

Fig. 3

The pseudotime analysis of fibroblasts. A The differentiation trajectory of fibroblasts, in which the variation in shade of color standed for the passage of time. B The three differentiation states (state 1, state 2 and state 3) distinguished from the differentiation trajectory. C The overview map of the six fibroblast subtypes’ distribution in the differentiation trajectory. D The individual distribution condition of the six fibroblast subtypes in the differentiation trajectory, labeled with different colors

Subsequently, we aimed to catch the differentiation-related fibroblast differential genes (FDGs) among the six subtypes. With the Monocle2 single-cell analysis toolset, 644 genes that were differentially expressed in the three differentiation states (p < 0.001) were extracted and shaped into the monocle2_sig gene set. On the basis of the bulk RNA-seq data and the clinical information of MESO patients recorded in the TCGA database, survival analysis was carried out using the Kaplan–Meier technique and the univariate Cox regression model to detect the survival-related genes and eventually form the KM_sig gene set (81 genes with p < 0.001) and the uniCox_sig gene set (150 genes with p < 0.001). Overlapping the three gene sets, a total of 39 genes in the intersection of the three gene sets were screened out and defined as fibroblast differentiation-related genes (FDGs) (Additional file 1: Table S1). According to the survival curves of the 39 FDGs (Additional file 1: Fig. S4), high expression of ADH1B, CFB, PRG4, PLAAT4, HP, and IFIT3 was favorable to survival, while the other 33 genes were considered risk factors. What’s more, five of the six favorable factors (CFB, PRG4, PLAAT4, HP, and IFIT3) were represented in state 3, as all 16 FDGs correlated with state 1 were marked as risk factors (Fig. 4A). Consequently, we deduced that in the MESO tumor microenvironment, fibroblasts with high degrees of differentiation might help to improve prognosis, whereas fibroblasts with low degrees of differentiation might have the opposite effect. Based on this inference, KRT19_FIB might be associated with an adverse prognosis.

Fig. 4

Construction of fibroblast differentiation-based classification of MESO patients in TCGA cohort. A The correlation network of 39 FDGs. The colors of the left half of 39 circles which correspond to the 39 FDGs, figured out the differentiation states of the genes, while the colors of the other side described the genetic effects on survival (Red: high expression was harmful for survival; Blue: high expression was favorable for survival). The size of circles was inversely proportional to the p value of Cox regression analysis. The links between circles displayed the relationship between genes where red referred to the positive correlation and blue indicated the negative correlation with remarkable statistical significance (p < 0.0001). B The Consensus matrix of k = 3. C The CDF curve and the delta area plot. D The cluster heatmap of FDGs expression in the patients’ MESO samples, based on bulk RNA-seq data from TCGA database, identifying three clusters. The clinical information of the patients in the cohort, including survival outcome, age, gender, race, TNM stage and metastatic condition were listed as well. The color in the heatmap from blue to red demonstrated the progression from low expression to high expression. E The comparison of survival conditions of the three clusters with the Kaplan–Meier curve. F The Kaplan–Meier survival analysis of the MESO patients in the high PCA score group and the low PCA score group. G The boxplot which displayed the distribution of the PCA scores of patients in the three clusters

Construction of a fibroblast differentiation-based classification (FDBC) of MESO

Depending on pseudotime analysis to identify the dynamics of gene expression over the course of consecutive cell development, we tried to seek out the different differentiation states of fibroblasts and differentiation characteristics of each subtype.

In order to explore the clinical value of FDGs, we obtained the bulk RNA-seq data and clinical information of MESO patients from the TCGA database and divided these patients into different classifications on the basis of their expression of FDGs, utilizing consensus clustering. Combining the heatmaps of the consensus clustering matrix, the CDF curve, and the delta area plot, three was considered the optimal k value, producing three clusters (Fig. 4B, C). As shown in Fig. 4D, cluster 1 highly expressed the genes in states 1 and 2, while cluster 2 had relatively high expression of state 3 genes (covering CFB, PRG4, PLAAT4, HP, and IFIT3), which suggested that cluster 1 might be associated with lowly differentiated fibroblasts and the fibroblasts in cluster 2 trended to be highly differentiated. The differential expression of all the FDGs was statistically significant (Additional file 1: Fig. S5C). Therefore, we designated cluster 1 “Lowly Differentiated Fibroblast-Related Mesothelioma (LDFM),” cluster 2 “Highly Differentiated Fibroblast-Related Mesothelioma (HDFM),” and cluster 3 “Moderately Differentiated Fibroblast-Related Mesothelioma (MDFM)”. Among the three clusters, the survival result of HDFM was remarkably better than that of LDFM and MDFM, while that of LDFM was the worst, embodying the prognostic value of the novel classification (Fig. 4E). We further calculated the PCA score of each patient based on 39 FDGs and divided patients into two PCA groups (low-score group and high-score group) (Additional file 1: Fig. S5A, B). Since the survival probability of patients in the high-score group was commonly higher than that of those in the low-score group, the PCA score probably had a positive correlation with survival (Fig. 4F). Meanwhile, we found that the PCA scores of patients in HDFM were generally higher than those in LDFM and MDFM (Fig. 4G). All the patients in LDFM flowed to the low-score group, and HDFM comprised a large proportion of the high-score group (Additional file 1: Fig. S5B). Accordingly, the novel classification might be an effective indicator for prognostic prediction to improve patient management.

The association between multi-omics alterations and FDGs

We performed multi-omics analysis in conjunction with the PCA score and the novel classification to better explain the prognostic difference caused by the differential expression of FDGs. This included the gene mutation landscape, tumor mutation burden (TMB), copy number variations (CNV), microsatellite instability (MSI), and immune infiltration. Although the major mutation genes were similar between the two PCA score groups, more gene types were discovered to mutate in the low-score group, such as FAT4, ALPK3, SDK1, SETDB1, and so on (Additional file 1: Fig. S6A), revealing the more aberrant biological activities of tumors in the low-score group to some extent. After that, we attempted to find out whether TMB and MSI changed along with the increase in PCA score. Unfortunately, the negative correlations between PCA score and both TMB (R = − 0.15, p = 0.18), and MSI (R = − 0.19, p = 0.084) were without statistical significance (Additional file 1: Fig. S6B, F). Uniting PCA score and TMB, the survival conditions of the low-TMB + high-PCA score group, the high-TMB + high-PCA score group, the low-TMB + low-PCA score group, and the high-TMB + low-PCA score group became worse in order (Additional file 1: Fig. S6C), which demonstrated the priority of PCA score over TMB in the survival assessment of MESO. Among the CNV alterations of 39 FDGs, the loss frequency of TMEM158 and COL7A1, which were respectively located on chromosomes 3 and 13, was arrestingly high (Additional file 1: Fig. S6D, F). Through comparing the immune infiltration of LDFM, HDFM, and MDFM, we noticed that the activated CD4+ T cells, CD56-dim natural killer cells, regulatory T cells, and type 2 T helper cells might have a tendency to gather in LDFM (Additional file 1: Fig. S7A), which is in line with the negative association between these cells and PCA score (Additional file 1: Fig. S7C). The expression of PD-L1 was not significantly different between the two groups (Additional file 1: Fig. S7B). On account of the vital role of immune cells in the anti-tumor and pro-tumor processes, the difference in the immune infiltration patterns of LDFM, HDFM, and MDFM may have contributed to the prognostic difference between them.

Investigation of the genes, transcription factors (TFs), and signaling pathways linked to Metastasis

Differential expression analysis was performed, contrasting RNA-seq data of primary tumors from the TCGA database with data of bone metastatic tumors from the MET500 database. The heatmap and volcano plot in Additional file 1: Fig. S8A visualized the expression patterns of every FDG in primary and metastatic tumors. There were 8 overexpressed FDGs, including PSAT1, FADS1, MFAP5, LDLR, COL4A1, COL4A2, MCAM, and KRT14 (four in state 1 and four in state 2), and 2 downregulated FDGs, including HP and CFB (all in state 3), in metastatic tumors (Additional file 1: Fig. S8B). Most of the upregulation and downregulation of transcription factors occurred without statistical significance (Additional file 1: Fig. S8C). Ultimately, we compared the enrichment degree of 50 typical pathways in primary and metastatic tumors and observed 49 pathways that were upregulated in metastasis, including the apical junction, epithelial-to-mesenchymal transition (EMT), IL2-STAT5 signaling, G2M checkpoint, E2F targets, TGF beta signaling, Hedgehog signaling, and so on (Additional file 1: Fig. S8D). These genes and pathways might play a role in the advancement of MESO as motivators or suppressors.

Establishment and internal verification of a novel prognostic prediction model

Nine MSigDB gene categories and the corresponding 42 pathways were recognized from FDGs in ORA analysis of the primary MESO and bone metastatic tumors, as shown in Fig. 5A. The hazard ratios for 39 FDGs were included in Fig. 5B based on a univariate Cox regression analysis between gene expression and related overall survival (OS), identifying six protective factors (ADH1B, CFB, PRG4, PLAAT4, HP, and IFIT3) and 33 hazardous factors with statistical significance (p < 0.001). In order to make the best use of the relationship between FDGs and survival, we applied Lasso regression analysis to select optimal genes from these FDGs for constructing a prognostic prediction model, decreasing bias, and finally including six key FDGs in the risk score formula (Additional file 1: Fig. S9A, B). The whole cohort of MESO was randomly divided into a train cohort and a test cohort at a ratio of 6:4 to validate the model internally. According to the individual risk score, patients in the three cohorts (all, train, and test cohorts) were segmented into a high-risk group and a low-risk group. The distribution of the two groups in each cohort was shown in Additional file 1: Fig. S9C. Among the six key FDGs, ADH1B and PLAAT4 (protective factors) were generally highly expressed in the low-risk group, while the expression of CDC20, CKS2, IPT1, and LDLR (hazardous factors) was relatively high in the high-risk group (Fig. 5C and Additional file 1: Fig. S9D). It was not surprising that the survival condition of low-risk was distinctly superior to the high-risk score in the three cohorts (Additional file 1: Fig. S9E), which was in accord with the phenomenon that the survival time showed a downward trend with the increase in risk score (Additional file 1: Fig. S9F). With ROC curves and the areas under the curves (AUC = 0.916 in the all cohort, AUC = 0.911 in the train cohorts, and AUC = 0.914 in the test cohort), the model also displayed excellent sensitivity and specificity (Additional file 1: Fig. S9G). Further, as depicted in Fig. 5D, we performed univariate and multivariate Cox regression analyses for age, gender, stage, bone metastasis, and risk score and proved that risk score could be regarded as a prognostic factor independent of these clinical measures (p < 0.001 in both Cox regression analyses; HR = 2828.449 in the univariate model and HR = 3201.584 in the multivariate model). Via a comparison of some clinical characteristics between the low-risk group and the high-risk group, we found that more patients were alive in the low-risk group (p = 0.003) despite the high mortality in both groups (Fig. 5E, F). In the GSEA analysis of both groups, pathways related to cell cycle and migration, such as the G2M checkpoint, E2F targets, and EMT, appeared to be active in the high-risk group, which might be the reason for the poor prognosis (Fig. 5G). Moreover, the contrast in immune infiltration between the two groups provided some information on the prognosis difference (Additional file 1: Fig. S10A–C). Survival analyses were carried out for immune cells and components that had statistically significant differences between the two groups (Additional file 1: Fig. S10D). The result demonstrated that the resting mast cells, plasma cells, neutrophils, and Treg T cells, which were highly infiltrative in the low-risk group, were connected to the higher survival probability with a p < 0.05. High aggregation of M1 and M2 macrophages, APC-co-inhibition, and APC-co-stimulation, on the other hand, were associated with a lower survival probability (p < 0.05).

Fig. 5

Construction of prognostic prediction model. A Over-representation analysis of 39 FDGs in the MSigDB gene sets, comparing the enrichment results of primary MESO with bone metastatic tumors. B The univariate Cox regression analysis between gene expression and overall survival (OS). C The expression levels of PLAAT4, ADH1B, CDC20, CKS2, JPT1 and LDLR, the genes selected to be included in the risk score calculation, in the low-risk group and high-risk group of the all, train and test cohorts. The all cohort referred to the entire TCGA cohort of MESO patients from which the train and test were randomly separated at the ratio of 6:4. D The forest plots showed the univariate and multivariate Cox regression analyses for clinical parameters (age, gender, stage, bone metastasis and distant metastasis) and risk score. E The clinical parameters’ differences between the low-risk group and high-risk group. F The comparison of survival outcomes (alive or dead) between the low-risk group and high-risk group with chi-square test. G The gene set enrichment analysis (GSEA) of the high-risk group and the low-risk group, based on the hallmark gene sets

The potential regulatory mechanism of LDFM, HDFM and MDFM

In order to investigate the possible functions and molecular mechanisms of FDGs in the development of MESO, integrated regulatory networks were created, identifying crucial transcription factors (TFs), pathways, immune cells, immune components, and RPPA that were associated with marker FDGs in the each MESO subtypes (Fig. 6A). Heatmaps were used to depict the correlation situations of regulatory network elements with one another (Fig. 6B).

Fig. 6

The prediction of regulatory mechanism in each subtype. A The regulatory networks of each MESO subtype, containing critical FDGs (blue circles), TFs (yellow circles), pathways (red circles), immune cells (blue-green circles), immune components (brown circles), and RPPA (green circles). The red lines between circles signaled positive correlation and the blue ones symbolized negative correlation with p < 0.0001. B The co-expression heatmaps reveal the interrelationship between components in each regulatory network. C The ATAC-seq results of the FDGs in the regulatory networks. D The violin plots demonstrated the IC50 distribution of the predicted 24 inhibitors among 3 different subtypes

MESO subtype 1: lowly differentiated fibroblast-related mesothelioma (LDFM)

SERPINE1 and INHBA, two risk genes related to differentiated state 1 (Fig. 4A), appeared in the regulatory network of LDFM (Fig. 6A), showing positive correlation with several signaling pathways, including apical junction, EMT, IL2_STAT5 signaling, TGF_beta signaling, and TNFA_signaling via NFκB. The expression level upregulation of the genes might be mediated by TFs, such as POLR3D, NOTCH1 and ETS1, which were widely expressed at high levels in bone metastatic tumor samples (Additional file 1: Fig. S11A). The ATAC-seq analysis revealed increased chromosome accessibility near SEPINE1 and INHBA, suggesting potential for binding of TFs or other regulatory elements to these regions (Fig. 6C). Also, the negative association of SERPINE1 and INHBA with resting mast cells that was shown to be favorable for survival in immune infiltration analysis (Additional file 1: Fig. S10D) and the positive association with APC_co_stimulation which were detrimental to survival, suggested the malignant regulatory mechanism of these two genes in LDFM. In addition, with the enrichment of some energy metabolic pathways (Additional file 1: Fig. S11B), the hypermetabolic trait of LDFM was unveiled, laterally reflecting the highly malignant nature of this subtype.

MESO subtype 2: highly differentiated fibroblast-related mesothelioma (HDFM)

In the regulatory network of HDFM (Fig. 6A), UBC2E and CDC20 had a positive relationship with the E2F target and G2M checkpoint, in agreement with the pathway enrichment result that HDFM was enriched in the cell cycle pathway (Additional file 1: Fig. S11B), hinting the highly proliferative ability. Another FDG in the network, COL4A2, was tied to EMT and TGF_beta signaling. The linkage between these three FDGs and immune components indicated that they might promote some immune infiltration that was inimical to survival. Moreover, chromosomal accessibility was also increased in regions close to these three genes (Fig. 6C), which was possibly in association with the three TFs (NCAPG, MYBL2 and NOTCH1) in the regulatory network.

MESO subtype 3: moderately differentiated fibroblast-related mesothelioma (MDFM)

Dissimilar to the genes identified in the LDFM and HDFM-related regulatory networks, both HP and CFB in the MDFM regulatory network were regarded as favorable factors for survival and associated with differentiation state 3 (Figs. 4A and 5B). Their negative connections with pathways, including EMT, hedgehog and TGF_beta signaling, also affirmed their advantageous role. At the same time, we found that their expression was also significantly lower in bone metastases than in primary tumors (Additional file 1: Fig. S11A), implying that the expression deficiency or downregulation of HP and CFB might lead to the augmented invasion of tumors. Actually, the expression of these two FDGs was higher in HDFM than in MDFM, as illustrated in Fig. 4D and Additional file 1: Fig. S5C, which partly explained the best prognosis of HDFM among the three subtypes.

Prediction of inhibitors

For increasing drug options to improve the survival of MESO patients with poor prognosis, we estimated the sensitivity of the three subtypes to various inhibitors by leveraging the pRRophetic R package. 24 inhibitors were finally filtered out, to which LDFM was more sensitive than HDFM (Fig. 6D), which were AG.014699, AP.24534, axitinib, AZD.0530, Bicalutamide, CCT018159, CCIR.99021, CMK, Embelin, FTI.277, Bexarotene, GDC0941, Imatinib, KU.55933, NU.7441, NVP.TAE684, OSI.906, BI.D1870, Parthenolide, RDEA119, Pazopanib, Pyrimethamine, PF.02341066, and PD.173074.

Clinical specimen validation

To further validate the clinical subtyping of MESO patients into LDFM, HDFM and MDFM, clinical specimens were enrolled for wet experiment validation. The immunohistochemical staining slides of the five different markers, INHBA (LDFM), UBE2C (HDFM), CDC20 (HDFM), HP (MDFM), and CFB (MDFM) in ×100 and ×400 fields of the light microscope were visualized in Fig. 7A. Subsequently, Pearson Chi-square tests were carried out to demonstrate significant differences in the differentiation (p = 0.02, Fig. 7B), immune infiltration (p = 0.02, Fig. 7C) and metastasis (p = 0.02, Fig. 7D) among LDFM, HDFM, and MDFM subtypes. The details were shown in Additional file 1: Table S2. From that, we could deduce that HDFM indicated good differentiation, high immune infiltration, and non-metastasis of MESO, while LDFM suggested poor differentiation, low immune infiltration, and metastasis of MESO, with MDFM in the middle, which was a decent validation for our clinical classification of MESO patients.

Fig. 7

Clinical specimen validation. A The immunohistochemical staining slides of the five different markers, INHBA, UBE2C, CDC20, HP and CFB in ×100 and ×400 fields of the light microscope in each MESO subtype. B A significant difference was revealed in the differentiation of LDFM, HDFM, and MDFM subtypes (p = 0.02). C A significant difference was obtained in the immune infiltration of LDFM, HDFM, and MDFM subtypes (p = 0.02). D A significant difference was shown in the metastasis of LDFM, HDFM, and MDFM subtypes (p = 0.02)