Deconstruction of a Metastatic Tumor Microenvironment Reveals a Common Matrix Response in Human Cancers

Abstract

We have profiled, for the first time, an evolving human metastatic microenvironment by measuring gene expression, matrisome proteomics, cytokine and chemokine levels, cellularity, extracellular matrix organization, and biomechanical properties, all on the same sample. Using biopsies of high-grade serous ovarian cancer metastases that ranged from minimal to extensive disease, we show how nonmalignant cell densities and cytokine networks evolve with disease progression. Multivariate integration of the different components allowed us to define, for the first time, gene and protein profiles that predict extent of disease and tissue stiffness, while also revealing the complexity and dynamic nature of matrisome remodeling during development of metastases. Although we studied a single metastatic site from one human malignancy, a pattern of expression of 22 matrisome genes distinguished patients with a shorter overall survival in ovarian and 12 other primary solid cancers, suggesting that there may be a common matrix response to human cancer.Significance: Conducting multilevel analysis with data integration on biopsies with a range of disease involvement identifies important features of the evolving tumor microenvironment. The data suggest that despite the large spectrum of genomic alterations, some human malignancies may have a common and potentially targetable matrix response that influences the course of disease. Cancer Discov; 8(3); 304-19. ©2017 AACR.This article is highlighted in the In This Issue feature, p. 253.

Figure 1. Study design and sample description

a) Overview of the samples and the analyses conducted on the same tissue specimen. b) Bar plot shows results from digital analysis of architecture of haematoxylin and eosin (H&E)-stained samples based on percentage of malignant cell area (tumor), stroma, and adipocyte area, coloured blue, green and red respectively. The combined percentage area occupied by tumor and stroma was used to determine the ‘disease score’ of each sample. Each G number represents one sample. Upper microscope images show H&E staining of a biopsy and the same biopsy pseudo-coloured as malignant cell area (tumor) blue, stroma green and adipocyte area, red. Bottom images show four different H&E-stained samples representative of sample range with increasing disease score. c) Schematic of the PLS regression method used to define higher-order features of the tumor microenvironment from molecular components.

Figure 2. The cells of the TME change with disease score

a) Adipocyte diameter negatively correlated with increasing disease score. Top panel, microscope images representative of low (left) and high (right) disease score tissue sections (stained for α-SMA by IHC) showing adipocytes. Scale-bar corresponds to 100µm. Bottom left panel, scatter plot illustrating mean ± sd of digitally quantified adipocyte diameter (linear regression, N = 16, R² = 0.66, p = 0.0001). b) Correlation of α-SMA positive cells against disease score. Top panel, representative low (left) and high (right) disease score tissue sections stained for α-SMA by IHC. Scale-bar corresponds to 100µm. Bottom panel, quantification of α-SMA+ area % against disease score (linear regression, N = 30, R² = 0.83, p < 0.0001). c) Correlation of α-FAP positive cells against disease score. Top panel, representative low (left) and high (right) disease score tissue sections stained for α-FAP by IHC. Scale-bar corresponds to 100µm. Bottom panel, quantification of α-FAP+ area % against disease score (power regression, N = 32, R² = 0.77, p < 0.0001). d) Cleveland plots of immune cell counts against disease score (Spearman’s correlation, N = 34). e-f) Heatmap of pairwise Pearson’s correlation coefficients of e) immune cell counts (N = 34), f) MSD-quantified cytokine/chemokine correlations against immune cell counts (N = 32). h) IHC of IL16 in HGSOC omental biopsies. Scale-bars correspond to 100µm. g) Heatmap of pairwise Pearson’s correlation coefficients of MSD-quantified cytokine/chemokine (N = 32).

Figure 3. Identification of matrisome proteins and genes that define tissue architecture

a) Matrisome data displayed as relative mass ratios. Top panels show individual matrisome proteins identified in low and high disease score tissue; bottom panels show the relative proportions of each of the major classes of matrisome proteins in lowest (N = 6) versus highest disease score (N = 10). b) Line graphs illustrating normalized protein abundance and local polynomial regression fitted trend lines of proteins that either decrease (top panel), or increase (bottom panel) with disease score. c) PLS-identified matrisome proteins and d) matrisome genes that define disease score. e) Scatter plot of gene and protein correlation with disease score, highlighted molecules denote significant correlations (Pearson’s correlation, N = 33, p < 0.05). f) IHC staining for four matrisome proteins, FN1, COMP, CTSB, COL11A1 identified from PLS analysis as highly significantly related to disease score. Scale-bars correspond to 200µm. g) Collagen fiber alignment; top panel shows representative images of high and low disease score tissue sections visualised using second harmonic generation, and bottom panel, semi-quantification of fiber alignment from images plotted as number of fiber occurrences per angle bin (predominant fibre direction normalized to 0 degrees) with local polynomial regression fitted lines and disease color-coding.

