Similarity and diversity of the tumor microenvironment in multiple metastases: critical implications for overall and progression-free survival of high-grade serous ovarian cancer

The tumor microenvironment is pivotal in influencing cancer progression and metastasis. Different cells co-exist with high spatial diversity within a patient, yet their combinatorial effects are poorly understood. We investigate the similarity of the tumor microenvironment of 192 local metastatic lesions in 61 ovarian cancer patients. An ecologically inspired measure of microenvironmental diversity derived from multiple metastasis sites is correlated with clinicopathological characteristics and prognostic outcome. We demonstrate a high accuracy of our automated analysis across multiple sites. A low level of similarity in microenvironmental composition is observed between ovary tumor and corresponding local metastases (stromal ratio r = 0.30, lymphocyte ratio r = 0.37). We identify a new measure of microenvironmental diversity derived from Shannon entropy that is highly predictive of poor overall (p = 0.002, HR = 3.18, 95% CI = 1.51-6.68) and progression-free survival (p = 0.0036, HR = 2.83, 95% CI = 1.41-5.7), independent of and stronger than clinical variables, subtype stratifications based on single cell types alone and number of sites. Although stromal influence in ovary tumors is known to have significant clinical implications, our findings reveal an even stronger impact orchestrated by diverse cell types. Quantitative histology-based measures can further enable objective selection of patients who are in urgent need of new therapeutic strategies such as combinatorial treatments targeting heterogeneous tumor microenvironment.


INTRODUCTION
Ovarian cancer is the most fatal gynecological malignancy [19,36]. Approximately 75% of patients have metastatic diseases at diagnosis [15], and survival rate for International Federation of Gynecology and Obstetrics (FIGO) stage III-IV ovarian cancer ranges from 25% to 37% [43]. The vast majority (70-80%) of ovarian cancer deaths are from high-grade serous carcinoma (HGSOC) with hardly any change in overall survival in the past years [6]. This is partly due to the distinctive biology of HGSOC: the lack of an anatomical barrier in the peritoneal cavity to cancer cell dissemination likely originated from the fallopian tube or precursor cells in the ovary [6,22,24,27,37,39].
To control this highly metastatic disease, concerted efforts have led to a new understanding of genetic heterogeneity in ovarian cancer and developments of

Research Paper
new targeted therapies and chemotherapy agents [2,4,5,18,30,33,48,49]. However, responses to therapies are often short-lived, and patients eventually relapse [51]. Understanding the contribution of the microenvironment to ovarian cancer progression will lead to the development of new biomarkers and provide the basis of improved therapeutic efficiency [4,20,29]. The desire to understand cancer-microenvironmental interactions has fuelled a developing interest in studying cancer from a novel perspective: ecology [3,17,23,32,40]. Seeing cancer cells as a 'species' disrupting the homeostasis of a complex ecosystems rather than just as a group of mutated cells [3] allows for drawing analogies to ecological hypothesis such as that complete eradication of a species is most likely impossible due to their ability to evolve adaptive strategies [17]. Another prominent example inspired by ecology is the "seed and soil" theory saying that in order to metastasize, the soil (microenvironment at the metastasis site) is as important as the seed (malignant cell) [34,41]. Understanding and early interventions of the primary "seeds" and metastatic "soil" are required for better clinical management of cancer. Within a Darwinian framework, such analyses can reveal distinct microenvironments, or, the habitats of cancer, that can employ a number of different ecological forces driving tumor development and spread [14]. Studies of these ecological forces occurring in tumors can benefit from application of statistical tools routinely used in ecological studies. Histology samples provide an abundance of data as input for these methods due to preserved spatial context. Thus, spatial analysis empowered by large-scale analysis of archival histology samples could facilitate studies of ecological interactions in human tumors with far-reaching implications [23,32].
Our previous work [25] demonstrated how fully automated image analysis using hematoxylin and eosin (H&E) slides of ovary tumors can enable the identification of ovarian cancer subtypes that were consistent with molecular subtyping previously reported [26,48]. These subtypes were defined based on proportions of cells in histological ovary tumor sections alone: a high lymphocyte ratio group with good prognosis and a high stromal ratio group with poor prognosis [25]. However, different types of cells interact and co-exist with high spatial heterogeneity within a patient, yet the collective effect of cell diversity has not been studied. Besides a recent paper [12] where the authors observed large variations in intracase comparisons of stromal characteristics between ovary and metastases, the microenvironment of tumors other than those from the ovary has received little investigation. In this paper, we apply ecologically inspired diversity measurements based on Shannon entropy to assess the heterogeneity of microenvironments in multiple local metastasis sites of HGSOC, including ovary, omentum, peritoneum, lymph node and appendix. Our aims are: 1) to establish the use of automated histology analysis for studying multiple metastasis sites besides ovary tumors, 2) to evaluate the similarity and ecological diversity of the microenvironment in these local metastases, 3) to understand the clinical implication of microenvironmental characteristics in HGSOC metastases.

