Neoadjuvant Chemotherapy Modulates the Tumor Microenvironment in High-Grade Serous Ovarian Cancer


 Background

While surgical reduction with adjuvant chemotherapy is the traditional treatment for high-grade serous ovarian cancer (HGSOC), neoadjuvant chemotherapy (NACT) has increasingly been applied. This work aims to investigate the expression profiles before and after NACT, explore changes in the tumor microenvironment, expand current treatments, and design a combination of treatment options for patients.
Methods

We downloaded 326 pre-NACT RNA sequencing data and 37 matched pre- and post-NACT samples from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Differentially expressed genes (DEGs) were determined with EdgeR, and Gene Ontology analysis was performed to identify the clusters responsible for the biological processes and pathways of HGSOC. Immune infiltration was analyzed using Single-sample Gene Set Enrichment Analysis (ssGSEA) and CIBERSORT. Kaplan-Meier (KM) survival analysis was performed to assess prognosis, and the potential correlations between modules and phenotypes were explored using weighted gene co-expression network analysis (WGCNA).
Results

After NACT, a total of 352 genes showed significant changes in RNA expression, among which 180 genes were up-regulated and 172 down-regulated. The most influential pathway was the positive regulation of mitogen-activated protein kinase (MAPK) cascade. Correlation analysis and KM survival analysis showed that overexpression of MAPK cascade genes correlated with shorter survival time in HGSOC patients. ssGSEA results showed that the expressions of anti-tumor cells (central memory CD4+ T cell and central memory CD8+ T cell) and pro-tumor cells (neutrophil and dendritic cells) were significantly increased after NACT. CIBERSORT showed that the abundances of memory B cells, NK cells, and monocytes were increased and the abundance of plasma cells was decreased after NACT. WGCNA and KM survival analysis showed that a lower abundance of Regulatory T cells (Tregs) was correlated with a better prognosis.
Conclusions

Gene expression of the MAPK pathway is up-regulated and the abundance of CD4+ T regulation cell decreases after NACT. Thus, the MAPK pathway may promote the differentiation of CD4+ T cells into Th17 cells while inhibiting Tregs development. The inhibited Tregs' development can lead to a better prognosis. Therefore, it is speculated that Tregs inhibitors combined with platinum-based NACT are potential treatment options for HGSOC.


Introduction
Ovarian cancer is the eighth leading cause of cancer deaths in women globally, with 295,414 new cases and 184,799 deaths in 2018 all over the world 1 , of which high-grade serous ovarian cancer (HGSOC) accounted for about 75% 2 . The standard treatment of HGSOC is debulking surgery followed by platinum/paclitaxel chemotherapy. After treatment, platinum-resistant cancer recurrence happens in approximately 25% of patients within 6 months 3 , and the 5-year overall survival rate is about 30% 4 . Recent treatment shifts towards increasing use of neoadjuvant chemotherapy (NACT) to reduce the burden of disease before surgical cytoreduction 56 . NACT is associated with an increased rate of optimal debulking and debulking to no-residual disease as well as a reduction in surgical morbidity and mortality 7 . Chemotherapy has been shown to affect gene expressions and cause molecular derangement over time 89 . The recent transition towards NACT for HGSOC provides an opportunity to assess the molecular response of individual patients to the treatment.
Although HGSOC is the most common histological subtype of epithelial ovarian cancer, many studies have shown that the tumor has signi cant heterogeneity 10 . However, it was still treated as a single entity in the previous research. TCGA analyzed the expression pro les of 489 HGSOC tumors and found that homologous recombination was defective in about half of the tumors 11 . Notch and FOXM1 signaling pathways are involved in the pathophysiology of serous ovarian cancer. Due to the high degree of heterogeneity and low overall operational mutation rate, it is almost impossible for a single treatment regimen being effective in all HGSOC patients. We need to improve the selection criteria for rst-line treatment and nd more successful treatments for cases. To this end, a closer link between molecular mapping and before and after NACT treatment is needed. Besides, because continuous tissue biopsy is often not feasible for patients, we need better strategies to assess genomic changes.
We hypothesized that changes in gene expressions in the tumor microenvironment (TME) before and after NACT may be decisive factors for relapse because they determine uncontrolled proliferation and resistance to treatment. Based on the molecular characteristics of before and after NACT, new therapeutic strategies for cancer progression can be developed. So far, before and after NACT HGSOC samples have not been systematically characterized. To investigate the molecular characteristics of recurrent HGSOC, TME, and its major counterparts during the treatment, we collected RNA sequencing data before and after NACT in HGSOC patients from TCGA and GEO 111213 and performed the differential expression analysis, Gene Ontology (GO), Single-sample Gene Set Enrichment Analysis (ssGSEA), CIBERSORT analysis, weighted gene co-expression network analysis (WGCNA), and Kaplan-Meier (KM) survival analysis. To explore the changes of TME before and after NACT and its relationship with prognosis, we explored the potential differences in the mitogen-activated protein kinase (MAPK) pathway and regulatory T cells (Tregs).

