Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis

doi:10.21203/rs.3.rs-4827560/v1

Download PDF

Research Article

Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis

https://doi.org/10.21203/rs.3.rs-4827560/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Tumor heterogeneity is associated with poor prognosis and drug resistance, leading to therapeutic failure. Here, we used tumor evolution analysis to determine the intra- and intertumoral heterogeneity of high-grade serous ovarian cancer (HGSOC) and analyze the correlation between tumor heterogeneity and prognosis, as well as chemotherapy response, through single-cell and spatial transcriptomic analysis. We collected and curated 28 HGSOC patients’ single-cell transcriptomic data from five datasets. Then, we developed a novel text-mining-based machine-learning approach to deconstruct the evolutionary patterns of tumor cell functions. We then identified key tumor-related genes within different evolutionary branches, characterized the microenvironmental cell compositions that various functional tumor cells depend on, and analyzed the intra- and intertumoral heterogeneity as well as the tumor microenvironments. These analyses were conducted in relation to the prognosis and chemotherapy response in HGSOC patients. We validated our findings in two spatial and seven bulk transcriptomic datasets (total: 1,030 patients). Using transcriptomic clusters as proxies for functional clonality, we identified a significant increase in tumor cell state heterogeneity that was strongly correlated with patient prognosis and treatment response. Furthermore, increased intra- and intertumoral functional clonality was associated with the characteristics of cancer-associated fibroblasts (CAFs). The spatial proximity between CXCL12-positive CAFs and tumor cells, mediated through the CXCL12/CXCR4 interaction, was highly positively correlated with poor prognosis and chemotherapy resistance in HGSOC. In this study, we constructed a panel of 24 genes through statistical modeling that correlate with CXCL12-positive fibroblasts and can predict both prognosis and the response to chemotherapy in HGSOC patients.

Tumor cell state

Functional clonality

Tumor evolution

HGSOC

Tumor transcriptomic heterogeneity

Microenvironment

Spatiotemporal transcriptome

CXCL12

Cancer-associated fibroblasts

Ovarian cancer (OC) ranks first in mortality among gynecologic malignancies [1]. High-grade serous ovarian cancer (HGSOC) is the most common and deadly subtype of OC [2]. Over 75% of HGSOC patients are diagnosed at an advanced stage with extensive malignant ascites and omental metastasis [3]. Complete staging surgery followed by platinum-based chemotherapy is the standard treatment for HGSOC [4]. However, over 25% of patients develop chemotherapy resistance within 6 months of initial treatment, and 70% experience recurrence within 2–3 years, ultimately succumbing to acquired drug resistance [5]. Therefore, platinum drug resistance is the primary cause of poor prognosis, recurrence, and death in HGSOC patients. Identifying clinical indicators that are closely associated with platinum chemotherapy resistance and HGSOC treatment prognosis is of great importance.

Inter- and intratumoral heterogeneity is a crucial factor influencing patient prognosis and treatment outcomes [6]. Intertumoral heterogeneity between patients is associated with diverse genetic mutations, epigenetic modifications, and transcriptional alterations. Intratumoral heterogeneity refers to the presence of different cell types within a tumor, resulting from interactions between tumor cells and the tumor microenvironment (TME) [6].

Previous research has focused on intertumoral heterogeneity by analyzing genomic characteristics to obtain molecular subtypes of tumors, identify different patient subgroups, and select targeted therapies. For instance, the composition types of cells within the tumor microenvironment can classify tumors as either cold tumors or hot tumors[3]. In HGSOC, it could be classified into four molecular subtypes: mesenchymal, immunoreactive, differentiated, and proliferative [7]. Among these, the immunoreactive subtype demonstrates a better prognosis and treatment response [7]. Identifying different mutational characteristics is crucial for treatment selection in patients. BRCA1/2, TP53, and PI3K/AKT/mTOR pathway genes, which are frequently altered in ovarian cancer, serve as prime targets for clinical diagnosis and therapy [8]. However, intratumoral heterogeneity causes drug resistance, leading to therapeutic failure [6]. Thus, we need to integrate molecular features of both inter- and intratumoral heterogeneity with functional heterogeneity to improve patient subclassification and the response to therapy.

Tumor heterogeneity is thought to adhere to the fundamental principles of Darwinian evolution, where individual cells harboring heritable mutations that enhance adaptability gain a survival advantage [9]. Natural selection drives clonal expansion, leading to the emergence of subclones with varying proliferative, migratory, and invasive capabilities [10]. The evolution of adaptive clones occurs within the dynamic tissue environment known as the TME, forming a complex local tumor ecosystem. Changes within the TME further influence genetic diversification and phenotypic outcomes, driving the compositional tumor heterogeneity [11–13]. Therefore, identifying the distinct tumor cell clone subtypes within HGSOC and the microenvironmental cell types they interact with, exploring the interactions between the microenvironment and tumor cells, and elucidating the patterns of tumor evolution are crucial for the treatment and prognosis prediction of HGSOC patients.

The TME, which includes stromal cells, immune cells, extracellular matrix, and a variety of soluble factors, can influence tumor progression and the response to therapy [13]. Single-cell and spatial transcriptomics technologies have significantly advanced the investigation of heterogeneity in HGSOC [7, 14]. For example, the presence of tumor-associated macrophages (TAMs) has been associated with poor prognosis in HGSOC, as these cells can promote angiogenesis, immunosuppression, and chemoresistance [15]. TGF-β (Transforming growth factor beta)driven cancer-associated fibroblasts (CAFs), mesothelial cells, and lymphatic endothelial cells are associated with a poor outcome, while plasma cells correlate with more favorable outcomes in HGSOC [7].

Understanding the evolutionary trajectories of tumor and TME cells is crucial to develope more effective and personalized treatment strategies [16]. The multifaceted intratumoral heterogeneity and the TME may shape the evolutionary dynamics of OC [17]. Understanding the tumoral heterogeneity and evolutionary biodiversity could provide conceptual knowledge about the evolving tumor cell community to improve patient outcomes.

2.1 Construction of a multi-omics atlas of HGSOC including chemotherapy response characteristics using single-cell and spatial transcriptomics

To enhance the understanding of the complex interplay between intra- and intertumoral heterogeneity in HGSOC, and to develop a comprehensive tumor evolutionary model that delineates the distinct tumor and microenvironmental profiles across diverse patient cohorts, we assembled a robust dataset (Fig. 1A). This dataset comprises four single-cell datasets and two spatial transcriptomics datasets from GEO datasets (Table S1), totaling 50 samples. Notably, 28 (57%) of these samples were obtained from patients prior to chemotherapy and who exhibited measurable responses to therapy (Fig. 1B). We also curated two bulk RNA-sequencing datasets that include chemotherapy treatment responses from GEO datasets, encompassing 24 samples (Table S1). Our analytical approach began with the application of advanced text-mining and machine-learning algorithms to three of the single-cell datasets. This enabled the construction of a sophisticated hierarchical evolutionary model specific to HGSOC. This model was then validated and further scrutinized for insights into the relationships between intra-tumoral microenvironmental heterogeneity and intertumoral population heterogeneity using the remaining two datasets. Subsequently, we focused on the spatial transcriptomics datasets to dissect the spatial distribution characteristics of microenvironmental heterogeneity In the merged test dataset, which consisted of 14 samples, we processed and analyzed a total of 87,288 cells after quality control (Fig. 1C). These cells were classified into several key subtypes: B cells (5,700, 6.53%), CAFs (17,664, 20.24%), endothelial cells (ECs; 1,174, 1.35%), epithelial cells (Epi; 20,273, 23.22%), monocytes (15,217, 17.43%), and T and NK cells (T_NK; 27,259, 31.23%). We used the CNV method and calculated the epithelial cell score to identify 18,273 individual tumor cells (Fig. 1C). Through this integrated approach, we aimed to not only elucidate the evolutionary dynamics of HGSOC but also provide a foundation for the development of personalized therapeutic strategies that account for the unique heterogeneity observed in each patient’s tumor.