Automated analysis of tumor microenvironment in local HGSOC metastases
We previously reported a high accuracy of our automated histological analysis system in identifying three major cell types (cancer cells, lymphocytes and stromal cells) in whole-section H&E histology samples of ovary tumors [25]. This system enabled the identification of ovarian cancer subtypes that were consistent with molecular subtyping previously reported [26,48]. These subtypes were defined based on proportions of cells alone in the ovary tumors: a high lymphocyte ratio group with good prognosis and a high stromal ratio group with poor prognosis [25]. Here, we evaluated the accuracy of our system in analyzing H&E sections of tumors from other local metastases including omentum, peritoneum, appendix, spleen, umbilicus and lymph node using two approaches (Figure 1, Methods). First, we tested the accuracy of single-cell classification using a collection of 4,633 single cells annotated by two pathologists that were not used to train the classifier (DNR and KK, Methods, Figure 2a). We found a good overall performance of our classifier compared with the pathologists' annotations with high mean balanced sensitivity and specificity scores across sites (cancer cells 0.90, lymphocytes 0.82 and stromal cells 0.89, Figure 2b). Stromal cell classification in peritoneum has the lowest balanced average measure due to a higher frequency of tissue contraction artifacts as tissue gaps originated from formalin fixation and paraffin embedding procedures. Secondly, we measured the correlation between automated and pathologist's cell proportion scores and observed a good level of correlation between automated and pathologist's scores of cell proportions (cancer: r = 0.8, lymphocytes: r = 0.7, stromal cells: r = 0.9, Figure 2c). The presence of a small amount of normal epithelial cells in the appendix was responsible for misclassification as cancer cells, resulting in a lower correlation with the pathologist's cancer cell proportion score. However, we observed in general a low amount of normal epithelial cells and the overall accuracy of cancer cell classification is high (median = 87%-96%, Figure 2b). We subsequently applied this automated image analysis tool on 192 local tumor lesions from 61 patients (Methods, Table 1). www.impactjournals.com/oncotarget

Comparison of microenvironmental composition between ovary tumors and other local metastasis sites
Based on the histology image analysis results, we first compared the cellular composition of ovary tumors by examining individual cell types. We found a low correlation between stromal ratio measured in the ovary and mean stromal ratio in other local metastasis sites of the same patient (Pearson correlation r = 0.302, p = 0.03, Figure 3a, Methods). Consistent with our previous study [25], the high stromal ratio group, defined by dividing patients into two equal-size groups based on the ratio of stromal cells to all cells in ovary tumors, was associated with poor overall survival (p = 0.03, HR = 1.99, 95% CI = 1.07-3.72, Figure 3b; Supplementary Figure 1b) However, stromal ratio defined in any local site other than the ovary, or the mean stromal ratio of those sites were not prognostic (p > 0.05, Table 3, Figure 3c). A moderate correlation between lymphocyte ratio measured in ovary tumor and mean lymphocyte ratio in other metastasis sites of the same patient was also observed (Pearson correlation r = 0.369, p = 0.003, Figure 3d). Low lymphocyte ratio was associated with poor overall survival (p = 0.04, HR = 0.52, 95% CI = 0.28-0.97, Figure 3e; Supplementary  Figure 1a). Neither lymphocyte ratio defined in any local site other than the ovary nor the mean lymphocyte ratio of these sites was prognostic (p > 0.05, Table 3, Figure 3f).