Data source
As of August 28, 2019, we downloaded all HGSOC RNA sequences and clinical data from https:// rebrowse.org. A total of 326 samples had RNA sequencing data, of which comprehensive survival data were available in 291 people (Supplemental 1). Thirty-seven matched before and after NACT samples were downloaded from GEO (GSE109934 and GSE71340). All of patients in the GEO database received 2 -6 cycles of neoadjuvant chemotherapy and underwent cytoreductive surgery.

Differential expression genes analysis
Three series matrix les were annotated with an o cial gene symbol using the data of RNA sequencing and microarray platform, and then gene expression matrix les were obtained. Gene expression matrix les were merged into one le, and the "sva" R package was used to conduct batch normalization of the expression data from the three different datasets. Finally, a normalized gene expression matrix le containing data from these three different datasets was obtained for differentially expressed gene (DEG) analysis. The "EdgeR" R package was used to conduct DEG analysis 14 . The threshold of DEGs was set as foldChange > 2 and adjusted P-value <0.05. To visualize the gene expression pattern, we generated a heat map using the ComplexHeatmap package in R v3.6.1.

GO enrichment analysis of DEGs
The "clusterPro ler" R package was used for GO enrichment analysis 15 . Enriched GO terms with adjusted P-value < 0.05 were considered as statistically signi cant. GO enrichment analysis used priori gene sets that have been grouped by their involvement in the same biological pathway or by proximal location on a chromosome. The potential statistically signi cant differences between the important pathways before and after NACT were analyzed. Histograms of enriched GO terms were implemented with the ggplot2 package (https://ggplot2.tidyverse.org/) in an R language environment.

CIBERSORT analysis
CIBERSORT was used to analyze the fractions of immune cells before and after NACT, during which an input matrix of gene expression signatures was used to calculate the relative proportion of each target cell type 16 . The software deconvoluted the mixture by linear support vector regression (SVR) and machine learning methods.

Co-expression network construction by WGCNA
Pearson's correlation coe cient (PCC) was used to assess the relationships between each pair of the 5000 genes by WGCNA 17 . These data were used to construct an unsupervised co-expression-based adjacency matrix with a soft threshold power of 4 based on the scale-free topology criterion to raise the matrix to simulate a realistic network structure. The intramodular connectivity (k.in) measured how connected, or co-expressed, a given gene was concerning the genes of a particular module. The connection strengths were assessed by calculating the topology overlap (TO), and the modules were de ned as sets of genes with a high TO. The topological overlap matrix (TOM) was the network distance measure for each gene pair from the adjacent matrix. A TOM-based dissimilarity measure (1-TOM) was utilized to achieve an average hierarchical linkage clustering. Gene modules were de ned using a dynamic hybrid branch-cutting algorithm with a cutoff of 0.95 and a minimum module size and cutoff of 30 in the hierarchical clustering dendrogram based on TOM dissimilarity. The module eigengene (ME) was calculated by a principal component analysis (PCA) by de ning the rst principal component of a given module. The module membership, also known as eigengene-based connectivity kME, related each gene expression pro le with the ME of a speci c module. The MEs of a summary pro le were used to assess the underlying correlation of gene modules with the clinical variables and survival.

Survival analysis
Survival analysis was performed using the survival R package with the hazard ratio (HR) and its corresponding 95% con dence interval (CI) determined by the Cox regression module and Kaplan-Meier survival analysis (http://cran.r-project.org/web/packages/survival/index.html). Overall survival (OS) or disease-free survival (DFS) were considered to be the survival endpoints. For a given ME/gene, the patients were split into high expression (≥median expression of the ME/gene) and low expression (<median expression of the ME/gene) groups.

Results
Screening of DEGs in pre-and post-NACT samples Gene expression was evaluated from RNA sequencing data and microarray using "EdgeR" R package.

GO and immune in ltration analysis of DEGs
Among the tumor-related pathways affected by the gene expression changes the positive regulation of MAPK cascade was the most-affected one, followed by regulation of protein serine/threonine kinase activity, peptidyl-tyrosine phosphorylation, peptidyl-tyrosine modi cation, and epithelial cell proliferation. Fig.2 shows a histogram of the displayed pathway for NACT effects.
GSEA was performed to compare changes in immune cells in pre-and post-NACT samples. The results showed that after NACT, the abundances of anti-tumor cells (central memory CD4 + T cell) and pro-tumor cells (neutrophil and plasmacytoid dendritic cell) were signi cantly increased (Fig.3A). Further analysis using CIBERSORT revealed that after NACT, the abundances of memory B cells (p = 0.099), NK cells (p = 0.001), and monocytes (p < 0.001) were signi cantly increased, the abundance of plasma cells was signi cantly decreased (p < 0.001) (Fig.3B).

MAPK pathways genes correlated with shorter survival time in HGSOC patients
According to the correlation coe cient, the nine most relevant genes in the MAPK pathway were selected: MAPK10, NTRK2, MAP3K5, PDGFRA, TGFBR2, FGF7, SOCS2, GADD45B, and MAP3K8 (Fig.4A). KM survival analysis was performed on these nine genes with TCGA HGSOC RNA sequencing data. The results showed that these nine genes signi cantly affected the prognosis of HGSOC patients (Fig.4B). In addition to MAP3K8, survival analyses of other genes indicated that the expressions of these genes were negatively correlated with the prognosis of HGSOC (Fig.4C-I). The survival analysis of total expression value also showed a negative correlation with prognosis.
A lower abundance of Tregs was correlated with a better prognosis WGCNA showed Tregs was highly correlated with the MEblue module (correlation coe cient> 0.4) (Fig.5A). In the MEblue module, each gene was highly correlated (cor = 0.8, P-value <0.05) (Fig.5B). We further screened the genes in the MEblue module and found the 12 most important genes: GFRA2, CLEC16A, RXRA, AVPR1A, MT3, REEP4, PLOD2, STAT1, NOL3, C1orf106, SIGLEC11, and GJB1. Three of these genes were related to the proliferation and differentiation of Tregs. We performed principal component analysis on the 12 most important genes and the 3 genes (GFRA2, RXRA, and STAT1), and calculated the expression values of these genes based on the contribution of each gene. Analysis of the expression values by the KM survival analysis showed that the reduction of these genes was related to the longer survival time of HGSOC patients (Fig.5C), which implied that a lower abundance of Tregs was correlated with a better prognosis.

Discussion
This study describes the change of gene expression and tumor microenvironment of pre-and post-NACT HGSOC patients. After NACT, altered expressions were detected in 352 genes. These genes were enriched for MAPK, serine/threonine kinase activity, and peptidyl tyrosine phosphorylation signaling pathways. Cisplatin used in NACT is a typical DNA damaging agent that can induce cross-linking between DNA and internal DNA strands and is known for its ability to induce apoptosis. Studies have con rmed that the addition of cisplatin to human cell lines activates the p38 MAPK pathway 21 . Inhibition of the MAPK pathway results in decreased expression of ATF3 messenger RNA and decreased cytotoxicity of cisplatin 22 . This explains why the gene expression involved in MAPK has changed. The top two upregulated genes were IL24 and NODAL. IL24 is known as an anti-cancer gene, and a previous study has shown that IL24 can be combined with cisplatin to enhance tumor cell death 23 . NODAL has been con rmed to induce apoptosis and inhibit its proliferation 24 . Thus, NACT has a positive effect on the antitumor treatment of HGSOC patients.
Analysis of immune in ltration found that after NACT, the expressions of anti-tumor cells (central memory CD4 + T cell, central memory CD8 + T cell) and pro-tumor cells (neutrophil and dendritic cell) were signi cantly increased. The abundances of memory B cells, NK cells, and monocytes were increased, and the abundance of plasma cells was decreased. GFRA2, RXRA, and STAT1, which were found by WGCNA analysis, were highly correlated with Tregs. In cancer, Tregs inhibits the anti-tumor immune response and contributes to the development of the TME, thereby promoting tumor invasion and cancer progression 2526 . Higher expressions of Tregs-related genes are associated with poor prognosis in many tumors, including ovarian cancer 2728 , pancreatic ductal adenocarcinoma 29 30 , lung cancer 31 , and glioblastoma 32 .
In our current study, CIBERSORT analysis showed that the abundance of Tregs decreased after NACT. The expression of the MAPK pathway was signi cantly increased in post-NACT patients, which was related to shortened survival time 33 . A possible explanation is that the MAPK pathway promotes the differentiation of CD4 + T cells into Th17 cells while inhibiting Tregs development. Survival analysis showed that inhibition of Tregs' development lead to a better prognosis. This may be because when the expression of Tregs is decreased, the inhibition of T cells is reduced, promoting the killing effect of T cells on tumors, thereby prolonging the overall survival of patients. Accordingly, targeted inhibition of Tregs development may promote the anti-tumor effect of the drug, thus improving the prognosis. Therefore, it is speculated that Tregs inhibitors combined with platinum-based neoadjuvant chemotherapy may be a potential treatment strategy for HGSOC.
Our current study was limited by few RNA-sequencing data in post-NACT patients, a lack of comprehensive analysis of genomic and proteomic data, and a lack of expression pro le data for relapsed patients after NACT. The above studies indicated that the NACT treatment activates the MAPK pathway, which had the effect of inhibiting the development of Tregs, and the decrease in the abundance of Tregs was associated with longer overall survival. Overall, we can draw two conclusions. 1. NACT treatment has a positive effect on anti-tumor therapy in patients with HGSOC. 2. Combined anti-Tregs therapy based on cisplatin chemotherapy may be bene cial to patients to prolong their overall survival.   Pathway changes in HGSOC before and after NACT A. GO enrichment analysis for the 352 genes expression changes after NACT identi ed multiple biological processes with FDR < 0.05. B. Clustering by the similarity of GO terms.