Single-Cell RNA-Sequencing Atlas Reveals the Tumor Microenvironment of Metastatic High-Grade Serous Ovarian Carcinoma

Ovarian cancer is the most common and lethal gynecological tumor in women worldwide. High-grade serous ovarian carcinoma (HGSOC) is one of the histological subtypes of epithelial ovarian cancer, accounting for 70%. It often occurs at later stages associated with a more fatal prognosis than endometrioid carcinomas (EC), another subtype of epithelial ovarian cancer. However, the molecular mechanism and biology underlying the metastatic HGSOC (HG_M) immunophenotype remain poorly elusive. Here, we performed single-cell RNA sequencing analyses of primary HGSOC (HG_P) samples, metastatic HGSOC (HG_M) samples, and endometrioid carcinomas (EC) samples. We found that ERBB2 and HOXB-AS3 genes were more amplified in metastasis tumors than in primary tumors. Notably, high-grade serous ovarian cancer metastases are accompanied by dysregulation of multiple pathways. Malignant cells with features of epithelial-mesenchymal transition (EMT) affiliated with poor overall survival were identified. In addition, cancer-associated fibroblasts with EMT-program were enriched in HG_M, participating in angiogenesis and immune regulation, such as IL6/STAT3 pathway activity. Compared with ECs, HGSOCs exhibited higher T cell infiltration. PRDM1 regulators may be involved in T cell exhaustion in ovarian cancer. The CX3CR1_macro subpopulation may play a role in promoting tumor progression in ovarian cancer with high expression of BAG3, IL1B, and VEGFA. The new targets we discovered in this study will be useful in the future, providing guidance on the treatment of ovarian cancer.

Introduction

Epithelial ovarian cancer usually occurs at an advanced stage and is the most common cause of death from gynecological cancer (1). High-grade serous ovarian carcinoma (HGSOC) is a common histological subtype of epithelial ovarian carcinoma, for 70 to 80%, while the endometrioid carcinoma (EC) subtype accounts for 10% (2). HGSOC is associated with a more fatal prognosis and frequent recurrences, with more than 85% of women with this type having a 10-year mortality rate of 70%, whereas endometrioid carcinoma is thought to originate in endometriosis, a histologic type that tends to have a better prognosis. Our knowledge of the molecular etiology and clinical pathology of HGSOC has greatly improved, and recent therapy has advanced (3–6). However, most patients are diagnosed at a late stage, when cancer has already metastasized, and the diagnosis results in 5-year survival of 30%. Accordingly, the development of effective therapies for metastatic ovarian cancer is urgently needed. To do so, we need to comprehensively characterize the cellular heterogeneity and define transcriptional features within the tumor microenvironment.

Genomic analysis of HGSOC revealed a mutation in the tumor suppressor gene TP53, which is also seen in endometrioid carcinoma (3), promoting ovarian cancer metastasis and chemoresistance, and defecting in homologous recombination (HR) DNA repair, which contributes to the somatic BRCA mutation (7). The Cancer Genome Atlas (TCGA) project has classified HGSOC into four transcriptional subtypes: ‘differentiated’, ‘immunoreactive’, ‘proliferative’, and ‘mesenchymal’ (8, 9).

Recently, a scRNA-seq study of six metastatic omental tumors that derived from primary HGSOCs unraveled the genetic signatures of immune cell subsets within the tumor microenvironment and identified NR1H2⁺ IRF8⁺ and CD274⁺ macrophage clusters, which were suggested with an anti-tumor response (10). Another scRNA-seq study revealed that the inhibition of the JAK/STAT pathway has potential anti-tumor activity (11). The heterogeneity of tumor cells and different immune cell types within the TME play a paramount role in shaping tumor behavior (12–14). Therefore, characterizing the complex interplay between tumor cells and immune cell phenotype within HGSOC will be beneficial to find the critical factors of ovarian carcinogenesis, metastasis, and targeted treatment.

In this study, we conducted single-cell RNA sequencing of five primary high-grade serous carcinomas samples (HG_P), three metastases from HGSOC to the peritoneum (HG_M), one normal ovarian sample, and two primary Endometrioid (EC_P) samples. By comparing HG_M with HG_P and EC_P, we comprehensively characterized the heterogeneity of tumor cells and immune cells in ovarian cancer lesions, as well as the dynamic changes in cell-type composition and intercellular interactions, providing new insights into the biological basis of the development of HGSOCs and ECs.

Results

Single-Cell Transcriptomic Profiling of the Cellular Heterogeneity of the HGSOCs

Droplet-based single RNA-seq (10X Genomic) was performed on a total of eight samples from five treatment-naive patients. For parallel analyses, the public scRNA-data of 1 HG_P and 2 HG_M samples (15) from the same patient were downloaded (Figure 1A, Tables S1, S2). After quality filtering, approximately 0.68 billion unique molecular identifiers (UMIs) were collected from 55802 cells with >250 genes detected. Of these cells, 28,571(51.2%) cells were from HG_P, 8925 (16%) cells were from HG_M, and 12751(22.9%) cells were from EC_P. All high-quality cells were used to perform canonical correlation analysis (CCA) and identify anchors or mutual nearest neighbors (MNNs). Then, we integrated all cells, conducting unsupervised graph-based clustering (16).

FIGURE 1

Figure 1 Overview of TME in primary HGSOCs, metastatic HGSOCs, and ECs. (A) Workflow of the samples collected and the data analysis strategy. (B) Cell populations identified. The UMAP projection of 55802 single cells from HG_P (n=5), HG_M (N=3), EC_P (n=2), HG_Nor (n=1) samples shows the 10 main clusters with annotation. Each dot corresponds to a single cell, colored according to cell type. (C) Canonical cell markers are used to identify the clusters. (D) Barplots of the cell type and cancer subtypes for all 11 tumors. (E) The cancer subtypes proportion for each pathological group.

By Uniform Manifold Approximation and Projection (UMAP) with the resolution of 1.1, we identified 10 major lineages (epithelial cells, B cells, NK cells, T cells, plasma cells, fibroblasts, mesenchymal stem cells (MSCs), endothelial cells, neutrophils, and myeloid cells) (Figures 1B, C, S1A). The cell types were mainly assigned based on canonical cell markers and functional categories according to significantly differential genes expressed from different clusters (17, 18). One remaining cluster was labeled as “N” because we could not confidently recognize this cell type (10). The respective proportion of each cell type was varied and significantly differed between tumors (Figure 1D, proportion test, df=10, p-value < 2.2e-16). Interestingly, the boxplot exhibited that the medians in B cells, NK/T cells, and myeloid cells were higher in HG_P than HG_M and EC_P, whereas plasma cells were more enriched in EC_P (Figure S1B).

Previously, The Cancer Genome Atlas (TCGA) had stratified HGSOC into four molecular subtypes. We wondered if the inter-patient variability among tumors were consistent with these subtypes. To this end, we assigned molecular subtypes to our samples with the consensusOV (9) classifier (Version 1.14.0) (Figure 1D, S2A). We found that all four subtypes were well presented in each ovarian lesion. EC_P expressed the DIF signature slightly over HG_P and HG_M, while HG_P was comparable with HG_M, and EC_P presented the lowest IMR signature, supporting low immune-cell infiltration (Figure 1E). What’s more, certain subtypes tended to be consistent with specific cell types (Figure S2B). Epithelial cells highly expressed the differential (DIF) signature and lesser expressed proliferative (PRO) signature. The mesenchymal (MES) signature was strongly expressed by the fibroblasts and MSCs cells, while the immunoreactive (IMR) signature mainly consisted of myeloid, T cells, and NK cells. Notably, more fibroblasts were classified as proliferative subtypes in our data, suggesting that there are more fibroblasts with relatively high tumor purity. Based on this result, we wanted to know which genes made the most contribution to this classification. Thus, we extracted markers that were used for classification and calculated their average expression in each cell type (Figure S3A). As expected, proliferative markers such as MARCKSL1, STMN1, UCHL1, MFAP2, TRO, etc. are indeed expressed at higher levels in fibroblasts than others. Likewise, we plotted these markers in the heatmap in each cell of fibroblasts (Figure S3B) and we found PRO-markers like MARCKSL1 and STMN1, especially for MARCKSL1, which are highly expressed in most cells. And these markers contribute the most to the classifications of proliferative-subtype, while other PRO-related markers, UCHL1 and MFAP2, are not expressed significantly per cell in fibroblasts though their average expression is higher than other cell types. These results illustrated the importance of fibroblasts in cancer progression and indicated that a subset of fibroblasts in our data is cancer-associated fibroblasts (CAFs).

In addition, we can assess TCGA-subtypes at the patient level from our single-cell data by calculating the average of each gene per sample. Validated with the dataset from a previous study (19) (Figure S4A-4C), we had a high degree of confidence to infer the subtypes of HG1_P, HG2_P, HG4_M2, and EC2_P (Figure S4D, 4E). Their subtypes are likely to be IMR, IMR, MES, and DIF. Overall, our findings fully illustrated the difference in subtype classification between bulk samples and single-cell data, where single cells can more accurately describe the TCGA-subtype and characterize the heterogeneity of ovarian cancer.