Tumor microenvironmental diversity of ovarian local metastases is associated with poor overall and progression-free survival
To measure the collective characteristics of the microenvironment, we next examined tumor diversity measures that consider all cell types simultaneously. The Shannon diversity index that measures species diversity in an ecosystem was calculated for each tumor (Methods). A high score indicates diverse cell types with similar proportions, whilst a low score suggests dominance of one cell type. We first compared the distribution of Shannon diversity for each tumor site. Appendix  Figure 2a). Analysis of averaged cell ratios in these sites showed that lymph node has in general a higher proportion of lymphocytes, as expected (Supplementary Figure 2b). The appendix, on the other hand, has higher proportions of stromal cells and lymphocytes, which explained the high diversity score.
To measure the level of microenvironmental diversity for each patient, we calculated the mean of diversity scores in all tumors other than those from the ovary (MetDiv). By relating it to survival, we found that a high MetDiv score was associated with poor OS (p = 0.002, hazard ratio (HR) = 3.18, 95% confidence interval (CI) = 1.51-6.68, Figure Kruskal-Wallis test was used for statistical testing.  Table 3; Supplementary Figure 1d). Patients with diverse microenvironments in metastases have 9% 5-year OS and 0% PFS compared with 42% and 20% respectively for patients with more homogenous microenvironments. The stability of these survival analyses was evaluated by random sampling of progressively fewer samples from 100% (n = 61) to 70% (n = 42). Mean metastasis diversity remained prognostic over 85% of the time for univariate analysis (OS; Figure 4c). This confirms the stability of this score as a predictor for survival. In contrast, measures of mean diversity scores in all tumors including those from the ovary, microenvironmental diversity in ovary tumors, or variability within local metastases (standard deviation) were not found to be prognostic (Table 3, p > 0.05). In addition, we compared our scores with another widely used ecological diversity index, the Simpson index (Methods). A high correlation between the two indices (r = 0.98, Supplementary Figure 3, Methods) was observed, which is consistent with previous studies that found very little qualitative differences between these scores [28,35]. We thus focused on the Shannon diversity index and MetDiv henceforth. Patient stratification based on MetDiv with respect to three cell types is illustrated in Figure 4de. A shift in microenvironmental features from the low diversity group to the high diversity group is evident. Metastasis sites with low MetDiv scores contained a high amount of cancer cells, whereas those with highly diverse microenvironment shifted towards the center of the triangle plot due to a more diverse microenvironmental makeup.

Metastasis diversity is not dominated by a cell type or metastatic site and remains stable with reduced amount of tissue
Subsequently, we tested whether the prognostic effect of our MetDiv score was driven by dominant cell types. When there is a dominant cell type, the Shannon diversity is low. Therefore a possible explanation of our observation is that a large amount of a single cell type in the metastases is responsible for a low diversity score and good prognosis. We then examined the distribution of cell types in the low and high MetDiv groups. We found no significant difference in the proportions of cancer cells, lymphocytes or stromal cells between the two groups (p > 0.1, Kruskal-Wallis test, Figure 4f). We next asked if MetDiv was driven by a specific metastasis site. MetDiv scores of these patients did not differ according to the presence of omentum, peritoneum, lymph node or appendix metastasis (p > 0.05, Supplementary Figure 4ad). We next examined the distribution of MetDiv groups according to clinicopathologic parameters. MetDiv was not associated with FIGO stage, debulking status and age (p > 0.05) but chemotherapy regimen (p = 0.002, Table 1).
A factor that could influence the reliability of MetDiv score is the amount of tissue used for our analysis. To test this, we simulated the scenario that only half of a tumor was available and computed the MetDiv score based on the new data (Methods). We found that the difference between the original and the new Shannon diversity score was very small (standard deviation = 0.0016-0.0006; Supplementary Figure 5a). The prognostic value of MetDiv for OS and PFS based on the new scores remained significant (p < 0.005). We further evaluated the robustness of the MetDiv score by randomly sampling 100% / 75% / 50% of cells from the whole-section images (Supplementary Figure 5b), revealing an even smaller standard deviation of the diversity score (0.0003), thereby supporting the stability of the MetDiv score against variable amount of tissue used for analysis.

