Single-cell landscapes demonstrate tumor heterogeneity due to CNVs
We performed a systemic pipeline to analyze the single-cell RNA-seq profiles of EOC (Fig. 1A). The 2794 single cells were classified into 11 clusters using UMAP (20) (Fig. 1B) and were annotated as seven cell types based on the cell markers (16) (Fig. 1C; see Supplementary Method). Contact the patient for clinical information, we found that benign, low-grade and high-grade serous EOC cells in cluster 1 were clearly distributed in three regions. Moreover, multiple cell clusters originated from several patients were found in malignant epithelial, fibroblasts and mesenchymal cells. All these suggest that the tumor heterogeneity exists in EOC, which was probably caused by copy number variations (CNVs) (33); this hypothesis was confirmed by InferCNV. Clusters 0, 6 and 9 defined as epithelial cells have a global accumulation of CNVs (Supplementary Figure S2A), suggesting the malignancy of these subsets. Despite tumor heterogeneity, high-grade serous ovarian cancers possessed deletions in chromosomes 13 and amplifications in chromosomes 8. Furthermore, immune cells were the most abundant component in addition to malignant epithelial cells (Fig. 1D), where the immune cells content of benign, low-grade and high-grade serous ovarian cancers increased sequentially. The overall proportion of immune cells in each clinical subtype of tumor may reflect the different efficiencies of immune depletion (16), which may be a plausible explanation for why high-grade serous ovarian cancers is more suitable for immune-targeted therapy (34). Immune cells and malignant epithelial cells originating from different patients were distributed in separate clusters, while this phenomenon was not observed in stromal cells (Supplementary Figure S2B and C). These results further confirmed that tumor heterogeneity and stromal cell unity.
Cell type-specific reprograming of metabolic pathway and hallmarks
Tumorigenesis is a progress of malignant transformation of cell physiological functions including metabolism, signal transduction and cell proliferation and apoptosis (35). Therefore, mining the physiological characteristics of TME components can reflect the reprogramming of tumor ecology. The pathway activity evaluation algorithm (36) (see Supplementary Method) was used to calculate the pathway score of each cell type. There were 69 metabolic pathways which were found to be significantly associated with different EOC cell clusters. More than 60 pathways were highly activated (activity score > 1 and p-value < 0.01) in at least one cell type (Supplementary Figure S3A). Similarly, there were 50 cancer hallmark genes sets which were highly activated in at least one cell type (Supplementary Figure S3B). As an example, the different activity of oxidative phosphorylation (OXPHOS), glycolysis, TCA cycle, and hypoxia in cell types was supported by direct comparison of the average expression of these pathways in each cell type (one-way ANOVA p < 0.01, Supplementary Figure S4). Compared with other cells, macrophages cells were significantly associating with more positive activation of metabolic pathways and cancer hallmarks (Fig. 1E). These up-regulated gene sets were included in several parts of physiological functions, such as OXPHOS, glycolysis and EMT. Macrophages cells have a global up-regulation of metabolic pathway activity compared to malignant epithelial cells, which is not observed in other cancer (36). Meanwhile, the heterogeneity of malignant epithelial cells reflected in metabolic activity and tumor proliferation and invasion mechanisms, where type-1 epithelial cell was hyper-metabolically active and higher EMT active, type-2 was meso-metabolically active, and type-3 was hypo-metabolically active (Supplementary Figure S3 and S5). Although the metabolic pathway activity was strongly influenced by cell subsets (Fig. 1F), there was a high correlation in oncogenic function between the three malignant epithelial cells (Fig. 1G). This pattern also was found in CAFs and mesenchymal cells. These results suggest that heterogeneity of tumor metabolic reprogramming does not largely influence oncogenic function.
Contrasting TME physiological features observed at single-cell and bulk resolution will contribute to the correction of EOC prior knowledge. The pathway activity scores were calculated based on microarray data for bulk tumor samples obtained from GEO (GSE26712) and was compared with the results detected by scRNA-seq profiles. We found that only 26 metabolic pathways and 14 hallmarks were up-regulated in tumor samples compared to normal samples (p-value < 0.01, Supplementary Figure S6 A and B). The correlation of physiological pathway activity between bulk tumor samples and malignant epithelial cells was not significant (Pearson’s R = 0.044, p-value = 0.718, Fig. 1H, Supplementary Figure S6 C and D), suggesting expression specificity between cell types was masked by the fact that bulk data measure the average expression levels over a mixture of multiple cell types. Moreover, the variation in the distribution of pathway activity for single-cell was greater than bulk data (average standard deviation of pathway activity = 0.075 for single cells compared to 0.015 for bulk tumor, Fig. 1I). Taken together, the programming of cellular-specific physiological functional activity can be better reflected at single-cell resolution than bulk resolution.
Intratumor heterogeneity of EOC depend on energy metabolism
For malignant epithelial cells, which have tumor heterogeneity reflected in molecular mechanisms, we performed re-cluster analysis and eight clusters were identified (Fig. 2A) that predominantly from seven patients, indicating that individual differences in tumor cells is evident. The phenomenon may be caused by the patient’s unique physiological environment including location-specific and nutrient supply (36). There is also an undefined cluster of cells overlapping with type 2 epithelial cells, suggesting that hypo-metabolically active malignant cells are widely distributed and physiologically similar between patients.
Reprogramming of cellular physiological mechanisms is largely determined by variation and space environment. To explore what aspects of cellular physiological mechanisms are affected by environmental factors. GSEA (23) was used to identify metabolic and cancer hallmarks enriched in genes that represents intratumor heterogeneity influenced by environment, mutation, etc. (Supplementary Figure S7). We found that the OXPHOS and NFKB signaling pathway had the highest enrichment scores in most tumor clusters (Fig. 2B and C). Similarly, glycolytic and hypoxic signaling pathways also have potential contribution to intratumor heterogeneity across these tumor clusters. These results suggest that energy metabolism and cellular resistance to apoptosis were major contributors to intratumor heterogeneity of malignant epithelial cells. The coefficient of variation (CV), standard deviation (SD), and information entropy was then used as weights for each gene accounting for variation between malignant cells to exclude potential prejudice for formula selection and found consistent result in identification of functional pathways (Supplementary Figure S8).
Furthermore, the association between energy metabolism and environmental factors such as oxygen supply was performed. The hypoxic signals (hypoxia inducible factor, HIF) that indirectly reflect the oxygen content of cells were evaluated (Fig. 2D). We found the high correlation between glycolysis and hypoxia (Pearson’ R = 0.51 for epithelial cells, Fig. 2E), which is consistent with previous studies demonstrating that hypoxia activates glycolysis (37, 38). Moreover, the activity of OXPHOS was significantly correlated with glycolysis (Pearson’ R = 0.62) and the response to hypoxia (Pearson’ R = 0.43). Notably, we did not detect strong correlation between OXPHOS and glycolysis (Pearson’ R = -0.13) and between OXPHOS and hypoxia (Pearson’ R = -0.13) using bulk RNA-seq data of TCGA (Fig. 2F). In addition to malignant epithelial cells, we had also identified the relationship between OXPHOS, glycolysis and hypoxia in stromal cells (Fig. 2G), indicating the particular relationship between energy metabolism and hypoxia at single-cell resolution. For cells with low activity in energy metabolism, the up-regulated genes were significantly enriched in microtubule-based movement (Fig. 2H), suggesting the subpopulation of malignant epithelial cells promote tumor progression (39). Taken together, energy metabolism major contributed to intratumor heterogeneity.
Reprogramming of energy metabolism during EMT
We revealed that the higher metabolism activity of malignant cell had higher EMT score above. However, the relationship between energy metabolism, an important contributor to tumor heterogeneity, and tumor metastasis remains unknown. Further, the pseudo-developmental trajectory of epithelial cells, mesenchymal were simulated using monocle (Fig. 3A). The EMT-cells were grouped into three branches (defined as “B1”, “B2” and “B3”) with different states. We found that the three epithelial cell subtypes were mainly concentrated in state1 and state3, while mesenchymal cells and CAFs were concentrated in state2 (Fig. 3B). The EMT pathway activity gradually increased with pseudo-time development and reached the highest at state 2 (Fig. 3C). The phenomenon was also observed in independent dataset (GSE147082) (Supplementary Figure S9A-E), indicating the pathology of epithelial cells transformation to mesenchymal cells. We found that the score of OXPHOS, glycolysis and hypoxia increased from branch 1 to branch 2 during EMT (Fig. 3D-F), which is consistent with previous studies suggesting that reduced energy metabolism provide environment for EMT and hypoxic conditions promote the EMT (40, 41). Several branching-dependent genes were identified from branch1 to branch2 and most of them related to EMT pathway and energy metabolism (Fig. 3G and Supplementary Figure S9I). Moreover, cell invasiveness and stemness follow similar trends to energy metabolism during EMT (Fig. 3H and I). The OXPHOS and glycolysis were higher correlated with cell invasiveness (Pearson’s = 0.59 for OXPHOS, Pearson’s = 0.56 for glycolysis) and cell stemness (Pearson’s = 0.37 for OXPHOS, Pearson’s = 0.26 for glycolysis, Fig. 3J). The high correlation between hypoxia and EMT (Pearson’s = 0.65) confirms the reliability of this study. Besides, 810 EMT-related genes (586 positively and 224 negatively; |R|>0.2 (25)) was identified based on pseudo time. GO functional enrichment results showed that positive EMT-related genes act mainly in the extracellular matrix and cell adhesion that may provide the requirements for cell transfer (Fig. 3K) and negative EMT-related genes are associated with the encoding of membrane proteins (Fig. 3L). These results indicate that hypoxia induces shift in the function of epithelial cells to promote tumor cell migration.
Prognostic marker EGR1 drives EMT in EOC
The driver genes of EMT related to pseudo-time were identified based on molecular regulatory network. In transcriptional regulatory network, we identified 41 positive-TFs and 10 negative-TFs from EMT-related genes (Fig. 4A), which targeted 50 positive-target genes and two negative-target genes (Fig. 4B, Supplementary Figure S9G). We found that FOSB, RUNX1, as important drivers of tumor proliferation and metastasis (42, 43), were significantly associated with OS of EOC patients (Fig. 4C and D) and up-regulated during EMT (Fig. 5F-H). It has been demonstrated that EGR1 can regulate angiogenesis to promote tumor cell metastasis (44), which up-regulated in the EMT process and regulated target-gene HSPG2 (Fig. 5F-J, Supplementary Figure S9F-H) to enhanced EMT in this study. Intriguingly, high expression of both EGR1 and HSPG2 were associated with poorer OS in EOC patients (Fig. 5E, Supplementary Figure S10A). Besides, the HMGA1 down-regulated expression of CDH1 encoding E-cadherin (Supplementary Figure S10B). Furthermore, the FOSB, RUNX1, HSPG2, and EGR1 were significantly different expressed between invasive epithelial (iE) and invasive mesenchymal (iM) samples defined from TCGA-OV data (45, 46) (Fig. 5K and Supplementary Figure S11), suggesting that several critical genes driving the EMT process can only be identified at single-cell resolution.
Knockdown EGR1 inhibit ovarian cancer cell migration and invasion
The TF EGR1 regulated the most EMT-related genes (Fig. 4B) and exhibited significant prognostic efficiency of patient survival (Fig. 4E). Positive correlation was also observed between EMT activity and EGR1 expression (Fig. 4C, Fig. 5A and Supplementary Figure S9J). Furthermore, experimental analysis by employing two ovarian cancer cell lines SKOV3 and A2780. Western blotting was used to show the Knockout efficiency of EGR1 in SKOV3 and A2780 (Fig. 5B). We evaluated the effect of EGR1 on the migration and invasion of EOC cells by transwell and wound healing (Fig. 5C-D). Migration and invasion capacity of SKOV3 and A2780 cells were significantly decreased compared to the control group, suggesting that the recovery of epithelial phenotype inhibits cell migration. On the other hand, we test the EMT-relative proteins in SKOV3 and A2780 after treatment with EGR1 inhibition (Fig. 5E). Our data showed that silencing of EGR1 significantly increased the expression of E-cadherin. While the protein of N-cadherin and Slug were down-regulated in EGR1 group. These results indicate that silence of EGR1 suppressed migration and invasion in ovarian cancer cells.
Suppressed EGR1 Decreases Expression of FAS and HSPG2
Next, the regulating mechanism of EGR1 was verified. The expression of two EGR1 target genes, FAS and HSPG2, were positively correlating with EGR1 based on the single cell expression profile and TCGA bulk dataset (Fig. 6A-D). Within the TME of EOC, cells with high FAS and HSPG2 expression level exhibiting increased EMT activity (Fig. 6E-F). To further confirm the transcriptional relationship between EGR1 and downstream target genes, we explored Chromatin Immunoprecipitation Sequencing (ChIP-seq) experiment dataset from ENCODE. Enriched EGR1 sequencing read peaks were found in the FAS and HSPG2 transcriptional factor binding region across different bio-samples, indicating directly transcriptional relationship between EGR1 and its targets (Fig. 6G-H). The experiment demonstrated that FAS and HSPG2 expression at the protein level was markedly down-regulated in EGR1-si group in SKOV3 and A2780 cells compared with NC-si group (Fig. 6I-J). These results indicated that EGR1 decreases the expression of FAS and HSPG2 and associating with EMT in the TME of EOC.
Cell-cell communication reveals the pro-metastatic effect of CAFs
The 42 signaling pathways was identified by CellChat (29) that enable cell-cell communication between nine cell subsets in tumor tissue, and the amount of signal emission in the three epithelial cell clusters correlated with metabolic pathway activity (Supplementary Figure S12A). Notably, CAFs cells exhibited most communication relationships with other cell types and have the strongest interactions with mesenchymal and macrophage cells (Fig. 7A). To explore the functions of these signaling pathways, we clustered the above 42 signaling pathways into three categories related to growth factors and metastatic determined by their functions (Fig. 7B). Further, we collected CCL and CXCL family ligand-receptor interaction pairs associated with tumor metastasis (Fig. 7C-E). CAFs highly expressed the ligand CXCL12 and their actionable receptors including CXCR4 and ACKR3 (Fig. 7C). We found that CXCR4 was highly expressed in immune cells such as dendritic, T cells, and macrophages and thus CAFs induced immune cells to control immune infiltration (47). CXCL12 on CAFs could also combined with CXCR4 in epithelial cells to induce EMT (8). CCL5 and CXCL2, secreted by CAFs, could interacted with ACKR1 that is highly expressed on the surface of endothelial cells (Fig. 7C and F). This signaling pathway is associated with tumor metastasis and invasion (48, 49), which also identified in EOC and we pinpointed their derived from CAFs. The role of CAFs in tumor proliferation and metastasis revealed that EOC patients have individual clinical phenotypes due to the cellular composition of the tumor tissue. The fraction of the nine cell subsets in the OV samples from TCGA was calculated using cibersortX (Fig. 7G). Correlating the cell fraction data with clinical information, we noted that part of stage III patients had higher fraction of CAFs but lower fraction of immune cells. The accumulation of CAFs was also associated with poorer OS, suggesting that CAFs plays an important role in tumor progression (Fig. 7H and I). However, the accumulation of mesenchymal and dendritic cells was however associated with better OS, suggesting that they may be protective factors for EOC (Fig. 7J, Supplementary Figure S12B). Taken together, CAFs were highlighted to promote EOC progression through the CXCL signaling pathway.