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
Introduction
Solid tumors consist of malignant cells surrounded and infiltrated by a variety of host cells that are recruited and “corrupted” by the cancer, aiding its growth and spread (1, 2). A dynamic network of soluble factors, cytokines, chemokines, growth factors, and adhesion molecules drives the interactions between malignant and nonmalignant cells to create this tumor microenvironment (TME; refs. 3, 4). The TME network stimulates extracellular matrix (ECM) remodeling, expansion of abnormal vascular and lymphatic networks, and migration of cells into and out of the tumor mass (5, 6). Solid tumors are also typically stiffer than the surrounding tissue due to aberrant ECM deposition and organization that has a major influence on cell and tissue mechanics (7, 8).
Although the TME is of critical importance during initiation and spread of cancer, relatively little is known about its evolution or the relationship between the molecular mechanisms of disease progression and higher-order features such as the extent of disease, nonmalignant cell density, and tissue stiffness. Studies on molecular mechanisms of human cancer have mainly focused on large-scale genomic and transcriptomic analysis of primary tumors (9) and the immune cell landscape (1). Human cancer evolution is also now being studied in multiple metastatic sites (10, 11), but mainly in terms of the genomics of the malignant cells. Also, most of these analyses focus on one stage of a cancer.
Here, for the first time, we have used multilayered TME profiling of a metastatic site, omental metastases of high-grade serous ovarian cancer (HGSOC), to identify molecular changes that predict the higher-order TME features. Our study differs from other genomic and transcriptomic studies in two important ways: First, we have integrated data from six different TME parameters from each metastatic sample studied; and second, we have studied the evolution of metastases by including samples that vary in the extent of disease.
HGSOC is one of the most lethal of the peritoneal cancers: Less than 30% of patients currently survive more than five years after diagnosis, with little improvement in overall survival in the past 40 years (12). As poor prognosis is mainly due to early dissemination into the peritoneal cavity (12, 13) and HGSOC metastases have a complex TME (14), there is a need for an integrated understanding of its different components (12). We chose to study the omental TME because it is the most frequent site for HGSOC metastases and is routinely resected during debulking surgery.
Using samples ranging from minimal to extensive disease, we conducted cellular, biomechanical, and molecular analyses on each biopsy. Integration of the different components using multivariate analyses allowed us to define, for the first time, gene and protein profiles that predicted extent of disease and tissue stiffness while also revealing how the ECM is remodeled during metastasis development. Of particular interest was an ECM-associated molecular signature, which we termed the matrix index, that predicted both extent of disease and tissue stiffness in our sample set. This novel signature distinguished patients with shorter overall survival not only in ovarian cancer, but also in 12 other cancer types irrespective of patient age, stage, or response to primary treatment, suggesting a common matrix response to human primary and metastatic cancers.
Results
Study Design
We analyzed six different parameters of omental biopsies from 36 patients with HGSOC: the extent of disease, densities of nonmalignant cells, tissue mechanics, cytokines, matrisome protein, and RNA profiles (Fig. 1A). The samples ranged from uninvolved or minimally diseased omentum to biopsies with extensive disease (Supplementary Table S1; Fig. 1B and Supplementary Fig. S1A).
The extent of disease in each biopsy was measured by digital histopathology on hematoxylin and eosin (H&E)–stained sections and was calculated as the percentage of tissue area occupied by malignant cells and stroma. We termed this the “disease score.” Remodeling of the omentum was extensive when malignant cells were present and the malignant cells comprised a minor proportion of the tissue (Fig. 1B). In order to monitor for any significant changes in sample architecture during tissue processing for the different analyses, we took serial sections for H&E staining. We did not observe any major changes in disease score between the different areas analyzed.
The density of the major nonmalignant cell populations was measured by immunohistochemistry (IHC) and digital histopathology in the same specimens. The biomechanical properties of the tissues were measured using a mechanical indentation methodology (15) that gave us the tissue modulus of each sample.
Twenty-nine cytokines and chemokines were measured in protein lysates using an electrochemiluminescence assay (Supplementary Table S2). For proteomic analysis of the same biopsies, we focused on the ECM and associated molecules using a method that enriches whole tissue lysates for the matrisome protein compartment (16). Using this technique, we detected 145 proteins (Supplementary Table S3). The term “matrisome” is defined as the ensemble of all core ECM proteins (collagens, proteoglycans, glycoproteins) and associated molecules (the secretome, ECM regulators, and ECM-affiliated molecules) of tissue ECMs (17). After alignment and filtering, RNA sequencing (RNA-seq) identified 15,441 protein-coding genes (Supplementary Table S4). We then used univariate analyses and a multivariate regression method—partial least squares (PLS; ref. 18)—to model the relationships between these different components of the metastases (Fig. 1C).
The Relationship between Cell Density and Disease Score
Using a tissue microarray constructed from the biopsies, we quantified the adipocytes, fibroblasts, and leukocytes that were the major nonmalignant components in the specimens and related this to disease score. The area occupied by adipocytes, the major cell type of the normal omentum, decreased with disease score (Fig. 1B), and this was in part due to a reduction in the diameter of the adipocytes (Fig. 2A), which may reflect research showing that adipocytes can provide energy for ovarian cancer cell growth (13). Using α-smooth muscle actin (α-SMA) and α-fibroblast activation protein (α-FAP) as markers of cancer-associated fibroblasts (CAF; ref. 19), we assessed the area of the tissue occupied by α-SMA+ and α-FAP+ cells and found a strong positive correlation with disease score for both markers (Fig. 2B and C).
We then stained and counted six major leukocyte subtypes and plotted cell density against disease score. In all cases there was a significant positive correlation between leukocyte density and disease score (P < 0.001; Fig. 2D; Supplementary Fig. S1B). Densities of T cells with the surface markers CD3, CD4, CD8, and CD45RO strongly correlated with each other (P < 0.001, r > 0.6), but CD68+ macrophage density only weakly correlated with the other leukocytes (P < 0.05, r < 0.5; Fig. 2E; Supplementary Table S5).
Therefore, as metastases developed in the omentum, the fatty tissue was replaced by a combination of fibroblasts, lymphocytes, and macrophages. The cellular composition that we observed illustrates the changes from a normal omental tissue, primarily composed of adipocytes with minimal immune cell infiltrate and little fibroblastic reaction, to heavily diseased tissues with profound tumor-associated inflammation and a large increase in all types of leukocytes and fibroblasts. However, although there was a general increase in stromal cell density across all of the markers we studied, we observed larger variance with high disease score samples for all immune cell counts. This is to be expected as it is already well documented that the TME of advanced HGSOC biopsies ranges from sparse to dense leukocyte infiltration.
The immune cell densities significantly correlated with their corresponding immune gene expression signatures extracted from the associated RNA-seq data for each sample (Supplementary Table S6), and levels of the adipogenic transcription factor PPARγ mRNA declined with disease score (Supplementary Fig. S1C). Thus, the cell density scores were validated by the gene transcription data.
Leukocyte Density and Cytokine Networks in the TME
Next, we correlated leukocyte density against levels of 29 different cytokines and chemokines in protein lysates of the metastases. There were eight significant correlations (Fig. 2F; Supplementary Table S7), the strongest of which was, unexpectedly, an association between IL16, a chemoattractant and modulator of T-cell function, and the density of CD3+, CD45RO+, and CD8+ cells. These correlations became stronger with the 10 samples with the highest disease score (Supplementary Fig. S1D). IHC revealed IL16 protein in both malignant and stromal areas, with a higher density in the former (Fig. 2G). There was also a high positive correlation between global cell proliferation assessed by Ki67 and Lymphotoxin A (LTA), IL17A, IL15, and CXCL10 (Fig. 2F).
As cytokine networks are major determinants of leukocyte density and phenotype in the TME (3, 20), we asked if the cytokine proteins and genes we detected in the tissue lysates could inform us about the networks that regulate omental metastases. We constructed heat maps showing pairwise comparisons of cytokine protein and gene transcription levels (Fig. 2H; Supplementary Fig. S1E; Supplementary Tables S8 and S9). Overall, the protein gene correlation was 30%, in line with other studies (e.g., 21). The heat maps show five significant coexpressions at both gene and protein levels: IL6 with IL1A, IL1B, and IL8; CSF2 with IL8; and CCL4 with CCL3. IL6 was of particular interest as we previously identified this as a major mediator of cytokine networks in ovarian cancer (20, 22). Finally, we asked if levels of any of the cytokines and chemokines associated with disease score. There were weak but significant associations with disease score with IL12B, IL13, IL16, VEGF, CCL11, CCL26, and CXCL10.
These results suggest that malignant cell–derived cytokine and chemokine networks in the omental metastases regulate leukocyte density and overall proliferative index. Unexpectedly, we identified the CD4 ligand IL16 as a potential major mediator of the leukocyte infiltrate. It is interesting that increased tissue and serum levels of IL16 have been reported during tumor development in laying hen models of ovarian cancer and in patients with ovarian cancer (23).
These cytokine and cellularity data confirm and extend previous research on the ovarian cancer and other TMEs and we believe validate our approach. Extensive study of another TME component, the matrisome, has recently become possible through proteomics methods that focus on these proteins (17). There is currently little information on how the matrisome evolves with disease progression. Therefore, our next aim was to study the ECM-associated proteins and genes, collectively termed the matrisome, in the same biopsies.
How the Matrisome Changes with Disease Progression
Using our matrisome-focused proteomic technique (16) and the RNA-seq data, we quantitatively assessed matrisome proteins and genes. In terms of relative mass ratios, the major matrix proteins in samples with the lowest disease score were collagen 1, 6, and 3, the glycoprotein fibrillin, the ECM regulator alpha-2-macroprotein, and the basement membrane proteoglycans lumican and heparin sulphate proteoglycan-2. In contrast, biopsies with the highest disease score had an expansion of ECM glycoproteins fibrinogen and fibronectin, as well as increases in proteoglycans, secreted factors, and affiliated proteins (FDR < 0.1; Fig. 3A).
Extending the analysis to the entire sample set, we found that as disease score increased, levels of some matrisome proteins decreased and others increased. Comparing the relative mass ratio of all matrisome proteins with disease score, we found that 18 proteins decreased and 49 proteins increased with disease progression (Fig. 3B; Supplementary Table S10). After these univariate analyses, we used the multivariate PLS regression method to rank genes and proteins according to their influence on disease score, and a permutation-derived threshold was applied to determine those that were most strongly associated with disease score (24, 25).
Of these, 58 proteins ranked top in PLS modeling of disease score (r2 = 0.70), defining a matrisome protein signature of disease score (Fig. 3C; Supplementary Table S11).
Four hundred twelve of the 764 matrisome genes detected in our transcriptomics dataset also predicted disease score (Supplementary Table S12). The top 60 genes are shown in Fig. 3D with 27 ECM-associated molecules predicting disease score at both the gene and protein levels (Fig. 3E; Supplementary Fig. S2A). We used IHC to detect four of these proteins, FN1, COMP, CTSB, and COL11A1, in HGSOC omentum, detecting all four within stromal regions (Fig. 3F). As collagen organization strongly influences cell behavior and tissue mechanics (26, 27), we utilized two-photon microscopy to visualize collagen fibers using second harmonic generation (SHG) label-free illumination (Fig. 3G). In low disease score tissues, collagen fibers were thin and arranged mostly around the adipocytes. In high disease score tissues, there were denser arrays of long collagen bundles with an apparent microscale orientation preference. Collagen orientation correlated strongly with disease score (Fig. 3G). It should be noted that Fig. 3A shows relative matrisome protein abundance as mass ratios, whereas Fig. 3G depicts the alignment of collagen in the tissues. Although the amount of collagen does increase in diseased tissue, its relative abundance goes down as other matrisome proteins are induced as the tissue becomes more diseased. Representative images for COL1A1 IHC staining are shown in Supplementary Fig. S2B.
Other Biological Processes Associated with Disease Score
We then analyzed the RNA-seq data to find other biological processes associated with disease score (Supplementary Table S13). Significantly associated pathways included cell metabolism, adhesion, and communication, as well as ECM organization and immune response pathways (Supplementary Fig. S2C; Supplementary Table S14).
We have described here, for the first time, how the matrisome changes with disease progression. As some of the strongest correlations with disease score were found with these ECM-associated proteins and genes, and increased tissue stiffness has been linked with tumor progression (28, 29), we next investigated how the changes in matrisome genes and proteins related to the biomechanical properties of the biopsies.
Relationships between Tissue Modulus (Stiffness) and Disease Score
We used mechanical indentation (15) to determine tissue modulus (which describes material stiffness independent of sample dimension) and the stress–relaxation behavior of the samples. We measured disease score from histologic sections of the area of the specimen that was indented (Fig. 4A; Supplementary Fig. S3). Biopsies with a high disease score displayed a nonlinear loading response and greater stress relaxation, whereas there was a relatively linear loading response in low disease score tissue (Fig. 4B; Supplementary Fig. S3C; Supplementary Table S15). Tissue modulus values in high disease score biopsies were one to two orders of magnitude higher than in low disease biopsies. There were significant positive correlations between tissue modulus and malignant cell area, the stromal area and disease score of each biopsy (Fig. 4C; Supplementary Fig. S3D). Thus, there was a significant log relationship between tissue modulus and disease score in the evolving TME, suggesting a close association of tissue stiffness with disease progression.
The Matrisome, Tissue Stiffness, and Disease Score
Using the PLS method, we identified 64 matrisome proteins, mainly glycoproteins, that accurately predicted tissue modulus (r2 = 0.69; Fig. 4D; Supplementary Fig. S4A; Supplementary Table S16). We then used 764 matrisome genes detected by RNA-seq and identified 405 that predicted tissue modulus (Fig. 4E; Supplementary Fig. S4B; Supplementary Table S17), of which 38 also featured as proteins in Fig. 4D. Thus, as with disease score, the tissue modulus could be predicted by a subset of ECM-associated genes and proteins of the matrisome.
We also modeled tissue modulus against the entire transcriptome of the metastases (Supplementary Fig. S4C; Supplementary Table S18). Genes associated with cell metabolism, cell communication, wound healing, ECM organization, as well as development, correlated with tissue modulus (Supplementary Fig. S4D; Supplementary Table S19). Figure 4F shows the PLS prediction plot and the top 50 genes from this signature. As expected, there was a strong overlap with disease score–associated genes and proteins (74% and 75%, respectively), and these were significantly associated with tissue modulus (Supplementary Table S20).
Collectively, the experiments described above demonstrated the complexity and dynamic nature of matrisome evolution during development of HGSOC metastases and the close relationship between tissue stiffness and extent of disease.
A Subset of Matrisome Molecules Models Both Disease Score and Tissue Modulus and Has Prognostic Significance in HGSOC
We next asked how many matrisome genes and proteins significantly defined both disease score and tissue modulus in our sample set (Fig. 5A; Supplementary Table S21). Twenty-two molecules were highly significant across all of our analyses with a gene–protein concordance of 68% (Fig. 5A; Supplementary Fig. S5A). Thirteen of the 22 proteins had documented protein–protein interactions (Fig. 5B). Using the ChEA database (30), we found that the 22 genes shared a range of common transcription factors, including RUNX2, STAT3, SMAD4, WT1, JUN, and TP53. These reflect pathways associated with the WNT signaling pathway, inflammation, and osteogenesis, whereas TP53 is of course the most frequently mutated gene in HGSOC, a major genetic driver of the disease (Supplementary Fig. S5B and Supplementary Table S22).
Using these 22 most significant molecules, we measured the ratio between the mean expression levels of the positively regulated genes and the mean expression levels of the negatively regulated genes. We termed this the matrix index because these molecules are all components of the matrisome (17). As would be expected, the matrix index of each sample significantly correlated with disease score and tissue modulus in our set of samples (P < 0.0001; Fig. 5C). There were also significant positive and negative correlations between matrix index and immune cell signatures in the corresponding RNA-seq data (Fig. 5D; Supplementary Table S23), notably regulatory T cell (Treg) and TH2 cell signatures, cell subtypes associated with tumor promotion and immune suppression (31), and a modest statistically significant relationship between disease score and entropy as a measure of clonal abundance for T and B cells (Supplementary Table S23). These data suggest that matrix index molecules may influence expansion of specific infiltrating cell populations. In support of these findings, there were significant linear correlations between the matrix index and the data in Fig. 2 in terms of CD4+ (Pearson r = 0.523, P = 0.001) and FOXP3+ (Pearson r = 0.52, P = 0.001) cells, but there was no correlation between matrix index and CD8+ cells (Pearson r = 0.29, P = 0.094).
Relevance of Matrix Index to Other Stages of HGSOC and Prognosis
As the matrix index positively correlated with disease score, tissue modulus and immune-suppressive signatures in our sample set, we wondered if it would distinguish patients with ovarian cancer with a poorer prognosis in transcriptomic data from untreated primary tumors. We extracted expression values from two publicly available HGSOC gene expression datasets and calculated the matrix index for each sample. The high and low matrix index groups were determined using a method described previously (32). High matrix index significantly correlated with shorter overall survival of HGSOC patients in both the International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA) gene expression datasets, as well as in our original sample set (Fig. 5E; Supplementary Fig. S5B–S5E). To test that the clinical outcome association of the matrix index was not a random finding, we conducted 200,000 simulations and found that the association was significantly above that expected from random signatures.
In order to account for the higher relative abundance of tumor cells compared with stroma present in the TCGA ovarian cancer samples compared with our samples, we plotted the correlation between matrix index and percentage tumor cell or percentage stroma in each TCGA sample. In both cases, we observed no association with matrix index (Supplementary Fig. S5F). Taken together with the immune cell correlations, this further suggests the matrix index is not only a measure of the tissue remodeling that accompanies HGSOC, but is also a measure of a matrisome composition that better supports tumor progression.
Interrogating the TCGA ovarian cancer dataset, we next evaluated the power of the matrix index against nine other well-known prognostic gene expression signatures in ovarian and other cancers, including signatures for stromal and immune responses (33–41). In terms of hazard-ratio scores, matrix index was in the top three after the 26-gene breast cancer stromal signature reported by Finak and colleagues (41) and the 193-transcriptional signature from TCGA (ref. 9; Fig. 5F, left). However, using multivariate analysis, matrix index was the single significant predictor of ovarian cancer survival independently of age, stage, grade, and treatment outcome (Fig. 5F, right; Supplementary Table S24).
At the protein level, we used matrix index to examine the recently released TCGA/Clinical Proteomic Tumor Antigen Consortium ovarian cancer proteomics dataset. Although the study was not focused on detecting ECM proteins, which requires matrisome protein enrichment prior to analysis, as described above, we found that there were 12 proteins from the matrix index with a significant association with survival (10 with P < 0.05 and a further 2 with P < 0.1).
Matrix Index in Other Human Cancers
ECM remodeling is a common feature of many human cancers, and significant desmoplasia and ECM deposition is found in other solid tumors. Because we hypothesize that the matrix index is a measure of a tumor-promoting matrisome in HGSOC, we wondered if it might also be a feature associated with poor outcome in other cancer types. We calculated matrix index values in 30 other publicly available gene expression datasets from epithelial, mesenchymal, and hematologic malignancies analyzing data from 9,215 human cancer biopsies. High matrix index was an indicator of poor prognosis in epithelial and mesenchymal cancers but not in hematologic cancers, melanoma, and glioblastoma (Fig. 6A and Supplementary Fig. S6A). Using univariate analysis, high matrix index predicted shorter overall patient survival in 15 datasets representing 13 major cancer types (P < 0.05; Supplementary Fig. S6B; Supplementary Table S25). The range of matrix index values across all these cancer databases had a median value close to 1.0 (Supplementary Fig. S6C). We believe this provides further evidence that the pattern of ECM-associated gene expression determined by the matrix index may be a common feature of human cancers. Remarkably, multivariate analysis showed that the prognostic value of the matrix index was independent of age, stage, grade, and response to primary treatment in 15 of the datasets representing 13 major cancer types (P < 0.05; Fig. 6B).
Using IHC, we confirmed the presence of four of the upregulated matrix index proteins, FN1, COL11A1, CTSB, and COMP, in three tissue microarrays from triple-negative breast cancer (TNBC), pancreatic ductal adenocarcinoma (PDAC), and diffuse large B-cell lymphoma (DLBCL; Fig. 6C). These cancers reflected the range of hazard ratios for high matrix index in Fig. 6B. Digital microscopy analysis showed the highest staining level in TNBC (Fig. 6D), in keeping with the matrix index score for this cancer (Supplementary Fig. S6C). FN1, COMP, and CTSB were present in stroma and fibroblastic cells of all tumors. COL11A1 was located within the malignant cells in all biopsies. FN1 was also found in malignant PDAC cells and in immune cells in DLBCL. CTSB was located in macrophages in TNBC and PDAC, and tumor cells in DLBCL.
Data Resource
All data in this paper are provided in a mineable web-based resource (http://www.canbuild.org.uk).
Discussion
In this paper, we have profiled, for the first time, an evolving human metastatic microenvironment, using analysis that includes gene expression, matrix proteomics, cytokine/chemokine expression, ECM organization, and biomechanical properties, all performed on the same sample. This gives a unique and informative snapshot of the evolving metastatic state of one type of ovarian cancer. Integration of the most significant features of this microenvironment may have identified a matrix response that is conserved in other cancers of epithelial or mesenchymal origins.
Our study has also shown that conducting multilevel analysis with data integration on well-characterized cancer biopsies with a range of disease involvement, and multiple analyses per sample, can identify important features representative of the evolving TME. This approach is complementary to “omic” molecular cancer datasets that have larger numbers of samples. The data presented here provide a unique resource regarding molecular, cellular, and mechanical regulation in the TME and a template for bioengineers who are building complex TME models.
Molecular genetics has revealed great intra- and intertumor heterogeneity. It is now accepted that malignant cell clones undergo Darwinian evolution, resulting in a high level of molecular heterogeneity. In contrast, this study shows that the interactions between malignant cells and the host to remodel the tissue matrix may be more consistent. It is already known that high lymphocyte density is a common indicator of good prognosis at different stages of disease in many malignancies, including HGSOC (14). We suggest that another common feature of TMEs may be patterns of matrisome genes and proteins and that these also have prognostic significance.
The upregulated genes that were most significantly related to disease score in our analysis, COL11A1, COMP, VCAN, FN1, COL1A1, and CTSB, have all been associated with cancer progression, poor prognosis, and malignant cell invasion in ovarian and/or other cancers (42–49). For example, fibronectin promotes ovarian cancer invasion and metastasis through an α5β1-integrin/c-MET/FAK/Src-dependent signaling pathway (44), and COL11A1 and VCAN feature in a 10-gene poor prognostic signature of collagen-remodeling genes regulated by TGFβ signaling in ovarian cancer (45). More recently, COL11A1 expression showed a positive association with poor prognosis in several epithelial cancers (49). Importantly, the matrix index appears to correlate with certain immune cell signatures that are also known to influence prognosis. For example, the matrix index does not positively correlate with CD8 molecular signatures, which are associated with good prognosis in HGSOC, but significantly correlates with Treg and TH2 signatures.
We have not focused on proteases in this study. However, in our samples, normal tissue ECM stained highly positive and evenly for COL1A1, but COL1A1 staining in diseased sections, although still strong, appeared more uneven. In particular, there was reduced staining around malignant cell areas. This may be due to expression of proteases such as MMP13, which degrade collagen structures (50), that we identified in our PLS analysis. Other matrix remodeling proteases that we identified such as MMP7, MMP11, CTSB, and ST14 are also able to degrade collagens, albeit to a lesser extent, but may also be capable of degrading matrix proteoglycans and glycoproteins (51, 52). In addition to protease activity, the apparent reduction in COL1A1 by mass ratio (Fig. 3A) was in part due to the relative increase in matrisome complexity in diseased tissues. Also, we used a modified matrisome analysis method, which increased our protein detection coverage, at the cost of underestimating the absolute amount of large fibrous core matrisome proteins such as COL1A1.
As we have shown, there was a core group of matrisome molecules that best predicted tissue modulus. This appears to be through an expansion of matrisome glycoproteins and proteoglycans (Fig. 3A) and a reorganization of fibrillar collagens (Fig. 3G). The expansion of the glycoprotein and proteoglycan compartment increases the potential for posttranslational modification within the extracellular space, which could significantly alter the mechanical properties of the tissue, particularly through glycosylation and cross-linking (53). For example, glycosaminoglycans on proteoglycans contribute to hydration, which contributes to tissue stiffness (54). Additionally, specific ECM cross-linking molecules, including the pro-lysyl hydroxylases which cross-link matrix proteins through collagen-like peptides, were associated with increasing disease in our samples (55). Molecules such as COL11A1 and COMP are normally present only in stiffer tissues such as cartilage and bone, whereas VCAN plays a role in the morphogenesis of these stiffer tissues. In addition, collagens 11 and 6 play roles in collagen fibril organization and matrix integrity.
Tumor mechanics have a profound effect on fibroblasts and cancer cells and can promote tumor progression and metastasis (56). Stiffening of ECM creates a feed-forward self-reinforcing loop that contributes to the activation state of the fibroblast (57). Elevated mechanosignaling in PDAC cells as a result of elevated tissue stiffening promotes tumor progression and aggression (8). In breast tumors, the stiffest regions are located at the invasive margins and tumors harboring the stiffest regions are the most aggressive (58).
There is a significant desmoplastic response in most solid tumors, but given that there is large intra- and intertumor heterogeneity, it is important to understand why our index of matrisome gene expression defines patients with poor prognosis in multiple human cancers. The matrix index does not seem to be simply a measure of the amount of desmoplasia or stromal component that accompanies cancer growth. We believe that it is a measure of a type of matrisome composition that is more able to promote tumor growth. We found a strong association between the density of α-SMA— and α-FAP–positive cells, two markers commonly associated with activation of cancer-associated fibroblasts, and disease score, and there are several examples in the literature of poor prognostic fibroblast, desmoplastic, wound healing, and stromal signatures in individual cancer types (39, 59). As fibroblasts are the predominant matrix producing cells in many tissues this may, at least partially, explain the commonality of the matrix index in different cancers. It is also interesting that some of the matrix index molecules we found to be downregulated as disease increases (LAMB1, LAMC1, LAMA4, COL15A, and HSPG2) are associated with the basement membrane, which is vital for maintaining tissue homeostasis.
Malignant cell response to tumor-associated fibrosis, and the stromal cell phenotypes that contribute to ECM deposition, can vary within and among major cancer types. This was shown in great detail recently in a study of experimental and human pancreatic cancers where a distinct malignant cell genotype modulated the fibrotic phenotype of the tissue and pathology (8). This does not argue against our finding because we have found the matrix index is variable among different cases of each cancer.
As the predictive power of the matrix index was independent of age, stage, and response to primary treatment, we suggest that the pattern of change in the matrisome may reflect increased propensity of the malignant cells to establish metastases. It is intriguing that five of the six upregulated ECM genes in our matrix index (COMP, VCAN, FN1, COL1A1, and CTSB) are typical of premetastatic niches (60). Another explanation for the association with poor prognosis could be that this configuration of ECM molecules prevents infiltration or effector function of host antitumor immune cells. A stiffened matrix can compromise T-cell antigen presentation and proliferation as well as TH1-cell differentiation (61). In addition, the ECM acts as a reservoir for angiogenic factors and is important for migration of endothelial cells during neoangiogenesis and vascular remodeling seen in cancer (62).
Many of the matrix index molecules described above circulate systemically as fragments from protease remodeling, sometimes as neoepitopes (63). Therefore, further investigation of the matrix index may have potential as a cancer diagnostic/prognostic blood test.
If we have identified a common and especially detrimental signature of the tumor-associated matrisome, then agents that could target or reconfigure the cancer matrisome could have wide applicability in solid cancers and may enhance the action of immunotherapies, especially given the association of high matrix index with immunosuppressive T-cell signatures.
Some molecules of the matrix index may also prove good targets for drug delivery to the tumor site. A recent study demonstrated collagen I targeting of an anti-EGFR mAb showed increased therapeutic efficacy (64). Targeting matrix index molecules that are not as ubiquitous as collagen I may provide a significant advancement to such strategies.
Methods
Ovarian Cancer Patient Samples
Patient samples were kindly donated by women with HGSOC undergoing surgery at Barts Health NHS Trust between 2010 and 2014. Blood and tissue that was deemed by a pathologist to be surplus to diagnostic and therapeutic requirement were collected together with associated clinical data under the terms of the Barts Gynae Tissue Bank (HTA license number 12199. REC no: 10/H0304/14). Each patient gave written informed consent and the study was approved by a UK national review board. The studies were conducted in accordance with the Declaration of Helsinki and International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS).
RNA Isolation and Sequencing
Total RNA was extracted from 10 × 50 μm cryosections from frozen tissue sections and placed directly into the RLT Plus buffer (Qiagen) and rigorously vortexed. Samples were then processed using Qiagen RNeasy Plus Micro kit according to the manufacturer's instructions. RNA quality was analyzed on Agilent bioanalyzer 2100 using RNA PicoChips according to the manufacturer's instructions. RNA integrity numbers (RIN) were between 8.1 and 9.9. RNA-seq was performed by Oxford Gene Technology to ∼42× mean depth on the Illumina HiSeq2500 platform, strand-specific, generating 101 bp paired-end reads, as previously described (65). The detailed methods regarding RNA-seq and bioinformatic analysis are provided in the Supplementary Methods.
Quantitative Proteomics
The ECM component was enriched from frozen whole tissue sections (20 × 30 μm sections, approximately 40–50 mg of tissue) as previously described (17). The extracted proteins were reduced, alkylated and digested with trypsin. Peptides were separated by nanoflow ultra-high pressure liquid chromatography (UPLC) and analyzed by mass spectrometry using a LTQ-Orbitrap XL mass spectrometer (ThermoFisher Scientific). The detailed methods of ECM component enrichment, peptide preparation, mass spectrometry, and bioinformatics analysis are provided in the Supplementary Methods.
Cytokine and Chemokine Analysis
Cytokines and chemokines were assayed using Mesoscale Discovery Platform (MSD SI2400) according to the manufacturer's instructions. Cytokine panel 1 (Human) K15050D, Proinflammatory panel 1 (Human) K0080087, and Chemokine panel 1 (Human) K0080125 were used. Samples used were lysates from the ECM-enrichment protocol (described above). The amount of total protein used from each sample was between 1 and 3 μg.
Mechanical Characterization
Mechanical characterization was performed on whole tissues using a flat-punch indentation methodology on an Instron ElectroPuls E1000 (Instron, UK) equipped with a 10 N load cell (resolution = 0.1 mN) in order to measure the modulus of the tissue samples (15). The detailed method of quantification is provided in the Supplementary Methods.
Histochemical Analysis
Frozen tissues that were subsequently used for RNA, proteomics, and cytokine analysis were cryosectioned to 8- to 10-μm slices. Sections were fixed in 4% paraformaldehyde (PFA) and stained with hematoxylin and eosin (H&E) using standard methods. Tissues used in mechanical characterization were cut in half at the center of the tissue dye marked area and perpendicular to the direction of indentation while still frozen. Tissue was then fixed in 4% PFA for 24 hours and paraffin embedded and sectioned (8 μm) using standard procedures followed by H&E staining. All tissue sections were scanned using a 3DHISTECH Panoramic 250 digital slide scanner (3DHISTECH), and the resulting scans were analyzed using Definiens software (Definiens AG). Disease scores were determined first by manually defining regions of interest in the tissue that represented tumor, stroma, fat (adipocytes), or other (lymphatic structure) and then training the software to recognize these regions of interest. Disease score was expressed as a percentage of the whole tissue area that contained tumor and/or stroma (Fig. 1B). Detailed methods of immunohistochemical analysis for quantification of immune cells, α-SMA– and α-FAP–positive cells, adipocyte diameters, ECM proteins, and second-harmonic generation microscopy are provided in the Supplementary Methods.
Statistical and Bioinformatics Analysis
All statistical analyses and graphics were performed in the statistical programming language R (version 3.1.3). Detailed methodology for PLS regression models, Matrix index, and its clinical association across cancer types is provided in the Supplementary Methods.
Accession Numbers
RNA-seq data have been deposited in Gene Expression Omnibus (GEO) under the accession number GSE71340. Proteomic data are available via the PRIDE database accession number PXD004060.
Data Availability
All of the primary data are deposited at http://www.canbuild.org.uk.
Disclosure of Potential Conflicts of Interest
J. Gribben reports receiving commercial research grants from Janssen and Acerta and is a consultant/advisory board member for AbbVie, Celgene, Kite, Roche, Novartis, and Janssen. H.M. Kocher reports receiving commercial research support from Celgene Inc. and is a consultant/advisory board member for Baxalta. J.S. Serody reports receiving commercial research support from Merck Inc. J.D. Brenton reports receiving a commercial research grant from Aprea Therapeutics, has ownership interest (including patents) in Inivata Ltd., and is a consultant/advisory board member for AstraZeneca. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: O.M.T. Pearce, R.M. Delaine-Smith, E. Maniati, J. Gribben, M. Lockley, M.M. Knight, F.R. Balkwill
Development of methodology: O.M.T. Pearce, R.M. Delaine-Smith, E. Maniati, S. Nichols, J. Wang, V. Rajeeve, A. Montfort, P.R. Cutillas, M. Lockley, C. Bessant, M.M. Knight, F.R. Balkwill
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O.M.T. Pearce, R.M. Delaine-Smith, E. Maniati, S. Nichols, J. Wang, S. Böhm, V. Rajeeve, T. Dowe, J. Gribben, H.M. Kocher, P.R. Cutillas, M. Lockley
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): O.M.T. Pearce, R.M. Delaine-Smith, E. Maniati, J. Wang, V. Rajeeve, P. Chakravarty, J.S. Serody, B.G. Vincent, J. Connelly, J.D. Brenton, C. Chelala, P.R. Cutillas, C. Bessant, M.M. Knight, F.R. Balkwill
Writing, review, and/or revision of the manuscript: O.M.T. Pearce, R.M. Delaine-Smith, E. Maniati, J. Wang, R.R. Jones, J. Gribben, H.M. Kocher, J.S. Serody, B.G. Vincent, J.D. Brenton, M. Lockley, C. Bessant, M.M. Knight, F.R. Balkwill
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R.M. Delaine-Smith, E. Maniati, J. Wang, S. Böhm, D. Ullah, J.L. Jones
Study supervision: M. Lockley, M.M. Knight, F.R. Balkwill
Other (provision of tissue): T. Dowe
Acknowledgments
The authors wish to acknowledge the role of the Breast Cancer Now Tissue Bank in collecting and making available the samples used in the generation of this publication. We thank Barts Trust Oncology Surgeons for sample provision and Prof. Kairbaan Hodivala-Dilke for useful discussion. We also thank Andrew Clear, Dr. Joanne ChinAleong, Dr. Prabhu Arumugam, and Dr. Sally Dreger for technical help with the tissue microarrays, George Elia and the BCI Pathology Core, Christof Smith and Dr. Dante Bortone for help with bioinformatics analysis of the immune cell signatures, and Dr. Jackie McDermott for histopathologic analysis of the TMA samples. Finally, we express our gratitude to the patients for donating the samples without which this work would not have been possible. This project was funded by the European Research Council (ERC322566) and Cancer Research UK (A16354, A13034, and A19694).