Figure 4. Identification of molecular components that define tissue modulus

a) Orientation of flat-punch indentation showing representative low and high disease score samples stained with H&E, dashed line indicates tissue area analysed for determining disease score. b) Representative load-displacement curve from loading phase obtained from high and low disease score samples. c) Optimal tissue modulus correlated against combined % tumor plus stroma (disease score) (N = 32, p < 0.05). d-f) Cross-validation plot of measured versus predicted tissue modulus values (diagonal line represents measured = predicted) and heatmap of PLS-identified d) matrisome proteins, e) matrisome genes, and f) all coding gene components that describe tissue modulus. Heatmap columns correspond to individual samples ordered by increasing tissue modulus. (N = 29, 30 and 30, respectively). Rows ordered by decreasing model weight values.

Figure 5. A matrix signature that predicts survival in ovarian cancer

a) Venn diagram showing the overlap of PLS-identified molecules associated to tissue modulus and disease score (DS) at both gene and protein level. A total of 22 ECM-associated molecules overlapped across all analyses, red colour denotes positive association and blue colour negative association of each molecule at gene (G) and protein (P) level with disease score and tissue modulus. b) Network of known protein:protein interactions from IntAct and BioGRID within the 22 ECM-associated. Visualisation was carried out using Cytoscape v.3.3.0. c) Based on gene expression levels of these molecules we calculated a matrix index as the ratio of average level of expression of genes positively associated to those negatively associated with disease score and tissue modulus. Scatter plots show the correlation of matrix index with tissue modulus (linear regression, N = 30, R² = 0.74, p < 0.0001) and disease score (linear regression, N = 35, R² = 0.76, p < 0.0001). d) Association of matrix index with immune gene signature expression. Barplot illustrates Spearman p-values, FDR corrected using the Benjamini & Hochberg method. Red denotes positive correlations, blue denotes negative and gray denotes insignificant associations. The dotted line specifies the significance cutoff p = 0.05. e) Kaplan-Meier survival curves with overall survival of TCGA and ICGC dataset for HGSOC divided by high or low matrix index. The x-axis is in the unit of years. f) Comparison of hazard ratio scores (HR, with 95% CI) derived from Cox proportional hazards model for matrix index and the indicated gene expression signatures extracted from literature on the ovarian TCGA dataset. Left panel corresponds to univariate analysis, right panel corresponds to multivariate analysis taking into account age, tumor stage, grade and treatment (i.e., primary therapy outcome success). The asterisks represent the significance in the KM analysis between the high- and low-index groups (***p < 0.001, **p < 0.01, *p < 0.05 and ^■0.05 < p < 0.1).

Figure 6. Matrix index reveals a common stromal reaction across cancers

a) Kaplan-Meier survival curves with overall survival from the indicated datasets divided by high or low matrix index. The x-axis is in the unit of years. b) Multivariate hazard ratio (HR, with 95% CI) derived from a Cox proportional hazards regression model across cancer types / datasets using the matrix index. In each cancer, patients were split into high and low index groups, and their association with the overall survival (OS) was tested taking into account age, stage, grade (T-factor), and treatment factors. Asterisks represent the significance in the KM analysis between the high- and low-index groups (***p < 0.001, **p < 0.01, *p < 0.05 and ■0.05 < p < 0.1). HR > 1 means that high index is inversely correlated with OS, while HR < 1 means high index positively correlated OS. c) Example IHC images from TNBC, PDAC and DLBCL biopsies digitally quantified using Definiens™ software on cancer tissue array cores for matrix index proteins FN1, COL11A1, CTSB, and COMP. High intensity staining = red, medium = orange, low = yellow. d) Quantification of IHC staining on tissue arrays from TNBC, PDAC and DLBCL biopsies using Definiens™ software. Box plots illustrate the percentage area of high intensity staining for each marker. Scale bar = 500µm. COL11A1 and FN1, N = 30, 36, 54; CTSB, N = 28, 35, 52; COMP, N = 29, 35, 54; for TNBC, PDAC and DLBCL respectively.