2.2 Inter- and intratumoral heterogeneity in high-grade serous ovarian cancer

Patient prognosis and treatment response are significantly affected by inter- and intratumoral heterogeneity. We used text mining and consensus clustering to dissect the clonal evolution within individual patients, followed by hierarchical clustering and bootstrap analysis to establish intertumoral clonal evolutionary relationships. We identified three evolutionary branches, designated as BR1, BR2, and BR3, across 69 clusters in 14 samples from GSE140819, GSE154600, and GSE184880 (Fig. 2A). Notably, individual tumor samples tended to cluster under a specific evolutionary branch. To validate the robustness of our method, we replicated the analysis in a subset of the data and confirmed consistent evolutionary relationships among cell populations in dataset GSE184880 (Figure S1).

Having established the clonal evolutionary relationships across different samples, we identified key genes within these branches and conducted functional analyses (Table S2). Key genes of BR1 were enriched in pathways related to myeloid leukocyte migration (GO:0097529), neutrophil degranulation (GO:0043312), and digestion (GO:0007586). BR2 genes were predominantly associated with cellular responses to chemokines (GO:1990869), T-cell migration (GO:0072678), regulation of response to type II interferon (GO:0060330), and autophagy (GO:0006914). BR3 genes were enriched in pathways such as the generation of amyloid fibrils (GO:1990000), cellular response to hypoxia (GO:0071456), positive regulation of growth (GO:0045927), and collagen fibril organization (GO:0030199) (Fig. 2B). These findings suggest that BR3 is characterized by highly fibrotic, immune desert-like tumors, BR2 by T-cell-infiltrated hot tumors, and BR1 by predominantly tumor-cell-dominated cold tumors [18].

Further analysis revealed significant heterogeneity in the distribution of different microenvironmental cells across the branches, with CAFs notably enriched in the BR3 branch (Fig. 2C). To delve deeper into intratumoral heterogeneity, we dissected the heterogeneity of epithelial, fibroblast, monocyte, and CD8 cells. We identified six epithelial cell subpopulations (E_HSPA1B, E_JUND, E_MALAT1, E_S100A4, E_SPP1, and E_TOP2A), with E_SPP1 and E_TOP2A predominantly found in the BR3 branch (Fig. 2D). Secreted phosphoprotein 1(SPP1) is highly expressed in multiple tumor types and interacts with fibroblasts or T cells via pathways such as SPP1-CD44, correlating with tumor malignancy[19]. Type IIA topoisomerase(TOP2A), a cell cycle gene, indicates high proliferative states when overexpressed. Cell velocity and pseudo-time analysis revealed that epithelial cells evolve toward the BR3 subpopulation, with terminal pathways enriched in the cell cycle, epithelial differentiation, and lipid metabolism pathways (Fig. 2D).

We identified five fibroblast subpopulations (F_ACTA2, F_COL1A1, F_JUNB, F_LGALS3, and F_MMP11), most of which were concentrated in the BR3 branch (Fig. 2D). Cell velocity analysis showed enrichment of cell adhesion-related pathways at the onset and metabolic, chemokine, and MAPK pathways at the end (Fig. 2D). Eight macrophage subpopulations were identified (M_A2M, M_CCL3, M_CXCL10, M_HLA-DPB1, M_HMGB2, M_HSPA1B, M_MMP9, and M_SPP1), with M_CCL3, M_CXCL10, M_HLA-DPB1, and M_SPP1 enriched in the BR3 branch (Fig. 2D). Temporal analysis revealed the enrichment of lipid metabolism, granulocyte differentiation, and c-Jun N-terminal kinase (JNK) pathway pathways in the terminal stage (Fig. 2D). Four CD8 cell subpopulations were identified, with CD8_GZMA and CD8_TIGIT enriched in BR1 and BR2, and BR3 dominated by CD8_CXCR4. Cell velocity analysis showed enrichment of T-cell activation, lymphocyte differentiation, and autophagy regulation pathways in the terminal stage of cell differentiation (Fig. 2D).

These results underscore significant inter- and intratumoral cellular heterogeneity in HGSOC, with epithelial cells in BR3 samples exhibiting high proliferation rates and infiltration by diverse fibroblast subtypes.

2.3 Tumor functional clonality is associated with patient prognosis and TME composition

In our previous analysis, we divided all single-cell samples into three branches and identified significant heterogeneity between these branches, along with 150 characteristic genes for each branch (Table S2). We further examined the prognostic clustering effect of these genes. Prognostic analysis across five datasets (TCGA, GSE14764, GSE26193, and GSE26712; Table S1) revealed strong prognostic clustering in these datasets, with the poorest prognosis samples exhibiting high levels of fibroblast infiltration (G1 in TCGA and G2 in GSE14764, GSE26193, and GSE26712; Fig. 3A, B). Correlation analysis also indicated a significant high correlation between fibroblast content and fibroblast-associated marker genes (COL3A1, COL1A2, COL1A1, FN1, and DCN; Figure S2). Subsequently, we used the obtained intertumoral heterogeneity branch genes as a reference gene set to perform a hypergeometric test in different bulk samples to calculate the distribution characteristics of samples in the single-cell BR branches. Bulk samples with the poorest prognosis and high fibroblast infiltration were significantly enriched in BR3 scores (Fig. 3C), further validating the reliability of using single-cell samples for studying intertumoral heterogeneity in tumors. Differential gene analysis in bulk samples also revealed the enrichment of numerous pathways related to epithelial cell migration, cell adhesion, and chemokines in samples with poor prognosis (Fig. 3D). These findings suggest that the fibroblast content in HGSOC samples may be highly correlated with tumor malignancy.

2.4 Identification of a subpopulation of CXCL12-positive fibroblasts associated with chemotherapy resistance and poor prognosis

