Unraveling tumor-immune heterogeneity in advanced ovarian cancer uncovers immunogenic effect of chemotherapy

Abstract

In metastatic cancer, the degree of heterogeneity of the tumor microenvironment (TME) and its molecular underpinnings remain largely unstudied. To characterize the tumor-immune interface at baseline and during neoadjuvant chemotherapy (NACT) in high-grade serous ovarian cancer (HGSOC), we performed immunogenomic analysis of treatment-naive and paired samples from before and after treatment with chemotherapy. In treatment-naive HGSOC, we found that immune-cell-excluded and inflammatory microenvironments coexist within the same individuals and within the same tumor sites, indicating ubiquitous variability in immune cell infiltration. Analysis of TME cell composition, DNA copy number, mutations and gene expression showed that immune cell exclusion was associated with amplification of Myc target genes and increased expression of canonical Wnt signaling in treatment-naive HGSOC. Following NACT, increased natural killer (NK) cell infiltration and oligoclonal expansion of T cells were detected. We demonstrate that the tumor-immune microenvironment of advanced HGSOC is intrinsically heterogeneous and that chemotherapy induces local immune activation, suggesting that chemotherapy can potentiate the immunogenicity of immune-excluded HGSOC tumors.

Conflict of interest statement

COMPETING INTERESTS

A.S. is a current employee of, and owns stock in, Merck. D.S.C. is on two medical advisory boards and invested in two surgical companies but none of those are related to this research. J.S.R.-F. is a paid consultant of Goldman Sachs Merchant Banking, Paige.AI and REPARE Therapeutics, and member of the scientific advisory board with paid honoraria of Paige.AI and Volition Rxand an ad-hoc member of the scientific advisory boards of Roche Tissue Diagnostics, Ventana, InVicro, Genentech, Novartis, GRAIL and Roche, outside the submitted work. D.Z. reports personal/consultancy fees from Merck, Synlogic Therapeutics, Biomed Valley Discoveries, Trieza Therapeutics, Tesaro and Agenus out of the scope of the submitted work. M.L.M. has received honoraria from GSK, but not related to this research. E.S. is a co-founder and shareholder of Cambridge AI Heath, consultant for Amazon and received honoraria from GSK, but none of these are related to this research. G.M. is founder and CTO of Pinpoint Oncology Ltd. A.J.-S. has owned and sold stocks while this work was in progress, none of which relates directly with this publication.

Figures

Fig. 1 |. Immune-related gene signatures contribute to the majority of the transcriptional variance observed across multiple tumor samples from treatment-naive HGSOC patients.

a, Presence, absence and replicate indication of samples and data types for treatment-naive samples (Cohort I). Metastases other than omentum were defined as “Other”. Samples from the same tumor are indicated with a connecting horizontal line (Supplementary Fig. 2a). Pseudoreplicates are samples from the same tumor and habitat, but from a different region within the habitat (see Methods). Age, age at diagnosis; BRCA, BRCA1/2 mutation status (Neg, negative; NA, data not available); WES, whole exome sequencing; IF, immunofluorescent staining. Extended clinical data can be found in Supplementary Table 1a,b. b, Flowchart of sample acquisition and analysis, c, t-SNE analysis of overall transcription profiles of multiple HGSOC tumor samples per patient, d, PCA of ssGSEA-based analysis of cancer hallmark gene sets, e, Principal component feature loadings (magnitude and direction) of c are shown in the variables factor map. Vectors are colored according to a major biological classification of cancer hallmark gene sets. Variation across classes in the PCs is shown in Supplementary Fig. 1c. Directionality of ESTIMATE’S tumor cellularity is represented with the map compass (n = 38 samples from n = 8 independent patients). Inset shows 95% non-parametric bootstrap confidence intervals for the means of loadings per hallmark gene set for PC1 (n = 36 samples from n = 8 independent patients). Approximated bootstrap P-values were calculated (see Methods) *P < 0.05, **P < 0.01 (immune vs. oncogenic P ≈ 0.01, stroma vs. oncogenic P ≈ 0.008). f, ESTIMATE immune score across patients and samples. The Case Study samples were taken from ref. 7. Box plots show median, interquartile range (25th and 75th percentiles) and 1.5x interquartile range. All samples are plotted. Abdomen image by Wenjing Wu/© 2018 Memorial Sloan Kettering Cancer Center.

