Skip to main content

ORIGINAL RESEARCH article

Front. Cell Dev. Biol., 21 December 2020
Sec. Molecular and Cellular Pathology
Volume 8 - 2020 | https://doi.org/10.3389/fcell.2020.614139

Immune Characterization of Ovarian Cancer Reveals New Cell Subtypes With Different Prognoses, Immune Risks, and Molecular Mechanisms

Shanshan Cong1 Qiuyan Guo1 Yan Cheng1 Yanan He1 Xibo Zhao1 Congcong Kong1 Shangwei Ning2* Guangmei Zhang1*
  • 1Department of Gynecology, The First Affiliated Hospital, Harbin Medical University, Harbin, China
  • 2College of Bioinformatics Science and Technology, Harbin Medical University, Harbin, China

Ovarian cancer (OV) is a considerable threat to the health of women due to its complex mechanisms and atypical symptoms. Various currently available treatments fail to substantially increase the survival rate of OV patients. The tumor microenvironment (TME) is gaining attention due to its role in tumorigenesis and tumor progression. This study mainly investigated the immune characteristics of OV by CIBERSORT and MCP-counter. We reclassified OV into four TME cell subtypes with different prognoses and evaluated the infiltration of the cells in each subtype. The immune risk of diverse subtypes was evaluated based on the immunoscore calculated by Cox regression analysis. The molecular mechanisms and hallmark pathways of the four subtypes were analyzed. The results indicate that the immune procancer cell subtype is associated with the worst prognosis, closely related to the high immune risk group, and characterized by low expression of checkpoints and MHC class I and II molecules, high expression of hypoxia-related genes, high enrichment of the EMT and hypoxia pathways, and low enrichment of the DNA repair and interferon α response pathways. This study contributes to the investigation of immune mechanisms and identifies more effective targets for immunotherapy of OV.

Introduction

Ovarian cancer is a highly lethal malignancy that ranks as the seventh leading cause of female cancer worldwide; there will be approximately 21,750 new cases of OV in 2020 (National Cancer Institute, 2020). Therefore, OV is a considerable threat to the health and lives of women. Current treatments for OV mainly involve a combination of surgery and chemotherapy based on platinum and paclitaxel. Studies of OV have validated a number of new therapies (targeted therapy and immunotherapy) (Lamichhane et al., 2017; Reislander et al., 2019). However, due to drug resistance and recurrence of OV, the mortality of OV remains higher than that of all gynecological malignancies, and the 5-year survival rate of OV has kept at approximately 45% (Ghisoni et al., 2019). Therefore, new perspectives and in-depth investigation of OV are needed.

The tumor microenvironment (TME) is the complex environment of a tumor that is mainly composed of the vasculature, cells, and extracellular matrix (Henke et al., 2019; Baghban et al., 2020). Fibroblasts, endothelial cells, and immune cells (lymphocytes, monocytes, macrophages, granulocytes, mast cells, etc.) are the major cellular components of the TME (Quail and Joyce, 2013; Guo and Deng, 2018; Peltanova et al., 2019); these cells are known as tumor microenvironment cells (TME cells). TME cells closely interact with the tumor cells that consequently influence the local immune response and progression of the tumors. Studies on TME in OV have identified multiple TME cells that are relevant to the prognosis of OV (Drakes and Stiff, 2018). However, the majority of publications on TME cells in OV are based on only a single or a few cell types and rarely integrate all TME cells to define their relationship with OV.

Therefore, we aimed to analyze the immune characteristics and roles of TME cells in OV in detail. An important step in our study was to link the RNA-seq data to TME cells by CIBERSORT (Drakes and Stiff, 2018) and MCP-counter (Becht et al., 2016). CIBERSORT is the most frequently cited immune cell infiltration analysis tool based on the principle of deconvolution linking the expression of 547 genes to 22 types of immune cells. MCP-counter includes 10 cell types, and fibroblasts and endothelial cells were used in our study because all other cell types in MCP-counter overlap with the cell types of CIBERSORT. Thus, our study used these two algorithms to estimate the proportion of TME cells to reclassify OV into four new subtypes; then, the TME cell infiltration pattern and immune risk of each subtype were assessed, and relevant molecular mechanisms were investigated.

Materials and Methods

Acquisition of Gene Profiles and Clinical Datasets of OV

The Cancer Genome Atlas (TCGA) RNA-seq data and clinical characteristics of OV were obtained from UCSC Xena1. The format of the TCGA RNA-seq data was fragments per kilobase million (FPKM). The primary solid tumor samples with the detailed FIGO stage, histological grade, and overall survival (OS) were selected. To validate the results, two additional raw datasets GSE26193 and GSE9891 were downloaded from the Gene Expression Omnibus (GEO)2.

OV Data Processing

The Ensembl IDs of the TCGA data were converted to gene names by the annotation dataset downloaded from the Ensembl genome browser3, and the probe IDs in the GEO data were processed using the relevant platform dataset. If multiple Ensembl or probe IDs matched a gene name, the average of all IDs was used to calculate the expression value of a gene. The FPKM OV data of TCGA were transformed to the transcripts per million (TPM) data format because TPM can be used to transform the RNA-seq expression data similar to the gene expression data (Vera Alvarez et al., 2019). The GEO OV data were normalized using the mas5 algorithm in the Affy R package (Gautier et al., 2004). All data used in the study were in a non-log format.

Evaluation of the Proportion of TME Cells in OV

