Immune profiling reveals prognostic genes in high-grade serous ovarian cancer

Study design and patient selection

The overall scheme is shown in Figure 1A. A total of 361 ovarian cancer samples were enrolled according to the criteria (eligible clinical information, overall survival > 1 month). Patient characteristics obtained from The Cancer Genome Atlas (TCGA) are detailed in Table 1. In general, we assessed the prognostic value of the TMEscore and then attempted to identify TME-related differentially expressed genes (DEGs) that contribute to HGSOC overall survival in the TCGA database. DEGs were further validated in the International Cancer Genome Consortium (ICGC), a separate HGSOC database.

Figure 1. Association between the TMEscore and prognosis in TCGA data analyzed by the CIBERSORT algorithm. (A) Workflow of the current study. (B) Based on the median TMEscore values, patients with HGSOC were divided into the high and low TMEscore groups. (C) As shown in the Kaplan-Meier plot, the median survival time in the high TMEscore group was longer than that in the low score group (p=0.0036). (D) Distribution of TMEscores by tumor stage for HGSOC patients. The boxplot shows that there was no association between tumor stage and TMEscore.

TMEscore is associated with overall survival in HGSOC

Of all HGSOC samples, 361 samples were evaluated. As calculated using the CIBERSORT algorithm, the TMEscore values ranged from -2.93883 to 2.27524 (Supplementary Table 1). Based on the median TMEscore values, 194 patients were in the high score group (53.7%), while 167 were in the low score group (46.3%, Figure 1B). The associations between the TMEscore and corres-ponding overall survival were then analyzed by generating Kaplan-Meier plots and evaluating the data with a log-rank test. The Kaplan-Meier plots demonstrated that a high TMEscore was positively correlated with favorable survival in HGSOC patients (Figure 1C, p=0.0036). Further evaluation revealed no correlation between TMEscore and tumor stage (p=0.54, Figure 1D).

DEGs between the TMEscore groups

We computed a heatmap to compare differential gene expression profiles in the high and low TMEscore groups (Figure 2A). We identified a total of 329 DEGs, including 311 upregulated and 18 downregulated genes (Supplementary Table 2). A volcano plot of gene expression profile data between patients with a high or low TMEscore is shown in Figure 2B. We carried out biological enrichment analyses, including Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, to outline the potential functions of DEGs. These genes were mainly involved in immune-related functions (Figure 2C–2F), such as T cell activation, cytokine receptor binding, T cell receptor complex and Th17 cell differentiation, suggesting that DEGs might exert functions in immune-related pathways in HGSOC.

Figure 2. DEG profiles by TMEscore in HGSOC. (A) Heatmap of the DEGs between the top half (high score) vs. bottom half (low score) of TMEscore values. A |log(fold change)|≥1 and an adjusted p-value < 0.05 were set as the cutoff criteria to screen for DEGs. (B) Volcano plot of gene expression profile data for patients with high and low TMEscores. (C–F) Functional enrichment analysis including Biological Process (BP), Cellular Components (CC), and Molecular Functions (MF) categories as well as KEGG pathways for 329 DEGs.

Survival analysis of DEGs in HGSOC

To explore the clinical role of individual DEGs in overall survival, we generated Kaplan-Meier survival curves from TCGA data. Based on the median value of each DEG, a total of 48 genes were shown to predict overall survival (Supplementary Table 3; p<0.05). The survival curves for the prognostic genes are shown in Figure 3A and Supplementary Figure 1. Consequently, enrichment analyses were performed with the prognostic DEGs. KEGG and GO analyses indicated that these genes were mainly involved in immune pathways such asTh17 cell differentiation, Cytokine−cytokine receptor interaction and Chemokine signaling pathways, similar to previous enrichment analysis of DEGs (Figure 3B–3E).

Figure 3. Discovery of prognostic TME-related DEGs with functional annotations in TCGA. (A) Kaplan-Meier survival curves were generated for selected DEGs with prognostic significance by a log-rank test. OS=overall survival time in months. (B–E) GO term and KEGG pathway analyses with the 48 prognostic DEGs.

Protein-protein interaction (PPI) network among prognostic genes