Fig. 2 |. T cell infiltrate variation across patients, within patients, and within tumors.

a, Multi-tumor sampling from 8 HGSOC patients are shown with each dot representing the percentage of T cell subsets in a quantified area within a given tumor section stained with multicolor IF for CD8, CD4 and FOXP3. Stromal areas were excluded based on H&E stains. Patient cases are indicated by different colors. Imaged-based phenotypic habitats are defined by the Greek letters α, β and γ. Habitats from the same tumor sample are indicated by connecting lines at the bottom (see Supplementary Fig. 2a for detailed examples). Box plots are sorted according to the median of CD8 T cell infiltration across patients, sites and habitats. Box plots show median, interquartile range (25th and 75th percentiles) and 1.5x interquartile range. All samples’ regions are plotted. The two pseudoreplicates in case 5 are indicated with a diagonal line in the tumor site symbol (n = 440 observations of n = 8 independent patients with 38 total samples, see Supplementary Table 2). b, Representative images of a.

Fig. 3 |. Unbiased analysis of tumor microenvironment heterogeneity in treatment-naive HGSOC tumors.

a, Treatment-naive cohort correlation between (i) total TME cell estimation scores and WES-derived tumor cellularity (TITAN), and (ii) percentage of CD8+, CD4+, and Tregs determined by immunofluorescent staining, b, TCGA ovarian cancer cohort correlation between total TME cell estimation scores for each method and WES-derived tumor cellularity (ABSOLUTE), and fitted multiple linear regression analysis using TCGA leukocyte methylation score as response variable and estimated immune cell types as explanatory variables (see Supplementary Note and Supplementary Table 3a). Adjusted R2, Akaike information criterion (AIC), and Bayesian information criterion (BIC) values were calculated to compare both goodness of fit and model simplicity. Kendall’s tau correlation coefficients and P-values were calculated for a and b with exact Kendall’s tau-b two-sided test. See below for significance levels. Normality and homoscedasticity assumptions were tested for all statistical comparisons, c, PCA of ssGSEA-based analysis using ConsensusTME estimations (n = 38 samples from n = 8 independent patients), d, Principal component feature loadings (magnitude and direction) of c. Vectors are colored according to cell-types as shown in Supplementary Fig. 3a, for example monocytes and macrophages M0, M1, M2 (orange), B cells and plasma cells (light blue), and CD8 and cytotoxic cells (yellow), e, Differential expression analysis of high and low tumor cellularity classified tumors using the WES-derived tumor cellularity (TITAN) score median of the cohort as a cutoff. Patient dependence was used as a covariate (n = 36 samples which have both WES and mRNA data, n = 8 independent patients). Limma moderated two-sided t-statistics by empirical Bayes moderation. Benjamini-Hochberg (BH) FDR corrected. FDR < 0.05 was considered differentially expressed, f, Gene ontology enrichment analysis of 28 significantly highly expressed genes in low tumor cellularity samples. Two-sided Fisher’s exact test, BH FDR corrected (Supplementary Table 3d). g,h, ssGSEA analysis of differentially expressed genes using hallmarks and ConsensusTME normalized enrichment scores (NES). Gene sets on the x-axes were ranked according to their NES (Supplementary Table 3b,c). High NES reflects high tumor cellularity. Dashed red lines indicate median and ±1.96 median absolute deviations (modified z-score). Kernel density plots of observed and fitted normal distribution are shown in the right margin. No significant difference between observed and fitted distribution was detected (Shapiro-Wilk test, D’Agostino’s K2 test, and Anderson-Darling test for normality distribution). ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Fig. 4 |. Mutation patterns in immune-excluded tumors associate with Wnt pathway genes and Myc target genes.