Previous research has demonstrated a strong correlation between the heterogeneity of TME composition, especially the infiltration of fibroblasts, and prognosis[2]. Given that chemotherapy is commonly used in the clinical treatment of ovarian cancer, we explored the relationship between intratumoral heterogeneity and chemotherapy resistance. We analyzed the differences in the cellular composition of the TME in tumors with varying responses to chemotherapy. Using the 150 characteristic genes obtained (Table S2), we conducted a tumor evolution analysis on five samples within the single-cell dataset (GSE154600). The chemotherapy-sensitive samples GSM4675276 and GSM4675276 belonged to the BR1 subgroup, while the chemotherapy-resistant samples GSM4675273 and GSM4675273 were classified under the BR3 subgroup (Fig. 4A; Table S3). Combining these five samples, a total of 50,571 cells were obtained, including B cells, CAF, monocytes, T cells, EPI_1 (epithelial cells with high levels of Matrix metallopeptidase 7 (MMP7) and E74 like ETS transcription factor 3 (ELF3)), and EPI_2 (epithelial cells with a high level of hes family bHLH transcription factor 1 (HES1)) (Figs. 4B and S3). Among these, B cells were significantly enriched in samples derived from BR1 and BR2, while CAFs were significantly enriched in BR3 samples, indicating that the degree of fibroblast infiltration not only affects patient prognosis but also correlates with the treatment response (Fig. 4C). To elucidate the relationship between CAFs and ovarian cancer treatment response, we further sub-classified CAFs and identified seven distinct subtypes (Fig. 4D). F_CCL21 and F_RAMP3 were significantly enriched in BR1, and F_PGF, F_IGFBP4, F_IGFBP7, and F_CXCL12 were significantly enriched in BR3, particularly the F_CXCL12 subtype, which was almost exclusively found in BR3 (Fig. 4E). We then used gene co-expression analysis to identify key characteristic genes in different CAF subpopulations, resulting in four gene modules, with the blue module showing the highest correlation with the BR3-specific F_CXCL12 subgroup (R = 0.88, P < 0.001, Fig. 4F) and the yellow module highly correlated with the BR1-specific F_RAMP3 module (R = 0.97, P < 0.001; Fig. 4F). Analysis of the hub genes in the blue module revealed that these genes were primarily enriched in pathways related to extracellular matrix remodeling and tube morphogenesis (Figs. 4G and S5). Given the significant enrichment of CAFs in BR3, we explored the communication between these fibroblasts and tumor cells and found high-intensity interactions between the fibroblast subpopulations F_CXCL12 and F_IGFBP7 and tumor cells in EPI_1 (Fig. 4H). Specifically, F_CXCL12 interacts with tumor epithelial cells via the ligand-receptor pair CXCL12/CXCR4 and through the secretion of placental growth factor (PGF) (Fig. 4I). These findings suggest that CXCL12-positive CAF cells may be associated with chemotherapy resistance.

To validate the relationship between fibroblasts and chemotherapy, we examined the relationship between fibroblasts and treatment response in the other dataset (GSE165897) that included chemotherapy response information (Table S4). Using the previously identified 150 characteristic genes, we performed a hierarchical evolutionary analysis between samples and conducted permutation tests. The platinum-free interval (PFI) values of samples of EOC3, EOC349, EOC540, and EOC87 in the BR3 subgroup were significantly lower than those in the groups BR1 and BR2, indicating significant chemotherapy resistance, while the samples EOC136 and EOC153 in BR1 showed chemotherapy sensitivity, again confirming the high correlation between intertumoral heterogeneity and treatment response (Fig. 5A; Table S4). From this dataset of nine samples, we obtained 17,898 cells, comprising 11 subpopulations, including two fibroblast subpopulations, F_CXCL12 and F_COL1A1, with significant F_CXCL12 infiltration in BR3 and significant F_COL1A1 enrichment in BR1 (Figs. 5B, C and S6). We then re-clustered the fibroblasts and identified five fibroblast subpopulations, including the F_CXCL12 subgroup (Figs. 5D and S6). Using gene co-expression analysis, we also obtained a gene module that was highly correlated with the F_CXCL12 subgroup (R = 0.85, P < 0.001; Fig. 5E), with characteristic genes primarily enriched in TGF-beta, MAPK, and cytokine interaction signaling pathways (Figure S7).

In both datasets, we identified a fibroblast subgroup with high CXCL12 expression, and a similarity correlation analysis between subgroups confirmed the high similarity of CXCL12-high fibroblast subgroups in both datasets (Fig. 5E).

2.5 Identification of a gene set affecting prognosis and drug resistance in HGSOC

Above, we identified a significant correlation between the highly invasive CXCL12-positive CAF subpopulation and both treatment response and prognosis in ovarian cancer. We next conducted a comparative analysis of the WGCNA hub gene results to consolidate the characteristic genes of the CXCL12-positive CAF subpopulation and identified 24 shared genes (Table S5). Subsequent Cox regression analysis using the GSE26193 dataset revealed that DCN, CXCL12, and TNFAIP6 were significantly positively associated with poor prognosis in ovarian cancer, with DCN specifically identified as a fibroblast-characteristic gene (Fig. 6A). Using these significantly differential genes and Cox regression coefficients, we performed a risk analysis across four additional datasets(GSE14764, GSE26712, TCGA and GSE9891, Table S1). The findings indicated that the risk coefficients for the high gene expression groups were notably higher than those for the low expression groups, and both CXCL12 and DCN showed significant positive correlations with poor overall patient prognosis (Fig. 6A, B). Through Friends analysis, we observed a high correlation of CXCL12 with other genes (Fig. 6C). Additionally, using these 24 genes, we applied enrichment analysis in three treatment datasets (GSE33482, GSE114206 [20], and GSE189843 [15]) to validate the association between CAFs and chemotherapy response (Fig. 6D; Table S6). Based on gene expression data from these samples, we categorized them according to their enrichment levels in the 24 genes. These 24 genes could effectively stratify samples into high and low groups, with the high CAF subpopulation exhibiting strong drug resistance. Moreover, the 24 genes showed excellent predictive efficacy for drug resistance, with AUC values of 1, 0.89, 0.83 in datasets GSE33482, GSE114206 and GSE189843 respectively (Fig. 6D). These findings underscore the critical role of fibroblasts in OC treatment and prognosis, and through further feature extraction, we have identified a reduced set of genes that are indicative of treatment response.

2.6 Spatial distribution characteristics of CXCL12-positive fibroblasts affect the clinical prognostic results of HGSOC

Given the strong correlation between CXCL12-positive fibroblasts and the treatment outcomes and prognoses of HGSOC, as well as their interaction with tumor cells via CXCR4, we elucidated the interaction mechanisms of CXCL12-positive fibroblasts in vitro and in spatial transcriptomic data. In vitro, we silenced the CXCL12 receptor gene CXCR4 in the ovarian cancer cell line SKOV3 (Fig. 7A). As shown by qPCR, all three siRNA knockdown fragments significantly suppressed CXCR4 mRNA expression levels in SKOV3 cells (P < 0.05; Fig. 7A). Subsequent addition of exogenous CXCL12 protein to both the control group and the CXCR4-si group demonstrated that CXCR4 knockdown significantly inhibited cell viability compared with the siNC control group (P < 0.05; Fig. 7B, C). Considering that fibroblasts and tumor cells interact through cellular communication, and the intensity of this communication is highly related to the spatial positioning of genes, we analyzed the spatial distribution of fibroblasts and tumor cells within a treatment response-associated spatial transcriptomic dataset. Using the inferCNV method and cell type identified, we identified tumor regions and the surrounding cellular environments. We analyzed the spatial distribution characteristics of cells in chemotherapy-resistant (GSM6506110) and chemotherapy-sensitive (GSM6506114) samples in GSE211956 dataset respectively (Fig. 7D, 7E, Table S7). In the chemotherapy-resistant samples, we identified 5 clusters (Fig. 7D(i)). Cell type analysis revealed that cluster 4 is enriched with fibroblasts with high expression level both CXCL12 and DCN (Fig. 7D(iii)). Furthermore, in the surrounding area of cluster 4 (Nbs_4), which are enriched in tumor epithelial cells and expressed high level of tumor cell marker KRT19 and CXCR4 (the receptor for CXCL12) (Fig. 7D(iii)). In contrast, in the chemotherapy-sensitive samples, we identified 4 clusters (Fig. 7E(i)). Among these, clusters 0, 1, and 3 are enriched with fibroblasts (Fig. 7E(ii)). Cluster 1 is specifically enriched with KRT19-positive tumor cells (Fig. 7E). However, in the regions dominated by KRT19-positive tumor cells, there is no high expression of CXCR4, forming a spatial interaction domain of CXCL12-positive CAFs. Finally, in our treatment response-associated clinical samples, we also validated the local spatial adjacency of CXCL12/DCN/KRT19 in chemotherapy-resistant samples (Fig. 7E).