CIBERSORT and MCP-counter were used to transform the RNA-seq data into the proportion of TME cells. The TCGA data were reformatted as required for CIBERSORT4. LM22 was selected as the signature file, and 1,000 permutations were used to acquire 22 immune cell fractions in 359 OV samples; the sum of the 22 immune cell fractions of each sample was assumed to be 1. The p-values of all TCGA samples were calculated, and P < 0.05 indicated the predominance of the immune cells in the sample. Screening of the results based on P < 0.05 yielded 172 TCGA samples with a certain fraction of 21 immune cell types (one cell type was removed because the fraction was zero in the results). The MCP-counter R package was used to evaluate the expression of two additional TME cell types, fibroblasts and endothelial cells. Then, the z-score was used to normalize the data acquired by two separate methods, and the results were combined to finally obtain the proportion of 23 TME cell types in 172 OV samples.

Clustering Methods

To determine a better classification method for investigating OV, consensus clustering of 172 OV data samples was performed by the ConsensusClusterPlus R package (Wilkerson and Hayes, 2010). The optimal inflection point (k = 4) of the elbow plot was used to define the best cluster number of OV samples. Finally, four new OV subtypes were defined based on TME cells; these subtypes were different from the stage or grade. Additionally, 23 TME cell types were also clustered into four groups by hierarchical clustering based on Euclidean distance and Ward’s linkage.

Cox Regression Analysis of TME Cells

Univariate Cox regression analysis of the proportion of 23 TME cell types and OS was performed for each sample. Six TME cell types were selected according to the hazard ratio and p-value of the Cox regression results to generate an immune risk model based on linear regression (Hijazi et al., 2016).

immunoscore=i=1nβjExp(celli)

In the model, Exp(celli) corresponds to the proportion of TME celli and βj corresponds to the regression coefficient of celli obtained by the Cox regression analysis. Each sample was assigned an immunoscore based on the model, and the median of all immunoscores was defined as a cutoff value to divide the OV data into the high and low immune risk groups.

Survival Analysis

Kaplan-Meier analysis was used to determine the relationship between the four subtypes and OS and to evaluate the prognostic significance of various subtypes by the survival R package (Lin and Zelterman, 2002). The survival analysis of the two immune risk groups was performed to determine the OS differences between the groups.

Gene Set Variation Analysis

Gene Set Variation Analysis (GSVA) was performed by the GSVA R package (Hanzelmann et al., 2013), which was used to calculate the GSVA score of each hallmark pathway in 172 TCGA OV samples. Each hallmark pathway was characterized by exclusive gene sets downloaded from GSEA5. GSVA converted the expression of the gene sets into the GSVA score of the corresponding hallmark pathway in each sample. A high GSVA score corresponds to the high enrichment of a hallmark pathway in a sample.

Statistical Analysis

Two groups were compared using Student’s t-test, and multivariate groups were compared by Kruskal-Wallis test. The receiver operating characteristic curve (ROC curve) was used to evaluate the sensitivity and specificity of the immunoscore-based prediction of OS. The area under the curve (AUC) was calculated by the timeROC R package (Blanche et al., 2013) to acquire the best predicted indicator. All analyses were implemented in R 3.6, and all p-values in the analyses were considered significant at P < 0.05.

Results

The Immune Landscape of OV Identified Four New Cell Subtypes

A total of 359 primary solid tumor samples of OV with detailed clinical information were downloaded from the TCGA dataset (Table 1). The gene expression data of these samples were converted to the ratio of 22 immune cell types by CIBERSORT filtered by the p-value to obtain the results (P < 0.05 corresponding to a high proportion of immune cells compared to non-immune cells in the tumor tissue; Newman et al., 2015). Two additional cell types (fibroblasts and endothelial cells) were acquired using the MCP-counter R package. As a result, the final data included the proportion of 23 TME cell types in 172 OV samples (CD4+ naïve T cells were removed because the proportion of this cell type was zero). Detailed clinical characteristics and cell proportions are listed in Supplementary Table S1.

TABLE 1
www.frontiersin.org

Table 1. The clinical features of 172 TCGA OV samples.

To analyze the immune landscape of OV, consistent clustering was used to divide OV samples on the basis of the proportion of the 23 TME cell types. As a result, four new cell subtypes of OV were identified and named according to the functions of the cells: the immune procancer cell subtype (IPCCS, n = 52) (Kalluri and Zeisberg, 2006; Yuan et al., 2017; Lee et al., 2019), immune killer cell subtype (IKCS, n = 62) (Uzhachenko and Shanker, 2019), immune proantibody cell subtype (IPACS, n = 24) (Oracki et al., 2010), and immune helper cell subtype (IHCS, n = 34) (Banchereau and Steinman, 1998; Crotty, 2019; Figure 1A and Supplementary Table S2). The OS analysis indicated that IPCCS was associated with the poorest prognosis, and the other three subtypes were associated with similar or better prognosis than IPCCS (P < 0.001, Figure 1B). The landscape of TME cell infiltration in OV was used to assess the proportion of each cell type and detailed clinical features of 172 TCGA OV samples (Figure 1C). Therefore, IPCCS may be a new OV type that can predict prognosis regardless of clinical characteristics. TME cell types included in IPCCS were different from those in three other subtypes (Figure 1C). Thus, additional investigation of IPCCS is needed to identify the mechanism of association with poor OS.

FIGURE 1
www.frontiersin.org

Figure 1. Division and analysis of four new cell subtypes in OV. (A) A clustergram of OV by the ConsensuSubtypePlus R package (k = 4). (B) The survival curves of four subtypes of OV. (C) The heatmap of 172 TCGA OV samples with four subtypes and clinical features.

Different Cell Subtypes Manifest Specific TME Cell Infiltration Patterns