a, Coding mutations and copy-number alterations (CNAs) in selected driver oncogenes and tumor suppressors based on known drivers in HGSOC. Tumor-normal DNA pairs from n = 10 independent patients, 47 independent tumor site samples (and 3 pseudoreplicates), were sequenced (whole-exome sequencing, WES) and analyzed for copy-number alterations using TITAN, resulting in estimates of tumor ploidy, cellularity, and subclonality. The top and side bars represent the summed column-wise and row-wise number of alterations, respectively. *, †, and ‡ indicate pairs of pseudoreplicates (see Fig. 1a). b, DNA copy-number signatures analysis showing tumor-specific exposure to all seven signatures ordered by tumor cellularity (TITAN) selecting samples with high confident purity and ploidy estimates (n = 42 samples, n = 10 independent patients), c, Copy-number signature 4 association with immune score. Using samples with both mRNA and WES data (n = 30 samples), ESTIMATE immune score was compared (Welch’s two-sided t-test) between samples with signature 4 exposure below (low s4) and higher than the median exposure (high s4). Boxplot with median midline, boxes representing the 1 to 3 quartiles and whiskers extending to extreme values at most 1.5 times the interquartile range. d, Functional mutation enrichment and tumor cellularity associations per hallmark gene set were tested through the Chi-squared likelihood ratio test using the tumor mutation load and the patient dependency as covariates. The difference in mean pathway mutation ratio between samples with high cellularity versus samples with low cellularity (WES-derived by TITAN, n = 50 samples) is plotted against BH FDR corrected P-values (Benjamini, Hochberg, and Yekutieli) from multiple linear regression (tumor cellularity~mutation ratio, n = 50 samples), e, Gene set enrichment analysis (GSEA) of gene-level correlation between absolute copy-numbers and gene expression (CNA~mRNA, per-gene median n = 36 samples), compared to enrichment in correlations between gene expression and tumor cellularity (mRNA~tumor cellularity, n = 36 samples), both estimated by means of Spearman’s rank correlation method (see Methods). BH FDR correction was performed, adjusting for the number of hallmark terms. Normality and homoscedasticity assumptions were tested for all statistical comparisons.

Fig. 5 |. Unbiased signaling pathway and TME cell decomposition analysis of chemotherapy treated HGSOC site-matched and unmatched tumor samples.

a, Presence, absence and replicate indication of samples and data types for the pre and post-chemotherapy matched and unmatched samples (Cohort II). Extended clinical data can be found in Supplementary Table 5. b, Flowchart of sample acquisition, clinical study design, and analysis, c, t-SNE analysis of overall transcription profiles of multiple HGSOC tumor samples per patient. d,e, PCA and principal component feature projections (magnitude and direction) of ssGSEA-based analysis of hallmark gene sets and ConsensusTME cells. Insets show paired comparison of pre- and post-treatment samples PC1s. Inset violin plots represent the full probability density of the data. Paired samples are connected with a line. Two-sided paired t-tests were conducted for hallmark, while Wilcoxon signed rank tests were conducted for ConsensusTME gene sets PC1 comparison. *P < 0.05. Normality and homoscedasticity assumptions were tested for all statistical comparisons. P-values were not corrected for multiple testing since the maximum number of tests within analysis was two (c). Pelvis image by Wenjing Wu/© 2018 Memorial Sloan Kettering Cancer Center.

Fig. 6 |. Chemotherapy-induced enrichment of NK cells evident in site-matched samples and is supported by preclinical data.