Metastasis diversity is prognostic independent of clinical parameters
We tested the prognostic value of clinicopathologic and prognostic scores including age, debulking status, chemotherapy treatment regime, staging, and our  Table 3). Random sampling was again performed for MetDiv in multivariate analysis of OS that include age, stromal and lymphocyte grouping (Methods). We found that with more than 80% (49) samples, more than 81.3% of the times MetDiv remained significant, supporting its stability as a prognostic marker (Figure 4g).

DISCUSSION
Influences from stromal endothelial cells, immune regulation and cancer-associated fibroblasts were known to be important factors for tumor progression with roles in angiogenesis, inflammation and phenotype switching [1,20,52]. Prognostic subtypes directly related to immune or stromal components in ovary tumors have been consistently reported in large-scale studies of genomics [16,33,48,49]. We have previously reported strong prognostic value of ovarian cancer subtypes based on the amount of immune or stromal infiltrate in ovarian tumor histological samples, where high amount of stromal cells was significantly associated with poor survival in contrast to lymphocytic infiltration which was found to correlate with favorable prognosis [25]. The most important finding in this study is that diverse microenvironment in multiple local metastases including omentum, peritoneal, lymph node and appendix is a novel feature of aggressive HGSOC. Its strong prognostic effect is independent of subtyping that consider a single cell type alone, microenvironmental diversity at a single tumor site, number of metastatic sites and standard clinical parameters in HGSOC. Therefore, our new study offers clear evidence that it is not just single cell type at play, but also that the collective contribution of diverse cells in the microenvironment parallels the acquisition of aggressive HGSOC phenotypes in locally advanced diseases.
Specifically, HGSOC patients with diverse microenvironment or high MetDiv score in their local metastases have a staggering 5-year overall survival of 9%, compared with 42% for patients with low MetDiv scores. Although there was a moderate degree of similarity in the amount of stromal cells between ovary tumors and other sites (Pearson correlation r = 0.302), only stromal cell abundance measured in ovary tumors but not in other sites was associated with poor overall survival in these locally advanced diseases (p = 0.03, HR = 1.99, CI = 1.07-3.72). Similar conclusions were drawn in our analysis comparing lymphocyte abundance in ovary and other sites. These results underscore the clinical significance of stromal and lymphocyte abundance in ovary tumors but not in other tumor sites. The MetDiv score is, however, independent of and stronger than stromal or lymphocyte abundance in multivariate survival analysis (MetDiv p = 0.002, HR = 3.40 95% CI = 1.59-7.27, age p = 0.01, HR = 2.72 95% CI = 1.27-5.83, stromal ratio p = 0.12, HR = 1.68, 95% CI = 0.88-3.21 and lymphocyte ratio p = 0.30, HR = 0.71, 95% CI = 0.38-1.34). Debulking status is known as a key clinical variable for advanced ovarian cancer [7,10,50], and was found to be prognostic in our previous study that included patients with only ovary tumor samples (OS and PFS p < 0.01) [25]. However, it was not found to be associated with OS or PFS in the patient subset with more than one local metastasis under study here (p > 0.05, Table 3). Reasons for the lack of significance can be either the small sample size or the now dated definition of optimal debulking (residual disease no larger than 1 cm in maximum diameter).
Note that the MetDiv score is different from the ecosystem diversity index in our recent study of primary breast tumors [31]. Ovarian cancer is unique in its high metastatic potential due the lack of an anatomical barrier to cancer cell dissemination in the peritoneal cavity [5]. MetDiv thus represents an ecologically inspired measure that was specifically designed for the spatial analysis of HGSOC microenvironments, using a collection of 192 local metastasis tumors spreading across 61 patients. Despite variations in background of different tumor sites and variable numbers of metastasis sites per patient, we showed that MetDiv provides a generalizable measure of cell diversity across multiple sites, emphasizing the importance of considering the cumulative effects of diverse cell types at multiple locations that can contribute to HGSOC progression. Further studies based on our findings may suggest novel solution to the problem of therapeutic resistance in HGSOC management [44,47] and lead to clinical innovations in treatment strategies. With recent advances in immunotherapy targeting immune checkpoints or cytokines, combinatorial therapies with these new agents and anti-VEGF drugs such as bevacizumab [9] to target multiple components may be carefully considered for patients with high MetDiv scores to mitigate the joint effect from multiple cell types. www.impactjournals.com/oncotarget A further study to illuminate the complex network of cytokine and chemokine signaling between malignant cells and normal cells intertwined with growth factors underlying the prognostic effect of metastasis diversity will shed further light on this aspect.
In summary, we demonstrated how automated histology analysis coupled with ecological measures is effective in characterizing the metastatic microenvironment in locally advanced HGSOC, leading to clinically significant findings. Advantages of our proposed method are that it utilizes objective automated histology analysis, widely available paraffin-embedded specimens and standard H&E staining methods, and we demonstrated that it is robust against variable amount of tissue ( > 50%) used for analysis. Limitations of our study include the lack of an independent test cohort, a twodimensional representation of complex tumors, no spatial association measured between cell types, two different treatment regimens, only one histologic subtype and only basic cell types (stromal cells, lymphocytes and cancer cells). Furthermore the definition of optimal surgery was refined over time and changed from < 1 cm to minimal/ none macroscopic residual disease [11]. Future efforts will focus on further validation of our automated system and the discrimination between different types of lymphocytes and stromal cells using immunohistochemistry approaches (e.g. fibroblasts, myofibroblasts, fibrocytes, vascular pericytes and endothelial cells) to further dissect the tumor microenvironment and consolidate their clinical relevance in HGSOC and other histologic subtypes as well as molecular profiling in independent cohorts with multiple sections per tumor.
Nevertheless, we expect that the need to understand the biological basis for coordinated influences from multiple microenvironmental cell types, and their concomitant effect on HGSOC treatment, will help shape future research on this devastating cancer type.