Significant discrepancies in IPCCS and the other three subtypes were identified in the case of 20 out of 23 TME cell types (Figure 2A and Supplementary Figure S1). M0 macrophages, M2 macrophages, activated mast cells, neutrophils, endothelial cells, and fibroblasts were considerably enriched in IPCCS, but a relatively low proportion of M1 macrophages and CD8+ T cells was found. IKCS was characterized by a high proportion of activated NK cells, CD8+ T cells, and regulatory T cells (Tregs). IPACS included a high number of naïve B cells and plasma cells. Activated dendritic cells, follicular helper T cells, and resting memory CD4+ T cells demonstrated high infiltration in IHCS, and Tregs were characterized by low infiltration.

FIGURE 2
www.frontiersin.org

Figure 2. The distribution of TME cells in four subtypes and different clinical features of OV. (A) The proportion of each TME cell type in four subtypes of OV. (B) The distribution of neutrophils and fibroblasts in stages II, III, and IV of OV. (C) The distribution of activated mast cells in the lymph node metastases of OV. (D) The distribution of naïve B cells and endothelial cells in various residual tumor diameter of OV (statistical comparison between the two groups was performed by Student’s t-test, and multivariate groups were compared by Kruskal-Wallis test. ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Not all p-values corresponding to M2 macrophages were significant; however, the expression of M2 macrophages was high in IPCCS similar to that of M0 macrophages (Figure 2A). The proportion of M2 and M1 macrophages demonstrated an opposite trend. This phenomenon confirmed that macrophage polarization to M1 and M2 may have different effects on the occurrence and development of cancer (Cheng et al., 2019). CD8+ T cells and NK cells destroy cancer cells, play a critical role in immune effector activity and are advantageous to the outcome of multiple tumors (McGrail et al., 2018). These observations are similar to our findings that CD8+ T cells and activated NK cells were low expression in IPCCS, which is associated with the worst prognosis of all four subtypes. Thus, CD8+ T cells and activated NK cells may kill tumor cells in OV, and a new cell subtype identified by us can predict the prognosis of OV based on the proportion of TME cell types. Cancer-associated fibroblasts (CAFs) are important TME cells that release factors and boost angiogenesis leading to poor tumor outcome and resistance to treatment (Gascard and Tlsty, 2016). Our data indicated that CAFs were considerably increased in IPCCS, which was associated with the worst prognosis.

The distribution of TME cell types in relation to various OV clinical features (FIGO stage, histological grade, lymph node metastasis, and residual tumor diameter) was investigated. The results indicated that neutrophils and fibroblasts were distributed in stages II, III, and IV of OV (Figure 2B). These two cell types have completely opposite trends. Lee et al. (2019) reported that an increase in neutrophils in the premetastatic omental niche could be stimulated in early stage of OV. Zhao et al. (2017) showed that fibroblasts are characterized by high expression of CXCL 14 that is overexpressed in advanced-stage OV. Activated mast cells were present in variable proportions in lymph node metastasis of OV (Figure 2C). Crivellato et al. (2008) showed that mast cells can stimulate tumor angiogenesis and lymphangiogenesis and that the inhibition of mast cells decreased the lymph node metastasis. At the same time, naïve B cells and endothelial cells had different distribution in various residual tumor diameters (diameter of residual tumor in the pelvic and abdominal cavity after surgical treatment of OV) (Figure 2D). Cecilia S. Leung et al. demonstrated that endothelial cells can promote tumor progression and chemoresistance (Leung et al., 2018) to subsequently increase the probability of residual tumor.

Certain Distinctive TME Cells Contribute to Diverse Immune Risk Grouping of Cell Subtypes

Given that TME cells are closely associated with the tumor, we generated a model that used TME cells to predict the immune risk of OV. Cox regression analysis identified six TME cell types (M1 macrophages, plasma cells, follicular helper T cells, neutrophils, fibroblasts, and M0 macrophages) with P < 0.05 according to the results of the analysis. These six TME cell types were used to establish an immune risk model of OV, and the results are shown as a forest plot (Figure 3A). The risk model included the results of the fraction levels multiplied by the regression coefficient of each of the six cell types. The immunoscore of 172 OV samples was calculated, and the median was used to separate the samples into the high and low immune risk groups (Supplementary Table S3). The Kaplan-Meier survival analysis of these two groups demonstrated considerable differences, and the high immune risk group was associated with poor OS (Figure 3B).

FIGURE 3
www.frontiersin.org

Figure 3. Calculation and application of the immunoscore. (A) The Cox regression analysis of 23 TME cell types of OV. (B) The survival curves of the immunoscore in OV. (C) The distribution of the immunoscore in four subtypes of OV. (D) The number of four subtypes of OV samples in the high and low immunoscore groups (statistical comparison between the multivariate groups was performed by Kruskal-Wallis test. *P < 0.05, ****P < 0.0001).

Comparison of the immunoscores of four subtypes indicated that IPCCS had the highest immunoscore and that the immunoscore of the four subtypes was gradually and slightly decreased (Figure 3C). This result is in agreement with the data of our previous survival analysis. IPCCS had the highest proportion in the high immune risk group and the lowest percentage in the low immune risk group, which was opposite to that of IPACS and IHCS (Figure 3D). The results indicated that IPCCS had a positive correlation with the high immune risk group. In contrast, IPACS and IHCS were positively correlated with the low immune risk group.

Different Subtypes Have Characteristics of Various Molecular Mechanisms and Pathways

To analyze the molecular mechanisms corresponding to various subtypes, current literatures were searched for data on chemokines and chemokine receptors (Nagarsheth et al., 2017), immune checkpoints (Huang et al., 2017), MHC class I and II molecules (Johnson et al., 2018), and hypoxia-related genes (Eustace et al., 2013). Initially, the expression of chemokines and chemokine receptors was compared in the four subtypes and two immune risk groups. The results showed that CXCL8, CXCL5, CXCL3, CXCL12, CXCR1, and CXCL1 were expressed at high levels in IPCCS and the high immune risk group (Figure 4A and Supplementary Figure S2A). Certain studies demonstrated that these chemokines or chemokine receptors are associated with the aggregation and activation of macrophages (Qin et al., 2017; Roca et al., 2018), neutrophils (De Filippo et al., 2013), and fibroblasts (Feig et al., 2013; Naito et al., 2019), which were also enriched in IPCCS; these factors may in part account for different TME cell infiltration patterns.

FIGURE 4
www.frontiersin.org

Figure 4. The expression of various molecules and hallmark pathways in the four subtypes. (A) The expression of chemokines and chemokine receptors in the four subtypes. (B) The expression of immune checkpoints in the four subtypes. (C) The expression of MHC class I and II molecules in the four subtypes. (D) The expression of hypoxia-related genes in the four subtypes. (E) The GSVA score of the hallmark pathways in the four subtypes (statistical comparison between the multivariate groups was performed by Kruskal-Wallis test. ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

The majority of immune checkpoints were expressed at low levels in IPCCS and the high immune risk group (Figure 4B and Supplementary Figure S2B), especially including well-known PD-1 (PDCD1), PD-L1 (CD274), LAG-3, and CTLA-4 molecules, similar to those described in a study of Xiao et al. (2019). Thus, checkpoints do not play an important role in the progression of OV in IPCCS. Hence, we hypothesized that the effect of immunotherapy with checkpoint blockers in IPCCS may not be as pronounced as that in the other subtypes. The expression of most MHC class I and II molecules was also low in IPCCS and the high immune risk group (Figure 4C and Supplementary Figure S2C). Low proportion of various T cells in IPCCS may suggest that T cells cannot proliferate and recognize tumor cells due to low expression of MHC class I and II molecules.

Most hypoxia-related genes were expressed at high levels in IPCCS and the high immune risk group indicated a certain degree of hypoxia in the environment of OV cells (Figure 4D and Supplementary Figure S2D). Bhandari et al. (2019) estimated that OV has a higher hypoxia score, and hypoxia contributes to deterioration and chemotherapy resistance of OV (Dorayappan et al., 2018). Therefore, hypoxia may influence the proportion of TME cells and lead to poor prognosis of OV.

To determine the differences in the enrichment of hallmark pathways between various subtypes and immune risk groups, GSVA of the hallmark pathways was performed using 172 TCGA OV data. The GSVA score of the hallmark pathways demonstrated that certain pathways (such as EMT, TNFα signaling via NF-kB, hypoxia, angiogenesis, and notch signaling) were highly enriched in IPCCS and the high immune risk group (Figure 4E and Supplementary Figure S2E). These well-known pathways are known to be related to cancer progression (Chi et al., 2006; Takebe et al., 2014; Ottevanger, 2017; Markowska et al., 2017; Yang et al., 2018). In contrast, some pathways (such as DNA repair, spermatogenesis, oxidative phosphorylation, interferon γ response, and interferon α response) had negative effects in IPCCS and the high immune risk group compared to those in other subtypes and the low immune risk group (Figure 4E and Supplementary Figure S2E). Differential enrichment of these pathways may be used as an indicator of the mechanism of differences in immune infiltration and to improve individualized treatment of the patients with various subtypes and immune risks.

Verification of Immune Risk Grouping in Independent Cohorts

To confirm the prognostic accuracy of immune risk grouping, the model was applied in two independent GEO cohorts (GSE63885 and GSE9891) and the entire TCGA cohort (n = 359). The same six TME cell types were used to compute the immunoscore of each sample in the GEO and TCGA cohorts, and each cohort was divided into the high and low immune risk groups according to the median immunoscore (Supplementary Table S4). We found significant differences in the OS of the high and low immune risk groups in the GEO and TCGA cohorts (Figures 5A–C). The survival analysis indicated that the high immune risk group had worse prognosis similar to our previous results.

FIGURE 5
www.frontiersin.org

Figure 5. Verification of the immunoscore in independent cohorts. (A) The survival curves of immunoscore in GSE9891. (B) The survival curves of immunoscore in GSE63885. (C) The survival curves of immunoscore in TCGA OV. (D) The ROC curves of the clinical features and immunoscore in TCGA OV (numbers after the clinical features and immunoscore represent the corresponding AUC).

The sensitivity and specificity of the clinical features and immunoscore in predicting OS were evaluated in 172 TCGA samples. The results of the ROC curves analysis indicated that the AUC of the immunoscore was gradually increased from 1 to 5 years (Supplementary Figure S3), and the immunoscore had the largest AUC compared to those of all clinical features (Figure 5D). Thus, the immunoscore had the highest associations with predicted OS in the TCGA OV samples. In summary, the immune risk grouping developed in the present study can be a powerful means to evaluate the prognosis of OV.

Discussion

Previous tumor studies have focused on aberrant genes, epigenetics, or non-coding RNA. Subsequent studies have gradually revealed associations between inflammation, immune cells, and the TME and tumorigenesis and tumor development. The effect of the TME on tumors is a complex and dynamic process. The proportion and activation status of TME cells differentially influence the proliferation, invasion, and metastasis of tumor cells (Quail and Joyce, 2013). Therefore, we investigated the effect of the TME on OV by a new approach designed in the present study.

Current studies on the TME of OV are based on a single or a few cell types; the expression of these cell types depends mainly on detecting several specific gene markers. Our study defined 21 TME cell types by CIBERSORT and 2 TME cell types by MCP-counter thus expanding the data coverage of our study of the TME of OV. Each cell type had dozens of representative genes that can increase the accuracy of the estimation of their expression. MCP-counter was selected mainly because it includes fibroblasts and endothelial cells, which have been shown to play key roles in the occurrence and development of cancer (Maishi and Hida, 2017; Thuwajit et al., 2018).

CIBERSORT and MCP-counter analyses of the TCGA OV data yielded the proportions of 23 TME cell types of OV. A classification diagram was used to divide all OV data into four subtypes based on the proportion of 23 TME cell types indicating that TME of OV has four main categories. These subtypes were different from the routine stage and grade characteristics, and each subtype had unique cell types that may account for different prognosis of OV. The Cox regression analysis of these cell types was used to select six TME cell types (P < 0.05) to construct an immune risk model. The model was used to calculate the immunoscore and evaluate the immune risk in each OV sample. The results indicated that different cell subtypes of OV are associated with variable immune risk. IPCCS was associated with high immune risk and the worst prognosis. In contrast, IPACS and IHCS were associated with low immune risk and better prognosis. The data indicated that the model can reliably distinguish differences in the survival time of the patients in the high and low immune risk groups; the high immune risk group was associated with poor prognosis. Consistent results were obtained by validating the model using the GEO database. The ROC curve analysis also showed that the immunoscore had higher sensitivity and specificity in the prediction of prognosis compared with that of other clinical features.

Previous studies demonstrated that chemokines and chemokine receptors play positive roles in inflammation and oncogenesis by regulating the trafficking of various inflammatory cells (Bian et al., 2019). Immune checkpoints play an important role in TME and can be potential targets for cancer treatment (Toor et al., 2020). MHC I and II molecules can present protein fragments to T cells that are essential for cell-mediated immunity and tumor immunity (Rock et al., 2016; Alspach et al., 2019). At the same time, hypoxia is frequent in TME and can change the components of TME leading to poor prognosis of cancer (Jing et al., 2019). Therefore, these genes are related to TME and have to be considered in the studies of immune characteristics of OV. The present study evaluated the expression of the chemokines, chemokine receptors, immune checkpoints, MHC class I and II molecules, hypoxia-related genes, and enrichment of the hallmark pathways in four subtypes to determine the mechanisms of different immune infiltration patterns. The results indicated that differences in TME cell infiltration in the subtypes were closely related to the expression of the chemokines and chemokine receptors because of the characteristics of these molecules that induced directional migration of immune cells (Sokol and Luster, 2015). High expression of the hypoxia-related genes in IPCCS was caused by hypoxic conditions in the TME (Jing et al., 2019). Low expression of immune checkpoints (such as PD-1, LAG-3, CTLA-4) which mainly exist on the surface of T cells (Dyck and Mills, 2017) in IPCCS mainly due to the lower proportion of various T cells in this subtype than the other three subtypes (IKCS, IPACS, and IHCS). Studies stated that checkpoints blockade might as new targets for cancer immunotherapy (Feng et al., 2019). Therefore, we could speculate that the therapy of checkpoints blockade might have better effects in the other three subtypes of OV. Numerous studies demonstrated that a number of pathways are linked to the progression and prognosis of OV (Wang et al., 2017; Garsed et al., 2018; Guo et al., 2018; Zeng et al., 2018). We compared differential enrichment of the hallmark pathways in four subtypes, and the results indicated that the mechanisms of the action of different cell subtypes are related to different pathways.

Although the OV dataset used in the present study was downloaded from TCGA, the number of OV samples was cut by almost 50% by CIBERSORT. In the future, the collection of additional clinical data on OV, corresponding clinical characteristics, and survival time is needed to increase the accuracy of our study. Our findings may require verification in additional investigations of TME cells and pathways. Differences in gene mutations, expression of ncRNAs, copy number variations, and other factors in the four subtypes require further studies.

Conclusion

In summary, our study reclassified OV into four subtypes (IPCCS, IKCS, IPACS, and IHCS) according to TME cell types and demonstrated infiltration of exclusive cell types in each subtype of OV. Then, six TME cell types were selected by Cox regression analysis to calculate the immunoscore that could assess the immune risk and predict the prognosis of OV. The results showed that IPCCS was associated with high immune risk and poor prognosis. Finally, the analysis of the mechanisms of various subtypes of OV was performed, and the results may assist in identifying effective therapeutic targets for OV.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: http://xena.ucsc.edu; https://www.ncbi.nlm.nih.gov/geo.

Author Contributions

GZ and SN designed and directed all research. SC, QG, YC, YH, XZ, and CK performed the data processing and experimental analyses. GZ, SN, SC, and QG drafted the manuscript. All authors reviewed and approved the final version of the manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Grant Nos. 81772780 and 81902646).

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.

Supplementary Material

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

Footnotes

  1. ^ http://xena.ucsc.edu
  2. ^ https://www.ncbi.nlm.nih.gov/geo
  3. ^ https://uswest.ensembl.org
  4. ^ https://cibersort.stanford.edu
  5. ^ https://www.gsea-msigdb.org

References

Alspach, E., Lussier, D. M., Miceli, A. P., Kizhvatov, I., DuPage, M., Luoma, A. M., et al. (2019). MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature 574, 696–701. doi: 10.1038/s41586-019-1671-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Baghban, R., Roshangar, L., Jahanban-Esfahlan, R., Seidi, K., Ebrahimi-Kalan, A., Jaymand, M., et al. (2020). Tumor microenvironment complexity and therapeutic implications at a glance. Cell Commun. Signal. 18:59. doi: 10.1186/s12964-020-0530-4

CrossRef Full Text | Google Scholar

Banchereau, J., and Steinman, R. M. (1998). Dendritic cells and the control of immunity. Nature 392, 245–252. doi: 10.1038/32588

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17:218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhandari, V., Hoey, C., Liu, L. Y., Lalonde, E., Ray, J., Livingstone, J., et al. (2019). Molecular landmarks of tumor hypoxia across cancer types. Nat. Genet. 51, 308–318. doi: 10.1038/s41588-018-0318-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bian, X., Xiao, Y. T., Wu, T., Yao, M., Du, L., Ren, S., et al. (2019). Microvesicles and chemokines in tumor microenvironment: mediators of intercellular communications in tumor progression. Mol. Cancer. 18:50. doi: 10.1186/s12943-019-0973-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanche, P., Dartigues, J. F., and Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat. Med. 32, 5381–5397. doi: 10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, H., Wang, Z., Fu, L., and Xu, T. (2019). Macrophage Polarization in the Development and Progression of Ovarian Cancers: An Overview. Front. Oncol. 9:421. doi: 10.3389/fonc.2019.00421

PubMed Abstract | CrossRef Full Text | Google Scholar

Chi, J. T., Wang, Z., Nuyten, D. S., Rodriguez, E. H., Schaner, M. E., Salim, A., et al. (2006). Gene expression programs in response to hypoxia: cell type specificity and prognostic significance in human cancers. PLoS Med. 3:e47. doi: 10.1371/journal.pmed.0030047

PubMed Abstract | CrossRef Full Text | Google Scholar

Crivellato, E., Nico, B., and Ribatti, D. (2008). Mast cells and tumour angiogenesis: new insight from experimental carcinogenesis. Cancer Lett. 269, 1–6. doi: 10.1016/j.canlet.2008.03.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Crotty, S. (2019). T Follicular Helper Cell Biology: A Decade of Discovery and Diseases. Immunity 50, 1132–1148. doi: 10.1016/j.immuni.2019.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

De Filippo, K., Dudeck, A., Hasenberg, M., Nye, E., van Rooijen, N., Hartmann, K., et al. (2013). Mast cell and macrophage chemokines CXCL1/CXCL2 control the early stage of neutrophil recruitment during tissue inflammation. Blood 121, 4930–4937. doi: 10.1182/blood-2013-02-486217

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorayappan, K. D. P., Wanner, R., Wallbillich, J. J., Saini, U., Zingarelli, R., Suarez, A. A., et al. (2018). Hypoxia-induced exosomes contribute to a more aggressive and chemoresistant ovarian cancer phenotype: a novel mechanism linking STAT3/Rab proteins. Oncogene 37, 3806–3821. doi: 10.1038/s41388-018-0189-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Drakes, M. L., and Stiff, P. J. (2018). Regulation of Ovarian Cancer Prognosis by Immune Cells in the Tumor Microenvironment. Cancers 10:302. doi: 10.3390/cancers10090302

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyck, L., and Mills, K. H. G. (2017). Immune checkpoints and their inhibition in cancer and infectious diseases. Eur. J. Immunol. 47, 765–779. doi: 10.1002/eji.201646875

PubMed Abstract | CrossRef Full Text | Google Scholar

Eustace, A., Mani, N., Span, P. N., Irlam, J. J., Taylor, J., Betts, G. N., et al. (2013). A 26-gene hypoxia signature predicts benefit from hypoxia-modifying therapy in laryngeal cancer but not bladder cancer. Clin. Cancer Res. 19, 4879–4888. doi: 10.1158/1078-0432.CCR-13-0542

PubMed Abstract | CrossRef Full Text | Google Scholar

Feig, C., Jones, J. O., Kraman, M., Wells, R. J., Deonarine, A., Chan, D. S., et al. (2013). Targeting CXCL12 from FAP-expressing carcinoma-associated fibroblasts synergizes with anti-PD-L1 immunotherapy in pancreatic cancer. Proc. Natl. Acad. Sci. USA. 110, 20212–20217. doi: 10.1073/pnas.1320318110

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, M., Jiang, W., Kim, B. Y. S., Zhang, C. C., Fu, Y. X., and Weissman, I. L. (2019). Phagocytosis checkpoints as new targets for cancer immunotherapy. Nat. Rev. Cancer 19, 568–586. doi: 10.1038/s41568-019-0183-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Garsed, D. W., Alsop, K., Fereday, S., Emmanuel, C., Kennedy, C. J., Etemadmoghadam, D., et al. (2018). Homologous Recombination DNA Repair Pathway Disruption and Retinoblastoma Protein Loss Are Associated with Exceptional Survival in High-Grade Serous Ovarian Cancer. Clin. Cancer Res. 24, 569–580. doi: 10.1158/1078-0432.CCR-17-1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Gascard, P., and Tlsty, T. D. (2016). Carcinoma-associated fibroblasts: orchestrating the composition of malignancy. Genes Dev. 30, 1002–1019. doi: 10.1101/gad.279737.116

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautier, L., Cope, L., Bolstad, B. M., and Irizarry, R. A. (2004). affy–analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 20, 307–315. doi: 10.1093/bioinformatics/btg405

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghisoni, E., Imbimbo, M., Zimmermann, S., and Valabrega, G. (2019). Ovarian Cancer Immunotherapy: Turning up the Heat. Int. J. Mol. Sci. 20:2927. doi: 10.3390/ijms20122927

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, C., Wang, X., Chen, L. P., Li, M., Li, M., Hu, Y. H., et al. (2018). Long non-coding RNA MALAT1 regulates ovarian cancer cell proliferation, migration and apoptosis through Wnt/beta-catenin signaling pathway. Eur. Rev. Med. Pharmacol. Sci. 22, 3703–3712. doi: 10.26355/eurrev_201806_15249

CrossRef Full Text | Google Scholar

Guo, S., and Deng, C. X. (2018). Effect of Stromal Cells in Tumor Microenvironment on Metastasis Initiation. Int. J. Biol. Sci. 14, 2083–2093. doi: 10.7150/ijbs.25720

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Henke, E., Nandigama, R., and Ergun, S. (2019). Extracellular Matrix in the Tumor Microenvironment and Its Impact on Cancer Therapy. Front. Mol. Biosci. 6:160. doi: 10.3389/fmolb.2019.00160

PubMed Abstract | CrossRef Full Text | Google Scholar

Hijazi, Z., Lindback, J., Alexander, J. H., Hanna, M., Held, C., Hylek, E. M., et al. (2016). The ABC (age, biomarkers, clinical history) stroke risk score: a biomarker-based risk score for predicting stroke in atrial fibrillation. Eur. Heart J. 37, 1582–1590. doi: 10.1093/eurheartj/ehw054

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, R. Y., Francois, A., McGray, A. R., Miliotto, A., and Odunsi, K. (2017). Compensatory upregulation of PD-1, LAG-3, and CTLA-4 limits the efficacy of single-agent checkpoint blockade in metastatic ovarian cancer. Oncoimmunology 6:e1249561. doi: 10.1080/2162402X.2016.1249561

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, X., Yang, F., Shao, C., Wei, K., Xie, M., Shen, H., et al. (2019). Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol. Cancer. 18:157. doi: 10.1186/s12943-019-1089-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D. B., Nixon, M. J., Wang, Y., Wang, D. Y., Castellanos, E., Estrada, M. V., et al. (2018). Tumor-specific MHC-II expression drives a unique pattern of resistance to immunotherapy via LAG-3/FCRL6 engagement. JCI Insight. 3:e120360. doi: 10.1172/jci.insight.120360

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalluri, R., and Zeisberg, M. (2006). Fibroblasts in cancer. Nat. Rev. Cancer 6, 392–401. doi: 10.1038/nrc1877

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamichhane, P., Karyampudi, L., Shreeder, B., Krempski, J., Bahr, D., Daum, J., et al. (2017). IL10 Release upon PD-1 Blockade Sustains Immunosuppression in Ovarian Cancer. Cancer Res. 77, 6667–6678. doi: 10.1158/0008-5472.CAN-17-0740

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, W., Ko, S. Y., Mohamed, M. S., Kenny, H. A., Lengyel, E., and Naora, H. (2019). Neutrophils facilitate ovarian cancer premetastatic niche formation in the omentum. J. Exp. Med. 216, 176–194. doi: 10.1084/jem.20181170

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, C. S., Yeung, T. L., Yip, K. P., Wong, K. K., Ho, S. Y., Mangala, L. S., et al. (2018). Cancer-associated fibroblasts regulate endothelial adhesion protein LPP to promote ovarian cancer chemoresistance. J. Clin. Invest. 128, 589–606. doi: 10.1172/JCI95200

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, H., and Zelterman, D. (2002). Modeling Survival Data: Extending the Cox Model. Technometrics 44, 85–86. doi: 10.1198/tech.2002.s656

PubMed Abstract | CrossRef Full Text | Google Scholar

Maishi, N., and Hida, K. (2017). Tumor endothelial cells accelerate tumor metastasis. Cancer Sci. 108, 1921–1926. doi: 10.1111/cas.13336

PubMed Abstract | CrossRef Full Text | Google Scholar

Markowska, A., Sajdak, S., Markowska, J., and Huczynski, A. (2017). Angiogenesis and cancer stem cells: New perspectives on therapy of ovarian cancer. Eur. J. Med. Chem. 142, 87–94. doi: 10.1016/j.ejmech.2017.06.030

PubMed Abstract | CrossRef Full Text | Google Scholar

McGrail, D. J., Federico, L., Li, Y., Dai, H., Lu, Y., Mills, G. B., et al. (2018). Multi-omics analysis reveals neoantigen-independent immune cell infiltration in copy-number driven cancers. Nat. Commun. 9:1317. doi: 10.1038/s41467-018-03730-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagarsheth, N., Wicha, M. S., and Zou, W. (2017). Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat. Rev. Immunol. 17, 559–572. doi: 10.1038/nri.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Naito, Y., Yamamoto, Y., Sakamoto, N., Shimomura, I., Kogure, A., Kumazaki, M., et al. (2019). Cancer extracellular vesicles contribute to stromal heterogeneity by inducing chemokines in cancer-associated fibroblasts. Oncogene 38, 5566–5579. doi: 10.1038/s41388-019-0832-4

PubMed Abstract | CrossRef Full Text | Google Scholar

National Cancer Institute (2020). Surveillance, Epidemiology, and End Results Program. Cancer stat facts: ovarian cancer. Available online at https://seer.cancer.gov/statfacts/html/ovary.html (accessed October 1, 2020).

Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Oracki, S. A., Walker, J. A., Hibbs, M. L., Corcoran, L. M., and Tarlinton, D. M. (2010). Plasma cell development and survival. Immunol. Rev. 237, 140–159. doi: 10.1111/j.1600-065X.2010.00940.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ottevanger, P. B. (2017). Ovarian cancer stem cells more questions than answers. Semin. Cancer Biol. 44, 67–71. doi: 10.1016/j.semcancer.2017.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Peltanova, B., Raudenska, M., and Masarik, M. (2019). Effect of tumor microenvironment on pathogenesis of the head and neck squamous cell carcinoma: a systematic review. Mol. Cancer 18:63. doi: 10.1186/s12943-019-0983-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, C. C., Liu, Y. N., Hu, Y., Yang, Y., and Chen, Z. (2017). Macrophage inflammatory protein-2 as mediator of inflammation in acute liver injury. World J. Gastroenterol. 23, 3043–3052. doi: 10.3748/wjg.v23.i17.3043

PubMed Abstract | CrossRef Full Text | Google Scholar

Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19, 1423–1437. doi: 10.1038/nm.3394

PubMed Abstract | CrossRef Full Text | Google Scholar

Reislander, T., Lombardi, E. P., Groelly, F. J., Miar, A., Porru, M., Di Vito, S., et al. (2019). BRCA2 abrogation triggers innate immune responses potentiated by treatment with PARP inhibitors. Nat. Commun. 10:3143. doi: 10.1038/s41467-019-11048-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Roca, H., Jones, J. D., Purica, M. C., Weidner, S., Koh, A. J., Kuo, R., et al. (2018). Apoptosis-induced CXCL5 accelerates inflammation and growth of prostate tumor metastases in bone. J. Clin. Invest. 128, 248–266. doi: 10.1172/JCI92466

PubMed Abstract | CrossRef Full Text | Google Scholar

Rock, K. L., Reits, E., and Neefjes, J. (2016). Present Yourself! By MHC Class I and MHC Class II Molecules. Trends Immunol. 37, 724–737. doi: 10.1016/j.it.2016.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sokol, C. L., and Luster, A. D. (2015). The chemokine system in innate immunity. Cold Spr. Harb Pers. Biol. 7:a016303. doi: 10.1101/cshperspect.a016303

PubMed Abstract | CrossRef Full Text | Google Scholar

Takebe, N., Nguyen, D., and Yang, S. X. (2014). Targeting notch signaling pathway in cancer: clinical development advances and challenges. Pharmacol. Ther. 141, 140–149. doi: 10.1016/j.pharmthera.2013.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuwajit, C., Ferraresi, A., Titone, R., Thuwajit, P., and Isidoro, C. (2018). The metabolic cross-talk between epithelial cancer cells and stromal fibroblasts in ovarian cancer progression: Autophagy plays a role. Med. Res. Rev. 38, 1235–1254. doi: 10.1002/med.21473

PubMed Abstract | CrossRef Full Text | Google Scholar

Toor, S. M., Sasidharan Nair, V., Decock, J., and Elkord, E. (2020). Immune checkpoints in the tumor microenvironment. Semin. Cancer Biol. 65, 1–12. doi: 10.1016/j.semcancer.2019.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Uzhachenko, R. V., and Shanker, A. (2019). CD8(+) T Lymphocyte and NK Cell Network: Circuitry in the Cytotoxic Domain of Immunity. Front. Immunol. 10:1906. doi: 10.3389/fimmu.2019.01906

PubMed Abstract | CrossRef Full Text | Google Scholar

Vera Alvarez, R., Pongor, L. S., Marino-Ramirez, L., and Landsman, D. (2019). TPMCalculator: one-step software to quantify mRNA abundance of genomic features. Bioinformatics 35, 1960–1962. doi: 10.1093/bioinformatics/bty896

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Zhang, X., Tang, W., Lin, Z., Xu, L., Dong, R., et al. (2017). miR-130a upregulates mTOR pathway by targeting TSC1 and is transactivated by NF-kappaB in high-grade serous ovarian carcinoma. Cell Death Diff. 24, 2089–2100. doi: 10.1038/cdd.2017.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatic 26, 1572–1573. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, Y., Ma, D., Zhao, S., Suo, C., Shi, J., Xue, M. Z., et al. (2019). Multi-Omics Profiling Reveals Distinct Microenvironment Characterization and Suggests Immune Escape Mechanisms of Triple-Negative Breast Cancer. Clin. Cancer Res. 25, 5002–5014. doi: 10.1158/1078-0432.CCR-18-3524

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Liu, L., Li, C., Luo, N., Chen, R., Li, L., et al. (2018). TRIM52 plays an oncogenic role in ovarian cancer associated with NF-kB pathway. Cell Death Dis. 9:908. doi: 10.1038/s41419-018-0881-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, X., Zhang, J., Li, D., Mao, Y., Mo, F., Du, W., et al. (2017). Prognostic significance of tumor-associated macrophages in ovarian cancer: A meta-analysis. Gynecol. Oncol. 147, 181–187. doi: 10.1016/j.ygyno.2017.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, B., Zhou, M., Wu, H., and Xiong, Z. (2018). SPP1 promotes ovarian cancer progression via Integrin beta1/FAK/AKT signaling pathway. Onco. Targets Ther. 11, 1333–1343. doi: 10.2147/OTT.S154215

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, L., Ji, G., Le, X., Wang, C., Xu, L., Feng, M., et al. (2017). Long Noncoding RNA LINC00092 Acts in Cancer-Associated Fibroblasts to Drive Glycolysis and Progression of Ovarian Cancer. Cancer Res. 77, 1369–1382. doi: 10.1158/0008-5472.CAN-16-1615

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: ovarian cancer, immune, tumor microenvironment, CIBERSORT, prognosis

Citation: Cong S, Guo Q, Cheng Y, He Y, Zhao X, Kong C, Ning S and Zhang G (2020) Immune Characterization of Ovarian Cancer Reveals New Cell Subtypes With Different Prognoses, Immune Risks, and Molecular Mechanisms. Front. Cell Dev. Biol. 8:614139. doi: 10.3389/fcell.2020.614139

Received: 05 October 2020; Accepted: 16 November 2020;
Published: 21 December 2020.

Edited by:

Lei Deng, Central South University, China

Reviewed by:

Guangchuang Yu, Southern Medical University, China
Chuan-Le Xiao, Sun Yat-sen University, China

Copyright © 2020 Cong, Guo, Cheng, He, Zhao, Kong, Ning and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Guangmei Zhang, guangmei_zhang@126.com; Shangwei Ning, ningsw@ems.hrbmu.edu.cn

Download