To explore the interplay among the 48 DEGs, we used the STRING database (https://string-db.org/) to analyze protein network interactions. We obtained a relationship network containing 222 edges and 43 nodes (Figure 4A). Moreover, in the PPI network, CXCL9, CXCL13, CCL5, GZMB and CD2 were the most remarkable nodes as they had the most connections with other nodes. Molecular COmplex DEtection (MCODE), a Cytoscape plugin, was used to identify the main coregulated modules. By applying the MCODE tool, we identified two closely related subgroups (Figure 4B). As shown in Figure 4B, we can see that all the module genes are up-regulated genes of DEGs, including CXCL9, CCL5, CXCL11, CD27, FOXP3, CD8A and other immune-related genes. The connected nodes for each gene intersection are listed in Supplementary Table 4.

Figure 4. Construction of the PPI network for the 48 prognostic DEGs. (A) The PPI network was constructed using the 48 prognostic DEGs with the R software package STRINGdb. (B) MCODE was used to identify the main coregulated modules. The most significant module is indicated in two closely related subgroups.

Prognostic validation in the ICGC database

To further validate the prognostic significance of 48 genes identified from the TCGA database, we downloaded and analyzed the gene expression data of a cohort of 81 HGSOC patients from the ICGC (Supplementary Table 5), an independent HGSOC database. Forty-seven of the 48 genes identified above had expression values, and further analyses validated seven genes (p<0.05, Figure 5A–5G) associated with prognosis (Supplementary Table 6). Notably, elevated mRNA levels of these genes correlated with better patient OS.

Figure 5. Validation of prognostic DEGs extracted from the TCGA database in the ICGC cohort. (A–G) Kaplan-Meier survival curves were generated for seven validated DEGs in an additional cohort of 81 HGSOC patients from the ICGC. p<0.05 by a log-rank test. OS=overall survival time in days.

GBP1, ETV7 and CXCL13 perform molecular functions in vitro

Crosstalk between tumor cells and their surrounding microenvironment is necessary for cell survival, growth, and proliferation and for epithelial-mesenchymal transition (EMT) and metastasis [16–18]. We sought to examine the molecular functions of these genes. The three genes with the lowest p-values validated from the ICGC database were selected for functional assays in HGSOC cells. A2780 cells were separately transduced with small interfering RNAs (siRNAs) targeting these three genes, and western blotting confirmed the efficiency of the siRNAs in A2780 cells (Figure 6A). The results showed that knockdown of GBP1 and ETV7 promoted the proliferation, colony formation, and migration of cells (Figure 6B–6D). In contrast, the opposite effects were observed when CXCL13 was downregulated (Figure 6B–6D).

Figure 6. In vitro assessment of GBP1, ETV7 and CXCL13. (A) The protein levels of GBP1, ETV7 and CXCL13 were measured by western blotting in A2780 cells transfected with siRNAs. (B) CCK-8 assays were performed to evaluate the proliferation of GBP1-, ETV7- and CXCL13-knockdown cells. (C) The colony-forming ability of A2780 cells was assessed to determine the effects of GBP1, ETV7 or CXCL13 downregulation on cell growth. (D) The invasion potential of cells was assessed using a Transwell assay. The scale bar represents 100 μm. NC: negative control. * indicates p< 0.05.

Confirmation of the expression levels of GBP1, ETV7 in patients

To further validate the clinical importance of GBP1, ETV7 and CXCL13, we aimed to perform immunostaining with antibodies on microarrays containing tissue from 165 HGSOC patients. Given that no antibodies specific for ETV7 and CXCL13 are commercially available, RT-PCR was adopted to confirm their expression levels and prognostic significance. We first assessed the protein expression of GBP1. Positive immunostaining signals of GBP1 were seen in the cytoplasm and membrane (Figure 7A). Of 165 patients with HGSOC, low expression (score 0-1) was observed in 94, and high expression (score 2-3) was observed in 71 (Figure 7B). The protein level of GBP1 was not correlated with any pathological characteristics of HGSOC but was potentially correlated with diaphragmatic metastasis (p=0.05; Table 2). We then plotted the overall survival curves according to GBP1 protein levels using the Kaplan-Meier method. As shown in Figure 7C, patients with high GBP1 protein expression exhibited more favorable overall survival (p=0.0003). Univariate analysis of overall survival revealed that GBP1 levels (p < 0.001), the residual tumor margins (p=0.003), diaphragmatic metastasis (p = 0.013) and mesenteric metastasis (p=0.021) were prognostic indicators in HGSOC (Figure 7D). Further multivariate analysis revealed that GBP1 expression (p < 0.001) and the residual tumor margins (p=0.003) were independent predictors for overall survival in HGSOC patients (Figure 7E). Moreover, we observed decreased ETV7 and CXCL13 mRNA expression in HGSOC tissues compared with normal ovary tissues (Figure 7F–7G). Kaplan-Meier analysis indicated that ETV7 and CXCL13 mRNA expression correlated with positive OS (p<0.05, Figure 7F–7G). Of note, the prognostic results of CXCL13 contradicted our in vitro experiments, perhaps because mRNA expression is not always consistent with protein levels. In addition, CXCL13, which is a cytokine, may exert its functions outside cells. In general, we demonstrated the prognostic value of GBP1 and ETV7 in HGSOC patients. Besides that, CXCL13 may require further large samples validation.

Figure 7. Levels of GBP1 and ETV7 expression correlated with overall survival in HGSOC. (A) Representative images of GBP1 expression in HGSOC tissues, visualized at 40× and 200× magnification. (B) Distribution of the immunoreactive score (IRS) in an HGSOC TMA. (C) Kaplan-Meier survival curve with log-rank analysis of overall survival showed statistical significance between the curves of patients with high GBP1 expression and those with low GBP1 expression (log-rank test, p=0.0003). (D) Univariate analysis was performed in 165 HGSOC patients. All bars correspond to 95% confidence intervals. (E) Multivariate analysis was performed in 165 HGSOC patients. All bars correspond to 95% confidence intervals. (F–G) Expression levels of ETV7 and CXCL13 were measured in HGSOC tissues compared with controls. Moreover, Kaplan-Meier method indicated the prognostic significance of ETV7 and CXCL13.

Ethics statement

The patients samples were used with approval from the Ethics Committee of Fudan University Shanghai Cancer Center and with consent from all patients. All procedures were performed in accordance with the Declaration of Helsinki and relevant policies in China.

Public data processing and TMEscore calculation

Level 3 data of HGSOC patients with available clinicopathological and survival data were downloaded from the TCGA data coordination center (https://tcga-data.nci.nih.gov/tcga/). For validation, gene expression profiles and clinical survival data for OC patients were obtained from the ICGC data portal (Supplementary Table 7, https://icgc.org/).

Based on the expression profile and clinical data downloaded from TCGA, heteroscedasticity was removed using the edgeR voom algorithm, and the proportion of LM22 cells was calculated with CIBERSORT, which allows sensitive and specific discrimination of 22 human immune cell phenotypes, including B cells, T cells, natural killer cells, macrophages, dendritic cells (DCs), and myeloid subsets [14, 15]. Then, according to the obtained cell proportions, three unsupervised clustering methods were used to explore sample classification, and the differences in expression were analyzed. According to the differential gene expression values obtained, consistency clustering was conducted to explore sample classification, and a chi-square test of independence was used to determine whether the two classifications before and after the three methods were consistent. Via a random forest model, signature genes were obtained by removing the redundancy of differentially-expressed genes (DEGs) with consistent classification. A Cox regression model was used to detect the relationship between signature genes and sample survival and to classify the genes. Finally, the TMEscores associated with every sample were calculated according to the gene classification. Data processing and TMEscore construction methods are detailed in the Supplementary Methods.

DEGs associated with the TMEscore

DEG analysis was performed with the Limma package [37] according to the median TMEscore value. The samples were separated into high and a low TMEscore groups. An empirical Bayesian approach was applied to estimate the gene expression changes using moderated t tests. The adjusted p-values for multiple testing were calculated using the Benjamini-Hochberg correction. A |log(fold change)|≥1 and an adjusted p-value < 0.05 were set as the cutoff criteria to screen for DEGs.

Module analysis and PPI network construction

The PPI network was retrieved from the STRING database [38]. Molecular interaction networks were visualized using Cytoscape software. Only individual networks with 10 or more nodes were included in further analyses. The connectivity degree of each node in the network was calculated. MCODE was then used to find clusters based on topology to locate densely connected regions [39].

Pathway enrichment analysis for molecular function

To further understand the functions of prognosis-related genes, we performed pathway enrichment analysis with the genes based on the KEGG and GO databases, including biological process (BP), molecular function (MF) and cellular component (CC) categories [40]. The top-ranked pathways with an adjusted p-value < 0.05 were identified as statistically significant. All analyses were performed using the R package ClusterProfiler [41].

Survival analysis

Kaplan-Meier plots were generated to illustrate the relationship between patient overall survival and the levels of DEGs. DEGs were identified as binary variables (high vs. low) using the median expression value as the cutoff value. A log-rank test was used to assess differences between survival curves.

Cell culture and treatments

The human OC cell line A2780 was cultured in Dulbecco’s modified Eagle’s medium (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco), 50 U/mL of penicillin and 50 μg/mL of streptomycin (Gibco). A2780 cell line was maintained at 37 °C and 5% CO2 in a humidified atmosphere.

Cell transfections and western blot analysis

SiRNAs targeting GBP1, ETV7, CXCL13 were purchased from GenePharma (Shanghai, China). Lipofectamine 3000 (Invitrogen, USA) was used for transfection experiments. All siRNA sequences are listed in Supplementary Table 8.

Transfected cells were cultured for 60 h and were then lysed in RIPA buffer (Sigma-Aldrich, USA) supplemented with a phosphatase inhibitor (Roche) and a protease inhibitor (Roche, Basel, Switzerland). A BCA protein assay kit (Thermo Scientific, USA) was used to measure protein concentration. The collected cell lysates were separated on 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) gels (EpiZyme, Shanghai, China) and were then transferred to polyvinylidene fluoride (PVDF) membranes (Millipore, Billerica, MA), which were blocked with 5% milk for 1 h at room temperature. PVDF membranes were incubated with primary antibodies (anti-β-actin, anti-GBP1, anti-ETV7, anti-CXCL13 antibodies, purchased from Proteintech) overnight at 4 °C, washed, and incubated with a secondary antibody for 1 h at room temperature. Signals were detected using chemiluminescence. β-Actin was used as the endogenous loading control.

Cell proliferation assays

A total of 2×10³ cells per well were seeded in 96-well plates 24 h before the experiment. A2780 cells were transfected with the targeting siRNAs or scrambled siRNA. Proliferation was measured using a CCK-8 kit (BOSTER, China) according to the manufacturer’s protocol. All experiments were performed in triplicate. Cell proliferation curves were plotted using the absorbance at each time point.

Transfected cells were digested with trypsin into single-cell suspensions 48 h later. For the colony formation assay, a total of 500 or 1,000 cells were plated in 12-well plates and incubated in corresponding medium containing 10% FBS at 37 °C. After being cultured for two weeks, the cells were fixed with anhydrous alcohol for 30 min and stained with 0.1% crystal violet for approximately 15 min. Then, visible colonies were manually counted. Triplicate wells were analyzed for each treatment group.

Migration assays

A cell invasion assay was performed using Transwell chamber inserts (8.0 mm, Corning, NY, USA) in a 24-well plate. Then, 4×10⁴ cells suspended in 200 μL of serum-free medium were added to the upper chambers. Culture medium containing 20% FBS was placed in the bottom chambers. The cells were incubated for 36 or 48 h at 37 °C. After incubation, the cells on the upper surface were scraped and washed away, whereas the cells on the lower surface were fixed with 20% methanol and stained with 0.1% crystal violet. The numbers of invaded cells in five randomly selected fields were counted using a microscope. The experiments were repeated independently in triplicate.

Immunohistochemistry

Immunohistochemistry was performed using the Envision System with diaminobenzidine (Gene Tech Co., Ltd., Shanghai) according to the manufacturer’s protocol. HGSOC tissue microarrays (TMAs) were constructed by Shanghai Outdo Biotech. In brief, specimens were incubated first with an anti-GBP1 antibody (1:200, Proteintech, China) overnight at 4 °C and then with a biotinylated secondary antibody (1:100, goat anti-rabbit IgG) for 30 min at 37 °C.

RNA isolation, reverse transcription and quantitative real-time PCR

Total RNA was extracted from the tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. Reverse transcription (RT) and quantitative real-time PCR (qRT-PCR) kits (Takara, Dalian, China) were utilized to evaluate the mRNA levels of the indicated genes. PCR primers were designed and synthesized using a primer design tool (Vector NTI®); The primers used in this study were listed as follows: ETV7: 5′-CTG CTG TGG GAT TAC GTG TAT C-3′ (Forward) and 5′-GTT CTT GTG ATT TCC CCA GAG TC-3′ (Reverse); CXCL13: 5′-GCT TGA GGT GTA GAT GTG TCC-3′ (Forward) and 5′-CCC ACG GGG CAA GAT TTG AA-3′ (Reverse). The relative quantification value for each target gene was expressed as 2^−ΔΔCT. β-Actin was used as an internal reference.

Statistical analysis

RStudio software (Version 1.2.1335) and GraphPad Prism 8 were used to construct the plots shown in the figures. Statistical analyses were conducted using RStudio software (Version 1.2.1335) and SPSS 22.0 (SPSS, Inc., Chicago, IL, USA). A value of P<0.05 indicated a statistically significant difference.