a, Pre/post NACT paired comparisons of hallmark gene sets and ConsensusTME inferred cells. Paired t-test, Welch’s t-test or Wilcoxon’s signed-rank test (all two-sided) were calculated according to the samples’ distribution and variance (see Methods). BH FDR P-value corrections were computed, b, Left histograms, multivariate analysis of cell type combinations associated with NACT (two-sided Hotelling’s T2 permutation test) and, right, notched box plots comparing pre- and post-NACT for sums (+) of cell scores: cytotoxic cells (i.e. signature of cytolytic activity from ConsensusTME), CD8+ T cells or NK cells (multiple linear regression without interaction, see Supplementary Fig. 6a). n = 8 independent patients; top row corresponds to n = 18 matched samples, bottom row to n = 38 unmatched samples. Notched box plots show median and interquartile range (25th and 75th percentiles), with paired data points indicated by a connecting line. Violin plots are shown in the background representing the probability density of the data. c, Flow cytometry analysis of homogenized tumors after UKP10 intraperitoneal inoculation of 8 weeks old C57/BL6 mice and treatment with 2 mg of cisplatin dissolved in 1 ml of PBS or 1 mL PBS as control as indicated. Analyzed immune cell types were CD8+ T-effector cells, NK1.1+ NK cells, CD4+ T-helper cells, FOXP3+ T-regulatory cells, CD19+ B cells, CD11b+ myeloid cells. Granzyme B (GrB) expression was used as a proxy for activity state of cytotoxic cells. Two-sided two-sample independent t-tests were conducted for CD8 T cells and two-sided Welch’s tests were conducted for NK cells. d, Similar to c, except peritoneal fluids of intraperitoneally inoculated ID8 cells were analyzed. CD4+ T-helper cells, FOXP3+ T-regulatory cells, CD19+ B cells, CD11b+ myeloid cells are shown in Supplementary Fig. 6f,g. Two-sided two-sample independent t-tests were conducted for CD8 T cells and NKGrB+ out of NK cells, while two-sided Welch’s tests were conducted for the rest of NK cell comparisons. Violin plots are shown in the background representing the full probability density of the data. All samples are plotted. Normality and homoscedasticity assumptions were tested for all statistical comparisons.

Fig. 7 |. Oligoclonal expansion of T cells and enrichment of shared TCRs after chemotherapy.

a, Comparisons of percentage of productive T cells (top), TCR productive clonality (middle), and maximum productive TCR frequencies (bottom) between pre- and post-NACT site-matched and site-unmatched samples. Paired (matched samples) parametric (TCR clonality and TCR freq.) and non-parametric Wilcoxon (% T cells) two-sided t-tests were computed. Unpaired (unmatched samples) non-parametric Wilcoxon (% T cells, TCR clonality and TCR freq.) two-sided t-tests were computed. TCR clonality is expressed as 1-entropy with values near 1 representing samples with one or a few predominant TCR rearrangements, while values near 0 represent more polyclonal samples. Notched box plots show median and interquartile range (25th and 75th percentiles), with paired data points indicated by a connecting line. Violin plots are shown in the background representing the full probability density of the data. b, Shared and unique TCR amino acid sequences between pre- and post-NACT site-matched and site-unmatched samples. Chi-squared test of independence of variables was conducted. c, Distributions of shared TCR amino acid sequences between patients pre- and post-NACT samples. One-way chi square test was conducted. d, Number of shared and unique TCRs pre- and post-NACT in site-matched and site-unmatched samples (top) and their productive frequencies (bottom). Notched box plots show median and interquartile range (25th and 75th percentiles), 1.5x interquartile range (top) and outliers. Violin plots are shown in the background representing the full probability density of the data. Widths of box and violin plots are proportional to the number of samples (top) and TCRs (bottom). Friedman ranking tests followed by Nemenyi post-hoc tests with associated adjusted p-values (q) were conducted. P-values were not corrected for multiple testing since the maximum number of tests within analysis was six (a). Normality and homoscedasticity assumptions were tested for all statistical comparisons.