Distinguish Worse Survival Cells from Cancer Epithelial Cells of HGSOCs

Based on the expression of PAX8 and CD24 (20, 21), we found that they were mainly expressed in epithelial cells as well as the subtype DIF (Figures S2A, S5A), which we termed “cancer epithelial cells”. In addition, we used T cells and myeloid cells from the normal sample (HG2_nor) as controls by inferring chromosomal copy number alterations (InferCNV, Version 1.2.1) (11) to confirm this (Figure S6). As shown below, we found that CNV trends in the same chromosomal region of different primary patients were distinct, whereas they were approximately consistent in the same patient (Figures 2A, B, S5B). Obviously, the genes generally mutated in ovarian cancer were more amplified in chr17 in HG3_M than in the primary tumors on both sides, such as the ERBB2 and HOXB-AS3 genes, which are generally mutated in ovarian cancer (8). These results demonstrated both intertumoral heterogeneity between patients and consistency within the same patient lesion.

FIGURE 2

Figure 2 Copy number profiles, intertumoral heterogeneity, and EMT signature subpopulations are identified. (A) The chromosomal landscape of copy number for 13,634 epithelial cells of seven primary tumors; amplification (red) and deletions (blue). (B) The chromosomal landscape of copy number for 2849 epithelial cells of metastatic tumors and primary tumors of the HG3 patient (L_HG3_P means the primary tumor from the left ovary in the HG3 patient; R_HG3_P means the primary tumor from the right ovary in HG3 patient). (C) The UMAP projection of 17,551 epithelial cells from 10 tumors of six patients (indicated by labels and colors) reveals tumor-specific clusters. (D) Differentially expressed genes of the top 10 genes (rows) that are differentially expressed in each cluster (columns). (E) Differentially expressed genes between Scissor⁺ cells and all other cells in HGSOCs, each point represents a gene. Red: significant genes; Black: NS genes. avg_logFC: log 2 fold-change of the average expression between the two groups. ((log-FC > 0.25, FDR <0.05) (F) Enrichment of significant genes related to Reactome and Hallmark pathways. (G) Kaplan-Meier plot shows that high expression of EMT signature has shorter overall survival in ovarian cancer. The high and low patients are split by the mean expression of the EMT-related gene set.

Re-running UMAP analysis on these cells, a total of eight sub-clusters were identified, of which revealed interpatient tumor-specific clusters (22, 23) (Figure 2C). Conversely, re-clustering the subsets like T cells and fibroblast cells without integration, we found that ovarian lesions from the same batch clustered together (Figure S5C). Differentially expressed genes in each cluster, interpatient GSVA, and cell cycle analysis were also shown (Figures 2D and S5D, 5E). We then used a Scissor algorithm (Version 2.0.0) to classify cells associated with worse survival (Scissor⁺ cells) from cancer epithelial cells, based on GDC TCGA bulk RNA-seq expression and clinical phenotype (24). As described previously, Scissor can quantify the similarity between single-cell data and bulk data through measurements, for example, Pearson correlations, and then it optimizes a regression model on the correlation matrix with the sample clinical phenotype. After that, it will provide feedback on three cell types, such as cells associated with worse prognosis (Scissor⁺ cells), cells related to better prognosis (Scissor^- cells), and cells that have no relationship with prognosis (Background cells).

In HGSOCs, we found that Scissor⁺ epithelial cells mainly accumulated in patients with metastasis (Figure S7A). To distinguish Scissor⁺ cells from Scissor^- cells and Background cells (All other cells), we compared their gene expression (Figure 2E, S7C, cut off: avg_log-FC > 0.25, FDR <0.05). Interestingly, high expression of EMT-related genes like GAS1, DCN, COL1A1, MGP, etc. in HGSOCs Scissor⁺ cells derived predominantly from patients with tumors metastasized (Figures 2E, S7B). Genes such as STMN1, CCND1, TUBA1B, and TUBB significantly expressed in Scissor⁺ cells of ECs were associated with the cell cycle (Figure S7E), which is a barometer of epithelial tumor cells proliferation (25). Furthermore, functional enrichment analysis of significantly expressed genes in Scissor⁺ cells also confirmed these (Table 1, Figures 2F, S7D), and survival analyses revealed that high levels of EMT and cell cycle signature were significantly related to poor overall survival in the Ovarian Cohort (Figures 2G, S7F). The observation of HGSOCs was consistent with the previous report that EMT is involved in increasing the invasion and metastasis of epithelial tumors (26–29). Therefore, targeting epigenetic regulation of EMT is a potentially powerful approach to inhibit the migration and invasiveness of HGSOCs.

TABLE 1

Table 1 Functional enrichment analysis based on the upregulated genes in HGSOCs Scissor⁺ cells versus All other cells (Scissor^- and Background cells).

Dynamic Trajectory Analysis During the Progression of HGSOCs

In the recent past, omentum metastasis has been reported (10). But in general, there are still few studies on the genetic dynamics of high-grade serous ovarian cancer metastasis, especially on the transcription factors involved in tumor progression. Based on the Monocle2 method (Version 2.21.1), pseudo-time reconstruction of epithelial cells was performed to infer the progression path of HGSOC (Figure 3A, and S8A). Macroscopically, the number of metastatic epithelial cells increased along the trajectory at the later stage. Besides, we also estimated the RNA velocities of every single cell by distinguishing un-spliced and spliced mRNAs, a function provided by velocyto (30) package. According to the direction of movement of each cell, the process of metastasis of HGSOCs can be clearly detected (Figure S8D).

FIGURE 3

Figure 3 Trajectory reconstruction during metastatic HGSOCs. (A) Monocle2 infers the development of epithelial cells along with pseudo-time (from patients HG3 and HG4 respectively, L_HG3_P means the primary tumor from the left ovary). Pseudo-time legend from dark to bright indicates cancer progression from the early to late stage. (B) Genswitches deduces the genes switch between cell states (left: L_HG3; right: HG4). (C, D) The heatmap displays the dynamic gene expression profiles during metastasis of ovarian cancer (from patients HG3 and HG4 respectively). The color key from blue to red indicates relative expression levels from low to light. The top annotated GO and KEGG terms in each cluster are shown. (E, F) Top 100 differentially expressed transcription factor genes (TFs; left) and the expression of specific TFs are on view along with the pseudo-time curve in (right). (G) Overexpression of proliferation and metastasis-related genes predicts poor prognosis in HGSOCs.

In particular, the dynamic gene expression profiles during the development of tumors were extracted. (Figures 3C, D). Interestingly, we found that the molecular mechanisms involved in metastatic HGSOCs were consistent whether the primary HGSOC was on the left or right side (Figures 3C, S8E). Genes related to the immune response were significantly decreased, whereas the genes related to DNA replication, cell cycle, epithelial cell proliferation, oxidative phosphorylation, and TCA cycle were significantly increased (Figures 3C, D). These results also suggested that inhibitors based on poly (ADP-ribose) polymerase (PARP), which aids in stopping the regeneration of cancer cells, such as Olaparib and Rucaparib, may be helpful for the treatment of metastatic HGSOCs (31–33). Meanwhile, the transcriptional factors (TFs) related to immune regulation, such as ARID5A, NFKB1, RORA, and ZNF683, were gradually downregulated along with the trajectory differentiation process (Figure 3E). And these immune-related TFs were scattered in the clusters of epithelial cells as well as other cell types (Figure S8B). Conversely, some well-known factors related to tumor growth promotion, such as HMGA1, GTF3A, PHF19, CENPX, and MBD2, were upregulated. Zingg emphasized that loss of cilia accelerates melanoma metastasis in benign cells by enhancing Wnt/β-Catenin Signaling (34). In our data, the expression of epithelial cilium movement markers, including FOXJ1, PRG, RFX2, and TMF1, although increased during the process, were mainly overexpressed in primary tumor cells (Figures 3A, F). This may indicate that disruption of cilia assembly leads to primary ovarian cancer that metastasizes to the peritoneum.

We also utilized GeneSwitches (35) (Version 0.1.0) to predict the genes that act as on/off switches between cell states in order during the tumor’s metastasis process (Figures 3B, S8C). Accordingly, overexpression of some acting on genes like FOSL2, NFIB, NFIC, and PROCR, which are associated with proliferation and metastasis, predicts poor prognosis in high-grade serous ovarian cancer (Figure 3G). Taken together, these results reveal dynamic gene expression profiles, highlighting several quintessential TFs and surface proteins that are dysregulated during ovarian cancer progression.

Cancer-Associated Fibroblasts with the EMT Program Enriched in Metastatic HGSOCs

Fibroblast is another vital biological cell type that synthesizes the extracellular matrix and collagen to maintain the structural integrity of connective tissue. In this study, 12,236 fibroblast cells were categorized into five distinct sub-clusters (Figures 4A, S9B, 9C). Fibro_1 cells were marked by STAR, an exclusive marker for ovarian stromal cells as previously reported (15). Fibro_2 cells were characterized by collagen (COL1A1, COL3A1) and cancer-associated fibroblast genes (CTHRC1, FAP). Fibro_3 cells expressed immunomodulatory (CFD, OGN) and tumor suppressor genes (CCDC80, PLA2G2A). Fibro_4 cells expressed growth factors (EGFR, IER2M KLF2). Fibro_5 cells were characterized by conserved and nuclear-enriched lncRNA (MALAT1, NEAT1), and MALAT1 modulates the expression of cell cycle-related genes in lung fibroblast and EMT-related genes in breast cancer (36, 37) (Figure 4C). The distribution of fibroblast sub-clusters in each tumor was varied (Figure 4B). In addition, comparisons based on hallmark gene sets of fibroblasts had been conducted. Fibroblasts from HG_M were more abundant in supporting tumor progress than HG_P, including angiogenesis, coagulation system, and EMT (Figure S9D). Fibroblasts from EC_P were more enriched in the structure formation of tissue than HG_P, such as myogenesis and adipogenesis, whereas inflammatory pathways including TNFA signaling via NFκB and inflammatory response more enriched in HG_P than EC_P (Figure S9E). Combined with the analysis of normal fallopian tubes (nFT) from previous studies (38, 39), we found that fibroblasts from HG_P were more abundant in supporting interferon-alpha/gamma response than HG_nor/nFT while estrogen responses early activity was more active in HG_nor/nFT than HG_P (Figure S9F).

FIGURE 4

Figure 4 Diversity of fibroblasts in HG_M. (A) The UMAP projection of 12,236 fibroblast cells of 11 samples from six patients (indicated by labels and colors). (B) Proportion and cell number of each fibroblast subtype in 11 samples. (C) Heatmap of marker genes expression. (D) Heatmap of functional gene sets. (E) GSVA analysis of differential pathways is scored per cell among five fibroblast subsets. (F) Active regulons in each fibroblast subsets.

Notably, Fibro_2 cells, the most enriched subtype in HG_M, accounting for 32%, expressed genes of cancer-associated fibroblasts (CAFs), angiogenesis, and collagen at a high level (Figures 4D, S9A). CAFs have been verified to promote tumor metastasis through upregulating genes like HSF1, which was involved in the pro-tumorigenic pathway (40). Furthermore, Hallmark pathway analysis also confirmed that Fibro_2 cells had more relevance with the pathways that sustain tumor growth, including angiogenesis, epithelial-mesenchymal transition (EMT), hypoxia, and PI3K/AKT/mTOR signaling (Figure 4E). Compared with HG_M and EC_P, Fibro_3 cells (19%) and Fibro_5 cells (3%) accounted for less proportion in HG_P. Intriguingly, Fibro_3 cells and Fibro_5 cells were both consistent with the characteristic of “antigen-presenting CAFs” as previously discerned (41), owning to express genes like CD74 and human leukocyte antigen (Figure 4D).

With the help of the SCENIC tool, we identified regulons unique to each fibroblast sub-clusters (Figure 4F). For instance, the transcription factors of SOX4 and SRF underpinned Fibro_2, while Fibro_5 cells were characterized by STAT1/STAT3 and NFKB1/NFKB2. Strikingly, SOX4 has been proven as an important co-factor of SMAD3, controlling pro-metastatic gene transcription and shaping the cell response to TGF-β in different scenarios, thereby promoting tumorigenesis (42). While STAT3 and NF-κB are pro-inflammatory regulators and they form transcriptional complexes that positively regulate gene expression in oncogenic pathways (43).

Heterogeneity of Tumor-Infiltrating Lymphocytes in HGSOCs

Infiltration of T cells into tumors modifies the natural course of the disease and plays a critical role in cancer immunotherapy (44, 45). From the ovarian cancer lesion, we classified a total of 7967 T and NK cells into eight subtypes: CD4⁺ T cells (CD4 IL7R; CD3D⁺ CD4⁺), regulatory CD4⁺ T cells (Tregs FOXP3; CD4, FOXP3), CD8⁺ T cells (CD8 GZMK, CD8 GZMH; CD3D⁺ CD8⁺), NK cells (NK CD56, NK IL7R; NCAM1, GNLY, TYROBP, NKG7), NKT cells (CD3D, CD8A, FCGR3A, GNLY), and Innate Lymphoid Cells (ILCs; CD3D) (Figures 5A, C and S10A). The number of cells and the proportion of each subtype in each tumor were shown (Figure 5B). Notably, 1665 T/NK cells were obtained from HG_M, while 5394 T/NK cells were from HG_P. Using a dendrogram to group the tumors based on the average expression of T cell markers, we found that HG_M showed a similar pattern to HG_P, whereas EC_P emerged with low expression in CD4⁺ T and CD8⁺ T cells (Figure 5D). This observation was consistent with the boxplot shown, the lower proportion of CD4⁺ T cells and CD8⁺ T cells was detected in EC_P than in HG_P and HG_M (Figure S10B). The low tumor-suppressive status in EC_P suggested that T cell-based immunotherapy may be inefficient in EC_P.

FIGURE 5

Figure 5 Subpopulations of tumor-infiltrating lymphocytes (TILs) in HG_M. (A) The UMAP projection of 7967 TILs of 11 samples from six patients (indicated by labels and colors). (B) Proportion and cell number of each subtype in 11 samples. (C) Dot plot (left) and UMAP-plot (right) display canonical cell markers. (D) Hierarchical clustering heatmap groups the tumors between HG_P, HG_M, and EC_P. (E) Reconstruction trajectory of CD8⁺ T cells inferred by Monocle2 (color by subtypes, expression of signature genes, and pseudotime). (F) Heatmap of the functional gene sets in TILs. (G) Active Regulons in each TILs. (H) Overexpression of the PRDM1 gene predicts a worse prognosis in ovarian cancer. (I) Cumulative distribution of cytotoxic CD8⁺ T cells between HG_P, HG_M, and EC_P. The cytotoxic score is calculated based on the average expression of cytotoxic markers. P-value was calculated by a two-sided unpaired Kruskal-Wallis rank-sum test.

Among CD4⁺ T cells, we identified naive (CD4 IL7R; TCF7, CCR7, SELL, LEF1) and regulatory (Tregs FOXP3; IL2RA, FOXP3, IKZF2). The Tregs FOXP3 cells highly expressed inhibitory genes, including TIGIT, CTLA4, and ENTPD1, and they also relatively expressed high levels of costimulatory molecules CD28, TNFRSF14, ICOS, TNFRSF9, which stimulate the inhibitory activities (Figure 5F). The tumor-suppressive microenvironment mediated by Tregs is a significant obstacle to successful immunotherapy, suggesting that depletion of Treg cells, like immune checkpoint blockade of CTLA-4 or PD1/PDL1, could be a potentially effective immunotherapy for ovarian cancer (46).

Among CD8⁺ T cells, CD8 GZMK and CD8 GZMH T cells were characterized by relatively high cytotoxic genes granzyme K (GZMK) and granzyme H (GZMH), respectively. Meanwhile, these cells also positively expressed T cell exhaustion markers, including LAG3 and PDCD1, indicating that the CD8⁺ T cells are exhausted after initial activation in ovarian cancer (Figure 5F). In addition, we inferred the gene regulatory networks across the TILs subtypes by SCENIC method (Figure 5G). The regulon of PRDM1 was upregulated in CD8 GZMK and CD8 GZMH T cells, which is connected with terminal T cell differentiation and contributes to the maintenance of an early memory phenotype and cytokine poly-functionality in TILs after knockout (47). Consequently, we concluded that PRDM1 might be one of the factors contributing to the exhaustion of CD8 GZMK and CD8 GZMH T cells. Indeed, verified by bulk RNA-seq from the TCGA dataset, it was a slightly significant Spearman correlation between PRDM1 expression and immune exhausted infiltrate in ovarian cancer (Cor:0.41, FDR: 2e⁻¹², Figure S10C). Beyond that, higher expression of PRDM1 predicted a worse prognosis in ovarian cancer (Figure 5H).

Next, we performed pseudo-time trajectory analysis to explore the dynamic states and cell transitions of CD4⁺ IL7R to CD8⁺ T cells via Monocle2. In the developmental trajectory, CD4⁺ IL7R started as a root, and gradually evolved into CD8 GZMK and CD8 GZMH, presenting a binary branched structure in which one side was the end of exhausted T cells, and the other side was the end of cytotoxic T cells (Figure 5E). In HG_M, the proportion of CD8⁺ cytotoxic T cells was higher than in HG_P (Figure 5I). On the contrary, the percentage of exhausted CD8⁺ T cells in HG_P was more than that in HG_M (Figure S10D).

We noticed three sub-clusters expressing the NK cells marker: TYROBP, GNLY, and NKG7. NK CD56 cells were characterized by NCAM1 (CD56), NK IL7R cells were characterized by IL7R, and NKT cells were identified by the specific T-cell markers including CD3D and CD8A (Figures 5C, D). NKT cells strongly expressed the GZMB, GZMA, GZMH, and PRF1 genes, indicating that they promoted tumor cytotoxicity in ovarian cancer (Figure 5F). Generally speaking, NK CD56 and NKT cells were more enriched in HG_P and HG_M than EC_P (Figure S10E).

Trajectory Reconstruction of HGSOCs Revealed Monocyte-to-Macrophage Differentiation

Tumor-infiltrating myeloid cells (TIMs) are critical regulators in tumor progression, playing essential roles in modulating tumor inflammation and angiogenesis (48, 49). Altogether, 7265 myeloid cells were collected, revealing 23 subsets through the ROGUE statistic (Version 1.0) (50) purified the cell population (Figures S11A, 11B). Then, we identified four common linages (cDCs, monocytes, or macrophage and monocytes derived DC) based on canonical cell markers and they were well presented in each sample (Figures 6A, B, D, S11D, 11E). Moreover, a subset of myeloid cells expressed myeloid/T-cells markers simultaneously (CD3D_undefined), which was not discussed below.

FIGURE 6

Figure 6 Subpopulations of myeloid cells in HG_M. (A) The UMAP projection of 7265 myeloid cells of 11 samples from six patients (indicated by labels and colors). (B) Proportion and cell number of each myeloid subtype in 11 samples. (C) RNA velocity of each myeloid subtype. (D) The dot plot displays canonical cell markers. (E) Trajectory reconstruction of monocyte evolved into macrophages. (F) Dynamics gene expression profile during monocyte-to-macrophage terminal differentiation. (G, H) Biological processes enrichment analysis of module 2 and module 4. (I) Heatmap of significant genes in each subtype. (J) Activate regulons in each myeloid subtype.

Monocytes are the progenitors of monocytes-derived macrophages and contribute to the overall coordination of immunity (51). Correspondingly, Monocytes can be separated from macrophages based on phylogenetic reconstruction (Figure S11F). CD14_mono cells were characterized by CD14, SELL, and S100A8/9, representing classical monocytes and being recruited during inflammation (Figure 6I). They also highly expressed FCN1, a complement system protein that defends against infectious agents (52). CD16_mono cells were less abundant and represented non-classical monocytes with high expression of FCGR3A (CD16), CDKN1C, LST1, and low expression of CD14. Similarly, CD16_mono cells expressed FCN1 at high levels but were more enriched in HG_M (Figure S11C).

DC cells were classified according to their origin and typical genes. CLEC10A_mono-DC with lower abundance was characterized by CLEC10A and CD14 representing monocyte-derived dendritic cells. CD1C_DC cells represented a classic cDCs subset, with high expression of CD1C, CD1E, and CLEC10A.

Macrophage cells were characterized by tissue-resident and their pro-inflammatory or anti-inflammatory function. CCL2_macro represented early-stage macrophage with the expression of CCL2. CX3CR1_macro expressed genes involved in immune modulation of chemokines, such as CCL3, CCL4, and CXCL8. What’s more, CX3CR1_macro significantly secreted BAG3 (Figure 6I), a multifunctional protein, which can combine with a specific receptor IFITM2 to induce the release of factors that sustain the growth and metastasis of tumor (53). MMP9_macro expressed genes related to inflammatory chemokines (CXCL2, CXCL3, CXCL8) and genes like MMPs (MMP19, MMP9), which play an important role in tumor tissue remodeling. CCL18_macro expressed both M1 marker (CD68) and M2 markers (CCL18, GPNMB). By the way, CCL18 played a key role in recruiting immunosuppressive myeloid cells (54). Finally, LYVE1_macro represented resident macrophages with a high level of LYVE1 and FOLR2.

Early works have supported that macrophages can either originate from monocyte cells or tissue-resident macrophages (55–57). Accordingly, we employed the RNA velocity to explore the ovarian cancer lineage trajectories and we found that a small number of CD14_mono evolved toward MMP9_macro, while LYVE1_macro evolved toward CCL18_macro as well as toward MMP9_macro (Figure 6C). To further investigate the dynamic change of genes during the differentiation of CD14_mono into macrophages in the ovarian lesions, we extracted classical monocytes (CD14_mono) and related macrophages (MMP9_macro, CCL18_macro, CX3CR1_macro) to delineate monocyte-to-macrophage differentiation by trajectory development analysis (Figure 6E). During the trajectory, CD14_mono were progenitor cells for MMP9_macrophage and then further separated into CX3CR1_macro and CCL18_macro. Profiling of gene expression dynamics along the trajectory had been divided into four modules (Figure 6F). Genes like CX3CR1, CXCL8, CXCR4, CCL3L1, VEGFA, and IL1B in Module 2 increased during the evolution of the branch of CX3CR1, whereas, they decreased during the evolution of the branch of CCL18_macro. Vice versa, Genes in Module 4 like CCL18, CCL7, and CCL8 were increased in the branch of CCL18_macro and reduced in the branch of CX3CR1. To reveal the biological characteristics of these two branches, we performed gene functional enrichment analysis. In the branch of CX3CR1_macro, the corresponding high expressed genes tend to be associated with inflammatory response and positive regulation of cell migration (Figure 6G), which is consistent with the suggestion that the high VEGF, CXCL8⁺ IL1β⁺ TAMs with the features of inflammatory could promote tumor progression in ovarian cancer (54, 58). Parallel results, in the branch which evolved into CCL18_macro, their corresponding high expressed genes were associated with metabolic precursor and energy production (Figure 6H).

We further characterized the functions of nine subtypes explained above by comparing pathway activities (Figure S11G). Pathways involved in angiogenesis, EMT, TNFA signaling via NFKB, and hypoxia were upregulated in CX3CR1_macro, MMP9_macro, and CD14_mono. These results indicate the potential tumor-promoting feature of CX3CR1_macro derived from CD14_mono. Finally, we applied SCENIC to identify TFs underlying each phenotype (Figure 6J). Interestingly, some recruited macrophage phenotypes shared similar TFs expression patterns with monocytes and tissue-resident macrophages.

Intercellular Communication Networks

As mentioned above, we have obtained cancer epithelial cells with poor prognosis (Scissor⁺ Epithelial cells), and we wondered what function they are involved in with the crosstalk of the tumor microenvironment in ovarian lesions. To this end, the CellphoneDB repository (Version 2.1.4) (59) was used to predict putative intercellular interactions between Scissor⁺ Epithelial cells and other cell types based on ligand-receptor signaling. Interestingly, many significantly overexpressed molecular pairs were associated with immunosuppression and HG_M showed a similar pattern to HG_P (Figures 7B, C). Compared with EC_P, Scissor⁺ Epithelial cells in HG_M and HG_P had more outgoing interactions with other cell types. Of note, macrophage and fibroblasts connected with Scissor⁺ Epithelial more frequently in HG_M (Figure 7A). When Scissor⁺ Epithelial cells expressed a relatively high level of EGFR as receptors, their corresponding ligands, such as AREG, COPA, GRN, MIF, and TGFB1, were widely expressed in other cells. When Scissor⁺ Epithelial cells expressed genes related to angiogenesis (VEGFA), the interactions (VEGFA_FLT1, VEGFA_KDR) were slightly abundant in HG_M. In addition, when Scissor⁺ Epithelial expressed ACKR2 acting as a receptor for chemokines including CCL3/CCL4/CCL5, their interactions were more abundant in HG_P than that of the other two groups (Figure 7B). It has been reported that ACKR2 is a scavenger receptor for chemokines and its deficiency against metastasis (60). And, CCL5 is important for the recruitment and activation of lymphocytes (61), so we proposed that ACKR2_CCL5 may weaken the recruitment and activation of lymphocytes, contributing to the metastasis of primary high-grade serous carcinoma. Among the three groups, the MDK_LRP1 molecule pair between Scissor⁺ Epithelial cells and myeloid cells expressed significantly, while their mean expression level in EC_P was higher than HG_P and HG_M (Figure 7B). MDK can combine with its receptor LRP1, which is beneficial to tumor-infiltrating macrophages, promoting myeloid inhibitory cell differentiation (MDSCs) (62, 63). Thus, MDK targeted therapy should suggest an effective treatment for ovarian cancer.

FIGURE 7

Figure 7 The intricate intercellular interplay in HG_P, HG_M, and EC_P. (A) Circos plot shows the intercellular interactions in HG_P, HG_M, and EC_P. Each line represents an interaction where one end represents a ligand that is expressed in one cell type and the other end represents a receptor that is expressed in another cell type. The thickness of each line corresponds to the number of distinct interacting pairs. (B, C) Dot plot shows the means of the average expression levels and the possibility of occurrence in selective interaction pairs.

In general, these results revealed that the crosstalk between Scissor⁺ Epithelial cell and other cell types via diverse receptor-ligand signals may profoundly affect ovarian cancer development and metastasis.

Discussion

Our current study has comprehensively characterized the dynamic variation of gene profiles during tumor progression in HGSOCs, as well as the heterogeneity of tumor cells, fibroblast cells, and immunophenotype, and the intricate intercellular interactions across HG_P, HG_M, and EC_P. We have identified unique subpopulations such as Scissor⁺ Epithelial cells, CAFs–Fibro_2, CX3CR1/CCL18 macrophages, and GZMK/GZMH CD8⁺ T cells. Those results provide a new perspective on the tumor microenvironment of ovarian cancer.

Epithelial cells were the largest cluster of cells, composing ~31% of the cells analyzed. However, current single-cell studies of ovarian cancer cannot explicitly distinguish cells with specific clinical phenotypes (24). To this end, Scissor⁺ cells with poor prognosis were identified by the Scissor algorithm. Consistent with the previous conclusion that the EMT signature is a potential factor for tumor invasion and metastasis (26, 27), the genes related to the EMT signature, including MGP, GAS1, and JUN, were found in Scissor⁺ cells of metastatic HGSOCs lesions. Moreover, during HGSOCs progression, several signaling pathways such as the cell cycle, tumor cell proliferation, oxidative phosphorylation, and TCA cycle, which needed energy metabolism, were markedly enhanced, suggesting that poly (ADP-ribose) polymerase (PARP) inhibition may be a targeted strategy for the treatment of metastatic HGSOCs (31–33). Strikingly, in our data, regardless of whether the primary tumor was from the left ovary or the right, the pathway activities generated by the tumor progression were basically consistent, as well as their copy number alterations in one segment of chromosomes.

Since the fibroblasts were the second-largest cluster cells analyzed, we observed Fibro_2, a subtype of CAF expressing the EMT program was specifically enriched in HG_M. This result supports that CAFs contribute to the EMT in HG_M and subsequently promote metastasis (64). Beyond that, the IL6_JAK_STAT3 signal pathway was also more enriched in Fibro_2. Consistently, JAK/STAT inhibitor JSI-124 has been proven to have an anti-tumor property in HGSOC cell lines (11). Combination therapies with the EMT or JAK/STAT inhibitor may help in the treatment of HGSOCs.

T cells are the crucial players in cancer immunotherapy (44, 45). Olalekan has revealed that ovarian cancer with high infiltration of CD8⁺ TOX⁺ and CD4⁺ GNLY T cells may be a good indication for patients (10). In our data, we found low T cell infiltration in EC_P compared with HG_P and HG_M. Furthermore, GZMH CD8 cells and GZMK CD8 cells simultaneously presented cytotoxicity and cell exhaustion. Notably, HG_M showed the highest cytotoxic of CD8⁺ T cells while the highest exhaustion of CD8⁺ T cells was in HG_P. Through the gene regulatory networks analysis, we speculated that PRDM1 may be involved in CD8⁺ T cells exhaustion and predicted that its high expression in ovarian cancer was associated with poor prognosis.

It has been manifested that TAM can promote the formation of niches before metastasis by secreting specific cytokines (65). They also can regulate the Tregs and Th17 cells to create immunosuppression, thereby promoting invasion and metastasis of ovarian cancer (66). On the basis of these theories, we have identified a subtype of TAM (CX3CR1_macro) with abundant production of BAG3, which can combine with IFITM2, leading to tumor metastasis (53). This TAM subtype also owns a high level of VEGFA, CXCL8, and IL1β, similar to previously reported TAM induced from monocyte with factor-1a stabilization in solid ovarian cancer that promoted tumor inflammation and metastasis (54).

By investigating the signaling network of Scissor⁺ Epithelial cells - other cells communication, we detected several receptor-ligand complexes that should be paramount for ovarian cancer development. Compared with HG_P and EC_P, stronger angiogenesis and tumor cell proliferation of intercellular interactions in HG_M. Metastasis of high-grade serous carcinoma may be related to the interaction between high expression of ACKR2 chemokine receptor and cytokines such as CCL5. The apparent pair of MDK_LRP1 among the three groups suggested that inhibition of the MDK_LRP1 pair might be an effective therapeutic target for ovarian cancer to reduce myeloid inhibitory cell differentiation (MDSCs).

We note that there are several limitations to our study. First, the number of patients with metastatic HGSOCs is small. Second, the clonal relationship of T cells was not investigated while T cell receptor therapy is an alternate therapy with great potential for ovarian cancer treatment (67). Hence, enlarging the cohort of metastatic HGSOCs and conducting immune profiling of T cell receptors and spatial transcriptomic may help unravel molecular mechanisms and elucidate the roles of different immune cells in HG_M.

In conclusion, our data have shed light on the tumor microenvironment of metastatic high-grade serous ovarian cancer at the single-cell level. Several novel markers and the potential therapeutic target detected in this study could provide valuable guidance for future clinical treatment.

Materials and Methods

Patients and Samples Collection

In this study, a total of eight samples, including one metastatic HGSOC matched two primary HGSOCs, one pair of matched primary HGSOC and normal ovary, one primary HGSOC and two ECs were collected from The Third Affiliated Hospital of Guangzhou Medical University, which have been subjected to pathological diagnosis. The clinical information of these patients is summarized in Table S1.

Preparation of Single-Cell Suspension

Specimens collected from patients with ovarian cancer were minced into fragments (< 1 mm3) and digested with 0.25% trypsin-EDTA (GIBCO) and DNase I (Roche Diagnostics) for 30 min at 37°C with agitation. The dissociated cell suspension was filtered through 70 μm strainer (BD Falcon), washed with cold PBS, and centrifuged at 4°C, x400g for 5 min. The cell pellet was resuspended in serum-free DMEM for further use.

Droplet-Based scRNA-Seq Data Preprocessing

The Cell Ranger (Version 6.0.2) pipeline generated raw gene expression matrices with human reference genome GRCh38. After removing doublets with DoubletFinder (Version 2.0.3) in each sample individually, the remaining cells were imported into R software (Version 4.1.0) for subsequent analysis by the Seurat package (Version 3.2.3). Cells with > 200 genes detected, genes expressed >5 cells, and genes expression >0 in all cells were selected for further analysis. Low-quality cells were removed according to the following criteria: unique molecular identifiers (UMIs) <500; genes <250 or genes >11000; UMIs derived from the mitochondrial genome >25%. After quality control, the gene expression was normalized by NormalizeData function, and cellular sequencing depth was adjusted by the SCTransform method.

Multiple Datasets Integration

To integrate multiple datasets across three conditions, we used the integration methods described at https://satijalab.org/seurat/v3.0/integration.html and https://hbctraining.github.io/scRNA-seq/lessons/06_SC_SCT_and_integration.html. The Seurat package (version 3.0.0) was used to assemble multiple distinct scRNA-seq datasets into an integrated and unbatched dataset. In brief, we used Sctransform to regress out confounding factors: number of UMIs, percentage of mitochondrial RNA, estimating the variance of the raw filtered data, and identifying the 3000 most variable genes. After that, we performed canonical correlation analysis (CCA) and then “integrated” the conditions to overlay cells that were similar or had a “common set of biological features” between groups.

Identification of Major Cell Types and Their Phenotype

Differential gene expression analysis was performed for clusters generated at various resolutions by both the Wilcoxon rank-sum test of Single-cell Transcriptomics (MAST) using the FindMarkers function. Annotation of the resulting clusters to cell types was based on the expression of marker genes along with the cell types assigned by SingleR packages (version ‘1.6.1’).

TCGA Subtype Classification

Single cell subtypes were classified by the consensusOV (version 1.14.0) package, with get.consensus.subtypes function, which also returns random forest probabilities for each subtype. The core principle of consensusOV is that it: 1. standardizes genes in each dataset to the same mean and variance, 2. computes binary gene pairs based on the standardized expression values.

CNV Estimation and Identification of Malignant Cells

We used an approach described previously to infer CNVs from the scRNA-seq data. Its R code was provided online. (https://github.com/broadinstitute/inferCNV) We set the cut off=0.1, denoise=TRUE, HSMM=TRUE, and hclust_method=‘ward.D2’. Immune cells from normal samples were considered as putative nonmalignant cells as control, and their CNV estimates were used to define a baseline. All epithelial cells from the ovarian tumor sample were used as input.

Distinguish Phenotype-Associated Cells

To link cells with specific phenotypes, we used the Scissor algorithm, a novel R package (Version 2.0.0) to identify the populations of the single-cell data associated with given phenotypes. (https://github.com/sunduanchen/Scissor)

Scissor integrates phenotype-associated bulk expression data and single-cell data by quantifying the similarity between every single cell and bulk sample. To identify relevant subpopulations, it then optimizes a regression model on the correlation matrix with the sample phenotype. The core formula of Scissor is as follows:

$$\frac{\text{min}}{\beta }-\frac{1}{n}l\left(\beta \right)+\lambda \left\{\alpha {\parallel \beta \parallel }_1+\frac{1-\alpha }{2}{\beta }^{T}L\beta \right\}$$

where L is a symmetric normalized Laplacian matrix; the tuning parameter λ controls the overall strength of the penalty, and α balances the amount of regularization between smoothing and sparsity. The phenotype-related cell subsets of interest are selected using the non-zero coefficient β solved for by the optimization model described above.

Statistical test: The scissors algorithm incorporates a reliability test to rule out false associations between identified cell subsets and bulk phenotypes. This statistical test can determine whether the inferred phenotype-cell association is reliable (P<0.05) or a false positive. First, it performs k-fold cross-validation (CV) on the correlation matrix S and estimates the cell coefficients in Scissor using only the training set. The predictive performance of the trained Scissor model is evaluated on the test set, and an averaged evaluation metric is obtained as the actual test statistic. Second, the bulk sample labels are randomized multiple times to break the original bulk phenotype-genotype relationship. Predictive performance at the random level was quantified by performing the same Scissor analysis and CV assessment using each permuted batch of data, obtaining the background distribution of the corresponding assessment measures. Finally, the actual test statistic calculated in the raw data is compared to the background distribution value. The reliability significance test p-value is the number of permutations-based test statistics above (or below) the actual test statistic divided by the number of permutations. In this study, the evaluation measures used in reliability significance tests were mean squared error (MSE) for linear regression (smaller is better), area under the ROC curve (AUC) for classification (higher is better), and agreement Sex index (C-index) for Cox regression (higher is better).

In actual operation, we set the family=“ cox”, alpha=0.077 for the HG group, and alpha=0.0265 for the EC group to select the phenotype-associated cell subpopulations by a Cox regression model. The number of the Scissor selected cells should not exceed 20% of the total cells in the single-cell data.

SCENIC Analysis

The regulons and TF activity (AUC) for each cell were calculated with the SCENIC (version 1.2.4) pipeline with motif collection version mc9nr, using per cell type with raw count matrices as input. We used GRNBoost (in Python) instead of GENIE3 to detect positive and negative associations for a bigger dataset. The function of exportsForGRNBoost was used to generate a gene expression matrix and TF list in special formats for GRNBoost to load.

Trajectory Inference Analysis

Trajectory analysis was performed using Monocle 2 (version 2.20.0), We assessed the significant gene by the differential gene expression analysis, and DEGs between the clusters were applied for dimension reduction using the reduceDimension function. Genes that changed along with the pseudotime were calculated and visualized with the plot_pseudotime_heatmap and the genes were clustered into subgroups according to the gene expression patterns. To identify the genes that separate cells into branches, the branch expression analysis modeling (BEAM) analysis was performed and genes resulting from the BEAM analysis with a q-value < 10−4 were separated into groups and visualized with the plot_genes_branched_heatmap function.

Estimations of RNA Velocities by Velocyto Package

In order to smoothly assess spliced and un-spliced mRNAs, we needed to convert the bam file to loom file by the function of run10x, provided by velocyto.py. Next, we merged multiple loom files by the function of loompy.combinem and then loaded the merged loom file into R software to combine analysis with the Seurat package.

Ordering the Gene Expression During Cell State Transitions

The genes, including the human surface proteins and transcription factors (TFs), act as on/off switches between cell states and are discovered by GenSwitches R packages (Version 0.1.0). The workflow of GeneSwitches is as follows:

1. Binarize gene expression data into 1(on) or 0(off).

2. Fit logistic regression on the binary states of gene expression and estimate switching time.

3. According to the default Settings, the poorly fitting genes are filtered and specific genes are extracted for plotting.

Gene Set Variation Analysis (GSVA)

Pathway analyses were performed on the 50 hallmark pathways described in the molecular signature database. We also evaluated the activity of 65 specific KEGG pathway activities from the Canonical pathway KEGG subset. To assign pathway activity estimates to individual cells, we applied GSVA (Version ‘1.34.0’) with standard settings.

Enrichment Analysis of Marker Genes

Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differential genes was implemented by the clusterProfiler (version 4.0.2) package. Reactome and Hallmark pathways analysis on differential genes was implemented by Metascape web-based portal (68).

Cell Cycle Scoring Assign

Each cell in the epithelial subpopulation was assigned a cell cycle-related score based on the gene expression of the G2/M and S phases. The “CellCycleScoring” function of the Seurat package was used to calculate the cell cycle score and store G2/M and S phase scores into data objects to predict the cell state.

TCGA Database

The cohort of TCGA ovarian cancer data was downloaded from https://xenabrowser.net/datapages/. The gene expression matrix and clinical phenotypes of HGSOCs and ECs were assessed using different datasets (GDC TCGA Ovarian Cancer (OV) & TCGA Endometrioid Cancer (UCEC)).

Survival Analysis

The Kaplan-Meier method evaluated ovarian cancer survival utilizing the KM plotter database (https://kmplot.com/analysis/).

We set split patients by ‘Auto select best cutoff’, which chooses the best performing threshold as cut-off among all possibilities between the lower and upper quartiles. All datasets provided by the KM plotter were taken into consideration for analysis.

Definition of Cytotoxicity and Exhaustion Scores

Cytotoxicity and exhaustion scores were calculated by the average expression of the genes from the predefined gene sets in CD8⁺ T cells of each group. To implement this method, the AddModuleScore function of the Seurat package was applied as previously described. We used eight cytotoxicity associated genes (GZMA, GZMB, GZMH, GZMH, IFNG, NKG7, PRF1, GNLY) and seven exhaustion associated genes (LAG3, TIGIT, PDCD1, HAVCR2, CTLA4, ENTPD1, BTLA) to define cytotoxicity and exhaustion scores

Correlation Analysis Between Immune Infiltration and Single Gene Expression

Spearman correlation between mRNA expression and exhausted infiltrate in ovarian cancer was performed by immune cell abundance pattern of Gene Set Cancer Analysis (GSCA) (69).(http://bioinfo.life.hust.edu.cn/GSCA/#/).

Statistics

Statistical analysis was carried out using R and Bioconductor.

Cell-Cell Communication Network

To investigate potential interactions between different cell types in the ovarian cancer tumor microenvironment, cell-cell communication analysis was performed as described previously by the CellPhoneDB Python package (Version 2.1.4), a publicly available repository of curated receptors and ligands and their interactions. Prediction of enriched receptor-ligand pairs between two cell types was derived from the expression of a receptor by one cell type and the expression of the corresponding ligand by another cell type. By default, only ligands and receptors expressed in at least 10% of cells in a given cell cluster were taken into consideration.

Pairwise comparisons between selective cell types were performed. We randomly permuted the cluster labels of all cells 1000 times to determine the mean of average ligand and receptor expression levels of the interactions, generating a null distribution for each receptor-ligand pair. A p-value for the likelihood of cell-type specificity of the corresponding receptor-ligand complex was obtained by calculating the proportion of the means as high as or above the actual mean. Then, we can select biologically relevant interactions.

Data Availability Statement

Single-cell RNA data reported in this paper have been deposited (PRJCA009148) in the Genome Sequence Archive in the BIG Data Center under accession number HRA002345 that are publicly accessible at https://ngdc.cncb.ac.cn/gsa-human/. All relevant data are available from the corresponding author upon reasonable request.

Ethics Statement

The studies involving human participants were reviewed and approved by The Third Affiliated Hospital of Guangzhou Medical University (No: 2017-075). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

XS and DS directed and oversaw the project. YT, DZ, and TC collected ovarian tissue samples. YD performed comprehensive bioinformatics analyses. YD, YB, and DS wrote the manuscript with all the other authors participating in the discussion, data interpretation, and manuscript editing. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Natural Science Foundation of Guangdong Province (no. 2020A1515010082), National Natural Science Foundation of China (NSFC, Project 81773012), National Science Foundation of China (grant 82170242 and 81570454), Science and Technology Planning Project of Guangzhou (202102010003), The Project for Key Medicine Discipline Construction of Guangzhou Municipality (no.2021-2023-17), and China Postdoctoral Science Foundation (2019M662858).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank all the participants and their families for supporting this study. We thank The Third Affiliated Hospital of Guangzhou Medical University in China for providing samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.923194/full#supplementary-material

References

1. Lheureux S, Gourley C, Vergote I, Oza AM. Epithelial Ovarian Cancer. Lancet (2019) 393(10177):1240–53. doi: 10.1016/s0140-6736(18)32552-2

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Stewart C, Ralyea C, Lockwood S. Ovarian Cancer: An Integrated Review. Semin Oncol Nurs (2019) 35(2):151–6. doi: 10.1016/j.soncn.2019.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Jayson GC, Kohn EC, Kitchener HC, Ledermann JA. Ovarian Cancer. Lancet (2014) 384(9951):1376–88. doi: 10.1016/s0140-6736(13)62146-7

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Bowtell DD, Bohm S, Ahmed AA, Aspuria PJ, Bast RC Jr., Beral V, et alRethinking Ovarian Cancer Ii: Reducing Mortality From High-Grade Serous Ovarian Cancer. Nat Rev Cancer (2015) 15(11):668–79. doi: 10.1038/nrc4019

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lheureux S, Cristea MC, Bruce JP, Garg S, Cabanero M, Mantia-Smaldone G, et al. Adavosertib Plus Gemcitabine for Platinum-Resistant or Platinum-Refractory Recurrent Ovarian Cancer: A Double-Blind, Randomised, Placebo-Controlled, Phase 2 Trial. Lancet (2021) 397(10271):281–92. doi: 10.1016/S0140-6736(20)32554-X

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Shih AJ, Menzin A, Whyte J, Lovecchio J, Liew A, Khalili H, et al. Identification of Grade and Origin Specific Cell Populations in Serous Epithelial Ovarian Cancer by Single Cell Rna-Seq. PloS One (2018) 13(11):e0206785. doi: 10.1371/journal.pone.0206785

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hao Q, Li J, Zhang Q, Xu F, Xie B, Lu H, et al. Single-Cell Transcriptomes Reveal Heterogeneity of High-Grade Serous Ovarian Carcinoma. Clin Transl Med (2021) 11(8):e500. doi: 10.1002/ctm2.500

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cancer Genome Atlas Research N. Integrated Genomic Analyses of Ovarian Carcinoma. Nature (2011) 474(7353):609–15. doi: 10.1038/nature10166

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Chen GM, Kannan L, Geistlinger L, Kofia V, Safikhani Z, Gendoo DMA, et al. Consensus on Molecular Subtypes of High-Grade Serous Ovarian Carcinoma. Clin Cancer Res (2018) 24(20):5037–47. doi: 10.1158/1078-0432.CCR-18-0784

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Olalekan S, Xie B, Back R, Eckart H, Basu A. Characterizing the Tumor Microenvironment of Metastatic Ovarian Cancer by Single-Cell Transcriptomics. Cell Rep (2021) 35(8):109165. doi: 10.1016/j.celrep.2021.109165

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Izar B, Tirosh I, Stover EH, Wakiro I, Cuoco MS, Alter I, et al. A Single-Cell Landscape of High-Grade Serous Ovarian Cancer. Nat Med (2020) 26(8):1271–9. doi: 10.1038/s41591-020-0926-0

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Sun Y, Wu L, Zhong Y, Zhou K, Hou Y, Wang Z, et al. Single-Cell Landscape of the Ecosystem in Early-Relapse Hepatocellular Carcinoma. Cell (2021) 184(2):404–21 e16. doi: 10.1016/j.cell.2020.11.041

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, et al. Single-Cell RNA Landscape of Intratumoral Heterogeneity and Immunosuppressive Microenvironment in Advanced Osteosarcoma. Nat Commun (2020) 11(1):6322. doi: 10.1038/s41467-020-20059-6

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ahmed N, Escalona R, Leung D, Chan E, Kannourakis G. Tumour Microenvironment and Metabolic Plasticity in Cancer and Cancer Stem Cells: Perspectives on Metabolic and Immune Regulatory Signatures in Chemoresistant Ovarian Cancer Stem Cells. Semin Cancer Biol (2018) 53:265–81. doi: 10.1016/j.semcancer.2018.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Qian J, Olbrecht S, Boeckx B, Vos H, Laoui D, Etlioglu E, et al. A Pan-Cancer Blueprint of the Heterogeneous Tumor Microenvironment Revealed by Single-Cell Profiling. Cell Res (2020) 30(9):745–62. doi: 10.1038/s41422-020-0355-0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive Integration of Single-Cell Data. Cell (2019) 177(7):1888–902 e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-Based Analysis of Lung Single-Cell Sequencing Reveals a Transitional Profibrotic Macrophage. Nat Immunol (2019) 20(2):163–72. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Geistlinger L, Oh S, Ramos M, Schiffer L, LaRue RS, Henzler CM, et al. Multiomic Analysis of Subtype Evolution and Heterogeneity in High-Grade Serous Ovarian Carcinoma. Cancer Res (2020) 80(20):4335–45. doi: 10.1158/0008-5472.CAN-20-0521

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Winterhoff BJ, Maile M, Mitra AK, Sebe A, Bazzaro M, Geller MA, et al. Single Cell Sequencing Reveals Heterogeneity Within Ovarian Cancer Epithelium and Cancer Associated Stromal Cells. Gynecol Oncol (2017) 144(3):598–606. doi: 10.1016/j.ygyno.2017.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kristiansen G, Denkert C, Schlüns K, Dahl E, Pilarsky C, Hauptmann S. CD24 is Expressed in Ovarian Cancer and is a New Independent Prognostic Marker of Patient Survival. Am J Pathol (2002) 161(4):1215–21. doi: 10.1016/s0002-9440(10)64398-2

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chen YP, Yin JH, Li WF, Li HJ, Chen DP, Zhang CJ, et al. Single-Cell Transcriptomics Reveals Regulators Underlying Immune Cell Diversity and Immune Subtypes Associated With Prognosis in Nasopharyngeal Carcinoma. Cell Res (2020) 30(11):1024–42. doi: 10.1038/s41422-020-0374-x

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hornburg M, Desbois M, Lu S, Guan Y, Lo AA, Kaufman S, et al. Single-Cell Dissection of Cellular Components and Interactions Shaping the Tumor Immune Phenotypes in Ovarian Cancer. Cancer Cell (2021) 39(7):928–44 e6. doi: 10.1016/j.ccell.2021.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sun D, Guan X, Moran AE, Wu LY, Qian DZ, Schedin P, et al. Identifying Phenotype-Associated Subpopulations by Integrating Bulk and Single-Cell Sequencing Data. Nat Biotechnol (2021) 40(4):527–538. doi: 10.1038/s41587-021-01091-3

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Tian T, Guo T, Zhen W, Zou J, Li F. Bet Degrader Inhibits Tumor Progression and Stem-Like Cell Growth via Wnt/B-Catenin Signaling Repression in Glioma Cells. Cell Death Dis (2020) 11(10):1–13. doi: 10.1038/s41419-020-03117-1

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Nieto MA, Huang RY, Jackson RA, Thiery JP. Emt: 2016. Cell (2016) 166(1):21–45. doi: 10.1016/j.cell.2016.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Bakir B, Chiarella AM, Pitarresi JR, Rustgi AK. EMT, MET, Plasticity, and Tumor Metastasis. Trends Cell Biol (2020) 30(10):764–76. doi: 10.1016/j.tcb.2020.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Xing X, Yang F, Huang Q, Guo H, Li J, Qiu M, et al. Decoding the Multicellular Ecosystem of Lung Adenocarcinoma Manifested as Pulmonary Subsolid Nodules by Single-Cell RNA Sequencing. Sci Adv (2021) 7(5):eabd9738. doi: 10.1126/sciadv.abd9738

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Majc B, Sever T, Zarić M, Breznik B, Turk B, Lah TT. Epithelial-To-Mesenchymal Transition as the Driver of Changing Carcinoma and Glioblastoma Microenvironment. Biochim Biophys Acta -Molecul Cell Res (2020) 1867(10):118782. doi: 10.1016/j.bbamcr.2020.118782

CrossRef Full Text | Google Scholar

30. La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, et al. RNA Velocity of Single Cells. Nature (2018) 560(7719):494–8. doi: 10.1038/s41586-018-0414-6

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Ledermann J, Harter P, Gourley C, Friedlander M, Vergote I, Rustin G, et al. Olaparib Maintenance Therapy in Platinum-Sensitive Relapsed Ovarian Cancer. New Engl J Med (2012) 366(15):1382–92. doi: 10.1056/NEJMoa1105535

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Moore K, Colombo N, Scambia G, Kim B-G, Oaknin A, Friedlander M, et al. Maintenance Olaparib in Patients With Newly Diagnosed Advanced Ovarian Cancer. New Engl J Med (2018) 379(26):2495–505. doi: 10.1056/NEJMoa1810858

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kristeleit R, Shapiro GI, Burris HA, Oza AM, LoRusso P, Patel MR, et al. A Phase I–II Study of the Oral PARP Inhibitor Rucaparib in Patients With Germline BRCA1/2-Mutated Ovarian Carcinoma or Other Solid Tumors. Clin Cancer Res (2017) 23(15):4095–106. doi: 10.1158/1078-0432.CCR-16-2796

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zingg D, Debbache J, Pena-Hernandez R, Antunes AT, Schaefer SM, Cheng PF, et al. Ezh2-Mediated Primary Cilium Deconstruction Drives Metastatic Melanoma Formation. Cancer Cell (2018) 34(1):69–84 e14. doi: 10.1016/j.ccell.2018.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Cao EY, Ouyang JF, Rackham OJJB. Geneswitches: Ordering Gene Expression and Functional Events in Single-Cell Experiments. Bioinformatics (2020) 36(10):3273–5. doi: 10.1093/bioinformatics/btaa099

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Tripathi V, Shen Z, Chakraborty A, Giri S, Freier SM, Wu X, et al. Long Noncoding RNA MALAT1 Controls Cell Cycle Progression by Regulating the Expression of Oncogenic Transcription Factor B-MYB. PloS Genet (2013) 9(3):e1003368. doi: 10.1371/journal.pgen.1003368

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Jadaliha M, Zong X, Malakar P, Ray T, Singh DK, Freier SM, et al. Functional and Prognostic Significance of Long non-Coding RNA MALAT1 as a Metastasis Driver in ER Negative Lymph Node Negative Breast Cancer. Oncotarget (2016) 7(26):40418. doi: 10.18632/oncotarget.9622

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Hu Z, Artibani M, Alsaadi A, Wietek N, Morotti M, Shi T, et al. The Repertoire of Serous Ovarian Cancer Non-Genetic Heterogeneity Revealed by Single-Cell Sequencing of Normal Fallopian Tube Epithelial Cells. Cancer Cell (2020) 37(2):226–42 e7. doi: 10.1016/j.ccell.2020.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ulrich ND, Shen YC, Ma Q, Yang K, Hannum DF, Jones A, et al. Cellular Heterogeneity of Human Fallopian Tubes in Normal and Hydrosalpinx Disease States Identified Using scRNA-Seq. Dev Cell (2022) 57(7):914–29 e7. doi: 10.1016/j.devcel.2022.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kalluri R. The Biology and Function of Fibroblasts in Cancer. Nat Rev Cancer (2016) 16(9):582–98. doi: 10.1038/nrc.2016.73

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Elyada E, Bolisetty M, Laise P, Flynn WF, Courtois ET, Burkhart RA, et al. Cross-Species Single-Cell Analysis of Pancreatic Ductal Adenocarcinoma Reveals Antigen-Presenting Cancer-Associated Fibroblasts. Cancer Discov (2019) 9(8):1102–23. doi: 10.1158/2159-8290.CD-19-0094

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Vervoort SJ, Lourenço AR, Tufegdzic Vidakovic A, Mocholi E, Sandoval JL, Rueda OM, et al. Sox4 can Redirect TGF-β-Mediated SMAD3-Transcriptional Output in a Context-Dependent Manner to Promote Tumorigenesis. Nucleic Acids Res (2018) 46(18):9578–90. doi: 10.1093/nar/gky755

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ji Z, He L, Regev A, Struhl K. Inflammatory Regulatory Network Mediated by the Joint Action of NF-KB, STAT3, and AP-1 Factors is Involved in Many Human Cancers. Proc Natl Acad Sci (2019) 116(19):9453–62. doi: 10.1073/pnas.1821068116

PubMed Abstract | CrossRef Full Text | Google Scholar

44. van der Leun AM, Thommen DS, Schumacher TN. CD8(+) T Cell States in Human Cancer: Insights From Single-Cell Analysis. Nat Rev Cancer (2020) 20(4):218–32. doi: 10.1038/s41568-019-0235-4

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Duong CP, Yong CS, Kershaw MH, Slaney CY, Darcy PK. Cancer Immunotherapy Utilizing Gene-Modified T Cells: From the Bench to the Clinic. Mol Immunol (2015) 67(2):46–57. doi: 10.1016/j.molimm.2014.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Liu X, Mo W, Ye J, Li L, Zhang Y, Hsueh EC, et al. Regulatory T Cells Trigger Effector T Cell DNA Damage and Senescence Caused by Metabolic Competition. Nat Commun (2018) 9(1):249. doi: 10.1038/s41467-017-02689-5

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Yoshikawa T, Wu Z, Inoue S, Kasuya H, Matsushita H, Takahashi Y, et al. Genetic Ablation of PRDM1 in Antitumor T Cells Enhances Therapeutic Efficacy of Adoptive Immunotherapy. Blood (2021). doi: 10.1182/blood.2021012714

CrossRef Full Text | Google Scholar

48. Cheng S, Li Z, Gao R, Xing B, Gao Y, Yang Y, et al. A Pan-Cancer Single-Cell Transcriptional Atlas of Tumor Infiltrating Myeloid Cells. Cell (2021) 184(3):792–809 e23. doi: 10.1016/j.cell.2021.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Engblom C, Pfirschke C, Pittet MJ. The Role of Myeloid Cells in Cancer Therapies. Nat Rev Cancer (2016) 16(7):447–62. doi: 10.1038/nrc.2016.54

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Liu B, Li C, Li Z, Wang D, Ren X, Zhang Z. An Entropy-Based Metric for Assessing the Purity of Single Cell Populations. Nat Commun (2020) 11(1):3155. doi: 10.1038/s41467-020-16904-3

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Murray PJ. Immune Regulation by Monocytes. Semin Immunol (2018) 35:12–8. doi: 10.1016/j.smim.2017.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Okuzaki D, Ota K, Takatsuki S-I, Akiyoshi Y, Naoi K, Yabuta N, et al. FCN1 (M-Ficolin), Which Directly Associates With Immunoglobulin G1, is a Molecular Target of Intravenous Immunoglobulin Therapy for Kawasaki Disease. Sci Rep (2017) 7(1):1–12. doi: 10.1038/s41598-017-11108-0

PubMed Abstract | CrossRef Full Text | Google Scholar

53. De Marco M, Basile A, Iorio V, Festa M, Falco A, Ranieri B, et al. Role of BAG3 in Cancer Progression: A Therapeutic Opportunity. Semin Cell Dev Biol (2018) 78:85–92. doi: 10.1016/j.semcdb.2017.08.049

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Paolini L, Adam C, Beauvillain C, Preisser L, Blanchard S, Pignon P, et al. Lactic Acidosis Together With GM-CSF and M-CSF Induces Human Macrophages Toward an Inflammatory Protumor Phenotype. Cancer Immunol Res (2020) 8(3):383–95. doi: 10.1158/2326-6066.CIR-18-0749

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BEP, Luu NT, et al. Resolving the Fibrotic Niche of Human Liver Cirrhosis at Single-Cell Level. Nature (2019) 575(7783):512–8. doi: 10.1038/s41586-019-1631-3

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Italiani P, Boraschi D. From Monocytes to M1/M2 Macrophages: Phenotypical vs. Functional Differentiation. Front Immunol (2014) 5:514. doi: 10.3389/fimmu.2014.00514

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Sharma A, Seow JJW, Dutertre CA, Pai R, Bleriot C, Mishra A, et al. Onco-Fetal Reprogramming of Endothelial Cells Drives Immunosuppressive Macrophages in Hepatocellular Carcinoma. Cell (2020) 183(2):377–94.e21. doi: 10.1016/j.cell.2020.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wang N, Wang S, Wang X, Zheng Y, Yang B, Zhang J, et al. Research Trends in Pharmacological Modulation of Tumor-Associated Macrophages. Clin Transl Med (2021) 11(1):e288. doi: 10.1002/ctm2.288

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhonedDB: Inferring Cell-Cell Communication From Combined Expression of Multi-Subunit Ligand-Receptor Complexes. Nat Protoc (2020) 15(4):1484–506. doi: 10.1038/s41596-020-0292-x

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Massara M, Bonavita O, Savino B, Caronni N, Mollica Poeta V, Sironi M, et al. ACKR2 in Hematopoietic Precursors as a Checkpoint of Neutrophil Release and Anti-Metastatic Activity. Nat Commun (2018) 9(1):676. doi: 10.1038/s41467-018-03080-8

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Tavares LP, Garcia CC, Goncalves APF, Kraemer LR, Melo EM, Oliveira FMS, et al. ACKR2 Contributes to Pulmonary Dysfunction by Shaping CCL5:CCR5-Dependent Recruitment of Lymphocytes During Influenza A Infection in Mice. Am J Physiol Lung Cell Mol Physiol (2020) 318(4):L655–L70. doi: 10.1152/ajplung.00134.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Filippou PS, Karagiannis GS, Constantinidou AJO. Midkine (MDK) Growth Factor: A Key Player in Cancer Progression and a Promising Therapeutic Target. Oncogene (2020) 39(10):2040–54. doi: 10.1038/s41388-019-1124-8

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Zhang Y, Zuo C, Liu L, Hu Y, Yang B, Qiu S, et al. Single-Cell RNA-Sequencing Atlas Reveals an MDK-Dependent Immunosuppressive Environment in ErbB Pathway-Mutated Gallbladder Cancer. J Hepatol (2021) 75(5):1128–41. doi: 10.1016/j.jhep.2021.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Erin N, Grahovac J, Brozovic A, Efferth TJDRU. Tumor Microenvironment and Epithelial Mesenchymal Transition as Targets to Overcome Tumor Multidrug Resistance. Drug Resist Update (2020) 53:100715. doi: 10.1016/j.drup.2020.100715

CrossRef Full Text | Google Scholar

65. Larionova I, Tuguzbaeva G, Ponomaryova A, Stakheyeva M, Cherdyntseva N, Pavlov V, et al. Tumor-Associated Macrophages in Human Breast, Colorectal, Lung, Ovarian and Prostate Cancers. Front Oncol (2020) 10:566511. doi: 10.3389/fonc.2020.566511

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Zhou J, Li X, Wu X, Zhang T, Zhu Q, Wang X, et al. Exosomes Released From Tumor-Associated Macrophages Transfer miRNAs That Induce a Treg/Th17 Cell Imbalance in Epithelial Ovarian Cancer. Cancer Immunol Res (2018) 6(12):1578–92. doi: 10.1158/2326-6066.CIR-17-0479

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Wu JW, Dand S, Doig L, Papenfuss AT, Scott CL, Ho G, et al. T-Cell Receptor Therapy in the Treatment of Ovarian Cancer: A Mini Review. Front Immunol (2021) 12:1141. doi: 10.3389/fimmu.2021.672502

CrossRef Full Text | Google Scholar

68. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat Commun (2019) 10(1):1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Liu C-J, Hu F-F, Xia M-X, Han L, Zhang Q, Guo A-YJB. GSCALite: A Web Server for Gene Set Cancer Analysis. Bioinformatics (2018) 34(21):3771–2. doi: 10.1093/bioinformatics/bty411

PubMed Abstract | CrossRef Full Text | Google Scholar