Solid-organ malignancies exhibit heterogeneity that enhances tumor cell survival and drug resistance through various molecular mechanisms [30]. Intratumoral heterogeneity evolves both spatially and temporally during tumor development, reprogramming the tumor microenvironment [31]. Understanding tumor clonality and its evolutionary path is crucial for comprehending tumor cell behavior and identifying driving factors [32, 33]. Tumors harbor numerous genetic and epigenetic alterations, many of which complicate the interpretation of their functional effects on tumor evolution [34]. Single-cell and spatial transcriptomics can delineate the roles of different cell types within tumors, providing insights into the functional diversity arising from genetic heterogeneity and adaptation [35–39].In this study, we integrated single-cell, bulk, and spatial transcriptomic data; elucidated the intratumoral heterogeneity (including tumor and TME cells of HGSOC); explored their evolutionary patterns of intra- and intertumoral heterogeneity; and identified critical factors influencing tumor prognosis and therapeutic responses.

Single-cell data processing faces challenges with large-scale datasets, necessitating efficient computational and statistical methods to handle data sparsity and ensure informative analysis [21]. Batch effects between datasets require robust correction techniques to guarantee accurate cross-dataset comparisons and integrations [21]. Determining the optimal number of clusters is essential for the analysis of intratumoral heterogeneity, yet this process is complicated by the lack of clear standards [22]. In this study, we initially identified intratumoral heterogeneity using results from Seurat’s initial clustering combined with CNV data. We extracted tumor-associated feature genes and identified “representative cells” that reflect the characteristics of different tumor clones from subpopulations using the TF-IDF method of text mining. We then analyzed the intratumoral heterogeneity of these cells. This approach ensures the extraction of the most accurate and essential tumor clone information while minimizing computational resource consumption. Additionally, it mitigates the effect of sparsity in single-cell transcriptome data on data analysis. Furthermore, we used the representative cells to determine the clustering numbers automatically and therefore confirm the intratumoral heterogeneity using SC3 consensus clustering and support vector machine methods. This strategy addresses the subjectivity associated with the manual definition of cluster numbers. Finally, by obtaining precise intratumoral heterogeneity and related feature genes, we constructed a tumor evolution model based on their expression levels and analyzed the characteristics of intra- and intertumoral heterogeneity in HGSOC. Through the integrated analysis of intra- and intertumoral heterogeneity, we can more comprehensively elucidate the roles of tumor and microenvironmental cells in the prediction of prognosis and chemoresistance of HGSOC. This provides a novel strategy for interpreting the correlation between intra- and intertumoral heterogeneity using expression data.

Using the aforementioned strategy, we analyzed and validated across multiple datasets that HGSOC patients can be classified into three evolutionary branches: BR1, BR2, and BR3. Notable heterogeneity was observed within each branch concerning the quantity and composition of tumor and microenvironmental cells. The BR1 branch was predominated by tumor cells and resembled cold tumors, BR2 was infiltrated with high level of B and T cells and resembled hot tumors, and BR3 was heavily infiltrated by fibroblasts, akin to immunologically excluded tumors [18]. By correlating tumor evolution with prognosis and resistance to therapy, we found that a high infiltration of fibroblasts was associated with a tumor cell subpopulation enriched in cell cycle pathways, low CD8 cell activity, increased cell proliferation, and malignant transformation in cell morphology, resulting in significantly worse prognosis and therapy resistance. Fibroblasts can be broadly classified into myofibroblasts, TGF-β-driven CAFs, inflammatory fibroblasts, and antigen-presenting fibroblasts and are highly associated with poor prognosis in various cancers [23–25]. In HGSOC, fibroblasts can form a dense physical barrier around tumor cells, impede the infiltration of immune cells, and promote tumor cell proliferation through the secretion of cytokines [26, 27]. Combining with WGCNA method, we identified a clinical predictive model comprising 24 genes and discovered that fibroblasts with high CXCL12 expression are significantly related to poor prognosis and chemotherapy resistance. The differences in the spatial distribution of cells, particularly the characteristics and functions of cells near the tumor boundary, play a crucial role in the development of tumor heterogeneity [28]. To this end, we developed a bioinformatics analytical workflow for tumor boundary identification and found that this type of fibroblast significantly infiltrates the margins of chemotherapy-resistant tumors. Spatial interactions between cell clusters may have a more profound effect on chemo responsiveness than the composition of the clusters alone. In HGSOC, fibroblasts can interact with other cell types through various pathways, ultimately contributing to chemotherapy resistance [25]. CXCL12-positive CAFs could promote cancer cell migration and invasion and upregulate the expression of PDL1 in bladder and pancreatic cancer cells [29, 30]. They can also attract monocytes through the CXCL12/CXCR4 pathway and induce their differentiation into M2 macrophages, which leads to enhanced tumor cell proliferation and reduced apoptosis in oral squamous cell carcinoma [31]. In this study, by analyzing the evolutionary pathways of tumor and microenvironmental cells, we determined that CXCL12-expressing fibroblasts, through spatially proximal interactions with tumor cells, influence patient prognosis and therapeutic outcomes, offering a new perspective for the prevention and treatment of HGSOC.

4.1 Preprocessing of single-cell RNA-seq data