Patient selection and characteristics
All patients with ovarian carcinoma, who received primary debulking surgery without neoadjuvant therapy between May 1999 and December 2010 at Sun Yat-sen University Cancer Centre, China, were reviewed. A total of 61 patients with International Federation of Gynecology and Obstetrics (FIGO) stage III-IV HGSOC ovarian cancer with at least one locally advanced resectable metastasis site were identified (Table 1). Resectable metastases were collected during debulking surgery. In total, 192 paraffin embedded blocks from 7 local metastasis sites (61 ovary, 51 omentum, 48 peritoneum, 9 appendix, 20 lymph node, 1 spleen, 2 umbilicus; Table 2) were obtained from 61 patients. For each patient 2-5 tumor sites were available (median n = 3). Patient consent and ethical approval were obtained by the institutional review board of Sun Yatsen University Cancer Centre. Paraffin-embedded tumor blocks were retrieved from the hospital archive. Clinical data including clinicopathologic characteristics, treatment regimen and follow-up 10-year overall survival (OS) and progression-free survival (PFS) were collected (Table  1). OS was censored at the date of death or, for living patients, the date of last contact. PFS was censored at the date of death or progression, whichever occurred first, or the date of last contact for the patients alive and without recurrent disease. Progression was defined as serially rising CA125 levels, or any clinical or radiographic evidence of new lesions as either local/regional relapse or distant metastasis. Median OS and PFS were 43.33 months (range 1.97-120 months) and 19.8 months (0.97-120), respectively. All patients underwent debulking surgery that mainly consisted of total hysterectomy, salpingo-oophorectomy, omentectomy, nodal dissection, or intestinal surgery followed by chemotherapy. Optimal debulking was considered as residual disease no larger than 1 cm in maximum diameter evaluated by surgeons at surgery [8,38]. The majority of patients (88.6%) had 6-8 cycles of paclitaxel/docetaxel plus carboplatin/cisplatin postoperatively, and 11.5% of patients received adjuvant chemotherapy of cyclophosphamide, bleomycin plus carboplatin before paclitaxel was available in China (Table  1).