We collected single-cell RNAseq matrix from five HGSOC datasets (GSE140819, GSE154600, GSE165897, GSE184880 and data in http://blueprint.lambrechtslab.org, Table S1). Raw gene expression matrices were processed using the Seurat R package (version 3.2.3) [32]. Cells of low quality were filtered out based on two criteria: 1) cells with fewer than 1000 unique molecular identifiers (UMIs) or with fewer than 100 genes detected; 2) cells with more than 30% of their UMIs originating from mitochondrial genes. For each matrix, gene expression matrices were normalized using the LogNormalize method via the NormalizeData function, and highly variable genes were identified using the scran package [33]. Dimensionality reduction was performed with principal component analysis (PCA), where the number of principal components (PCs) was selected based on the JackStraw function(Seurat). Cell clusters were identified using the FindClusters function(Seurat)and visualized through uniform manifold approximation and projection (UMAP). Differential genes between clusters were calculated by the Wilcoxon rank-sum test with a family-wise error rate set at 5%.

4.2 Identification of tumor cells

Copy number variations (CNVs) were estimated in each cell by analyzing the averaged expression profiles across chromosomal intervals. The initial estimation of CNVs for each chromosomal region was conducted using the infercnv R package (version 1.0.4; inferCNV of the Trinity CTAT Project. https://github.com/broadinstitute/inferCNV). In epithelial cells, CNVs were determined by comparing the expression levels derived from the single-cell RNA-sequencing (scRNA-seq) dataset, with a cutoff value of 1 and a noise filter set at 0.2. For each sample, gene expression data were re-standardized, and the values were constrained within the range of − 1 to 1. The CNV for each cell was calculated as the sum of the squares of the CNV values across the respective regions.

4.3 Calculation of epithelial cell scores

We used the expression levels of epithelial scores based on the mean expression of five well-established markers: KRT19, KRT7, KRT18, KRT8, and EPCAM, as reported in previous studies [24].

4.4 Tumor evolution analysis

4.4.1 Selection of representative tumor cells

We input the clustering information from Seurat and simultaneously identified the tumor cell subgroups using the CNV method. We used the normalized matrix to calculate the PC scores for each cell. Subsequently, we used the Wilcoxon rank-sum test to obtain PCs that showed significant differences between classes (FDR < 0.05, false discovery rate). Finally, we selected the top 10 genes with the highest absolute coefficients of the differentially expressed PCs as the characteristic genes for the corresponding PCs.

We used the characteristic genes obtained from the tumor subgroups and applied the term frequency-inverse document frequency (TF-IDF) method. This technique is typically used in information retrieval and text mining to evaluate how important a word is to a document in a collection or corpus. In the context of gene expression analysis, TF-IDF (Term Frequency-Inverse Document Frequency) can be used to identify the most relevant genes in the tumor subgroups. TF-IDF measures the importance of a gene based on how frequently it appears in a particular subgroup compared to its occurrence across all subgroups in the dataset, thus highlighting genes that are both frequent and unique within specific tumor subgroups.

$$\:{S}_{i,j}=\frac{{N}_{i,j}}{\sum\:_{i=1}^{m}\:{N}_{i,j}}*log(1+\frac{n}{\sum\:_{j=1}^{n}\:{N}_{i,j}})$$

Where m and n represent the number of genes and cells in the single-cell expression matrix, respectively. $\:{N}_{i,j}$ is the $\:m$ by $\:n$ normalized gene expression matrix defined above. Within the TF-IDF method, we set the lower quartile of a gene's expression across all cells as the threshold. If a gene's expression was below this threshold, we set its expression to 0. This approach ensures that only genes with a certain level of expression are considered significant, effectively filtering out genes that are expressed at very low levels across the cells. After applying the TF-IDF method to obtain the gene weight matrix $\:{S}_{i,j}$, we calculated the 95th percentile (A) of each gene's weight across all cells and set the lower threshold for gene weights at 0.25 × A. Any gene weight below this threshold in $\:{S}_{i,j}$ was set to 0, resulting in a corrected $\:{S}_{i,j}$ matrix. Subsequently, for each cell group, we summed the weights of the characteristic genes to derive a total weight and set the cutoff value $\:\text{T}$ to establish the cell weight threshold utilizing the formula:

$$\:\text{T}=\frac{1}{2}{Q}_{2}+\frac{1}{4}({Q}_{1}+{Q}_{3})$$

where Q1, Q2, and Q3 represent lower, median, and upper quartile values, respectively. Cells with a score higher than T were identified as representative tumor cells.

4.4.2 Clustering and refinement

Using the high-purity tumor cells obtained in the first step for each tumor subgroup, we employed consensus clustering method (SC3[34]) for consensus clustering to obtain subclassification labels for each subgroup. Subsequently, we merged the subcategories of each subgroup to identify those with statistically significant differences. Specifically, we first extracted the top 1000 highly variable genes from the PCA of high-purity malignant cells. For each subgroup’s subclassification obtained through the consensus clustering method SC3 [34], we trained a support vector machine model using the first 10 PCs and validated the classification accuracy through 10-fold cross-validation. Next, we randomly shuffled the labels 100 times to perform a permutation test to calculate the classification significance P-values between each pair of clusters. These P-values were then sorted in descending order, and an iterative merging process was conducted until the largest P-value was less than 0.05, at which point the iteration ceased. Through this method, we identified the subclusters within each tumor subgroup that were of the highest purity and had statistically meaningful classifications.

4.4.3 Hierarchical evolutionary analysis

We intersected the features derived from TF-IDF with those obtained through PCA and obtained tumor-associated feature genes. Using these genes, high-purity tumor cells, and the purified subpopulations of each subgroup, we computed the mean expression values of the corresponding subgroup’s feature genes to generate a pseudo-bulk dataset. We then applied the removeBatchEffect function in the Limma package (version 3.58.1) to remove batch effects across various datasets. Subsequently, we performed hierarchical evolutionary analysis between subgroups using the pvclust function from the pvclust package (version 2.2, nboot = 1000).

Following the aforementioned step, we obtained highly purified subpopulations within tumors, facilitating the analysis of intratumoral heterogeneity.

4.4 Cell-type enrichment analysis

Each cell-subtype analysis encompassed a range of tumor stages, for which we computed enrichment scores (EScores) to quantify the cell-subtype prevalence at different stages. These EScores reflect the ratio of cell subtype numbers at specific stages to their overall distribution and highlight when a subtype is predominantly enriched, with values greater than 1 signifying enrichment at that stage [24].

4.5 Weighted correlation network analysis in cell subclusters and acquisition of fibroblast-related Genes

The normalized expression matrix was used to construct a weighted gene co-expression network via the weighted correlation network analysis (WGCNA) R package (version 1.69). To mitigate the effect of noise and outliers, analysis was performed on “pseudo cells,” which represent the average gene expression of 10 randomly selected cells within each distinct cell type [35]. Network construction was achieved using the “blockwiseModules” function, applying the default settings. For each identified module, a PCA was performed using the module eigengenes. The correlation between module eigengenes and cell-type metadata was calculated to evaluate the relevance of each module using Pearson’s correlation test. Subsequently, hub genes within significant modules were identified based on their modular connectivity, which refers to the absolute value of Pearson’s correlation between genes (module membership) and their relationship with clinical traits, defined as the absolute value of Pearson’s correlation between individual gene expression and cell type. We performed WGCNA in the fibroblast subpopulations of datasets GSE154600 and GSE165897 from GEO database. [36] and took the intersection of characteristic genes of cell subpopulations significantly infiltrating the BR3 branch as the final fibroblast-related characteristic genes, which were used for the construction of the risk prognosis and drug resistance analysis model below.

4.6 Similarity analysis among cell subpopulations

We used the single-cell subpopulations and the characteristic genes obtained through WGCNA to calculate the mean expression levels of these genes in each subpopulation and sample phenotype (the evolutionary branch to which the sample belongs). We used the R package ggcor (https://github.com/hannet91/ggcor) to compute the correlations between cell subpopulations and between cell subpopulations and phenotypes and visualize these correlations.

4.7 Clustering of bulk RNA data samples

Bulk RNA data were retrieved from the Cancer Genome Atlas (TCGA; https://www.cancer.gov/tcga) and GEO datasets. Subsequently, non-negative matrix factorization (NMF) clustering methods were performed on the normalized expression data using the NMF R package (version 0.23).

4.8 Mapping of single-cell and Bulk RNA subpopulations

To establish the correspondence between the intertumoral clusters obtained from single-cell evolutionary analysis and the clusters identified through clustering methods in bulk, we first extracted the characteristic gene sets from the single-cell BR1, BR2, and BR3 subpopulations. Subsequently, we used a hypergeometric test to calculate the enrichment scores of the bulk RNA-seq samples within these three gene sets. We then selected the connections with the highest enrichment scores that were significantly enriched as the mapping relationship between single-cell and bulk RNA-seq samples.

4.9 Friends analysis

The Friends analysis approach assesses the functional correlation among various genes within a pathway, suggesting that the interaction of a gene with others in the same pathway enhances its likelihood of expression. Using the R package GOSemSim [37], we calculated the functional correlations among genes linked to the prognosis of HGSOC and drug resistance.

4.10 Spatial transcriptome data analysis

We collected two HGSCO spatial transcriptome datasets including chemotherapy treatment response (GSE189843 and GSE211956, Table S1) from GEO database. Both datasets are spatial transcriptomics obtained using the 10x Genomics platform, and they include matrices of spatial gene expression information where rows represent genes and columns represent spatial spots. Raw gene-spot matrices were analyzed with the Seurat package (version 3.2.3) in R. Spatial transcriptome data were qualitatively controlled using parameters including total spots, media UMIs/spot, median genes/spot, and median mitochondrial genes/spot. Spots used in the subsequent analysis were filtered for a minimum detected gene count of 200 genes while genes expressed in fewer than three spots were removed. Normalization across spots was performed with the SCTransform function. Dimensionality reduction and clustering were performed with PCA at a resolution of 1 with the first 30 PCs. We conducted cluster analysis using FindClusters and then used the standardized expression matrix to calculate the average expression levels of immune-related genes (PTPRC, CD2, CD3D, CD3E, CD3G, CD5, CD7, CD79A, MS4A1, and CD19) [28]. The subpopulation with the highest expression levels was selected as the normal control. We used the InferCNV method to identify subpopulations of tumor cells. Subsequently, using the RegionNeighbours function from the R package STutility (version 1.1.1; https://ludvigla.github.io/STUtility_web_site/), we determined the cells at the tumor edge. We then used the FindMarkers function to calculate the differentially expressed genes in the region adjacent to the tumor edge.

4.11 Immune infiltration analysis

We assessed the immune score of various immune cells in HGSOC patients by employing xCell (R package, version 1.1) on RNA-seq datasets and microarray datasets [38]. The microarray datasets underwent quantile normalization, while the RNA-seq dataset was quantified in terms of fragments per kilobase million(FPKM)

4.12 Cell transfection

The human ovarian cancer cells SKOV3 were purchased from the ATCC cell bank.

were digested with trypsin and resuspended. After cell counting, cells were plated at a density of 4 × 104 per well in a 24-well plate and incubated overnight in a cell culture incubator. siRNA transfection was performed according to the Lipofectamine 3000 (L3000015, Thermo Fisher Scientific (Waltham, MA, USA)) protocol and using 15 pmol siRNA and 1.5 µL lipofectamine per well. After 48 hours of incubation, cells were collected for qRT-PCR and Western blot. The sequences of siRNA for the knockdown and control groups were as follows:

siCXCR4-1 GGCAAUGGAUUGGUCAUCCUGGUCA;

siCXCR4-2 UGGUUGGCCUUAUCCUGCCUGGUAU;

siCXCR4-3 UGUUUCCACUGAGUCUGAGUCUUCA;

siNC UUCUCCGAACGUGUCACGUTT.

4.13 qRT-PCR

qRT-PCR was used to detect the transcriptional level of the CXCR4 gene and the knockdown effect of different siRNA sequences. RNAiso Plus(9109, TaKaRa, Tokyo, Japan)was added to the cell samples for lysis, and the samples were processed according to the Chloroform-Trizol RNA Extraction protocol After air-drying at room temperature, RNA pellets were dissolved in 20 µL Rnase-free H₂O. The concentration and A260/280 of RNA samples were measured using a Microplate UV-Vis Spectrophotometer. Reverse transcription was performed according to the Goldenstar® RT6 cDNA Synthesis Kit Ver.2 ༈TSK302S, Tsingke, Beijing, China)instructions, and 2 µg of total RNA was used as the template. The reaction mixture was incubated at 50°C for 5 minutes and 85°C for 2 minutes. The qPCR system was prepared according to the ArtiCanATM SYBR qPCR Mix༈TSE501, instructions, Tsingke, Beijing, China༉ and the real-time PCR program was set as follows: holding stage, 95°C 30 s; cycling stage (40 cycles), step 1 95°C 15 s, step 2 60°C 30 s; and melt curve stage: 95°C 15 s, 60°C 60 s, 95°C 15 sec. The primer sequences used were as follows:

CXCR4-h-F ACTACACCGAGGAAATGGGCT;

CXCR4-h-R CCCACAATGCCAGTTAAGAAGA;

GAPDH-h-F TGACAACTTTGGTATCGTGGAAGG;

GAPDH-h-R AGGCAGGGATGATGTTCTGGAGAG.

4.14 Western blot

Western blotting was used to detect the protein level of CXCR4. Cell samples were treated with RIPA(Radio-Immunoprecipitation Assay )lysis buffer (R0278, Sigma-Aldrich, St. Louis, MO, USA) containing PMSF༈Phenylmethanesulfonyl fluoride༉and protease inhibitors; sonicated on ice; and centrifuged at 4°C, 12,000 g for 10 minutes. Then, the supernatant was taken for protein quantification using the BCA (Bicinchoninic Acid Assay )method. After mixing with loading buffer and denaturing at 100°C for 5 minutes, proteins were separated using 10% SDS-PAGE gel and then transferred to a PVDF(Polyvinylidene fluoride) membrane. After blocking with 5% skimmed milk for 1 hour at room temperature, the membrane was incubated withCXCR4 (E3Q4B) Rabbit mAb (, Cell Signaling Technology, Danvers, MA, USA) diluted at a ratio of 1:1000 in TBST(Tris-Buffered Saline with Tween) containing 5% BSA(Bovine Serum Albumin) ,.After washing with TBST three times for 5 minutes each, incubate the membrane with Anti-rabbit IgG HRP-linked Antibody(7074S, Cell Signaling Technology, Danvers, MA, USA ) at a ratio of 1:2000 in TBST .Visualize the protein bands using ECL reagent (E-IR-R307, Elabscience Biotechnology Wuhan, China)according to the manufacturer's instructions.

4.15 CCK8 assay

The CCK8 assay was used to detect the effect of CXCL2/CXCR4 on cell viability. When cells reached 80–90% confluence, cells were digested with trypsin and resuspended. After cell counting, cells were plated at a density of 3000 cells/well in 96-well plates. After incubation in a 37°C incubator for 24 hours, the medium was replaced with complete medium containing 0, 50, 100, 200, or 300 ng/mL CXCL12 and cultured for 48 hours. CCK8 reagent (10 µL) was added per well, and the absorbance at 450 nm was measured using a microplate reader after a 1-hour incubation at 37°C. In separate experiments, siCXCR4 and siNC transfections were performed 24 hours after plating. Six hours later, CXCL12 (200 ng/mL) was added, and the CCK8 assay was conducted using the same method.

In this study, we aimed to dissect the complex tumor heterogeneity of high-grade serous ovarian cancer (HGSOC) and its implications on prognosis and chemotherapy response using advanced single-cell and spatial transcriptomics. By analyzing data from 28 HGSOC patients across five datasets, we developed an innovative text mining-based machine learning method to unravel the evolutionary dynamics of tumor cell functions. This approach revealed critical tumor-related genes and the varied microenvironmental compositions on which different tumor cell functions rely. Our findings highlighted a significant connection between increased tumor cell state heterogeneity and worse patient outcomes, including prognosis and treatment resistance. We validated these insights using additional spatial and bulk transcriptomic datasets, encompassing a total of 1,030 patients. We further identified that heightened intra- and inter-tumoral functional clonality is closely linked with the characteristics of cancer-associated fibroblasts (CAFs). Notably, the spatial proximity between CXCL12-positive CAFs and tumor cells, facilitated by the CXCL12/CXCR4 axis, emerged as a strong predictor of poor prognosis and chemotherapy resistance. Moreover, we developed a panel of 24 genes that are highly correlated with CXCL12-positive fibroblasts. This gene panel effectively predicts both prognosis and chemotherapy response in HGSOC patients. Our study underscores the critical role of tumor heterogeneity in therapeutic outcomes and offers new insights into potential biomarkers for personalized treatment strategies.

Supplementary Materials

The following supporting information can be downloaded online.

Conflicts of Interest

The authors have declared that no competing interest exists.

Funding

This study was supported by the Henan Province and Ministry of Health of Medical Science and Technology Program (SBGJ202302028 for Tingjie Wang and SBGJ202101009 for Yongjun Guo), Dalian Science and Technology Innovation Fund (2022JJ12SN049 for Jun Yang), and the Fundamental Research Funds for the Central Universities.

Author contributions

The study design and supervision were conducted by Yongjun Guo, Jun Yang, Jun Li, and Tingjie Wang. Data analysis, in vitro experiments, and manuscript writing were performed by Tingjie Wang, Lingxi Tian, and Ruitao Long. Clinical data collection, sorting, analysiss and multicolor immunofluorescence staining analysis were performed by Bing Wei, Cuiyun Zhang, Bo Wang, and Yougai Zhang. The patients’ clinical data and their therapeutic responses were verified and evaluated by Yougai Zhang and Xiaofei Zhu. The manuscript has been reviewed and approved by all authors.

Acknowledgments

The authors thank Dr. Mengyun Ke for the sample collection and storage.

Data Availability Statement

The datasets analyzed in this study have been deposited in the Gene Expression Omnibus (GEO) repository under the accession numbers in supplementary tables.

Coburn SB, Bray F, Sherman ME, Trabert B (2017) International patterns and trends in ovarian cancer incidence, overall and by histologic subtype. Int J Cancer 140(11):2451–2460
Chowdhury S, Kennedy JJ, Ivey RG, Murillo OD, Hosseini N, Song X, Petralia F, Calinawan A, Savage SR, Berry AB, Reva B, Ozbek U, Krek A, Ma W, da Veiga Leprevost F, Ji J, Yoo S, Lin C, Voytovich UJ, Huang Y, Lee SH, Bergan L, Lorentzen TD, Mesri M, Rodriguez H, Hoofnagle AN, Herbert ZT, Nesvizhskii AI, Zhang B, Whiteaker JR, Fenyo D, McKerrow W, Wang J, Schurer SC, Stathias V, Chen XS, Barcellos-Hoff MH, Starr TK, Winterhoff BJ, Nelson AC, Mok SC, Kaufmann SH, Drescher C, Cieslik M, Wang P, Birrer MJ, Paulovich (2023) A. G., Proteogenomic analysis of chemo-refractory high-grade serous ovarian cancer. Cell 186(16):3476–3498e35
Zheng X, Wang X, Cheng X, Liu Z, Yin Y, Li X, Huang Z, Wang Z, Guo W, Ginhoux F, Li Z, Zhang Z, Wang X (2023) Single-cell analyses implicate ascites in remodeling the ecosystems of primary and metastatic tumors in ovarian cancer. Nat Cancer 4(8):1138–1156
Matthews BG, Bowden NA, Wong-Brown MW (2021) Epigenetic Mechanisms and Therapeutic Targets in Chemoresistant High-Grade Serous Ovarian Cancer. Cancers (Basel) 13, (23)
Silva R, Glennon K, Metoudi M, Moran B, Salta S, Slattery K, Treacy A, Martin T, Shaw J, Doran P, Lynch L, Jeronimo C, Perry AS, Brennan DJ (2023) Unveiling the epigenomic mechanisms of acquired platinum-resistance in high-grade serous ovarian cancer. Int J Cancer 153(1):120–132
Liu J, Dang H, Wang XW (2018) The significance of intertumor and intratumor heterogeneity in liver cancer. Exp Mol Med 50, (1), e416
Olbrecht S, Busschaert P, Qian J, Vanderstichele A, Loverix L, Van Gorp T, Van Nieuwenhuysen E, Han S, Van den Broeck A, Coosemans A, Van Rompuy AS, Lambrechts D, Vergote I (2021) High-grade serous tubo-ovarian cancer refined with single-cell RNA sequencing: specific cell subtypes influence survival and determine molecular subtype classification. Genome Med 13(1):111
Zhu JW, Charkhchi P, Akbari MR (2022) Potential clinical utility of liquid biopsies in ovarian cancer. Mol Cancer 21(1):114
Zellmer VR (2014) Evolving concepts of tumor heterogeneity. Cell bioscience 4:1–8
Polyak K (2014) Tumor heterogeneity confounds and illuminates: a case for Darwinian tumor evolution. Nat Med 20(4):344–346
Cordani M, Dando I, Ambrosini G, Gonzalez-Menendez P (2024) Signaling, cancer cell plasticity, and intratumor heterogeneity. Cell Commun Signal 22(1):255
Ciriello G, Magnani L, Aitken SJ, Akkari L, Behjati S, Hanahan D, Landau DA, Lopez-Bigas N, Lupianez DG, Marine JC, Martin-Villalba A, Natoli G, Obenauf AC, Oricchio E, Scaffidi P, Sottoriva A, Swarbrick A, Tonon G, Vanharanta S, Zuber J (2024) Cancer Evolution: A Multifaceted Affair. Cancer Discov 14(1):36–48
de Visser KE, Joyce JA (2023) The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 41(3):374–403
Ferri-Borgogno S, Zhu Y, Sheng J, Burks JK, Gomez JA, Wong KK, Wong STC, Mok SC (2023) Spatial Transcriptomics Depict Ligand-Receptor Cross-talk Heterogeneity at the Tumor-Stroma Interface in Long-Term Ovarian Cancer Survivors. Cancer Res 83(9):1503–1516
Stur E, Corvigno S, Xu M, Chen K, Tan Y, Lee S, Liu J, Ricco E, Kraushaar D, Castro P, Zhang J, Sood AK (2022) Spatially resolved transcriptomics of high-grade serous ovarian carcinoma. iScience 25(3):103923
Khatib S, Pomyen Y, Dang H, Wang XW (2020) Understanding the Cause and Consequence of Tumor Heterogeneity. Trends Cancer 6(4):267–271
Iacobuzio-Donahue CA, Litchfield K, Swanton C (2020) Intratumor heterogeneity reflects clinical disease course. Nat Cancer 1(1):3–6
Hegde PS, Chen DS (2020) Top 10 Challenges in Cancer Immunotherapy. Immunity 52(1):17–35
Ma L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DA, Forgues M, Candia J, Hernandez MO, Kelly M, Zhao Y, Tran B, Hernandez JM, Davis JL, Kleiner DE, Wood BJ, Greten TF, Wang XW (2021) Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol 75(6):1397–1408
Veskimae K, Scaravilli M, Niininen W, Karvonen H, Jaatinen S, Nykter M, Visakorpi T, Maenpaa J, Ungureanu D, Staff S (2018) Expression Analysis of Platinum Sensitive and Resistant Epithelial Ovarian Cancer Patient Samples Reveals New Candidates for Targeted Therapies. Transl Oncol 11(5):1160–1170
Lahnemann D, Koster J, Szczurek E, McCarthy DJ, Hicks SC, Robinson MD, Vallejos CA, Campbell KR, Beerenwinkel N, Mahfouz A, Pinello L, Skums P, Stamatakis A, Attolini CS, Aparicio S, Baaijens J, Balvert M, Barbanson B, Cappuccio A, Corleone G, Dutilh BE, Florescu M, Guryev V, Holmer R, Jahn K, Lobo TJ, Keizer EM, Khatri I, Kielbasa SM, Korbel JO, Kozlov AM, Kuo TH, Lelieveldt BPF, Mandoiu II, Marioni JC, Marschall T, Molder F, Niknejad A, Raczkowska A, Reinders M, Ridder J, Saliba AE, Somarakis A, Stegle O, Theis FJ, Yang H, Zelikovsky A, McHardy AC, Raphael BJ, Shah SP, Schonhuth (2020) A., Eleven grand challenges in single-cell data science. Genome Biol 21, (1), 31
Yu L, Cao Y, Yang JYH, Yang P (2022) Benchmarking clustering algorithms on estimating the number of cell types from single-cell RNA-sequencing data. Genome Biol 23(1):49
Zhang L, Cascio S, Mellors JW, Buckanovich RJ, Osmanbeyoglu HU (2024) Single-cell analysis reveals the stromal dynamics and tumor-specific characteristics in the microenvironment of ovarian cancer. Commun Biol 7(1):20
Wang T, Xu C, Zhang Z, Wu H, Li X, Zhang Y, Deng N, Dang N, Tang G, Yang X, Shi B, Li Z, Li L, Ye K (2022) Cellular heterogeneity and transcriptomic profiles during intrahepatic cholangiocarcinoma initiation and progression. Hepatology 76(5):1302–1317
Hussain A, Voisin V, Poon S, Karamboulas C, Bui NHB, Meens J, Dmytryshyn J, Ho VW, Tang KH, Paterson J, Clarke BA, Bernardini MQ, Bader GD, Neel BG, Ailles LE (2020) Distinct fibroblast functional states drive clinical outcomes in ovarian cancer and are regulated by TCF21. J Exp Med 217, (8).
Mao X, Xu J, Wang W, Liang C, Hua J, Liu J, Zhang B, Meng Q, Yu X, Shi S (2021) Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer 20(1):131
Cai J, Tang H, Xu L, Wang X, Yang C, Ruan S, Guo J, Hu S, Wang Z (2012) Fibroblasts in omentum activated by tumor cells promote ovarian cancer growth, adhesion and invasiveness. Carcinogenesis 33(1):20–29
Xun Z, Ding X, Zhang Y, Zhang B, Lai S, Zou D, Zheng J, Chen G, Su B, Han L, Ye Y (2023) Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis. Nat Commun 14(1):933
Feig C, Jones JO, Kraman M, Wells RJ, Deonarine A, Chan DS, Connell CM, Roberts EW, Zhao Q, Caballero OL, Teichmann SA, Janowitz T, Jodrell DI, Tuveson DA, Fearon DT (2013) Targeting CXCL12 from FAP-expressing carcinoma-associated fibroblasts synergizes with anti-PD-L1 immunotherapy in pancreatic cancer. Proc Natl Acad Sci U S A 110(50):20212–20217
Zhang Z, Yu Y, Zhang Z, Li D, Liang Z, Wang L, Chen Y, Liang Y, Niu H (2023) Cancer-associated fibroblasts-derived CXCL12 enhances immune escape of bladder cancer through inhibiting P62-mediated autophagic degradation of PDL1. J Exp Clin Cancer Res 42(1):316
Li X, Bu W, Meng L, Liu X, Wang S, Jiang L, Ren M, Fan Y, Sun H (2019) CXCL12/CXCR4 pathway orchestrates CSC-like properties by CAF recruited tumor associated macrophage in OSCC. Exp Cell Res 378(2):131–138
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd;, Hao Y, Stoeckius M, Smibert P, Satija R Comprehensive Integration of Single-Cell Data. Cell 2019, 177, (7), 1888–1902 e21.
Lun AT, McCarthy DJ, Marioni JC (2016) A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Res 5:2122
Kiselev VY, Kirschner K, Schaub MT, Andrews T, Yiu A, Chandra T, Natarajan KN, Reik W, Barahona M, Green AR, Hemberg M (2017) SC3: consensus clustering of single-cell RNA-seq data. Nat Methods 14(5):483–486
Tosches MA, Naumann YTM (2018) Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles. Science 360(6391):881–888
Zhang K, Erkan EP, Jamalzadeh S, Dai J, Andersson N, Kaipio K, Lamminen T, Mansuri N, Huhtinen K, Carpen O, Hietanen S, Oikkonen J, Hynninen J, Virtanen A, Hakkinen A, Hautaniemi S, Vaharautio A (2022) Longitudinal single-cell RNA-seq analysis reveals stress-promoted chemoresistance in metastatic ovarian cancer. Sci Adv 8(8):eabm1831
Yu G, Li F, Qin Y, Bo X, Wu Y, Wang S (2010) GOSemSim: an R package for measuring semantic similarity among GO terms and gene products. Bioinformatics 26(7):976–978
Aran D, Hu Z, Butte AJ (2017) xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 18(1):220
Disclaimer/Publisher’s Note The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content

The authors declare no competing interests.

SupplementaryMaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Decoding the effect of fibroblast heterogeneity on prognosis and drug resistance in high-grade serous ovarian cancer through tumor evolution analysis

Status:

Version 1

Abstract

Figures

1. Introduction

2. Results

2.2 Inter- and intratumoral heterogeneity in high-grade serous ovarian cancer

2.3 Tumor functional clonality is associated with patient prognosis and TME composition

2.4 Identification of a subpopulation of CXCL12-positive fibroblasts associated with chemotherapy resistance and poor prognosis

2.5 Identification of a gene set affecting prognosis and drug resistance in HGSOC

2.6 Spatial distribution characteristics of CXCL12-positive fibroblasts affect the clinical prognostic results of HGSOC

3. Discussion

4. Materials and Methods

4.1 Preprocessing of single-cell RNA-seq data

4.2 Identification of tumor cells

4.3 Calculation of epithelial cell scores

4.4 Tumor evolution analysis

4.4.1 Selection of representative tumor cells

4.4.2 Clustering and refinement

4.4.3 Hierarchical evolutionary analysis

4.4 Cell-type enrichment analysis

4.5 Weighted correlation network analysis in cell subclusters and acquisition of fibroblast-related Genes

4.6 Similarity analysis among cell subpopulations

4.7 Clustering of bulk RNA data samples

4.8 Mapping of single-cell and Bulk RNA subpopulations

4.9 Friends analysis

4.10 Spatial transcriptome data analysis

4.11 Immune infiltration analysis

4.12 Cell transfection

4.13 qRT-PCR

4.14 Western blot

4.15 CCK8 assay

5. Conclusions

Declarations

Supplementary Materials

Conflicts of Interest

Funding

Author contributions

Acknowledgments

Data Availability Statement

References

Additional Declarations

Supplementary Files

Status:

Version 1