Histology image analysis and validation
H&E histological sections of multiple metastases from paraffin-embedded tumor blocks were generated, digitalized (Aperio, 20x, 0.5µm/pixel) and analyzed using open source R package CRImage [53] (Figure 1). The cell classifier was trained based on 100 features comprised of mostly morphological and some textural measurements as described in [25]. Cells were classified based on their morphological differences of the nucleus positive for haematoxylin stain without using immunohistochemical target stains. Immune cells typically display small, round and homogeneously basophilic nuclei; cancer cells in general have nuclei of larger size and greater variability in texture and shape. Both cell types can be differentiated from stromal cells such as fibroblasts, myofibroblasts, fibrocytes, vascular pericytes and endothelial cells that contain more elongated nuclei. The average amount of cells classified in each tumor section was 362,417 (interquartile range 208,129-502,789) for cancer cells, 92,861 (interquartile range 39,075-106,220) for stromal cells and 67,433 (interquartile range 27,944-90,852) for lymphocytes.
Validation of the classifier on tumors other than those from the ovary was performed using two approaches. First, a total of 4,633 single-cell annotations from two pathologists were collected from all six sites (cancer = 2,668, lymphocyte = 1,120, stromal = 845; Figure 2a, Supplementary Figure 6). Note that this data was only used for testing and not for training the classifier. Subsequently, sensitivity and specificity as well as the balanced average of both were computed (Figure 2b, Supplementary  Figure 7 and 8). The second approach was performed by comparing proportions of cancer cells, lymphocytes and stromal cells to total amount of cells in 24 randomly selected fields of view. The fields of view were 2000 x 2000 square pixel images at a resolution of 0.5µm per pixel, representing 6 different metastasis sites each with four fields of view. Automated image analysis results were benchmarked against scores from the pathologists (DNR, KK) using Pearson's correlation test.

Tumor microenvironmental diversity measures
We quantified the ecosystem diversity for a tumor j using the Shannon diversity index [45]: Eq.1 where p ij is the proportion of the i th cell type in tumor j, and R is the number of different cell types. A high value of the Shannon diversity index represents a heterogeneous environment populated by similar proportions of different cell types, whilst a low value indicates a homogeneous ecological composition. A known disadvantage of Shannon diversity is that it only accounts for species abundance and not absence [42]. We therefore also compared it with the Simpson diversity index [46]. The Simpson index of tumor j was defined as follows: Eq.2 where p ij is the proportion of the i th cell type in tumor j and R the number of different cell types.

Other statistical methods
Survival analysis was carried out using R. Survival curves were plotted according to Kaplan-Meier [21] method and the difference in survival was assessed using the log-rank test. Cox proportional hazards model [13] was used to identify prognostic factors in univariate and multivariate analyses. To be considered for multivariate analysis, a factor had to be statistically significant (p < 0.05) in univariate analysis. Continuous variables were dichotomized due to their skewed distribution (skewness ranging between -0.4 and 1) by employing either an iterative thresholding approach or if no optimal threshold could be found by dividing patients into two equal-size groups based on the median. Lymphocyte ratio was dichotomized into a group of lower 33% and a group of higher 67%, stromal ratio into a group of lower 50% and a group of higher 50%, both in line with our previous study [25]. To evaluate the proposed uni/multivariate survival models, progressively decreasing number of samples from 100% to 70% were randomly drawn 1000 times from the cohort without replacement, survival analysis was performed and log-rank test p-value was recorded. This was used to test the stability of the uni/multivariate survival model by computing the percentage of times where it remained statistically significant. Kruskal-Wallis test was used for testing association between a continuous variable and a categorical variable. Fisher's exact test was used to test for association between categorical variables. Cell ratio was measured as the ratio of a cell type to all cells in a tumor. Variability of scores was measured using standard deviation. For testing the amount of tissue required for our analysis, we divided each section image into two halves horizontally or vertically, thereby obtaining equal sized tissue areas. New scores were computed using half of the tissue and compared with the original scores using absolute difference.

Data and software availability
R code (Sweave) and data files to reproduce the reported results are available as supplementary data on www.yuanlab.org.