# METASTASIS: FROM CELL ADHESION AND BEYOND

EDITED BY : Vasiliki Gkretsi and Triantafyllos Stylianopoulos PUBLISHED IN : Frontiers in Oncology

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-953-7 DOI 10.3389/978-2-88945-953-7

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# METASTASIS: FROM CELL ADHESION AND BEYOND

Topic Editors: Vasiliki Gkretsi, European University Cyprus, Cyprus Triantafyllos Stylianopoulos, University of Cyprus, Cyprus

This eBook consists of a collection of articles focused on fundamental processes of cancer cell metastasis, such as cell-Extracellular matrix adhesions, epithelial-to-mesenchymal transition (EMT) and lymph node metastasis as well as on upcoming research fields including the effects of biomechanical factors, the use of analytical and statistical tools and experimental techniques to further understand and characterize the invasive and metastatic potential of tumors.

Citation: Gkretsi, V., Stylianopoulos, T., eds. (2019). Metastasis: From Cell Adhesion and Beyond. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-953-7

# Table of Contents

*04 Editorial: Metastasis: From Cell Adhesion and Beyond* Vasiliki Gkretsi and Triantafyllos Stylianopoulos

#### SECTION 1

#### CELL ADHESION AND MATRIX STIFFNESS IN CANCER CELL METASTASIS

*07 Cell Adhesion and Matrix Stiffness: Coordinating Cancer Cell Invasion and Metastasis*

Vasiliki Gkretsi and Triantafyllos Stylianopoulos

*14 Defining the Role of Solid Stress and Matrix Stiffness in Cancer Cell Proliferation and Metastasis*

Maria Kalli and Triantafyllos Stylianopoulos

#### SECTION 2

#### CANCER CELL INVASION AND METASTASIS TO THE LYMPH NODES

*21 Molecular Mechanisms and Emerging Therapeutic Targets of Triple-Negative Breast Cancer Metastasis*

Christiana Neophytou, Panagiotis Boutsikos and Panagiotis Papageorgis


Kyo Suk Goo, Seaok Kim, Sunho Moon, Hayong Jung, Qifa Zhou, Robert H. Chow and K. Kirk Shung

#### SECTION 3

#### INHIBITION OF IMMUNE SUPPRESSION AS A MEANS TO DEAL WITH CANCER

*74 The Expression and Prognostic Impact of Immune Cytolytic Activity-Related Markers in Human Malignancies: A Comprehensive Meta-analysis*

Constantinos Roufas, Dimitrios Chasiotis, Anestis Makris, Christodoulos Efstathiades, Christos Dimopoulos and Apostolos Zaravinos

## Editorial: Metastasis: From Cell Adhesion and Beyond

Vasiliki Gkretsi 1,2 \* and Triantafyllos Stylianopoulos <sup>1</sup> \*

*<sup>1</sup> Cancer Biophysics Laboratory, Department of Mechanical and Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus, <sup>2</sup> Biomedical Sciences Program, Department of Life Sciences, School of Sciences, European University Cyprus, Nicosia, Cyprus*

Keywords: metastasis, cell adhesion, lymph nodes, epithelial to mesenchymal transition, triple negative breast cancer, cluster of circulating tumor cells, solid stress, stiffness

**Editorial on the Research Topic**

**Metastasis: From Cell Adhesion and Beyond**

#### METASTASIS

Metastasis is a complex multistep process during which cancer cells within a tumor dissociate from one another, migrate, and invade through surrounding tissues to finally enter the circulation or the lymphatic system, being thus transported to other sites of the body where they establish a new metastatic tumor. Many different approaches have been followed so far to study this process and find ways to prevent it. The collection of articles in this Frontiers Research Topic depicts exactly that.

## CELL ADHESION

Firstly, a pivotal role in the metastatic process is played by the cell-extracellular matrix (ECM) adhesion proteins as well as their interaction with actin cytoskeleton. Gkretsi and Stylianopoulos provide a concise review of the recent literature on important determinants of the cell's adhesome at cell–ECM adhesion sites that affect its invasive properties. Multiple protein-protein interactions define this adhesome linking the ECM directly or indirectly with the actin cytoskeleton (1, 2) and downstream effectors such as RhoGTPases (3) that collectively coordinate metastasis-related cellular processes. Furthermore, ECM accumulation within the tumor often leads to desmoplasia, an intense fibrotic response, causing tumor stiffening. Stiffening in turn, adds a biomechanical perspective to the whole concept of tumor growth and metastasis (4, 5). In that regard, Gkretsi and Stylianopoulos also emphasize the importance of keeping stiffness in mind when developing in vitro model systems.

### THE MECHANICAL COMPONENT

Adding to the biomechanical aspect of metastasis and tumor growth, Kalli and Stylianopoulos, define the concepts of stiffness and solid stress in tumors, as so far it is not clear whether matrix stiffness and solid stress are interrelated or if they have distinct roles in tumor progression. Pointing out that increased solid stress and stiffness are two distinct biomechanical abnormalities of the

Edited and reviewed by: *Paolo Pinton, University of Ferrara, Italy*

#### \*Correspondence:

*Vasiliki Gkretsi vasso.gkretsi@gmail.com Triantafyllos Stylianopoulos tstylian@ucy.ac.cy*

#### Specialty section:

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

> Received: *19 October 2018* Accepted: *11 March 2019* Published: *02 April 2019*

#### Citation:

*Gkretsi V and Stylianopoulos T (2019) Editorial: Metastasis: From Cell Adhesion and Beyond. Front. Oncol. 9:214. doi: 10.3389/fonc.2019.00214*

**4**

tumor microenvironment, they present a review of the different effects of these two parameters on the behavior of cancer and stromal cells. They also review and compare the in vitro experimental approaches that have been employed so far to analyze the effect of stiffness and solid stress providing a useful guide for similar studies.

#### TRIPLE-NEGATIVE BREAST CANCER METASTASIS AND METASTASIS TO THE LYMPH NODES (LN)

Along the same lines, Neophytou et al. focus on one of the most desmoplastic types of cancer, breast cancer, and triple negative breast cancer (TNBC), in particular, and provide a thorough analysis of the molecular mechanism involved during epithelial-to-mesenchymal transition (EMT), as well as the genes activated in this aggressive cancer type during the different stages of metastasis (metastasis promoting genes and metastasis suppressors). Moreover, they discuss recent advances on TNBC treatment, at the preclinical level, using agents that remodel tumor microenvironment and enhance the effects of chemotherapy delivery as well as advances emerging from novel molecular targets.

As with TNBC, in many cancer types, the first sites of metastasis for the original tumor are lymph nodes (LN). In fact, LN metastasis has been associated with worse prognosis although the mechanism is still vague. Jones et al. provide herein, an overview of the seeding, growth, and dissemination of LN metastases based on recent literature. Emphasis is given on how tumor cells and their secreted molecules decrease anti-tumor immunity and promote tumor growth in the LN.

#### CLUSTERS OF CIRCULATING TUMOR CELLS

Tripathi et al. focus on another aggressive type of breast cancer, which is also desmoplastic, inflammatory breast cancer (IBC). In this type of cancer, metastasis occurs not only through circulating tumor cells (CTCs) but rather via the generation of CTC clusters. CTC clusters may be rare and are thought to retain some epithelial characteristics, as they do not undergo a complete EMT, but account for more than 90% of metastases. Tripathi et al. based their work on a theory suggesting that the more hierarchically organized a physical system is, the more adaptable it can become. Thus, in the research article presented in this special issue, Tripathi et al. use the cophenetic correlation coefficient (CCC) to quantify the hierarchical organization in terms of gene expression of two different gene sets. They show that indeed high CCC, of both collective dissemination-associated genes and the IBC-associated genes, is associated with higher metastatic relapse rate in breast cancer patients.

#### A HIGH-THROUGHPUT, FUNCTIONAL TECHNIQUE FOR ASSESSING CANCER CELL INVASION POTENTIAL

Interestingly, in a more applicable point of view, the research article by Weitz et al. introduces a novel high-throughput, functional method for assessing cancer cell invasion potential. This method takes advantage of the biophysical changes occurring during metastasis that enable a cancer cell to invade the surrounding tissue. Using this technique, prostate, and bladder cancer cells are labeled with a fluorescent calcium dye and imaged during stimulation with low-intensity focused ultrasound; invasive cell lines exhibit calcium elevation which is not true for non-invasive cells (Weitz et al.). Thus, this method provides a means of assessing tumor invasion potential which could prove useful in cytology studies and ultimately improve clinical management (Weitz et al.).

### INTRATUMORAL IMMUNE CYTOLYTIC ACTIVITY

Last but not least, Roufas et al. provide us with a different view of dealing with metastasis focusing on immune checkpoint blockade therapy. Contrary to the approach taken by most anticancer immunotherapies, immune checkpoint blockade aims at blocking immune responses by inhibiting immune suppressor molecules, thus awakening the cytotoxic T lymphocytes from dormancy and enabling them to kill the cancer cells they infiltrate (6). Here, Roufas et al. conduct a comprehensive metaanalysis to evaluate the intratumoral immune cytolytic activity (CYT) in different cancer types, as judged by the expression of toxins granzyme A (GZMA) and perforin 1, and investigate differences between primary and metastatic tumors (data obtained from The Cancer Genome Atlas and Genotype-Tissue Expression project databases). They show that the cytolytic index among other associations with tumor-infiltrated immune cells promotes evasion from immunosurveillance in certain cancers Roufas et al..

## CONCLUDING REMARKS

In this research topic, we presented a collection of articles focused on fundamental processes of cancer cell metastasis, such as cell-ECM adhesions, EMT and LN metastasis as well as on upcoming research fields including the effects of biomechanical factors, the use of analytical and statistical tools and experimental techniques to further understand and characterize the invasive and metastatic potential of tumors.

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## REFERENCES


metastasis. Cold Spring Harb Symp Quant Biol. (2016) 81:189–205. doi: 10.1101/sqb.2016.81.030817

6. Wieder T, Eigentler T, Brenner E, Rocken M. Immune checkpoint blockade therapy. J Allergy Clin Immunol. (2018) 142:1403–14. doi: 10.1016/j.jaci.2018.02.042

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Gkretsi and Stylianopoulos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## Cell Adhesion and Matrix Stiffness: Coordinating Cancer Cell invasion and Metastasis

#### *Vasiliki Gkretsi1 \* and Triantafyllos Stylianopoulos <sup>2</sup> \**

*1Department of Life Sciences, Biomedical Sciences Program, School of Sciences, European University Cyprus, Nicosia, Cyprus, 2Cancer Biophysics Laboratory, Department of Mechanical and Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus*

Metastasis is a multistep process in which tumor extracellular matrix (ECM) and cancer cell cytoskeleton interactions are pivotal. ECM is connected, through integrins, to the cell's adhesome at cell–ECM adhesion sites and through them to the actin cytoskeleton and various downstream signaling pathways that enable the cell to respond to external stimuli in a coordinated manner. Cues from cell-adhesion proteins are fundamental for defining the invasive potential of cancer cells, and many of these proteins have been proposed as potent targets for inhibiting cancer cell invasion and thus, metastasis. In addition, ECM accumulation is quite frequent within the tumor microenvironment leading in many cases to an intense fibrotic response, known as desmoplasia, and tumor stiffening. Stiffening is not only required for the tumor to be able to displace the host tissue and grow in size but also contributes to cell–ECM interactions and can promote cancer cell invasion to surrounding tissues. Here, we review the role of cell adhesion and matrix stiffness in cancer cell invasion and metastasis.

#### Keywords: extracellular matrix, cell–extracellular matrix adhesion, actin cytoskeleton, cell invasion, metastasis, stiffness, solid stress, desmoplasia

Cancer cells undergo certain fundamental changes in terms of cell physiology to attain a malignant phenotype. They acquire self-sufficiency in growth signals, insensitivity to growth-inhibitory signals, limitless replicative potential, evasion of apoptosis, sustained angiogenesis, and tissue invasion capacity that enables them to metastasize to distant sites of the body (1, 2). In fact, the latter is the unique "hallmark of cancer" that differentiates benign and malignant tumors and truly defines cancer (3).

Metastasis is a complex process in which cancer cells spread from a primary site to other organs in the body. It consists of several steps and the involvement of the extracellular matrix (ECM), and the cytoskeleton is indisputable. During this process, malignant cells dissociate from the original tumor mass, reorganize their attachment to the ECM though alterations in cell–ECM adhesion dynamics, and start degrading surrounding ECM to eventually invade through adjacent

#### *Edited by:*

*Michelle Matter, University of Hawaii Cancer Center, United States*

#### *Reviewed by:*

*Santos Mañes, Consejo Superior de Investigaciones Científicas (CSIC), Spain Leonardo Freire-de-Lima, Universidade Federal do Rio de Janeiro, Brazil*

#### *\*Correspondence:*

*Vasiliki Gkretsi v.gkretsi@euc.ac.cy; Triantafyllos Stylianopoulos tstylian@ucy.ac.cy*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

*Received: 16 December 2017 Accepted: 20 April 2018 Published: 04 May 2018*

#### *Citation:*

*Gkretsi V and Stylianopoulos T (2018) Cell Adhesion and Matrix Stiffness: Coordinating Cancer Cell Invasion and Metastasis. Front. Oncol. 8:145. doi: 10.3389/fonc.2018.00145*

**7**

**Abbreviations:** 3D, three dimensional; CAF, cancer-associated fibroblasts; CTC, circulating tumor cell; ECM, extracellular matrix; EMT, epithelial to mesenchymal transition; ERK, extracellular signal-regulated kinase; FAK, focal adhesion kinase; FLNA, filamin A; FOXC2, forkhead box protein C2; ILK, integrin-linked kinase; LOX, lysyl oxidase; MDSC, myeloid-derived suppressor cell; MMP, matrix metalloproteinase; MSC, mesenchymal stem cell; PARVA, parvin alpha; PARVB, parvin beta; PARVG, parvin gamma; PI3K, phosphatidylinositide 3-kinase; PINCH, particularly interesting new cysteine-histidine rich protein; RSU-1, Ras suppressor-1; ROCK, Rho-associated protein kinase; SMA, smooth muscle actin; STAT, signal transduced and activator of transcription; TGF-β, transforming growth factor-β; TNF-α, tumor necrosis factor-α; VASP, vasodilatorstimulated phosphoprotein; ZEB1/2, zinc finger E-box-binding homeobox.

tissues and/or intravasate into blood vessels and travel through the circulation to distant sites of the body (4). The establishment of a metastatic tumor at the new site is not random but rather seems to follow a pattern known as "metastatic tropism." Cancer cells that have managed to survive in the circulation find a metastatic niche that, based on the "seed and soil" theory, is suitable for their growth (5–7). Hence, some cancer types metastasize according to circulation patterns or based on the anatomical proximity of neighboring organs or the host–organ microenvironment. For instance, prostate cancer shows a preference toward the bone, pancreatic cancer forms metastases to the lung and liver, and breast cancer metastasizes to the bone, liver, lung, and the brain (6–8). Notably, biophysical and biochemical cues from the tumor ECM affect each one of the "hallmarks of cancer" (9) and control cell–cell and cell–ECM adhesions, which in turn determine cancer cell invasion and metastasis (7). Thus, integrins, ECM-related adhesion proteins and cell–cell adhesion proteins play a vital role in regulating the various stages of metastasis and defining the aggressiveness of cancer cells (10).

#### CELL–CELL AND CELL–ECM ADHESION PROTEINS IN CANCER CELL METASTASIS

Cancer cells are able to invade the surrounding ECM in the form of single cells or as collective groups of cells moving together, depending on whether cell–cell adhesion proteins, such as E-cadherin, are completely or partially lost in the original tumor, respectively (11). Although integrin-independent migration has also been described (12), both modes of invasion are considered to be heavily dependent on integrin-mediated adhesion to the ECM, whereas collective invasion also requires dynamic cell–cell adhesions so that loosening of cell junctions becomes sufficient for invasion. Thus, E-cadherin expression or its localization in cell–cell junctions is often lost in advanced cancers and has been linked to higher incidence of metastasis (11).

However, the actual outcome in terms of invasion is ultimately dependent upon the balance between E-cadherin-mediated adhesions and integrin cell–ECM adhesions (11). Integrins connect the ECM with the interior of the cell transmitting extracellular signals through the assembly of multiple protein complexes that act as adaptor proteins and also bear strong attachments to actin cytoskeleton (10). There are more than 180 cell–ECM proteins forming networks of protein–protein interactions at cell–ECM adhesion sites, which altogether comprise what is known as cell's *adhesome* (13). Critical determinants in cell–ECM adhesions that also link the ECM directly or indirectly with the actin cytoskeleton include talin, paxillin, kindlins, vinculin, integrin-linked kinase (ILK), parvins [parvin alpha (PARVA), parvin beta, and parvin gamma], particularly interesting new cysteine–histidine rich protein (PINCH)-1, Ras suppressor-1 (RSU-1), vasodilator-stimulated phosphoprotein (VASP) and its interactor Migfilin (14), and α-actinin (15–17). Upon integrin, activation protein tyrosine kinases Src and focal adhesion kinase (FAK) are also activated promoting further cytoskeletal changes as well as activation of downstream signaling pathways vital for cell adhesion, proliferation, survival, migration, and invasion (**Figure 1**) (18). Small Rho GTPases, Rho, Rac, and Cdc42, as well as Rho-associated protein kinase (ROCK) are such downstream effectors known to coordinate cytoskeletal reorganization and cell migration. Interestingly, most of these components of the cell–ECM adhesions have been found to be significantly deregulated in most cancer types with their expression being associated with higher metastatic potential or lower survival rates (19–26). Moreover, increased levels of RhoA, RhoB, RhoC, Rac1, Cdc42, and ROCK, have been found in latestage tumors and metastases with prognostic relevance in breast cancer (27, 28). This suggests a strong involvement of cell–ECM adhesion molecules in cancer cell metastasis, although the exact molecular mechanisms involved can be different depending on cell type, tumor location, or grade. In fact, research has shown that cancer cells can have different modes of invasion, and thus a different molecular mechanism activated every time (29, 30). For instance, Rho signaling through ROCK promotes a rounded bleb-associated mode of motility, whereas elongated cell motility is associated with Rac-dependent F-actin-rich protrusions and does not require Rho or ROCK (30).

### DISRUPTION OF CELL–CELL ADHESION AND EPITHELIAL TO MESENCHYMAL TRANSITION (EMT)

All the above-described changes in cell–ECM and cell–cell adhesion components are important for the detachment of cancer cells from the original tumor mass and their invasion through adjacent tissues and contribute to the "epithelial to mesenchymal transition," also termed EMT. EMT refers to a transition of polarized epithelial cells toward cells exhibiting mesenchymal properties that enables them to metastasize. Thus, during EMT, epithelial cells reorganize their cytoskeleton, dissociate from one another, and begin expressing mesenchymal genes. These genes may vary significantly in different cells and tissues but there are certain transcription factors, such as TWIST1/2, SNAIL1/2, zinc finger E-box-binding homeobox, and forkhead box protein C2 that are indispensable for EMT in all cases (37–40). In fact, EMT-activating transcription factors have been proposed to have pleiotropic functions acting on all stages of cancer progression from initiation to metastasis (41). Also, several cytokines such as transforming growth factor-β (TGF-β), tumor necrosis factor-α, and interleukin-6, as well as ECM proteins such as collagen I, fibronectin, and hyaluronan are crucially involved in EMT in various tumors (37). Notably, several types of cancer cells have been found to acquire a more mesenchymal-like phenotype which also correlates with their resistance to cytotoxic drugs (38, 42), providing a link between EMT and cancer therapy. Moreover, expression of EMT markers has been also found in circulating tumor cells (CTCs) that are fundamental in the metastatic process. These markers facilitate detection of CTCs while also giving more insights into tumor diagnosis, treatment, and prognosis (43).

All in all, current studies have demonstrated the complexity of the EMT process which raises important and exciting questions for future investigation (41).

#### TUMOR MICROENVIRONMENT AND DESMOPLASIA

Apart from cancer cells and the ECM, tumors exhibit an additional aspect of complexity that accounts for the heterogeneity attributed to them and plays an important role in metastasis. They contain a number of allegedly normal cells that comprise the "tumor microenvironment" (1, 44). Hence, structural components of the tumor microenvironment are the tumor blood and lymphatic vessels, and the stromal cell constituents of the tumor that can be subdivided into three categories: (a) angiogenic vascular cells, which include endothelial cells and pericytes, (b) infiltrating immune cells, which include platelets, mast cells, neutrophils, inflammatory monocytes, myeloid-derived suppressor cells (45), macrophages (46), CD8<sup>+</sup> T-cells, NK T-cells, CD4<sup>+</sup> T-cells (47), and B cells, and (c) cancer-associated fibroblasts (CAF) cells, which include activated tissue fibroblasts, activated adipocytes, a-smooth muscle actin (α-SMA) positive myofibroblasts, and mesenchymal stem cells (48). As expected, the exact composition of a tumor's microenvironment varies depending on the tumor type and its location, which justifies the observed heterogeneity among tumors, rendering every tumor unique.

The ECM is a fundamental constituent of the tumor microenvironment that closely interacts with cancer cells for the transmission of signals in and out of the cell through integrins (10), while also providing the necessary growth factors for tumor growth (49). Moreover, upregulation of ECM remodeling molecules, such as TGF-β, are considered to be responsible for the development of *desmoplasia* in tumors (50). Desmoplasia is an intense fibrotic response characterized by the formation of dense ECM consisting of increased levels of total fibrillar collagen, fibronectin, proteoglycans, and tenascin C that accumulates within the tumor. It is associated with increased production and secretion of inflammatory and tumorigenic growth factors, and it is also characterized by an abnormally large population of stromal cells. Moreover, a large percentage of tissue fibroblasts are transformed to CAFs that contain high levels of α-SMA. Therefore, it is proposed that TGF-β activates fibroblasts to become CAFs, which in turn produce more ECM fibers leading to desmoplasia (50). Apart from that, molecules that remodel the ECM, such as matrix metalloproteinases and lysyl oxidase, are also critical for desmoplasia development (51). Collectively, desmoplastic tumors are considered to be more aggressive and are, in fact, associated with worse prognosis in several cancer types (52, 53).

#### ROLE OF TUMOR STIFFNESS IN CANCER CELL INVASION AND METASTASIS

desmoplasia is highly related to tumor *stiffening*, which is perhaps the only mechanical property of tumors that clinicians can really appreciate. Stiffness, which defines how rigid a material is or the extent to which a material resists deformation in response to an applied force (54), depends on the composition and organization of the structural components of a material and describes the extent to which it deforms in response to an applied force or the magnitude of the developed force when the material is subject to a specific strain. Therefore, the stiffer a material is, the more resistant to deformations and more prone to develop higher stresses (i.e., force per unit area) becomes. In tumors, in particular, which are known to grow at the expense of the host tissue, the stress exerted from the tumor on the host should balance the reciprocal stress applied from the host to the tumor. Therefore, the developed stresses within a tumor depend on the relative stiffness between the two tissues and from a biomechanical point of view, stiffening is required for a tumor to be able to displace the host tissue and grow in size (55, 56). Using mathematical modeling, we have previously estimated that tumors should be at least 1.5 times stiffer than their surrounding normal tissue, otherwise confinement by the host prevails to tumor expansion (57).

As mentioned earlier, tumor stiffness is mainly determined by the amount of ECM, particularly collagen and hyaluronan contained in the tumor. Given the fact that the interior of the tumor is subject to compression (58), its stiffness is mainly determined by hyaluronan, which owing to its fixed negative charges creates hydrated, gel-like regions within the tumor capable of resisting compressive stresses (59–62). At the tumor periphery, tumor growth can remodel the collagen fibers and change their orientation toward the tumor circumference. As a result, collagen fibers can be stretched and develop tensile stresses. Therefore, stiffness at the periphery should also depend on the amount of collagen (63, 64).

For the study of ECM stiffness, cancer cells usually grow in three dimensions (3D) within a collagen, hyaluronan, or similar gel that mimics the ECM and parameters that most often vary to modulate stiffness are either the gel's concentration or the degree of collagen crosslinking for gels that contain collagen. Cancer cell spheroids are also employed for the study of cell invasion through the matrix (65–70). Increasing ECM stiffness has been shown to induce malignant phenotype (71–73) characterized by Rhodependent cytoskeletal tension that leads to enhanced cell–ECM adhesions, disruption of cell–cell junctions and increased growth (69) (**Figure 1**) and is actually associated with activated FAK and extracellular signal-regulated kinase signaling (69). Finally, another proof that stiffness is crucially involved in cancer cell metastasis comes from preclinical studies showing that disruption of tumor ECM integrity halts metastasis (74).

Cells can sense ECM stiffening through integrins by cytoskeletal filaments that coordinate cell migration and induce changes within the cell. As that, a stiffer ECM can induce production of fibronectin, a glycoprotein of the ECM that binds from one side to extracellular collagen, fibrin, and heparan sulfate proteoglycans and from the other side to integrins. ECM stiffening can also enhance cell–ECM adhesions that connect the ECM to the cytoskeleton through local adhesion proteins, and increase cytoskeletal tension by Rho/ROCK signaling activation (69, 75). Therefore, integrin clustering can initiate the recruitment of focal adhesion signaling molecules such as FAK, ILK, PARVA, Src, paxillin, as well as Rac, Rho, and Ras that cause cell contractility and can promote tumor progression (**Figure 1**) (76, 77). In addition, stiffening of the ECM can enhance phosphatidylinositide 3-kinases activity, which regulates tumor invasion (78–80). Furthermore, the cell–ECM adhesion protein RSU-1 was found to be significantly upregulated in increased stiffness conditions in a 3D collagen-based *in vitro* culture system, while tumor spheroids made of cells lacking RSU-1 lost their invasive capacity through the 3D matrix in all stiffness conditions (65). Moreover, lack of the actin polymerization regulator VASP, also inhibited tumor spheroid invasion through matrix of increasing stiffness indicating that both actin cytoskeleton and cell–ECM adhesions play pivotal role in tumor spheroid invasion through 3D matrix (81), an *in vitro* property that mimics tumor invasion in a real tumor setting.

In addition, in pancreatic tumors with mutant SMAD4, matrix stiffening was associated with elevated ROCK activity that in turn stimulated increased production of ECM, assembly of focal adhesions and signal transducer, and activator of transcription-3 (STAT-3) signaling driving tumor progression (82). Matrix stiffening can also induce EMT, leading to the acquisition of a more aggressive phenotype that promotes cancer cell invasion owing to a loss of intercellular adhesions (83), and it is hypothesized to contribute to the transformation of cancer cells to stem cell-like cancer cells that can survive under the harsh hypoxic conditions of the tumor microenvironment, are more resistant to cytotoxic drugs, and can migrate and invade through surrounding tissues (84).

#### EFFECTS OF SOLID STRESS ON CANCER CELL BEHAVIOR

It should be noted, however, that even though ECM stiffness can be related to the magnitude of solid stress, the two quantities are distinct and thus, one should not be used to replace the other (85). Solid stress is defined as the force per unit area of the structural components of a tissue, which can cause either compaction (compression) or expansion (tension) of the material, whereas stiffness refers to the extent to which the tissue can resist deformations or external forces (54).

Different experimental procedures have been also developed to study the effects of solid stress and ECM stiffness on cancer cell behavior. For the study of solid stress, transmembrane pressure devices, cancer cell spheroids, or modifications of these are most often used. In transmembrane pressure devices, cells grow as single cells embedded in a matrix or as a monolayer on a transwell insert membrane, and a piston with adjustable weight is placed on the top to apply a predefined stress (86–89). This method has been used to study stress induced changes in gene expression, invasion, and migration. In the tumor spheroid model, cancer cell aggregates form spheroids that are embedded in a matrix that mimics tumor ECM, such as agarose, collagen, or matrigel. The matrix exerts an external stress to the cells and pertinent studies focus on the effect of solid stress on cancer cell proliferation and apoptosis (87, 90–93). This method is limited, however, in that the applied by the matrix stress cannot be directly quantified. When applied to compress cancer cells, solid stress can inhibit proliferation, induce apoptosis, and increase cancer cell invasive and metastatic potential (87, 88, 90). Compressive solid stress can also activate fibroblasts to become CAFs (similarly to TGF-β), which in turn can facilitate not only development of desmoplasia but also cancer cell invasion to the surrounding, normal tissues (89, 94).

#### CONCLUDING REMARKS AND PERSPECTIVES

Cell–ECM adhesion proteins, actin cytoskeleton, and ECM stiffness evidently play a major role in driving cancer cell invasion and metastasis being involved in virtually all steps of the metastatic process from cell dissociation from the original tumor, to invasion through surrounding ECM until the final step of cancer cell homing in the new metastasis site. For this to happen, ECM is connected through integrins to the cell's adhesome at cell–ECM adhesion sites where multiple protein–protein interactions take place connecting the ECM to the actin cytoskeleton so that response to external stimuli is well coordinated. In fact, signals from these cell-adhesion proteins appear to be crucial for defining the invasive potential of cancer cells, while evidence shows that they may also prove potent targets for inhibiting cancer cell invasion and thus, metastasis (65, 81). Moreover, as ECM stiffness is also a driving force in metastasis (72, 73), it also needs to be taken into account when studying cancer cell metastasis both *in vitro* and *in vivo* in an attempt to better recapitulate tumor

#### REFERENCES


microenvironment in a physiologically relevant manner. Thus, the development of appropriate and physiologically relevant *in vitro* systems is needed to define the molecular determinants in the process and open new avenues in the discovery of novel therapeutic candidates to block metastasis.

From another point of view, solid stress is a distinct parameter that affects cancer cell behavior and should be also considered in *in vitro* tumor models. Furthermore, solid stress is exerted not only on cancer cells but also on the endothelial cells that form the tumor micro-vessels. As a result, blood vessels can be compressed or totally collapsed, creating large avascular regions within the tumor thus causing hypo-perfusion and hypoxia (58, 95) which ultimately inhibit systemic administration of drugs to the tumors (55, 60) and can promote tumor progression in multiple ways (96). Notably, recent *in vivo* evidence has shown that modulating the tumor microenvironment through administration of drugs that alleviate intratumoral solid stress (such as anti-fibrotic agents) reduces mechanical stresses, decompresses tumor vessels, and improves tumor drug delivery (60–62, 97), once again suggesting that modulation of ECM is of fundamental significance for tumor biology and cancer therapeutics.

#### AUTHOR CONTRIBUTIONS

Both authors contributed equally to writing, editing, and approving the final version of this article.

#### ACKNOWLEDGMENTS

The authors are grateful to Dr. Andreas Stylianou and Ms. Maria Kalli for insightful suggestions on the design of the figure used in this manuscript.

#### FUNDING

This work has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC Grant Agreement No. 336839-ReEngineeringCancer.


adhesions and migration. *J Biol Chem* (2006) 281:12397–407. doi:10.1074/jbc. M512107200


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Gkretsi and Stylianopoulos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Defining the Role of Solid Stress and Matrix Stiffness in Cancer Cell Proliferation and Metastasis

*Maria Kalli and Triantafyllos Stylianopoulos\**

*Cancer Biophysics Laboratory, Department of Mechanical and Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus*

Solid tumors are characterized by an abnormal stroma that contributes to the development of biomechanical abnormalities in the tumor microenvironment. In particular, these abnormalities include an increase in matrix stiffness and an accumulation of solid stress in the tumor interior. So far, it is not clearly defined whether matrix stiffness and solid stress are strongly related to each other or they have distinct roles in tumor progression. Moreover, while the effects of stiffness on tumor progression are extensively studied compared to the contribution of solid stress, it is important to ascertain the biological outcomes of both abnormalities in tumorigenesis and metastasis. In this review, we discuss how each of these parameters is evolved during tumor growth and how these parameters are influenced by each other. We further review the effects of matrix stiffness and solid stress on the proliferative and metastatic potential of cancer and stromal cells and summarize the *in vitro* experimental setups that have been designed to study the individual contribution of these parameters.

#### *Edited by:*

*Barbara Zavan, Università degli Studi di Padova, Italy*

#### *Reviewed by:*

*Federica Barbieri, Università di Genova, Italy Elisabetta Rovida, Università degli Studi di Firenze, Italy Deepak Gurbani, University of Texas Southwestern Medical Center, United States*

> *\*Correspondence: Triantafyllos Stylianopoulos tstylian@ucy.ac.cy*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

*Received: 01 December 2017 Accepted: 20 February 2018 Published: 12 March 2018*

#### *Citation:*

*Kalli M and Stylianopoulos T (2018) Defining the Role of Solid Stress and Matrix Stiffness in Cancer Cell Proliferation and Metastasis. Front. Oncol. 8:55. doi: 10.3389/fonc.2018.00055*

Keywords: extracellular matrix, fibroblasts, externally applied stress, growth-induced stress, *in vitro* models

### UNRAVELING THE TUMOR MICROENVIRONMENT

Tumor stroma and biomechanical abnormalities developed during tumor growth comprise dominant regulators of cancer progression (1–3). The tumor stroma is composed of an extracellular matrix (ECM), which consists of immune cells, fibroblasts, capillaries, and fibrillar proteins, such as collagen I, elastin, and fibronectin, as well as hyaluronan and other sulfated glycosaminoglycans (4). Fibroblasts are key regulators of ECM composition and organization and physiologically remain in the quiescent state with negligible metabolic and transcriptomic activities (5, 6). In response to tissue damage, fibroblasts become activated and are characterized by the expression of alpha-smooth muscle actin (α-SMA). In this activated state, fibroblasts overproduce ECM proteins, mainly collagen I and fibronectin, secrete cytokines and growth factors, and exert contractile forces modifying tissue architecture (5, 6).

In tumors, fibroblasts tend to acquire a constantly activated phenotype as a response to several growth factors secreted from the highly proliferative cancer cells, including transforming growth factor-β (TGFβ), epidermal growth factors (EGFs), and bone morphogenetic proteins (BMPs) (5, 6). Activated fibroblasts, which are commonly known as cancer-associated fibroblasts (CAFs), start a chronic wound healing-like response toward cancer cells, leading to an excessive accumulation of fibrillar ECM proteins, a condition known as desmoplasia (5). Under this desmoplastic reaction, CAFs continuously produce and remodel the tumor ECM increasing tumor stiffness (1, 5). Desmoplasia and ECM stiffening characterize many tumor types, especially breast and pancreatic cancers, and it usually promotes tumor progression (1, 7, 8).

As the density of cancer cells, stromal cells, and ECM constituents increase within the restricted environment of the host tissue, it leads to the development of mechanical stress (i.e., force per unit area) within the tumor (1, 3, 9–11). This stress, derived from the structural components of a tumor, is known as solid stress and can be divided into two parts. A part of it, known as growthinduced (or residual) stress, is accumulated during tumor growth by microscopic interactions among structural components of the tumor microenvironment, and it remains within the tumor even if the tumor is removed (3). These interactions might include collagen stretching by cancer cells and CAFs, and hyaluronan and cancer cell swelling to resist compression (12–15). Moreover, as tumors grow and exert forces on the adjacent host tissue, a reciprocal compressive stress is applied from the host tissue to the tumor, to resist tumor expansion (1). This stress is known as externally applied stress, and it diminishes after tumor excision (1). The total solid stress in a tumor interior is compressive (i.e., tends to reduce the size of an object), while near the interface between the tumor and normal tissue, the stress is tensile (i.e., tends to increase the size of an object) (16, 17).

#### THE DEFINITION OF ECM STIFFNESS AND SOLID STRESS

It is not clearly defined in the pertinent literature whether matrix stiffness and solid stress refer to the same term or they are two distinct biomechanical abnormalities of a tumor that are related to each other. By definition, stiffness is a material property, which describes the extent to which a material resists deformation in response to an applied force, while solid stress is a force per unit area, which can cause either compaction (compression) or expansion (tension) of a material. In solid tumors, the stiffness is mainly determined by ECM composition and organization, while solid stress arises by the sum of the physical forces exerted during tumor growth. These forces can be generated in the subcellular level by cytoskeletal filaments that control cellular processes such as filopodia formation and extension. At the cellular level, forces are exerted by cell contractions (such as in CAFs) and cell–ECM interactions during migration of cancer and stromal cells, while at the tissue level, forces are exerted between the tumor and the host tissue (18–21).

The relationship between tumor stiffness and solid stress can be described using the analogy of a spring of specific elastic modulus (E) that obeys Hooke's law (**Figure 1**). According to the equation of Hooke's law for linear elastic materials, σ ε = ⋅ *E* , when a tumor of elastic modulus E grows and pushes the surrounding host tissue of elastic modulus E′, it causes a deformation ε1 and a subsequent stress σ1. As a consequence, the host tissue returns an equal and opposite stress σ1′, the so-called externally applied solid stress. At the same time, growth-induced solid stress is accumulated in the tumor interior owing to interactions among tumor components (**Figure 1A**). Thus, the total solid stress accumulated intratumorally is the sum of the externally applied and the growth-induced solid stress. In the case that the stiffness of the tumor E2 is greater than E1, then the tumor can displace the host tissue with a greater deformation and the externally applied solid stress σ2 can be greater than σ1 (**Figure 1B**). Therefore, in

this case, a solid tumor creates a stiffer matrix to push against the normal tissue and grow in size. Indeed, it has been demonstrated using mathematical modeling that the stiffness of a solid tumor should be at least 1.5 times greater than that of the host tissue, in order for the tumor to displace the tissue and grow (14).

As for the growth-induced solid stress, it increases during tumor growth, while the matrix stiffness might stop changing (16, 17). In this case, the further increase in total solid stress accumulated intratumorally can become less depended on matrix stiffness (**Figure 1C**). This hypothesis was confirmed by an elegant study by Nia et al. (16), showing that the total solid stress transmitted into the cells can depend only in part on tumor stiffness, and thus, the two terms should not be used without a distinction. Specifically, Nia et al. found that primary pancreatic tumors exhibited larger stresses compared to those in metastatic sites, while the opposite effect was observed for colon tumors (16). Interestingly, tumor stiffness was similar in the primary and metastatic tumor for both the pancreatic and colon cancer models, showing that tumor stiffness and solid stress are not necessarily coupled to each other. In addition, they found that solid stress increased in breast tumors of larger size despite the fact that stiffness did not change with tumor size. In line with our analysis, these observations can be explained by the fact that growth-induced solid stress generated owing to microscopic interactions among structural components in the tumor interior contributes to the accumulation of an additional to the externally applied solid stress. Therefore, the effects of matrix stiffness and solid stress on tumorigenesis and metastasis should be studied separately (22). Following, we provide a summary of these effects on cancer and stromal cell behavior, elaborating on the less studied contribution of solid stress and the pertinent experimental setups.

### EFFECTS OF MATRIX STIFFNESS ON CANCER AND STROMAL CELLS

The effect of ECM stiffness on cancer and stromal cells has been studied using *in vitro* two-dimensional substrates (2D) and threedimensional tumor analogs (3D). In 2D models, cells are seeded on coating substrates such as collagen or fibronectin (23–26), while the 3D models include single cells or tumor spheroids embedded in gels composed of collagen or matrigel (27–34) (**Figure 2A**). In both cases, stiffness is increased by changing the protein density or the degree of crosslinking of the matrix to study the effects of ECM-originating mechanical cues on cancer and stromal cells.

Matrix stiffness can activate intracellular signaling pathways to regulate cellular behavior. Cancer cells recognize the increase in ECM stiffness and respond by generating increased traction forces on their surroundings through actomyosin and cytoskeleton contractility (9, 35, 36). Moreover, the changes in matrix rigidity are sensed and transmitted intracellularly through mechanosensors such as p130 CRK-associated proteins, growth factor receptors, or integrin-ECM adhesion plaques (9, 10, 23, 35, 37–40). These mechanosensors can recruit focal adhesion molecules such as FAK, SRC, paxillin, RAC, RHO/RAS GTPases, and Rho-associated kinase to trigger signaling cascades and cytoskeleton organization (9, 10, 35, 36, 39, 41–44). These cascades finally regulate gene expression and induce quantifiable changes in cell shape, survival, migration, and invasion (9, 35, 39, 42). For example, it has been shown that tissue stiffness activates the nuclear translocation of the transcription factor TWIST1 in breast cancer cells, which inhibits the expression of E-cadherin and promotes cell invasion (35, 45). Furthermore, in a 3D model consisting of breast tumor spheroids growing in collagen matrix, the Ras suppressor-1 (RSU-1), a cell-ECM adhesion protein, was shown to be upregulated as a response to increasing stiffness. Interestingly, tumor spheroids knockdown for RSU-1 or actin polymerization regulator (VASP) lost their invasiveness through the 3D matrix (46, 47). Matrix stiffening is also shown to induce fibroblast activation and migration, which leads to a fibrotic response setting a positive feedback to matrix stiffness (13, 15, 35, 48, 49). However, in these studies, it cannot be distinguished explicitly whether the observed effects are emerged by increased cell-ECM adhesion sites owing to increased ECM density or by stiffness-induced solid stress generation.

#### EFFECTS OF SOLID STRESS ON CANCER AND STROMAL CELLS

While the role of ECM stiffness in cancer and stromal cells is actively studied, data regarding the effect of solid stress on tumor progression are elusive. There are several experimental setups mimicking the solid stress developed in the tumor microenvironment. These setups include models consisting of tumor spheroids growing in a confined environment that causes the development of solid stress (50–57) and models employing a transmembrane pressure device that applies a mechanical compression on a cell monolayer or on single cells embedded in a matrix (51, 58–60) (**Figure 2B**).

Regarding the first method, cancer cells are grown as spheroids in a polymer gel (e.g., agarose), which leads to the development of solid stress that resists to spheroid expansion (**Figure 2B**, i). Helmlinger et al. (55) using spheroids of colon adenocarcinoma cells estimated that the accumulated solid stress was in the range of 45–120 mmHg (6–16 kPa), depending on the concentration of the agarose gel and the size of the spheroid. In an analogous study, Cheng et al. (51) estimated the solid stress to be ~28 mmHg (~3.73 kPa) when mammary carcinoma cell spheroids were growing in a 0.5% agarose matrix. Recent *in vivo* measurements of breast, colon, pancreatic, and brain tumors estimated that the growth-induced stress is in the range of 1.56–142.4 mmHg (0.21–20 kPa) (3, 11, 16, 54). Differences in the magnitude of solid stress among *in vitro* studies and between *in vitro* and *in vivo* methods should depend on the tumor model and the experimental procedure used in each study. However, the conclusion that increasing compressive stress inhibits tumor growth is common (51, 52, 55, 57), while this effect was reversed when loads were removed (51, 55). It was also observed that solid stress can regulate tumor morphology since mechanical loads can induce apoptotic cell death in regions with high compressive stress and allow proliferation in low-stress regions of the tumor spheroid, suggesting that anisotropic loads result in anisotropic tumor growth (51).

Figure 2 | Experimental methods employed to analyze the effects of stiffness and solid stress on cancer and stromal cells *in vitro*. (A) Experimental setups studying the effect of ECM stiffness on cancer and stromal cells. There are two-dimensional models (2D) consisting of (i) a cell monolayer seeded on coating substrates (e.g., collagen type I or fibronectin) and three-dimensional models (3D) consisting of (ii) tumor spheroids or (iii) single cells embedded in a matrix (e.g., collagen type I, matrigel). Both models were aimed to investigate the effect of changes in extracellular rigidity on the transduction of mechanical signals into the cells as well as on the migration, invasion, proliferation and gene expression of cancer and stromal cells (B) Experimental setups studying the effect of solid stress on cancer and stromal cells. Setups include tumor spheroids that grow within (i) a polymer matrix, (ii) within elastic capsules, or (iii) in a confined polymer device. (iv,v) The setups are composed of cells seeded on the inner chamber of a transwell insert on the top of which an agarose cushion is placed or are embedded in a polymer matrix. A piston with adjustable weight applies a predefined and measurable compressive solid stress on the cells. These models provided useful information about the direct effect of solid stress on tumor growth and morphology as well as on cancer cell proliferation, migration, and gene expression. (C) A summary of *in vitro* and *in vivo* studies for the effect of solid stress in tumor progression.

More recent studies developed novel techniques to mimic solid stress during tumor growth in the absence of a matrix. Alessandri et al. (50) employed a microfluidic method based on the encapsulation and growth of cells inside permeable, elastic, and hollow microspheres (**Figure 2B**, ii). This approach offered the ability to produce size-controlled multicellular spheroids growing in confined conditions. They found that the confined spheroids exhibited a necrotic core compared with the unconfined spheroids. In contrast, peripheral cells were more proliferative and migratory, suggesting that mechanical cues from the surrounding microenvironment may trigger cell invasion from a growing tumor (50). Desmaison et al. (53) designed polymer polydimethylsiloxane microdevices to restrict the growth of spheroids and subsequently to induce the development of mechanical stress (**Figure 2B**, iii). They showed that the mitosis of mechanically confined spheroids was suppressed compared to spheroids grown in suspension (53). Furthermore, it was demonstrated that a population of cells within the confined tumor spheroids was arrested at mitosis, which was due to the inhibition of bipolar spindle assembly (53). Later, Fernández-Sánchez et al. (54) developed a method that allows the delivery of a defined mechanical pressure *in vivo* by subcutaneously inserting a magnet close to the mouse colon. The implanted magnet generates a magnetic force on ultramagnetic liposomes stabilized in the mesenchymal cells of the connective tissue surrounding colonic crypts after intravenous injection (54). The induced pressure was similar in magnitude to the endogenous stress (54), in the order of 9.0 mmHg (1.2 kPa), without affecting tissue stiffness, as monitored by ultrasound strain imaging and shear wave elastography (54). The magnetic pressure stimulated Ret activation and the subsequent β-catenin phosphorylation, impairing its interaction with E-cadherin in adherens junctions (54). These data suggested that tumor progression could be driven by signaling pathways that are directly activated by mechanical pressure.

To study the effect of a predefined solid stress on cancer cells, the transmembrane pressure device has been introduced (**Figure 2B**, iv). Setups employed consist of a transwell insert that fits in a well of a 6-well culture plate. The insert is separated in the lower chamber containing culture medium and the upper chamber containing the cell monolayer. A piston of a preferable weight is applied on the cell monolayer, while water, nutrients, and oxygen from the culture media are diffused through the pores of the transmembrane. This device provides a tool to mimic solid stress in a predefined manner according to the stress magnitudes measured in native tumor tissues.

Cheng et al. (51) used this device to study the effect of solid stress on murine mammary carcinoma cells. In this study, they applied a stress ranging from 0 to 60 mmHg (0–8 kPa), and they observed increased apoptosis with increased stress levels. In a following study, they used the same experimental setup to study the migration of cancer cells using a scratch wound assay (60). They applied a stress of 5.8 mmHg (0.77 kPa), and concluded that in these levels of compression, cancer cells stopped proliferating and started to create a leader cell formation, which allowed them to move toward the scratch having an invasive phenotype. Mitsui et al. (59) used a similar device for bone osteosarcoma cells to identify the effect of compressive stress on the expression of matrix metalloproteinases and plasminogen activators. They observed enhanced protein and mRNA levels of these molecules under low mechanical compression of bone cells (0–2.20 mmHg/0–0.29 kPa) (59). Recently, Chen et al. (61) observed increased migration and mesenchymal-like phenotype of renal carcinoma cells that were compressed by 0–5.0 mmHg (0–0.66 kPa), while Kalli et al. (62) found that normal fibroblasts become activated as a response to solid stress to promote pancreatic cancer cell migration.

Another device that was developed to study the effect of solid stress in a more realistic way involved the use of single cancer cells growing in an agarose matrix (**Figure 2B**, v). This device was composed of two custom-made parts, the well pressor and the optic pressor (58). Both devices consisted of a chamber containing a 3D gel with single cells embedded, a screw and a nut for pressure application, and their housing support. Specifically, the well pressor applied a strain that compressed the cell-contained agarose gels to 50% of their original volume. This stress was estimated to be ~0.37 mmHg (~0.05 kPa), much smaller than loads measured by other studies (3, 51, 55, 58). However, this stress was sufficient to cause differential expression profiles of metastasis-associated genes in glioblastoma and breast cancer cells. In addition, the optic pressor provided quantifiable changes in cell circularity and orientation with respect to the direction of the applied force (58).

Collectively, these *in vitro* studies suggest that mechanical forces can regulate tumor morphology, tumor growth, and metastatic potential of cancer cells in the absence of matrix stiffness. However, as indicated in **Figure 2C**, there is a discrepancy among the levels of solid stress applied or estimated in the pertinent studies due to the variability of the experimental procedures and the cancer models used. Therefore, it should be given special attention when performing experiments to study the effect of solid stress on tumor progression, taking into account the estimations derived from *in vivo* studies.

#### CONCLUSION AND FUTURE PERSPECTIVES

In light of recent studies showing that increased matrix stiffness and elevated solid stress are two distinct tumor abnormalities, and given the fact that most pertinent studies are focused on the effects of stiffness, it becomes clear that scientific efforts need to focus on the implications of solid stress in tumor progression and metastasis (16, 22).

Regarding the implications of tumor stiffness in tumor progression, most pertinent *in vitro* models include only cancer cells and ECM matrix. However, tumor stiffness might also depend on the presence of stromal CAFs that continuously interact with the fibrillar proteins. CAFs-ECM interactions remodel the ECM organization and fibers orientation for cancer cells to migrate and invade into the matrix (1, 63, 64). Regarding the effects of solid stress on tumor progression, further studies are required to shed light upon the mechanisms by which solid stress is transmitted and guides cellular behavior of both cancer cells and CAFs. Moreover, CAFs exert contractile forces that contribute to the accumulation of solid stress in the tumor interior. Therefore, it is necessary to include both cell types when solid stress and ECM stiffness are being studied.

It has been also shown that CAFs dynamically interact with cancer cells to promote tumor progression (62, 64). In fact, CAFs mediate the invasiveness of colon, pancreatic, and breast cancer cells when co-injected into mice (64–68), while breast and prostate tumors containing CAFs grew faster than tumors injected with normal fibroblasts (69, 70). Nevertheless, there is no pertinent study taking into account the effect of ECM stiffness and solid stress on the interaction of cancer cells and CAFs and *vice versa* the implication of tumor-stromal interactions in ECM stiffening and solid stress accumulation.

Concerning the complexity of the tumor microenvironment, new experimental setups consisting of cancer cells, CAFs, and changes in matrix stiffness and solid stress, in combination or separately, should be introduced to broaden our knowledge about the role of each component in the evolution and malignancy of cancer.

#### AUTHOR CONTRIBUTIONS

MK and TS planned and wrote the manuscript of this mini review.

### FUNDING

This research has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC Grant Agreement No. 336839-ReEngineeringCancer.

### REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any potential conflict of interest.

*Copyright © 2018 Kalli and Stylianopoulos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Molecular Mechanisms and emerging Therapeutic Targets of Triple-negative Breast Cancer Metastasis

*Christiana Neophytou1 , Panagiotis Boutsikos2 and Panagiotis Papageorgis2 \**

*1Department of Biological Sciences, School of Pure and Applied Sciences, University of Cyprus, Nicosia, Cyprus, 2Department of Life Sciences, European University Cyprus, Nicosia, Cyprus*

Breast cancer represents a highly heterogeneous disease comprised by several subtypes with distinct histological features, underlying molecular etiology and clinical behaviors. It is widely accepted that triple-negative breast cancer (TNBC) is one of the most aggressive subtypes, often associated with poor patient outcome due to the development of metastases in secondary organs, such as the lungs, brain, and bone. The molecular complexity of the metastatic process in combination with the lack of effective targeted therapies for TNBC metastasis have fostered significant research efforts during the past few years to identify molecular "drivers" of this lethal cascade. In this review, the most current and important findings on TNBC metastasis, as well as its closely associated basal-like subtype, including metastasis-promoting or suppressor genes and aberrantly regulated signaling pathways at specific stages of the metastatic cascade are being discussed. Finally, the most promising therapeutic approaches and novel strategies emerging from these molecular targets that could potentially be clinically applied in the near future are being highlighted.

#### *Edited by:*

*Scott Bultman, University of North Carolina at Chapel Hill, United States*

#### *Reviewed by:*

*Saraswati Sukumar, Johns Hopkins University, United States Daniele Vergara, University of Salento, Italy*

#### *\*Correspondence:*

*Panagiotis Papageorgis p.papageorgis@euc.ac.cy*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

*Received: 28 November 2017 Accepted: 30 January 2018 Published: 22 February 2018*

#### *Citation:*

*Neophytou C, Boutsikos P and Papageorgis P (2018) Molecular Mechanisms and Emerging Therapeutic Targets of Triple-Negative Breast Cancer Metastasis. Front. Oncol. 8:31. doi: 10.3389/fonc.2018.00031*

Keywords: triple-negative breast cancer, metastasis, targeted therapy, tumor microenvironment, dormancy

#### INTRODUCTION: TUMOR HETEROGENEITY AND CURRENT CHALLENGES IN TRIPLE-NEGATIVE BREAST CANCER (TNBC) TREATMENT

Breast cancer is the most frequently diagnosed cancer among women in the United States and Europe (1, 2). Despite the relative improvement in patient survival rates, breast cancer remains the most commonly diagnosed cancer and the second leading cause of cancer deaths in women worldwide. One of major challenges for the effective treatment of breast cancer is its intertumoral and intratumoral heterogeneity (3). Breast cancer can be initially classified into three different types based on the presence or absence of estrogen receptors (ERs), progesterone receptors (PRs), and the human epidermal growth factor receptor 2 (Her2/neu) (4). Hormone receptor-positive breast cancers that express ER and/or PR constitute approximately 60% of all breast cancers (5). The Her2/neu receptor is overexpressed in approximately 20% of all breast cancer cases; while TNBC constitute approximately 20% of breast cancer cases and are negative for the expression of ER, PR, and Her2/neu (6, 7).

Based on their molecular profile, breast cancers may also be clustered into basal-like and luminal subsets. Luminal breast cancers are more heterogeneous compared to basal cancers in terms of gene expression, mutation spectrum, copy number changes, and patient outcomes and can be further subdivided into luminal A and B subtypes (8, 9). The luminal A subtype represents 50–60% of breast cancer cases and is characterized by low histological grade and good prognosis. Luminal A cancers express ER and PR and have a low frequency of P53 mutations (9). Luminal B represents 10–20% of all breast cancers; compared with the luminal A subtype, these cancers are more aggressive; they have a higher grade, worse prognosis, and worse proliferative index. Luminal B display an increased expression of proliferation genes; they are ER+, PR+/−, Her-2+/−, and EGFR+ and have a higher frequency of P53 mutation (9). Because luminal cancers have a high frequency of PIK3CA mutations, the gene that encodes the p110α catalytic subunit of the phosphatidylinositol 3-kinase (PI3K), agents targeting the PI3K/AKT/mammalian target of rapamycin pathway may be useful for their treatment (10).

The basal-like subtype represents 10–20% of breast cancer cases. They are characterized by high proliferation, high histological grade, and poor prognosis. Basal-like cancers can be triple negative and have a high frequency of P53 mutations combined with loss of Rb1 (9, 11). However, not all basal-like cancers are triple negative; studies have shown that 5–45% of basal-like cancers express ER while 14% express Her2/Neu (12, 13). TNBC is a diverse group of malignancies and can be further categorized to different subtypes. An analysis of 21 breast cancer data sets containing 587 TNBC cases identified seven subtypes based on differential expression of a set of 2,188 genes: two basal like (BL1 and BL2), a mesenchymal (M), a mesenchymal-stem cell-like, an immunomodulatory, a luminal androgen receptor/luminal-like, and an unclassified type (14).

The deregulation of adult mammary stem cells (aMaSC) during tumorigenesis is believed to contribute to the development of TNBC. aMaSCs give rise to common progenitor cells that can differentiate either to basal progenitors that develop mature basal cells, or luminal progenitors. Disruption in the homeostasis of luminal progenitor cells may lead to the development of TNBC. Contributors in the development of TNBC include aberrantly activated signaling pathways, such as Wnt/β-catenin and Notch, transcriptional factors, like Snail, and embryonic stem cell markers including Sox2, Nanog, and Oct4. These alterations allow the restoration of proliferation capacity as well as the de-differentiation of these progenitor cells, leading to the accumulation of mutations that give rise to TNBC (15).

Traditionally, due to the lack of ER, PR, and Her2/Neu expression, the ineffectiveness of current breast cancer targeted therapies as well as due to the challenges in identifying key molecular drivers of TNBC progression, chemotherapy has been the foundation of treatment for patients with this disease over the last decades. Despite its sensitivity to chemotherapy, TNBC is associated with a higher risk of distant recurrence, high rates of metastases, higher probability of relapse and worse overall survival (OS) compared to other subtypes (16, 17).

#### COMPLEXITY OF TNBC METASTASIS

The dissemination of breast cancer cells and eventual metastatic growth to distant organs—predominantly the bone, lungs, and brain—represents a significant clinical problem, as metastatic disease is incurable and is the primary cause of death for the vast majority of TNBC patients. Metastatic spread of tumor cells is a highly complex, yet poorly understood process, and consists of multiple steps, including acquisition of invasive properties through genetic and epigenetic alterations, angiogenesis, tumor–stroma interactions, intravasation through the basement membrane, survival in the circulation, and extravasation of some cancer cells to distal tissues (18). However, disseminated cells that survive pro-apoptotic signals in their new environment often remain quiescent in secondary organs undergoing long periods of latency, also known as the dormancy period (19). It is well established that the outgrowth of metastatic cells in a foreign tissue microenvironment is a highly inefficient process and is considered as the rate-limiting step of breast cancer metastasis (20) (**Figure 1**). During this stage, breast cancer cells are usually difficult to detect and exhibit resistance to chemotherapy due to lack of proliferation (19). This remains a major clinical problem since patients, often considered as "survivors," can develop metastatic disease years later. Disseminated tumor cells (DTCs) can enter a state of dormancy in secondary organs by exiting the proliferative cycle for an indefinite period or by achieving a balanced state of proliferation and apoptosis. Successful emergence from dormancy is the result of further evolution of surviving DTCs, by accumulating molecular alterations as well as *via* permissive interactions with the tumor microenvironment (19). By acquiring these characteristics, metastatic populations can optimally adapt to the host microenvironment and initiate colonization. While significant progress has been made to highlight some of the specific processes required for the breast tumor initiation, efforts have recently been focused on elucidating the roles of critical genes, the underlying molecular mechanisms and signaling pathways involved in the fatal late stages of metastatic dissemination. These studies are of outmost importance for the development of novel effective treatments against metastasis of TNBC.

#### GENES IMPLICATED IN MULTISTEP TNBC METASTASIS

#### Local Invasion/Intravasation

Upon accumulation of genetic and/or epigenetic alterations, breast cancer cells at the primary tumor initially acquire properties, such as self-renewal, ability to migrate, and invade the surrounding normal tissues. During local invasion, breast cancer cells undergo epithelial-to-mesenchymal transition (EMT), a highly orchestrated transcriptional program, initially described during embryonic development, associated with dramatic remodeling of cytoskeleton, loss of apico-basolateral polarity, dissolution of cell–cell junctions, concomitant with downregulation of epithelial markers and upregulation of mesenchymal genes (21). This process is triggered by EMT-master regulators,

endothelial blood vessel wall to a secondary organ where they enter a prolonged dormant state by forming micrometastases. Finally, the activation of metastasiscolonizing genes and the interaction with the local microenvironment create permissive conditions for macrometastatic outgrowth. Red: metastasis promoters, green: metastasis suppressors.

such as the transcription factors Slug, Snail, and Twist to promote TNBC cell migration and intravasation in the circulation (22–24). The TGFβ pathway plays a critical role in regulating this early metastatic event. During intravasation, TGFβ promotes overexpression of musculoaponeurotic fibrosarcoma oncogene family protein K (MAFK) to induce EMT and enhance tumor formation and invasion *in vivo* (25). The TGFβ-Smad signaling axis controls the EMT step in the malignant progression of breast cancer cells either by inducing the expression of master transcriptional regulators of EMT, as described above, or by epigenetic silencing of epithelial genes, including CDH1 (26). The EMT program regulated by TGFβ/Smad signaling also involves WAVE3, a WASP/WAVE family actin-binding protein. In TNBC cells, depletion of WAVE3 expression prevented TGFβ-induced EMT phenotype (27). However, despite numerous studies using cell lines and animal models suggesting a functional role of EMT and EMT-inducing transcription factors in promoting breast cancer metastasis, the *in vivo* role and clinical relevance of this process remains controversial (28–31).

Moreover, the majority of genes implicated in TNBC metastasis have been reported to play a major role at the initial stages of cancer cell dissemination which include migration, invasion, and intravasation. This is not surprising given the fact that cancer cell dissemination is thought to be an early event during breast cancer evolution and that primary and metastatic tumor growth is likely to progress in parallel (32). For example, activation of CXCR4 receptor *via* its ligand CXCL12 or ANGPTL2 was found to induce MLK3 and Erk1/2 signaling and promote intravasation which leads to the development of lung and bone metastases (33–39). This hyperactive signaling axis may also function in multiple stages of the metastatic cascade, including angiogenesis, extravasation, and osteolysis at the secondary organ. At the same time, it is becoming increasingly clear that trans-endothelial migration and invasion of breast cancer cells in the vasculature is inhibited by metastasis suppressors, including TP63, LIFR, lysyl oxidase-like 4 (LOXL4), FOXF2, SSBP1, RAB1B, and TIEG1 (25, 40–47), suggesting that the migratory and invasive potential of breast cancer cells is ultimately determined by the balance in the activity of these molecules. The identification of numerous genes implicated in the initial stages of TNBC metastasis highlights the significant challenges for early molecular diagnosis and therapy.

#### Survival in Circulation

Upon entering the blood vessels, circulating tumor cells express proteins that have antiapoptotic and pro-survival functions which allow them to attach to and infiltrate specific secondary sites. Neurotrophic tyrosine kinase receptor TRKB was shown to inhibit anoikis, a form of cell death caused by lack of adhesion, *via* the PI3K/Akt pathway. These studies indicated that TRKB induces survival and proliferation of breast cancer cells to promote infiltration in the lymphatic and blood vessels and colonization in distant organs (48). In TNBC cells, brain-derived neurotrophic factor (BDNF) binds and activates TRKB receptor to regulate a network consisting of metalloproteases and calmodulin and thus modulate cancer–endothelial cells interaction. Importantly, Erk1/2 inhibitors were able to block the BDNFinduced phenotype, suggesting that blocking this pathway may be explored for therapeutic purposes against TNBC metastasis (49). In addition, the binding of platelets with circulating breast cancer cells has been shown to essential for their survival, evasion of pro-apoptotic signals, whereas interfering with this interaction inhibits the development of lung metastasis in TNBC mouse models (50, 51).

#### Extravasation in Distal Sites

Many of the genetic alterations found to be involved in intravasation are also implicated in extravasation (**Table 1**) since, in large part, these two processes are considered "mirrored" to each other. The TGFβ pathway plays an important role in regulating both these metastatic steps. More specifically, TGFβ induces the assembly of a mutant-p53/Smad protein complex to inhibit the function of the metastasis suppressor TP63 and promote cell migration and invasion (40). During extravasation, TGFβ induces angiopoietin-like 4 (ANGPTL4) expression *via* the Smad signaling pathway; the increased levels of ANGPTL4 enhance the retention of cancer cells in the lungs by disrupting vascular endothelial cell–cell junctions, thus increasing the permeability of lung capillaries to facilitate trans-endothelial passage of breast cancer cells (52). Moreover, targeting the decoy interleukin-13 receptor alpha 2 (IL13Ra2) upregulates the metastasis suppressor TP63 in an IL13-mediated, STAT6 dependent manner and impairs extravasation of basal-like breast cancer cells to the lungs (41). Several reports also highlight the importance of the synergistic effects of genes in promoting metastasis by regulating specific stages of the process. For example, EREG, COX2, MMP1, and MMP2 can collectively promote metastatic extravasation to the lungs. These four genes were found to be overexpressed in TNBC cells independently of VEGF. Individual reduction of each gene or their silencing in different combinations produced limited effects on tumor growth *in vivo* while concurrent silencing of all four achieved nearly complete growth abrogation (53).

#### Metastatic Colonization

Following extravasation and infiltration at the secondary site, a genetic program is initiated so that cancer cells can escape dormancy and form micro and macrometastatic tumors. Initially, EMT plasticity and the reversal to MET phenotype have been shown to be important for metastatic colonization (113). During this process, epithelial phenotype becomes re-established through miR-200-mediated downregulation of ZEB1, SIP1 to promote metastatic colonization (114, 115). Also, breast DTCs in the bone marrow gain the ability to form typical osteolytic metastases by producing parathyroid hormone-related protein (PTHLH), tumor necrosis factor-α (TNFα), interleukin-6 and/or interleukin-11. These factors stimulate the release of receptor activator of nuclear factor-κB ligand (RANKL) from osteoblasts which induces osteoclast formation (33, 58, 83, 116). Furthermore, inflammation in the lung microenvironment could also be responsible for triggering the escape of metastatic breast cancer cells from latency leading to metastatic colonization (117). A subset of genes contributing to primary tumor growth can also promote survival and growth at the secondary site. Chemokines CXCL1/2 mediate chemoresistance and lung metastasis by attracting myeloid cells into the tumor, which produce low molecular weight calcium-binding proteins S100A8/9 that enhance cancer cell survival by binding to the receptor for advanced glycation end products (RAGE) (59). Another calcium binding protein, S100A7 has been found to enhance tumor growth and metastasis, by binding to RAGE and activating Erk and NFκB signaling (88, 90). Furthermore, fibroblast growth factor receptor (FGFR) was shown to trigger pro-survival signals through PI3K/Akt signaling and promote outgrowth of metastatic breast cancer cells to the lungs (62). However, it needs to be highlighted that cellular and genetic context among cancers influences whether proteins act as tumor suppressors or metastasis promoters. One controversial example is LOXL4 which has been shown to recruit bone marrow-derived cells and facilitate colonization of TNBC to the lungs *via* a HIF1α-dependent mechanism (118). However, in another study, knockdown of LOXL4 expression in TNBC cells promoted primary tumor growth and lung metastasis which was associated with thickening of collagen bundles and remodeling of the extracellular matrix (ECM) within tumors (25). Overall, it is noteworthy that while some genes have been associated only with TNBC metastasis so far (i.e., TIEG1, MAFK, MLK3, SDPR), the majority is also involved in other tumor types, suggesting a more fundamental role in cancer progression.

### CONCLUDING REMARKS ON CURRENT AND FUTURE PERSPECTIVES ON TNBC METASTASIS THERAPY

Due to their molecular heterogeneity, there are no drugs that can target the entire spectrum of TNBC tumors and each subtype is vulnerable to specific therapeutic approaches. Despite the lack of FDA-approved targeted therapies for TNBC to date, ongoing clinical trials are assessing the efficacy of single or combinatorial approaches that tackle different TNBC molecular alterations. Up to 20% of TNBC have been associated with germ-line mutations in BRCA1 (119). TNBC tumors with loss of function of BRCA1 or BRCA2 are sensitive to poly(ADP-ribose) polymerase inhibitors and alkylating agents that induce DNA double-strand breaks (120). Olaparib has been the most successful PARP inhibitor


(*Continued*)

Mechanisms and Therapy of TNBC Metastasis

#### TABLE 1 | Continued


(*Continued*)

TABLE 1 | Continued

Metastasis-promoting genes


**27**

Mechanisms and Therapy of TNBC Metastasis

#### TABLE 1 | Continued


*A comprehensive list of genes implicated in various stages of the metastatic cascade, their reported functions, upstream or downstream regulatory signaling pathways involved, gene ontology, as well as the secondary organs which become affected.*

against BRCA-mutated TNBC, inducing partial responses in 54% of patients when administered as a single agent (121) and an overall response rate of 88% when combined with carboplatin (122). Anti-androgens as well as FGFR inhibitors have been tested in clinical trials against TNBCs that are androgen receptor-positive or harbor FGFR amplification, respectively (123, 124). Gammasecretase inhibitors that block the NOTCH pathway are currently in clinical trials for TNBC patients with upregulated NOTCH signaling (125). All together clinical trials have shown that each agent alone provides small or no benefit in TNBC patients suggesting that further effort is needed to discover novel targets of TNBC and to identify each patient's molecular profile that will lead to a more individualized treatment.

Toward this goal, some of the metastasis-promoting genes reported here could be further exploited for the future development of promising targeted therapies. Since local invasion, intravasation and possibly extravasation are thought to occur relatively early in the metastatic process (32), a plausible strategy would be to target dormancy and the outgrowth of macrometastatic tumors in distal organs. Since this final stage is considered the critical "rate-limiting" step of the "invasion-metastasis" cascade requiring even years to be completed, it provides a window of opportunity for effective therapy. Therefore, different approaches could aim against "druggable" molecules that facilitate metastatic colonization, such as overexpressed receptors or secreted molecules (i.e., CXCL1/2, FGFR, TGFβ1, WNT1, ANGPTL2, CSF2, RANKL), which target commonly deregulated signaling networks at this late-stage (**Table 1**). Ongoing clinical trials are evaluating the efficacy of the TGFβR1 inhibitor LY2157299 with paclitaxel (NCT02672475), whereas the FGFR inhibitor Lucitanib is also under testing (NCT02202746) for patients with metastatic TNBC. The ultimate goal would be, if not to completely eliminate dormant metastatic breast cancer cells, to prolong dormancy period and hopefully transform this stage into a chronic inactive cancer cell state.

Importantly, recent studies have shown that tumor cells are able to evade immune responses by activating negative regulatory pathways, also known as immune checkpoints, that block T-cell activation through cytotoxic T-lymphocyte protein 4 (CTLA4) or *via* binding of the programmed cell death protein 1 (PD1) receptor expressed on T-cell surface to the PDL1 ligand expressed by cancer cells in response to various cytokines (126). The recent development and FDA approval of anti-CTLA4, anti-PDL1, and anti-PDL1 monoclonal antibodies that elicit

#### REFERENCES


antitumor clinical responses in a variety of solid cancers created enthusiasm for cancer therapy (127). Currently, several clinical trials are underway to evaluate the efficacy of this approach in TNBC as well (128).

However, a major clinical problem is that breast cancer is considered one of the most desmoplastic tumor types due to the production of excessive amounts of ECM components, such as collagen and hyaluronan, which generate mechanical stresses within the growing tumor (129). This results in blood vessel compression, hypoperfusion, and hypoxia which promote cancer progression and metastasis as well as hinder drug delivery (130). Therefore, targeting components of the tumor microenvironment has also been recently proposed as another promising strategy for TNBC therapy by improving tumor penetration and delivery of cytotoxic drugs (131). For example, targeting of cancer-associated fibroblasts using pirfenidone, an FDA-approved drug for idiopathic pulmonary fibrosis, has been shown to suppress metastasis of TNBC in combination with doxorubicin (132). This effect is likely to be mediated through remodeling of tumor microenvironment which reduces ECM components through suppression of TGFβ signaling, improves perfusion and delivery of chemotherapy (133). Similar effects have also been demonstrated using the anti-fibrotic drug Tranilast or the anti-hypertensive drug Losartan in combination with chemotherapy or nanotherapy in mouse models for TNBC (134–136).

In conclusion, this evidence suggests that efforts in the near future should be focused toward the development and testing of novel anti-metastatic targeted therapies for late-stage TNBC that could be used in combination with existing chemotherapies, immunotherapies as well as with microenvironment-remodeling agents that can improve drug penetration and overall therapeutic efficacy.

#### AUTHOR CONTRIBUTIONS

CN and PB wrote the paper and helped with illustrations. PP conceived the theme, wrote the paper, and prepared illustrations.

### FUNDING

This work was supported by the Cyprus Research Promotion Foundation (KOULTOURA/ΒΡ-ΝΕ/0415/03).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Neophytou, Boutsikos and Papageorgis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Growth and immune evasion of Lymph node Metastasis

*Dennis Jones1,2, Ethel R. Pereira1,2 and Timothy P. Padera1,2\**

*1Edwin L. Steele Laboratories for Tumor Biology, Department of Radiation Oncology, MGH Cancer Center, Massachusetts General Hospital, Boston, MA, United States, 2Harvard Medical School, Boston, MA, United States*

Cancer patients with lymph node (LN) metastases have a worse prognosis than those without nodal disease. However, why LN metastases correlate with reduced patient survival is poorly understood. Recent findings provide insight into mechanisms underlying tumor growth in LNs. Tumor cells and their secreted molecules engage stromal, myeloid, and lymphoid cells within primary tumors and in the lymphatic system, decreasing antitumor immunity and promoting tumor growth. Understanding the mechanisms of cancer survival and growth in LNs is key to designing effective therapy for the eradication of LN metastases. In addition, uncovering the implications of LN metastasis for systemic tumor burden will inform treatment decisions. In this review, we discuss the current knowledge of the seeding, growth, and further dissemination of LN metastases.

Keywords: lymphatics, lymph node, tumor, metastasis, immunity

INTRODUCTION

#### *Edited by:*

*Triantafyllos Stylianopoulos, University of Cyprus, Cyprus*

#### *Reviewed by:*

*Marc Achen, Peter MacCallum Cancer Centre, Australia Nancy H. Ruddle, Yale University, United States*

*\*Correspondence:*

*Timothy P. Padera tpadera@steele.mgh.harvard.edu*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

*Received: 18 November 2017 Accepted: 01 February 2018 Published: 21 February 2018*

#### *Citation:*

*Jones D, Pereira ER and Padera TP (2018) Growth and Immune Evasion of Lymph Node Metastasis. Front. Oncol. 8:36. doi: 10.3389/fonc.2018.00036*

Metastasis is the leading cause of death from cancer (1) and represents a challenging clinical problem. Lymph nodes (LNs) are common sites of metastasis and nodal disease predicts increased mortality in many cancer types. Meanwhile, LNs are critical for initiating antitumor immune responses. Thus, cancer cells that have metastasized to LNs must escape immune detection to avoid destruction. The process of lymphatic metastasis is regulated at several steps and by several different molecules (**Figures 1** and **2**, respectively), beginning with the orchestration of lymphangiogenesis and preparation of a LN microenvironment favorable for tumor growth (premetastatic niche). Cancer cells then invade tumor-associated lymphatic vessels at the primary site *en route* to tumor-draining LNs (TDLNs), where they survive and grow. In a metastatic node, immunological destruction of cancer cells depends on the degree of cancer cell immunogenicity and the extent of nodal immunosuppression. Similar to primary tumors, cancer cells in LNs shape their interactions with the host immune system by controlling the infiltration and reactivity of immune cells. The local microenvironment of the LN also dictates the growth and response of LN metastases to therapeutic intervention. For example, only a small fraction of drugs delivered systemically accumulate in LNs (2). Identifying effective therapy for LN metastases takes on new urgency as cancer cells in LNs have also been proposed to disseminate to other metastatic sites by lymphatic or hematogenous routes. In this review, we summarize recent progress in the understanding of lymphatic metastasis and metastatic outgrowth. We also discuss the consequences of lymphatic metastasis and therapeutic efforts to target LN lesions in experimental mouse models and humans.

### LYMPHATIC ENDOTHELIAL CELLS (LECs) AND TUMOR IMMUNITY

#### Mediators of Immunosuppression

Recent studies suggest that in addition to serving as a portal for tumor dissemination, lymphatic vessels facilitate tumor growth through immune suppression (3). To generate an antitumor T cell response, migratory dendritic cells (DCs) from primary tumors cross-prime naïve T cells in TDLNs (4).

Figure 1 | Progression of lymphatic metastasis from primary tumor to tumor-draining LN (TDLN). Primary tumors induce lymphangiogenesis to facilitate lymphatic metastasis and release of immunomodulatory molecules, including exosomes, which lead to immunosuppression of TDLNs. Lymph node (LN) lymphatic endothelial cells (LECs) capture tumor antigen and tolerize T cells *via* programmed death-ligand 1 expression. Tumor-associated lymphatic vessels and tertiary lymphoid organs have been implicated in immune suppression and immune activation. High endothelial venules found in primary tumors can allow infiltration of naive T cells that may further differentiate into effector T cells. Tumor-associated lymphatic vessels recruit both cancer cells and immune cells by releasing chemoattractants (see Figure 2). Cancer cells, T cells, and dendritic cells enter lymphatic capillaries and migrate through collecting lymphatic vessels to LNs. Cancer cells in lymphatic vessels can attach to the lymphatic endothelium en route to LNs. Active mechanisms, such as CCL1/CCR8 signaling, control cancer cell entry into the LN. Polyclonal cancer cells proliferate to form a metastatic lesion that invades deeper into the LN parenchyma, where it can grow and replace LN tissue in the absence of new blood vessel growth. The immune response to a growing metastatic lesion is limited; some immune cells are excluded from LN lesions, while other immune cells are present, but unable to eliminate cancer cells (not shown). Some cancer cells may exit through the efferent lymphatic vessel and seed secondary draining LNs. Recent evidence suggests LEC sphingosine-1-phosphate (S1P) helps shape the antitumor immune response.

The adhesion ligand Mac-1 on DCs can bind to the adhesion molecule intercellular adhesion molecule-1 (ICAM-1), which is upregulated on endothelial cells of collecting lymphatic vessels in response to inflammation (5), including inflammation generated by the tumor microenvironment. This Mac-1/ ICAM-1 interaction inhibits DC maturation (5) and may blunt the ability of DCs in an inflamed tumor microenvironment to prime antitumor T cells. LECs further inhibit antitumor T cell responses by inducing tolerance to tumor antigens. LECs can scavenge tumor and other peripheral antigens and cross-present them to CD8 T cells, but LECs lack co-stimulatory molecules needed for full activation of CD8 T cells (6). Programmed death-ligand 1 (PD-L1), a ligand for the T cell inhibitory receptor PD-1, is expressed on tumor-associated lymphatic vessels (7) and LN LECs (8). The engagement of PD-L1 on LECs with T cell PD-1 induces CD8 T cell tolerance to tumor antigens (7, 8). Given the paucity of antitumor T cells and functional lymphatic vessels within some tumors (9, 10), the degree of CD8 T cell interaction with tumor-associated lymphatic vessels and their degree of inhibiting antitumor immunity is unclear. Recent studies show that the presence of LECs inside tumors make the tumors more responsive to anti-PD-1 therapy, suggesting lymphatics can have a potent inhibitory effect on T cell function (11).

In normal physiology, LECs produce sphingosine-1-phosphate (S1P), which is secreted into lymph and controls lymphocyte exit from LNs (12). Neutralization of systemic S1P with a therapeutic antibody suppresses lung metastasis (13). More recently, a genome-wide functional screen identified the S1P transporter spinster homolog 2 (*Spns2*) as a regulator of metastatic colonization in animals with experimental lung metastases (14). Spns2 is expressed on LECs and is critical for LEC release of S1P (15). Global and lymphatic-specific deletion of *Spns2* decreased pulmonary metastases following intravenous tumor cell injection, with fewer total T cells present in lungs relative to WT mice (14). However, a higher proportion of effector memory T cells to regulatory T cells (Tregs) were found in the lungs of *Spns2* deficient animals (14). Coupled with enhanced KLRG1<sup>+</sup>, CD69<sup>+</sup>CXCR3<sup>+</sup> T cell activation, these findings were suggestive of an adaptive immune response against lung metastases (14). In addition, the natural killer (NK) cell population in the lungs of *Spns2* deficient animals was increased and limited the growth of lung metastases, even after CD8 T cell depletion (16). By contrast, S1P signaling has the potential to promote antitumor T cell responses. LECproduced S1P appears to function not only as a regulator of lymphocyte circulation, but also supports naïve CD4 T cell survival by maintaining their mitochondrial content through an S1P1 receptor-dependent mechanism (16). Targeting S1P signaling to

decrease metastatic burden requires a better understanding of its temporal and spatial role in shaping antitumor T cell responses.

The chemokine CCL21 is produced by LECs and mobilizes DCs to LNs (17). Although CCL21 can also establish an immunosuppressive tumor microenvironment (18), it was recently shown that tumor-associated lymphatic vessels also facilitate naïve T cell recruitment into melanoma tumors through CCL21 production (11). The presence of T cells allows local immune priming and the ability to unleash potent antitumor activity in response to vaccination or immune checkpoint blockade (11). So while LECs and CCL21 themselves help suppress antitumor immune responses, their ability to recruit T cells sensitizes lymphangiogenic tumors to immunotherapy. Cancer patients with elevated VEGF-C (indicative of increased lymphangiogenesis and lymphatic metastasis) typically have poor prognosis. Based on recent data, VEGF-C [which upregulates CCL21 (19)] may now be used as a biomarker in these patients to predict their response to immunotherapy. The ability of T cells with effector function to leave primary tumors and travel to TDLNs through afferent lymphatic vessels also suggests that tumor-associated lymphatic vessels may assist in dampening tumor growth (20).

The above examples illustrate the potential for contextdependent benefit of inhibiting, altering or utilizing intrinsic properties of LECs to maximize effective antitumor immune responses.

### Lymphangiogenesis

Lymphangiogenesis is a hallmark of many solid tumors. The expansion of the lymphatic network is primarily mediated by VEGF-C and its receptors VEGFR-2/3 (21). VEGF-D also binds VEGFR-2/3 and potently induces lymphangiogenesis (22, 23). Both VEGF-C and VEGF-D expression correlate with increased LN metastasis (24). Furthermore, *Vegf-d* deficient mice displayed less lymphatic metastasis relative to tumor-bearing wild type mice (25). Many preclinical studies have shown prevention of lymphatic metastasis by blocking VEGF-C or VEGF-D-mediated lymphangiogenesis (26–28). Clinically, this point of intervention is challenging as lymphangiogenesis is an early event in the natural history of cancer progression and many patients will already have LN metastases on initial presentation. However, additional opportunities may exist to target lymphangiogenesis. For example, lymphangiogenesis was identified as a mechanism of tumor resistance to paclitaxel chemotherapy in mice (29). In response to paclitaxel in mouse models, tumor-infiltrating macrophages secrete cathepsin, which activates heparanase. Active heparanase, by unknown mechanisms, increases both VEGF-C transcription and tumor invasiveness (30). In another study, VEGFR-3 reporter mice were used to image lymphangiogenesis in distant LNs, liver, lungs, and spleens in tumor-bearing mice. Following tumor resection, VEGFR-3 levels declined but reemerged before tumor relapse, suggesting a defined window of opportunity to inhibit lymphangiogenesis in distal premetastatic organs (31).

Although there are a plethora of targets that exist to inhibit lymphangiogenesis (32), many studies find that VEGF-C/-D-VEGFR-2/3 signaling directly or indirectly promotes lymphangiogenesis in response to a wide range of stimuli (22, 24, 25). Surprisingly, few clinical trials targeting the lymphatic endothelium in cancer are ongoing, although several small molecules that non-selectively target VEGFR-3-mediated lymphangiogenesis are approved for cancer indications (33). A Phase I study was recently completed that assessed the effect of LY3022856 (IMC-3C5), a monoclonal antibody targeting human VEGFR-3, on colorectal cancer (CRC) (34). While LY3022856 was well tolerated at the given dose, minimal antitumor benefit was noted in patients with CRC. The impact of LY3022856 on tumor lymphangiogenesis and lymphatic metastasis was not assessed. As mentioned earlier, inhibiting lymphangiogenesis through targeting VEGF-C/D or VEGFR-3 is also complicated by the uncertainty of the effect that lymphatic vessels have on antitumor T cell responses (11).

Independent of lymphatic vessel growth and lymphangiogenesis, VEGF-C can promote cancer metastasis by disruption of the vascular endothelial cadherin/β-catenin complex at intercellular junctions of LECs (35). The authors concluded that enhanced permeability of intestinal lymphatic vessels caused by VEGF-C can increase CRC transmigration and metastasis. Thus, how therapies targeting VEGF-C/D signaling will impact cancer progression will depend on the specific context of the disease as well as other therapies being used in conjunction.

#### ESTABLISHMENT OF PREMETASTATIC NICHE IN TDLNs

#### Extracellular Vesicles (EVs)

In addition to the delivery of cells, lymphatic vessels deliver primary tumor-derived soluble and vesicle-associated factors to condition TDLNs before the arrival of cancer cells. Exosomes, a type of EV, were shown to modulate the immune and stromal response in TDLNs (36, 37). Melanoma cells injected into the footpad and taken up by local lymphatic vessels had a similar distribution pattern as premetastatic melanoma-derived exosomes previously injected into the footpad, suggesting exosomes influence the recruitment of cancer cells to the LN (38). Mechanistically, exosomes upregulated host genes that promoted the retention, recruitment, and progression of LN metastases (38). It is unclear what components of exosomal cargo (e.g., mRNA, miRNA, or proteins) were necessary for changes in nodal gene expression. Melanoma exosomes have also been shown to enhance metastasis by "educating" and mobilizing bone marrow-derived cells to metastatic sites (36), including LNs, where they facilitated cancer cell invasion. Melanoma-derived EVs were identified in afferent lymphatic vessels of patients (39). Cocultures of EVs from human melanoma cells with DCs resulted in inhibition of DC maturation (39). In premetastatic LNs, CD169<sup>+</sup> subcapsular sinus (SCS) macrophages capture tumor-derived extracellular vesicles (TeVs) (40) and protect host LNs from TeV-mediated immunosuppression. TeV release was required to accelerate tumor progression after the investigators depleted host macrophages. It is, however, unclear how TeVs can escape capture by SCS macrophages to deliver their payload under normal conditions. In contrast to the pro-tumor effects of EVs, migratory DCs acquire tumor-secreted vesicles released by circulating tumor cells in the lung (41). From here, the vesicleloaded DCs migrate to mediastinal LNs to interact with and potentially activate antitumor T cells to limit metastatic growth. Taken together, these data suggest that TeVs have immune regulatory functions as well as help initiate and support the growth of LN metastases.

#### Lymphatics in a Premetastatic LN

It is known that lymphangiogenesis occurs in premetastatic LNs (42, 43). Nodal lymphangiogenesis has been shown to be tumor antigen independent and B cell dependent (42, 44) through production of VEGF-A and VEGF-C (43–46). Recently, midkine, a heparin-binding factor produced by tumor cells, was identified as a critical factor for mTOR-dependent lymphangiogenesis in premetastatic sites including skin, LN, spleen, and lung (31). Furthermore, midkine mediated tumor cell adhesion to LECs and promoted tumor colonization in distant organs.

It is unclear how LN lymphangiogenesis results in metastasis. One hypothesis is that LN lymphangiogenesis may lead to more efficient delivery of cancer cells to LNs and distant organs (47). This may be facilitated by the increased lymph flow that accompanies increased LN lymphangiogenesis (42). Increased lymphatic drainage from primary tumors was also associated with LN enlargement (48) and nodal remodeling, which may alter the distribution of antigen and soluble factors. Increased lymphatic drainage also coincided with collagen and hyaluronic acid deposition in premetastatic TDLNs of B16F10 tumors. Parental B16 melanoma cells failed to increase TDLN matrix remodeling and were inefficient at metastasis, suggesting that in addition to lymphangiogenesis, an increase in TDLN matrix may be a prerequisite for formation of LN metastases (48, 49).

### Fibroblast Reticular Cells

The LN contains an array of stromal cells, including fibroblastic reticular cells (FRCs). Much is known about the tumor-promoting effects of cancer-associated fibroblasts, but few reports have characterized the FRC response to cancer cells. Transcriptional profiling of FRCs from non-tumor-bearing animals revealed abundant expression of chemokines critical for lymphocyte recruitment, including *CCL19*, *CCL21*, *CXCL12*, and *CXCL13* (50). FRCs also produce several forms of collagen (50), indicative of their role in forming the conduit system that delivers small antigens deep into the LN for antigen presentation (51). FRCs express genes necessary for MHC class 1/2 presentation (50) and can present peripheral antigens to T cells (52). Similar to LECs, FRCs contribute to peripheral tolerance by facilitating deletional tolerance (52, 53) and dampening effector T cell proliferation (54, 55). A recent transcriptional analysis revealed FRCs in premetastatic LNs are "reprogrammed" to favor tumor growth (56). In spontaneous and orthotopic models of melanoma, TDLN FRCs proliferated, but produced less *IL-7* and *CCL21*, which are critical for T cell survival and guidance, respectively. The reduction in *IL-7* and *CCL21* resulted in disruption of the TDLN architecture, with loss of clear delineation between B and T cell zones. In a separate study, the loss of FRC *CCL21* in the TDLN was associated with disorganized T cell and B cell zones in premetastatic LNs (57). The perturbation of LN architecture due to altered FRC signaling molecules suggests altered immune responses to tumors. Since LNs are priming sites for adaptive immune responses, the disordered LN architecture may fail to elicit systemic protection from subsequent heterogeneous cancer cell clones that arrive in the TDLN (56). In metastatic LNs, collagen production was increased relative to tumor-free LNs (58). Although unclear whether recruited fibroblasts, FRCs, or cancer cells are the source of additional collagen, the investigators speculate that the increased density of collagen fibers may allow cancer cells to adhere and migrate within metastatic LNs. It is unknown how tumor cells influence FRC transcriptional status.

#### TUMOR CELL MIGRATION TO LNs

Cancer cells enter lymphatic vessels and travel with the lymph to establish LN metastasis (59). Cancer cells may actively migrate into lymphatic capillaries in response to molecular cues (19, 60) or they may passively enter into lymphatic capillaries (19, 60). Metastasis to the LN likely depends on a combination of intrinsic cancer cell properties and signals in the tumor microenvironment. VEGF-C and lymphatic flow both upregulate CCL21 in lymphatic endothelium (19, 61), attracting CCR7<sup>+</sup> tumor cells (62). In a triple-negative breast cancer model, CCL21 was sufficient to recruit RORγt <sup>+</sup> innate lymphoid cells (ILCs) into the primary tumor and promote metastasis to LNs (63). Furthermore, CXCL13 was required for clustering of ILCs and induction of epithelial–mesenchymal transition, likely driving invasion of cancer cells. In breast cancer patients, the presence of ILCs was significantly associated with lymphatic invasion at the primary tumor.

Several studies have shown that another chemokine, CXCL12, facilitates lymphatic metastasis of CXCR4<sup>+</sup> tumor cells (64–66). CXCL12 expression is found on lymphatic vessels within primary tumors and guides CXCR4*<sup>+</sup>* melanoma cells toward lymphatic vessels. Migration and invasion of CXCR4<sup>+</sup> papillary thyroid carcinoma cells are dependent on CXCL12, which was produced by senescent cancer cells at the invasive border (67). These senescent cells invaded lymphatic vessels and persisted in metastatic foci, suggesting that they may promote lymphatic metastases. CXCR4 is also expressed on the surface of LECs (68) and is critical for lymphangiogenesis through CXCL12 stimulation, independent of the VEGFR-3 pathway (68). Thus, targeting the CXCR4/CXCL12 may provide a dual benefit of inhibiting cancer cell migration and lymphangiogenesis to curb lymphatic metastasis.

After entry of cancer cells into lymphatic vessels, it is thought that lymph flow allows cancer cells to traverse the collecting lymphatic vessel network until they reach TDLNs (59). Based on 3D modeling, it was predicted that smaller breast cancer cells may have a survival advantage over larger breast cancer cells in the lymphatic circulation because of the lower wall shear stress that they encounter (69). Several studies have shown that inflammation causes dilation and inhibits contractile ability of collecting lymphatic vessels (70, 71). More work needs to be done to determine if tumor-induced collecting lymphatic dilation (10, 22, 59) or reduced contraction (72) enhances tumor cell dissemination by decreasing the shear stress on cancer cells. It is known that tumor cells can arrest within lymphatic vessels while "in-transit" to LNs (73). Compromised barrier integrity of lymphatic vessels may allow arrested cancer cells to escape lymphatic vessels and form metastases (74, 75). Additional characterization of the mechanism of how tumor cells attach to lymphatic endothelium and grow within lymphatic vessels is needed to treat in-transit metastases.

Recently, the chemokine CCL1 and its receptor CCR8 were demonstrated to be important for melanoma cell entry into TDLNs. CCL1 is produced by SCS LECs and mediated entry of CCR8<sup>+</sup> melanoma cells into LNs (60). Tumor cells in the SCS can also bypass the LN parenchyma and drain through cortical and medullary sinuses to exit LNs *via* efferent lymphatic vessels (76). The enzyme lipoxygenase 15 (ALOX15) metabolizes arachidonic acid to 12(S)-hydroxyeicosatetraenoic acid [12(S)-HETE] and 15(S)-hydroxyeicosatetraenoic acid [15(S)-HETE]. Cancer cell-derived 12(S)-HETE forms discontinuities in the walls of lymphatic vessels, allowing LN metastases to invade nodal lymphatic vessels (77). The fate of these cancer cells is unclear, although TDLN lymphangiogenesis has been reported to be involved in further lymphatic spread of human breast cancer (78) and the presence of lymphatic vessel invasion by LN metastases is associated with worse survival (79). It is possible that cancer cells circulate to additional nodes through lymphatic vessels and eventually enter the systemic circulation through the thoracic duct.

#### IMMUNE EVASION IN TDLNs

#### Macrophages

Lymph node SCS macrophages are the first line of defense against tumor cells entering the LN. SCS macrophages capture microbes, antigen–antibody complexes and dead cancer cells for delivery of these antigens to nearby immune cells (80, 81). In premetastatic LNs, an experimental antigen (a fluorescent protein overexpressed in tumor cells) from the primary tumor was captured by SCS macrophages and distributed to follicular DCs, resulting in antibody production against the antigen (82). SCS macrophages can also directly cross-present tumor antigens to CD8 T cells (81). Sinus macrophages in regional LNs of CRC patients made direct contact with CD8 T cells and a high density of sinus macrophages is associated with increased overall survival (83). On the other hand, tumor-associated macrophages are often associated with poor prognosis and promotion of tumor growth (84). Strategies to deplete TAMs include targeting colony-stimulating factor 1-receptor (CSF1-R) (85), which controls macrophage chemotaxis. Interestingly, an increase in the burden of LN metastases was found following treatment with an anti-CSF1-R antibody (86). This increase in metastatic burden was associated with the loss of SCS macrophages due to anti-CSF1-R therapy (86). Tumor-promoting (M2) macrophage depletion strategies should examine the effect on SCS macrophages to avoid unintended growth of LN metastases.

#### Neutrophils

Neutrophils, such as macrophages, are heterogeneous and have been reported to have either pro-tumor or antitumor phenotypes in primary tumors (87). High levels of tumor-associated neutrophils are associated with LN metastasis and poor prognosis (88). Granulocyte colony-stimulating factor was necessary to expand and polarize neutrophils to an immunosuppressive phenotype in mice bearing mammary tumors (89). The immunosuppressive neutrophils, whose expansion was also driven by IL-17, were able to suppress cytotoxic T cells and facilitate the establishment of LN metastases. Neutrophils can also secrete pro-inflammatory leukotrienes and initiate LN metastases *via* leukotriene receptors (LTR) on cancer cells (90). LTR expression was found in a cohort of primary and LN tumors from breast cancer patients. More research is needed to characterize the phenotype of neutrophils found in LN metastases and their role in metastatic progression.

### T Cells

The TDLN often fails to produce effective antitumor immunity and instead tolerizes the patient to tumor antigens (91). The mechanisms that induce systemic tolerance include cross-presentation of tumor antigen by tolerogenic antigen-presenting cells, apoptosis of antigen-presenting cells (92), and suppression of antitumor T cells by an expanded pool of Tregs. Experimentally, subcutaneous injection of B16 melanoma resulted in tolerized CD8 T cells and lethal metastatic outgrowth (93). However, B16 cells implanted directly into LNs—without a primary tumor were rejected, supporting previous evidence (94) that showed the primary tumor exerts a tolerogenic effect on the TDLN (94). However, initial metastatic deposits in lymphoid organs are important for the induction of antitumor CD8 T cells and tumor rejection (94). Notably, tumor cells injected into LNs using different cancer models have been shown to persist and disseminate to distant organs (95). These differences may be explained by factors such as the immunogenicity and antigen presentation capabilities of different cancer cells. Increased LN metastasis was found in cancer patients with tumor downregulation of MHC I (96).

### NK Cells

Natural killer cells are cytotoxic lymphocytes of the innate immune system that are often recruited to tumors including prostate (97), melanoma, kidney, liver, and breast (98). NK cells are able to recognize and eliminate cells with aberrant ligand or altered/absent MHC expression (99, 100). However, NK cell infiltration into primary tumors is limited; NK cells that enter tumors are often found in primary tumor stroma and lack direct contact with cancer cells (98). Likewise, NK cells in metastatic LNs were adjacent to metastatic melanoma lesions (101). Moreover, NK cells isolated from metastatic human LNs showed significantly reduced cytotoxicity (101, 102). NK cells that were isolated from metastatic LN lesions and stimulated with IL-2 or IL-15 displayed more efficient lysis of cancer cells (101), suggesting that the LN tumor microenvironment suppresses NK cell function. Thus, despite their presence in the TDLN, immunosuppressed NK cells may lack the ability to eliminate cancer cells.

### B Cells

The number of B cells in premetastatic TDLNs is significantly increased (42). EVs from tumor cells increased immunosuppressive B cells in premetastatic LNs (40). Depletion of SCS macrophages in tumor-bearing animals permitted interactions of TeVs with B cells and led to increased antibody production. Although the tumor antigen specificity of the antibodies are unknown, transfer of these antibodies to wild-type recipient mice correlated with enhanced tumor growth compared with antibody transfer from tumor-bearing animals without disruption of SCS macrophages. Regulatory B cells (Bregs) were recently identified in mouse and human blood and secondary lymphoid organs (103). Bregs can secrete immunosuppressive cytokines, such as TGF-β and IL-10 (104), to dampen the effector activity of antitumor T cells. However, the data conflict on whether Bregs support or suppress tumor growth (105). The phenotypic markers that identify Bregs also remain unclear.

Together, these data suggest that multiple immune cell types in premetastatic and metastatic LNs have suboptimal killing activity for cancer cells. Identifying molecular targets to reverse immune suppression of several cell types will be of therapeutic benefit in treating metastatic cancer.

### ECTOPIC LN IMMUNITY

High endothelial venules (HEVs) and their homeostatically associated chemokines are crucial for entry of naïve lymphocytes into LNs (106). Ectopic HEVs in primary human breast cancer and melanoma tumors allow lymphocytic intravasation and predict a favorable prognosis (106, 107). LTα3–TNFR signaling, not LTβR, was critical for the generation of HEV-like vasculature expressing peripheral node addressin (PNAd) in models of lung cancer and melanoma (108) studied by Peske et al. The PNAd<sup>+</sup> vasculature was critical for infiltration of naïve T cell into tumors. Growth of HEV-containing tumors was delayed although treatment with the S1P antagonist FTY720 retained tumor-specific T cells in secondary lymphoid organs (108). These data suggest that naïve T cells can differentiate into effector T cells within primary tumors.

Tertiary lymphoid organs (TLOs) are aggregates of immune cells that mimic the structure and function of LNs. They are formed in several diseases associated with chronic inflammation, including cancer (109). TLOs include lymphocytes, lymphatic vessels, and HEVs (110). TLOs within tumors are associated with improved outcomes for patients and function as sites of immune priming for the generation of antitumor lymphocytes. In animal models of pancreatic and breast cancer, tumor-associated blood vessels developed HEV markers and formed TLO-like structures in response to a combination of antiangiogenic (DC101) and immune checkpoint (PD-L1) therapy (111). Formation of LTβR signaling-dependent HEVs resulted in enhanced infiltration and activation of antitumor T cells, and better antitumor responses (111). By contrast, the presence of TLOs and ectopic HEVs has been shown to promote tumor growth in other experimental conditions. Finkin et al. found that the cytokine-rich environment of liver TLOs promoted the growth of hepatocellular carcinoma progenitor cells in a hepatocellular carcinoma model (112). From studies using a murine model of lung cancer, Joshi et al., found that TLOs in tumor-bearing lungs were a site of antigen presentation (113). However, the abundant Tregs within the TLOs presumably suppressed antitumor T cell responses. The use of immune checkpoint inhibitors, such as the anti-CTLA-4 antibody ipilimumab, can lead to the depletion of Tregs (114) and may enhance the antitumor function of effector lymphocytes in tumors and TLOs with a high Treg: effector T cell ratio. In a separate study, depletion of Tregs resulted in intratumoral HEV formation and higher numbers of activated CD4 and CD8 T cells in carcinogen-induced fibrosarcomas, resulting in reduced tumor burden (115, 116). TNF-producing T cells that signaled through TNFR were critical for intratumoral HEV neogenesis (116).

Given the context-dependent benefit of TLOs, it is unclear whether induction of TLOs is a viable therapeutic strategy to limit tumor progression. TLO formation might need to be accompanied by a strategy to prevent Treg formation to generate an effective antitumor response. Moreover, the signaling mechanisms to induce TLO formation appear to be dependent on various cytokines and receptors of the lymphotoxin/TNF family as well as the local microenvironment, requiring further research to build our understanding.

#### GROWTH OF CANCER CELLS IN LNs

#### LN Vasculature

In primary tumors, hypoxia is associated with expression of VEGF, which in turn leads to the sprouting of nascent blood vessels. Primary tumor hypoxia has been shown to increase the frequency of LN metastasis (117) by upregulating the integrin α<sup>5</sup> subunit, which is required for 3D cell migration *in vitro* (117). We, and others, have found hypoxic tumor cells in the LN (118, 119). It is unclear whether these hypoxic cancer cells in the LN sinus maintain this status from their state in the primary tumor or if tumor cells become hypoxic on arrival and proliferation in the avascular LN sinus.

Angiogenic and non-angiogenic-dependent metastases have been found in LNs and other metastatic tissues, such as the lung and liver (120–123). The presence of hypoxic cancer cells correlates with endothelial cell proliferation in some LN metastases, a pattern that could be predicted by the characteristics of the primary tumor (124). We found elevated VEGF and angiogenesis in primary tumors but did not find the same in metastatic LNs, despite the presence of hypoxic cancer cells (118). In agreement with the findings of the lack of angiogenesis in LN metastasis from our laboratory, other studies have shown that the vascular density of metastatic LNs is lower than that of non-metastatic nodes (121, 125, 126).

Although overexpression of VEGF leads to the expansion of the LN lymphatic vessel network (43–45, 127), a limited number of studies suggest VEGF, or other growth factors, have an effect on the number of blood vessels within the LN. Overexpression of VEGF has been reported to increase the number of HEVs (45). By contrast, other reports demonstrate that VEGF only increases the diameter of LN blood vessels, perhaps through proliferation of endothelial cells of existing vessels (43, 127). The scarcity of evidence concerning VEGF and inflammationinduced sprouting angiogenesis relative to lymphangiogenesis raises questions about the mechanistic control of the LN vasculature. During the progression of LN metastases, the existing LN vasculature may be sufficient to support tumor growth. It has been proposed that angiogenesis is redundant for tumor growth in the LN due to the rich native vascularity of LNs; the vessel density of the LN is comparable to that of the primary tumor (125). It has also been proposed that remodeled HEVs in TDLNs can nurture established metastatic lesions in LNs (128). Although studies of VEGF in other metastatic organs require investigation, the unresponsiveness of LN metastases to antiangiogenic therapy (59, 118) adds another explanation for the poor outcomes of anti-VEGF therapy in adjuvant settings.

Recent studies have investigated mechanisms of resistance to antiangiogenesis therapy in metastatic disease. The tyrosine kinase inhibitor sunitinib, whose targets include VEGF receptors, stimulated transcription and mRNA stabilization of *VEGF-C* in a xenograft model of renal cell carcinoma (129). As a result, lymphangiogenesis and lymphatic metastasis were increased. A recent study identified vascular cooption in breast cancer liver metastases as another resistance mechanism to anti-VEGF therapy (123). In a "replacement" pattern of growth, liver metastases replaced hepatocytes and were physically associated with liver sinusoidal blood vessels. Silencing of Actin Related Proteins 2/3 (ARP 2/3), which mediate breast cancer cell motility, effectively inhibited vascular cooption of these liver metastases. It remains to be determined if this mechanism of vascular cooption is active in LN metastases. As multiple modes of growth were seen in liver metastases, metastatic growth and vascularization may depend on multiple factors, including the tumor type. Further investigation is needed to tailor treatments targeting the growth of metastases, including LN metastases.

#### CLINICAL MANAGEMENT AND TREATMENT OF LN METASTASIS

Lymph node metastasis is a critical prognostic indicator for patients with solid tumors including melanoma, breast, and gastric cancers. However, the role LN metastasis plays in cancer progression has been debated in the clinic for decades (130). The guidelines and treatment strategies for patients with nodal disease are evolving as clinical trials are conducted to improve the ability to contain and treat metastases. Until recently, axillary LN dissection (ALND) had been the standard of care for breast cancer patients with sentinel-node involvement. However, recent clinical trials designed to define the benefit of ALND has changed the way breast cancer patients are being treated. The International Breast Cancer Study Group Trial 23-01 (IBCSG 23-01) clinical trial was conducted to determine the benefit of ALND in patients with limited sentinel-node involvement (1–2 micrometastatic nodes) and tumors less than 5 cm (131). Five-year follow-up showed no disease-free survival benefit in patients that underwent ALND compared with those that did not. These findings are consistent with the ACOSOG Z0011 trial involving patients with limited sentinel-node involvement undergoing breast-conserving surgery, in which these patients were assigned randomly to receive either ALND or no further axillary surgery (132, 133). Both trials suggest that ALND in early-stage breast cancer patients with limited nodal involvement do not have a survival benefit compared with patients that do not undergo ALND. These trials have led to a reduction in the number of breast cancer patients undergoing ALND, which has also reduced the morbidities associated with ALND. However, all of these patients received traditional systemic adjuvant therapy and radiation that potentially eliminated any residual disease in the LNs (132). Thus, radiation and systemic therapies may be sufficient to control nodal disease, making ALND unnecessary for breast cancer patients with limited LN involvement.

To test this hypothesis, recent clinical trials have assessed the benefit of axillary radiation therapy (ART) in early-stage breast cancer patients as an alternate treatment strategy to ALND. Although ALND provides excellent regional control of the disease, patients experience debilitating side effects such as lymphedema. The After Mapping of the Axilla: Radiotherapy or Surgery (AMAROS) randomized phase III trial compared sentinel-node-positive T1-2 breast cancer patients who received either ART or ALND as adjuvant treatment after systemic therapy (134). The results from this trial showed that ART provided excellent local control of disease and the outcomes were comparable to ALND. Furthermore, ART patients had fewer complications compared with those receiving ALND.

For melanoma patients, the Multicenter Selective Lymphadenectomy Trial-I (MSLT-I) provided evidence that patients who undergo sentinel LN (SLN) biopsy have an increased survival rate (135). In this trial, patients were randomized to either (i) wide excision of the melanoma with SLN procedure, followed by complete nodal dissection in patients with a positive SLN or (ii) wide excision and nodal observation, with lymphadenectomy at the time of LN recurrence. However, the importance of complete LN dissection remained controversial, as the main difference between the treatments for patients with disease in LNs was the timing of when lymphadenectomy occurred, with patients having disease removed earlier (SLN patients) having better outcomes. More recently, results of the MSLT-II, a randomized, multicenter, phase-3 clinical trial conducted on 1,934 melanoma patients were reported (136). The MSLT-II trial evaluated the importance of complete LN dissection by randomizing patients with sentinel-node metastases to either immediate complete LN dissection or nodal observation with ultrasonography. Results from this trial showed that immediate complete LN dissection in patients with sentinel-node metastases was not associated with increased 3-year melanoma-specific survival. However, patients that underwent complete LN dissection had increased rates of local-disease control compared with the observation group at 3 years (92 + 1.0 vs. 77 + 1.5%) and a slightly higher rate of disease-free survival compared with the observation group at 3 years (68 + 1.7 and 63 + 1.7%, respectively). The results of the MSLT-I and MSLT-II trials show that removing positive SLNs improves outcomes, but that further LN removal does not improve 3-year overall survival. However, complete nodal dissection in positive SLN patients did improve disease recurrence. The importance of recurrent disease in cancer progression and ultimate patient survival is thus being called into question. It will be critical to address whether there is a difference in 5- and 10-year overall survival in the MSLT-II trial, as these data may be able to account for systemic progression from the recurrent disease. The risk of systemic progression from recurrent disease must be weighed against the very clear reduction in the rate of lymphedema in patients that did not undergo complete LN dissections.

Taken together, the results from these clinical trials have revolutionized the way clinicians manage cancer treatment. Current clinical practice has adopted the theory that less extensive LN dissections for patients with early-stage disease reduce complications without changing overall survival. In these cases, systemic adjuvant therapy and radiotherapy are able to manage any residual disease. Treating local disease in the LN that has not yet progressed to other locations with systemic therapies and radiation could decrease the long-term mortality rate in patients (137), potentially mitigating the impact of less extensive surgeries that leave some LNs with metastatic cancer cells in patients. However, fundamental questions in the biology of LN metastases still remain unanswered, including: What is the fate of cancer cells once they metastasize to the LN? Do LN metastases spread beyond the node and contribute to disease progression? Answers to these questions through continued research will give us better insights into the most effective strategy to manage the progression of solid tumors.

Although there is no direct experimental evidence to show that cancer cells can escape the LN, a few genetic studies provide evidence that nodal disease can colonize distant organs. A recent study used phylogenetic tracing to analyze tumor cell evolution in CRC patients. The analyses from this study revealed that 35% of distant liver metastasis has LN metastasis as the closest phylogenetic neighbor (138). An earlier study characterized the somatic evolution of mutations in cancer cells from primary and metastatic tumors by genome sequencing in a genetically modified mouse model of small cell lung carcinoma (SCLC) (139). This analysis suggested that multiple related primary tumor subclones can seed the LN and a single clone can then spread further from the node to the liver. From the experiments in the SCLC model, the authors were unable to conclude if LN metastases preceded systemic dissemination and if the LN microenvironment altered the genetic and epigenetic makeup of cancer cells before distant dissemination.

Even with radiation and systemic therapies, some patients still have LN recurrences. Systemic therapies that minimize toxicity while targeting disease in LNs represent a promising approach. Only a small fraction of systemically delivered chemotherapy drug accumulates in LNs (2). We have found that increased vascular permeability of LN blood vessels did not increase chemotherapy penetration into the LN parenchyma (140). Many current strategies focus on optimizing lymphatic delivery and retention properties of therapeutics that target tumor-associated lymphatic vessels (141). Targeted delivery of immunomodulatory agents or cancer cell-specific cytotoxic drugs into TDLNs can improve cancer vaccination strategies (142) and eradicate disease from LNs, respectively (141).

### CONCLUSION

A growing interest and ability to study LEC biology has provided enormous insights into the role lymphatic vessels play in tumor metastasis and cancer progression. While it is known that tumorassociated lymphatic vessels are a route of metastasis, it is now appreciated that tumor cells also use these vessels to establish immunological tolerance in TDLNs. Although antitumor immune responses can be generated locally in primary tumors, TDLNs serve as critical sites for antitumor immune responses. It remains unknown how the systemic adaptive immune response to cancer is shaped by immunosuppressed LNs. Several additional important areas of exploration remain, including understanding the influence of the LN microenvironment on cancer cell behavior and determining the contribution of LN metastases to distant organ metastasis, which at present remains controversial. Therapies targeting LN metastases must also consider enhancing antigen presentation to tumor-specific T cells. Moreover, therapies to activate tumor-specific T cells should be considered in parallel with strategies to break tolerance and other immunosuppression mechanisms. With continued research focus on the LN, we will gain deeper insights into mechanisms of immune evasion

### REFERENCES


by cancer cells. A more thorough understanding of the molecular signals that facilitate tumor metastasis to TDLNs and beyond may also provide therapeutic targets to control the further dissemination of lymphatic metastases. With these continued advances, patient survival from metastatic cancer will continue to improve.

### AUTHOR CONTRIBUTIONS

DJ, EP, and TP contributed to the writing of the manuscript.

#### ACKNOWLEDGMENTS

This work was supported by the NIH DP2OD008780 (TP), R01CA214913 (TP), R01HL128168 (TP), NCI Federal Share/ MGH Proton Beam Income on C06 CA059267 (TP), and the Massachusetts General Hospital Executive Committee on Research ISF (TP), and Fund for Medical Discovery Postdoctoral Fellowship Award (EP), and Burroughs Wellcome Postdoctoral Enrichment Program Award (DJ). The authors would like to acknowledge Dr. Lance Munn and Sonia Pereira for their assistance with the graphical illustration.

egress and lymphatic patterning. *J Exp Med* (2010) 207:17–27. doi:10.1084/ jem.20091619


142. Maisel K, Sasso MS, Potin L, Swartz MA. Exploiting lymphatic vessels for immunomodulation: rationale, opportunities, and challenges. *Adv Drug Deliv Rev* (2017) 114:43–59. doi:10.1016/j.addr.2017.07.005

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Jones, Pereira and Padera. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

#### *Shubham Tripathi 1,2, Mohit Kumar Jolly <sup>2</sup> , Wendy A. Woodward3,4, Herbert Levine1,2,5,6 and Michael W. Deem1,2,5,6\**

*1PhD Program in Systems, Synthetic, and Physical Biology, Rice University, Houston, TX, United States, 2Center for Theoretical Biological Physics, Rice University, Houston, TX, United States, 3Department of Radiation Oncology, The University of Texas MD Anderson Cancer Center, Houston, TX, United States, 4MD Anderson Morgan Welch Inflammatory Breast Cancer Research Program and Clinic, The University of Texas MD Anderson Cancer Center, Houston, TX, United States, 5Department of Bioengineering, Rice University, Houston, TX, United States, 6Department of Physics and Astronomy, Rice University, Houston, TX, United States*

#### *Edited by:*

*Triantafyllos Stylianopoulos, University of Cyprus, Cyprus*

#### *Reviewed by:*

*Alessandro Giuliani, Istituto Superiore di Sanità, Italy Ovidiu Radulescu, Université de Montpellier, France*

> *\*Correspondence: Michael W. Deem mwdeem@rice.edu*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

> *Received: 02 February 2018 Accepted: 18 June 2018 Published: 04 July 2018*

#### *Citation:*

*Tripathi S, Jolly MK, Woodward WA, Levine H and Deem MW (2018) Analysis of Hierarchical Organization in Gene Expression Networks Reveals Underlying Principles of Collective Tumor Cell Dissemination and Metastatic Aggressiveness of Inflammatory Breast Cancer. Front. Oncol. 8:244. doi: 10.3389/fonc.2018.00244*

Clusters of circulating tumor cells (CTCs), despite being rare, may account for more than 90% of metastases. Cells in these clusters do not undergo a complete epithelial-tomesenchymal transition (EMT), but retain some epithelial traits as compared to individually disseminating tumor cells. Determinants of single cell dissemination versus collective dissemination remain elusive. Inflammatory breast cancer (IBC), a highly aggressive breast cancer subtype that chiefly metastasizes *via* CTC clusters, is a promising model for studying mechanisms of collective tumor cell dissemination. Previous studies, motivated by a theory that suggests physical systems with hierarchical organization tend to be more adaptable, have found that the expression of metastasis-associated genes is more hierarchically organized in cases of successful metastases. Here, we used the cophenetic correlation coefficient (CCC) to quantify the hierarchical organization in the expression of two distinct gene sets, collective dissemination-associated genes and IBC-associated genes, in cancer cell lines and in tumor samples from breast cancer patients. Hypothesizing that a higher CCC for collective dissemination-associated genes and for IBC-associated genes would be associated with retention of epithelial traits enabling collective dissemination and with worse disease progression in breast cancer patients, we evaluated the correlation of CCC with different phenotypic groups. The CCC of both the abovementioned gene sets, the collective dissemination-associated genes and the IBC-associated genes, was higher in (a) epithelial cell lines as compared to mesenchymal cell lines and (b) tumor samples from IBC patients as compared to samples from non-IBC breast cancer patients. A higher CCC of both gene sets was also correlated with a higher rate of metastatic relapse in breast cancer patients. In contrast, neither the levels of *CDH1* gene expression nor gene set enrichment analysis (GSEA) of the abovementioned gene sets could provide similar insights. These results suggest

**47**

that retention of some epithelial traits in disseminating tumor cells as IBC progresses promotes successful breast cancer metastasis. The CCC provides additional information regarding the organizational complexity of gene expression in comparison to GSEA. We have shown that the CCC may be a useful metric for investigating the collective dissemination phenotype and a prognostic factor for IBC.

Keywords: collective dissemination, inflammatory breast cancer, epithelial-to-mesenchymal transition, hierarchy, hybrid E/M, cophenetic correlation coefficient

#### INTRODUCTION

Metastasis is responsible for 90% of deaths from solid tumors (1). It involves the escape of cancer cells from the site of the primary tumor, their entry into the circulatory system, and finally, colonization of and proliferation at a distant organ. However, this process is highly inefficient. Only an estimated 0.2% of the disseminated tumor cells are able to form a lesion at distant organ sites (2, 3). A well-studied mechanism of metastasis is single cell dissemination where carcinoma cells acquire migratory and invasive traits *via* an epithelial-to-mesenchymal transition (EMT) (4). These cells can then utilize blood or lymph circulation to reach distant organ sites, where they reacquire epithelial traits of cell–cell adhesion and apico-basal polarity *via* a mesenchymalto-epithelial transition (MET) to establish metastases (4).

Recent studies have highlighted that EMT is not a binary process. Rather, cells *en route* to a mesenchymal phenotype can acquire a stable hybrid epithelial–mesenchymal (hybrid E/M) phenotype (5, 6). These observations have called into question the indispensability of a complete EMT followed by MET in metastasis (7). Instead, collective migration of tumor cells *via* clusters of circulating tumor cells (CTCs) has been suggested as an alternate mechanism of metastasis (8). Clusters of tumor cells had been detected in the bloodstream of cancer patients even before the characterization of EMT as a driver of cancer metastasis (9, 10). These clusters of tumor cells can efficiently seed secondary tumors, exhibiting up to 50 times the metastatic potential of individually migrating tumor cells (11). Tumor cell clusters accounted for >90% of metastases in a mouse model of breast cancer (12). Abundance of CTC clusters in the bloodstream has been associated with significantly poor prognosis in breast cancer and in small cell lung cancer (SCLC) (11, 13).

Multiple factors are believed to be responsible for the heightened metastatic potential of these CTC clusters. These include effective response to mechanical signals and chemical gradients by cells in CTC clusters as compared to migrating single tumor cells (14, 15), better evasion of the host immune system (16), and potential cooperation among heterogeneous cell types in CTC clusters (17, 18). Studies have shown that collectively invading tumor cells from the primary lesion often co-express epithelial and mesenchymal markers (19–21). Thus, cells in CTC clusters tend to manifest a hybrid epithelial–mesenchymal (hybrid E/M) phenotype and to retain cell–cell adhesion characteristics (8).

Inflammatory breast cancer (IBC) is a highly aggressive breast cancer subtype that has been reported to predominantly metastasize *via* CTC clusters (22). Characterized by breast erythema, edema, and *peau d'orange* presenting with or without a noticeable tumoral mass (23, 24), IBC involves tumoral infiltrate in the dermal lymphatics and about 30% of IBC patients have distant metastases at the time of diagnosis as compared to only 5% of non-IBC type breast cancer patients (25). Though only 2–4% of breast cancer cases each year in the United States are of the IBC type, IBC patients account for 10% of the annual breast cancer-related mortalities. A hallmark of IBC is the presence of cohesive clusters of tumor cells in the local lymph nodes (26) and IBC patients have larger and a higher frequency of CTC clusters as compared to non-IBC breast cancer patients (27). Abundance of CTC clusters has been shown to be associated with poor progression-free survival in IBC patients (27). Despite their great propensity to metastasize, tumor cells in the primary lesion and in metastatic lesions of IBC maintain a high expression E-cadherin, a hallmark of epithelial cells (26). IBC thus presents an exciting model for the study of collective dissemination of tumor cells *via* CTC clusters and of the prognostic potential of these clusters of migrating tumor cells. The results presented here strengthen the argument for investigating IBC to elucidate the mechanisms underlying collective dissemination of tumor cells.

Here, we invoke concepts from theoretical models of evolution to investigate cluster-based dissemination of tumor cells and analogous IBC characteristics. Theoretical studies suggest that systems with a more hierarchical structure are more adaptable (28–30) due to their ability to efficiently span the space of possible states. Hierarchical systems are also more robust to perturbations because a hierarchical network architecture has a buffering effect that hinders the propagation of local perturbations to a majority of nodes (30, 31). Hierarchical organization, thus, emerges over time in physical systems that are evolving in a changing environment with a rugged fitness landscape exhibiting numerous peaks and valleys (29). Given that tumor cells involved in metastasis and invasion progress through many different microenvironments (32–34), one can expect the expression of genes associated with a metastatic phenotype to be more hierarchically organized in instances of successful macrometastases as compared to instances with no metastasis.

We quantified the hierarchical organization in the expression of two distinct sets of genes, one associated with collective dissemination of tumor cells and the other related to IBC, in cancer cell lines and in breast cancer patients. For this purpose, we used the cophenetic correlation coefficient (CCC) metric. The CCC for a set of genes takes into consideration the collective expression of all genes within the given set and the correlations between the expression levels of different genes. It captures the level of hierarchical organization in the collective expression of genes in the given set. A higher CCC indicates greater hierarchical organization in the expression of genes. The CCC was first used for comparing tree-like relationships represented by different dendrograms (35). It has been used previously to quantify the differences in expression of metastasis-associated genes in breast cancer patients with different clinical outcomes (36) and to quantify the differences in expression of genes predictive of clinical outcome in adult acute myeloid leukemia in patients belonging to different risk categories (37).

The goal of the present study was to determine whether the hierarchical organization in the expression of two sets of genes of interest is different in cell lines exhibiting different EMTassociated phenotypes and in tumor samples from breast cancer patients exhibiting features of IBC and non-IBC type disease. The first set of genes investigated here includes 87 genes reported to be associated with collective dissemination of tumor cells as CTC clusters: genes differentially expressed in cells forming CTC clusters as compared to individual CTCs (12). The second gene set includes 78 genes reported to be differentially expressed in tumor samples from IBC patients in comparison to tumor samples from non-IBC breast cancer patients (38). We observed that the CCC for both of these gene sets was higher in (a) epithelial cell lines as compared to mesenchymal cell lines and (b) tumor samples from IBC patients as compared to tumor samples from non-IBC breast cancer patients. A higher CCC further correlated with worse disease progression in breast cancer patients. In light of these observations, we propose that the metastatic aggressiveness of IBC potentially derives from the hierarchical organization in the expression of collective dissemination-associated genes in metastasizing tumor cells.

### MATERIALS AND METHODS

#### Genes Associated With Collective Dissemination of Tumor Cells

Using multicolor lineage tracking, Cheung et al. showed that polyclonal seeding by disseminated clusters of tumor cells is the dominant mechanism for metastasis in a mouse model of breast cancer (12). These clusters accounted for more than 90% of distant organ metastases in mice. Circulating tumor cell clusters were observed to be enriched in expression of the epithelial protein keratin 14 (K14), and 87 genes with enriched or depleted expression in K14<sup>+</sup> primary tumor cells as compared to K14<sup>−</sup> primary tumor cells were identified. Broadly, expression of adhesion complex-associated genes was enriched and that of MHC Class II genes was depleted in K14<sup>+</sup> cells. We used this set of genes as a signature of the collective dissemination phenotype.

### Genes Associated With the IBC Phenotype

Van Laere et al. obtained tumor samples from patients with breast adenocarcinoma: 137 samples from IBC patients and 252 samples from patients with non-IBC type breast cancer (non-IBC) (38). IBC patients were selected in accordance with the consensus diagnostic criteria described by Dawood et al. (23). RNA from the tumor samples was hybridized onto Affymetrix GeneChips (HGU133-series) to obtain the corresponding mRNA expression profiles. Linear regression models were employed to identify a set of 78 IBC specific genes, which were differentially expressed in IBC tumor samples as compared to non-IBC tumor samples, independent of the molecular subtype of the tumor (38). We used this set of genes as a signature of the IBC phenotype in breast cancer patients. There were no genes common between this set of IBC-associated genes and the set of collective disseminationassociated genes described above. Both gene sets are available as Supplementary Material. The statistical methods used previously to obtain these gene sets are summarized in the Supplementary Material.

#### Gene Expression Data From Different Cell Lines

We used two different datasets of gene expression in cell lines, each cell line classified as epithelial (E), mesenchymal (M), or hybrid epithelial–mesenchymal (hybrid E/M). The first dataset was from the study by Grosse-Wilde et al. (39), Gene Expression Omnibus (GEO) accession number GSE66527. A total of 24 clones established from HMLER cell lines [normal human mammary epithelial cells immortalized and transformed with hTERT and the oncogenes *SV40LT* and *RAS* (40)] were sorted into 13 *CD24*<sup>+</sup>/*CD44*<sup>−</sup> E clones and 11 *CD24*<sup>−</sup>/*CD44*<sup>+</sup> M clones. The E clones and the M clones displayed cobble-stone like morphology and dispersed, fibroblast morphology, respectively.

The second dataset included gene expression from the National Cancer Institute 60 anticancer drug screen (NCI60), which includes panels of cell lines representing nine distinct types of cancer: leukemia, colon, lung, central nervous system, renal, melanoma, ovarian, breast, and prostate (41). The 60 cell lines have been classified into epithelial (E) (*n* = 11), mesenchymal (M) (*n* = 36), and hybrid epithelial-mesenchymal (hybrid E/M) (*n* = 11) categories on the basis of protein levels of E-cadherin and Vimentin (42). The gene expression data for these cell lines obtained using the Affymetrix Human Genome U133A array platform were downloaded from the CellMiner database (43, 44).

#### Gene Expression Data From Tumor Samples From IBC and Non-IBC Breast Cancer Patients

We used three different datasets of gene expression in tumor samples obtained from breast cancer patients. Each patient in the three datasets was diagnosed with either IBC or non-IBC type breast cancer (non-IBC). Iwamoto et al., GEO accession number GSE22597, collected tumor biopsies prospectively from 82 patients with locally advanced disease. A clinical diagnosis of IBC was made in 25 of these patients (45). Boersma et al., GEO accession number GSE5847, examined primary breast tumor samples from 50 patients, 15 of whom were diagnosed with IBC on the basis of the pathology and medical reports (46). Finally, Woodward et al., GEO accession number GSE45584, obtained tissue samples from core biopsies of breast tissue in 40 breast cancer patients, 20 IBC and 20 non-IBC (24).

In Iwamoto et al. and Woodward et al., IBC diagnosis was made in patients with clinical presentation of breast erythema and edema over more than one-third of the breast. In Boersma et al., nine IBC patients presented with erythema and edema, while six IBC patients exhibited pathology indicating dermal lymphatic invasion and tumor emboli.

The microarray platforms and the normalization techniques used previously to obtain the gene expression profiles for different cell lines and for tumor samples from cancer patients have been outlined in the Supplementary Material.

#### Definition of Gene Network for Different Phenotypic Groups

For each phenotypic group, e.g., NCI60 cell lines labeled as epithelial (E) or patients in the Iwamoto et al. (45) dataset diagnosed with IBC, and gene set, e.g., the set of genes associated with IBC or the set of collective dissemination-associated genes, we defined a network with the genes as nodes and weighted edges between these nodes. The weight of the edge between gene *i* and gene *j* in the phenotypic group *G* was defined as

$$I\_{ij}^G = \left| \sum\_{k \in G} \frac{(e\_i^k - \mu\_i^G)(e\_j^k - \mu\_j^G)}{\sigma\_i^G \sigma\_j^G} \right| \tag{1}$$

Here, *em k* is the expression of gene *m* in the sample *k* (patient/ cell line), µ*<sup>m</sup> <sup>G</sup>* and σ*<sup>m</sup> <sup>G</sup>* are the mean and SD of the expression of gene *m* in the phenotypic group *G*, respectively, and the summation is over all patients or cell lines belonging to the group *G*. This definition resulted in a fully connected network for each phenotypic group and gene set. Since Eq. 1 is symmetric in *i* and *j*, the networks obtained were undirected.

We constructed such networks for the epithelial and mesenchymal cell lines in the Grosse-Wilde et al. (39) dataset and for the epithelial, mesenchymal, and hybrid epithelial–mesenchymal cell lines in the NCI60 dataset. Such networks were also constructed for IBC and non-IBC breast cancer patients in the three breast cancer datasets, Iwamoto et al. (45), Woodward et al. (24), and Boersma et al. (46), using each of the two gene sets described above, genes associated with collective dissemination of tumor cells and genes associated with the IBC phenotype.

#### Calculation of the CCC

To quantify the hierarchy in the expression of a set of genes in different groups of patients and cell lines, we used a metric called the CCC (35). The CCC is a measure of how well a hierarchical clustering of nodes in a network reproduces the distances between nodes in the original network. Intuitively, the CCC is a measure of how tree-like a network is. Since a tree topology is a prototypical hierarchical structure, a measure of the tree-like characteristic of a network allows us to aptly quantify the underlying hierarchy in the structure of the network.

For calculating the CCC of a given network, we defined the distance between nodes *i* and *j* as the Euclidean commute time distance (ECTD) between the two nodes, which is given by the square root of the mean first passage time taken by a random walker to travel from node *i* to node *j* and then back to node *i*. The ECTD between nodes *i* and *j* depends not only on the weight of the edge between nodes *i* and *j* but also on the number of different possible paths between the two nodes. The ECTD decreases as the number of possible paths between the two nodes increases, and increases if any path between the two nodes becomes longer (47). This makes the ECTD suitable for clustering tasks. As described before, the network obtained for each phenotypic group and gene set was undirected and fully connected. This ensures that the ECTD between any pair of nodes will be finite. For a network with *N* nodes, we generated a *N* × *N* matrix *D* such that *Dij* is the ECTD between nodes *i* and *j* (48). The matrix *D* was then used as an input to the average linkage hierarchical clustering algorithm (49), which generates a tree topology (*T*), i.e., a dendrogram, which best approximates the distances between the nodes of the network given by the matrix *D*. We then calculated the CCC as the correlation between the original pairwise distances and the corresponding distances in the tree topology:

$$\text{CCCC} = \frac{\sum\_{i$$

Here, *d Dij* = is the mean of the original pairwise distances and *t Tij* = is the mean of the pairwise distance in the tree topology. If the original network is hierarchical, the distances between nodes in the tree topology obtained *via* hierarchical clustering (*T*) will be highly correlated with the distances between nodes in the original network (*D*). Hence, the CCC will be high. However, if the original network lacks any hierarchical organization, this correlation will be weak, and the CCC will be low.

To test the sensitivity of the calculated CCC to the choice of ECTD as the network distance metric for hierarchical clustering, we alternatively defined the distance between node *i* and node *j* in the network as the resistance distance (50) between the two nodes. The resistance distance between any two nodes is given by the effective electrical resistance when a battery is connected across the two nodes. Like the ECTD, the resistance distance depends on all possible paths between nodes *i* and *j* and is, therefore, suited for clustering tasks. Using the resistance distance to create the matrix *D* where *Dij* is the resistance distance between the nodes *i* and *j*, we calculated the CCC as described above.

The CCC calculated for a network was normalized with respect to the CCC of random networks with the same set of nodes but re-distributed edge weights. For this, we generated 10 such random networks by shuffling entries in the matrix *D* and then calculated the average of the CCCs of these random networks (CCCrand). The normalized CCC was then defined as

$$\text{CCC}\_{\text{norm}} = \frac{\text{CCC} - \text{CCC}\_{\text{rand}}}{1 - \text{CCC}\_{\text{rand}}} \tag{3}$$

Finally, to obtain the error in the estimate of CCCnorm, we used the bootstrap method (51). The method assumes that the distribution of gene expression in a patient or cell line group is the empirical distribution function of the observed expression in samples within the group. For a patient or cell line group with size *n*, we drew *n* samples from the group with replacement and calculated CCCnorm for the sampled group. This sampling process was repeated 100 times to obtain 100 CCCnorm values. The SE in the estimate of the CCCnorm for the group was then given as the sample SD of the 100 sampled CCCnorm values. All *p*-values were also calculated using the bootstrap method.

The MATLAB code used for calculating the CCC is available at https://github.com/st35/gene-network-CCC.

#### RESULTS

#### Higher CCC for the Collective Dissemination-Associated Gene Network in Epithelial Cell Lines and in IBC Patients

We constructed networks with genes associated with collective dissemination of tumor cells (12), hereafter referred to as "collective dissemination-associated" genes, as nodes and the weight of the edge between nodes in a pair defined according to Eq. 1. Such networks were constructed for the E and M cell lines from the gene expression data from Grosse-Wilde et al. (39) and for the cell lines in the NCI60 anticancer drug screen (41) that have been categorized into E, M, and hybrid E/M classes (42). Representative networks for E and M cell lines from Grosse-Wilde et al. (39) are shown in Figures S1A,B in Supplementary Material. The normalized CCC for these networks was calculated using the method described above, and the results are shown in **Figure 1**. E cell lines exhibited a significantly higher CCC as compared to M cell lines (*p*-value = 0.01) for the collective disseminationassociated gene network in the dataset from Grosse-Wilde et al. (39), **Figure 1A**. In the NCI60 dataset, the normalized CCC of the collective dissemination-associated gene network was higher for E cell lines as compared to the pooled M and hybrid E/M cell lines, **Figure 1B**. The bootstrap distribution of normalized CCC values for E cell lines was distinct from the distribution for M cell lines in the dataset from Grosse-Wilde et al. (39) and from the distribution for pooled M and hybrid E/M cell lines in the NCI60 dataset (Kolmogorov–Smirnov test, *p*-value < 0.01). We further calculated the CCC using the resistance distance instead of the ECTD and observed a similar trend in CCC values in the two datasets, Figures S2A,C in Supplementary Material.

We constructed similar networks for IBC and non-IBC patients using Affymetrix U133A profiles obtained by Iwamoto et al. (45). Representative networks for IBC and non-IBC breast cancer patients are shown in Figures S1C,D in Supplementary Material. Normalized CCC values for patients in the two groups are shown in **Figure 2A**. IBC patients exhibited a higher CCC for the network associated with collective dissemination of tumor cell clusters as compared to non-IBC breast cancer patients. The difference between the two groups in the dataset was significant (*p*-value < 0.02). Further, bootstrap distributions of the normalized CCC values for the two groups were

Figure 1 | Normalized cophenetic correlation coefficient of the collective dissemination-associated gene network for (A) 13 epithelial (E) and 11 mesenchymal cell lines from the study by Grosse-Wilde et al. (39), and (B) 11 epithelial (E) and 47 hybrid epithelial-mesenchymal (E/M) + mesenchymal (M) cell lines from the NCI60 dataset (41, 42). Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \**p*-Value < 0.05.

statistically distinct with *p*-value < 0.01 for the Kolmogorov– Smirnov test. The same trend in CCC values was observed from calculation of CCC using the resistance distance, Figure S2E in Supplementary Material. However, we did not observe a significant trend for the breast cancer samples characterized by Boersma et al. (46) and for the samples characterized by Woodward et al. (24), **Figures 2B,C**.

#### Higher CCC for the IBC-Associated Gene Network in Epithelial Cell Lines and in IBC Patients

We constructed networks with genes differentially expressed in tumor samples from IBC patients as compared to tumor samples from non-IBC breast cancer patients, hereafter referred to as "IBC-associated" genes, as nodes. The weight of the edge between nodes in a pair was defined using Eq. 1. Such networks were constructed for the E and M cell lines in the dataset from Gross-Wilde et al. (39) and for the E and pooled M + hybrid E/M cell lines in the NCI60 dataset. Normalized CCC values for these groups of cell lines calculated using the method described above are shown in **Figure 3**. E cell lines displayed a higher CCC for the IBC-associated gene network as compared to the other cell lines in both datasets [*p*-value = 0.03 for the cell lines in the study by Grosse-Wilde et al. (39) and *p*-value = 0.02 for the cell lines in the NCI60 dataset]. The bootstrap distributions of normalized CCC values for the two groups of cell lines were statistically distinct for both datasets (*p*-value < 0.01 for the Kolmogorov–Smirnov test in each case). A higher CCC for the epithelial cell lines was also observed on using the resistance distance to calculate the CCC, Figures S2B,D in Supplementary Material.

Using Affymetrix U133A profiles from Iwamoto et al. (45), we constructed similar networks with IBC-associated genes as nodes for both IBC and non-IBC breast cancer patients. Normalized CCC values for the two breast cancer patient groups are shown in **Figure 4A**. The IBC group exhibited a significantly higher CCC for the IBC-associated gene network as compared to the non-IBC patient group (*p*-value = 0.01). Bootstrap distributions for the two groups were again statistically distinct (*p*-value < 0.01 for the Kolmogorov–Smirnov test). This trend in CCC values was also observed on using the resistance distance to calculate the CCC, Figure S2F in Supplementary Material. A similar trend in the CCC values for IBC and non-IBC patient groups was observed for breast cancer patients in the two other independent breast cancer datasets, Boersma et al. (46) (*p*-value = 0.02) and Woodward et al. (24) (*p*-value = 0.06), **Figures 4B,C**.

Saunders and McClay had used a well-understood gene regulatory network in the sea urchin embryo to identify transcription factors that control cell changes during EMT by perturbing individual transcription factors (52). They further determined 30

of CCCnorm calculated using the bootstrap method. \**p*-Value < 0.05.

human transcription factors homologous to those identified in sea urchins. We calculated the CCC of a network with these transcription factors, hereafter referred to as "canonical drivers of EMT," as nodes for the IBC and non-IBC samples from each of the three breast cancer datasets, Iwamoto et al. (45), Boersma et al. (46), and Woodward et al. (24). The weight of the edge between any two transcription factors was defined using Eq. 1. We observed that the IBC patient group exhibited a lower CCC for the network composed of canonical EMT drivers as compared to the non-IBC patient group in data from each of the three studies, **Figure 5**.

#### Higher CCC for the Two Networks Correlates With Early Metastasis Posttreatment

Having analyzed the differences in CCCs of collective dissemination-associated and IBC-associated gene sets in epithelial and mesenchymal cell lines and in tumor samples from IBC and non-IBC breast cancer patients, we next investigated if the CCC of these gene sets could provide insights into the metastatic propensity of tumors. We constructed networks with the two sets of genes, collective dissemination-associated and IBC-associated, as nodes for breast cancer patients who exhibited metastatic relapse within 5 years posttreatment (53). These patients were classified into two groups, those with metastatic relapse within 30 months posttreatment and those with metastasis between 30 and 60 months posttreatment, **Figures 6A,B**. Edge weights were defined, once again, using Eq. 1. For both collective dissemination-associated and IBC-associated gene sets, the CCC was significantly higher for the patient group with early metastatic relapse of breast cancer, i.e., relapse within 30 months of treatment, as compared to the patient group with relatively late relapse, i.e., metastatic relapse after 30 months posttreatment, **Figures 6A,B**. The *p*-values were 0.02 and 0.01 for the collective dissemination-associated gene network and the IBC-associated gene network, respectively. The same trend was observed upon considering only estrogen-receptor-positive patients, Figure S3 in Supplementary Material. There were too few samples from estrogen-receptor-negative patients for a similar analysis.

To investigate if the observation that a more hierarchical expression of collective dissemination-associated genes correlates with early relapse posttreatment can be generalized to other cancer types, we calculated the CCC of collective disseminationassociated and IBC-associated genes for SCLC patients. SCLC is a highly aggressive cancer subtype that is known to form tumor emboli and metastasize quickly, predominately *via* clusters of CTCs (55–57). SCLC patients with fewer than 10 months of disease-free survival posttreatment exhibited a higher CCC for both collective dissemination-associated and IBC-associated gene sets as compared to patients with greater than 10 months of disease-free survival posttreatment as computed from the data in the study by Rousseaux et al. (54), **Figures 6C,D**.

The metastasis of cancer to different organs is characterized by organ-specific bottlenecks (58). While tumor cells from the site of the primary lesion can easily migrate to the local lymph nodes by moving passively with the lymph flow, migration to other organs such as skin or liver is much more challenging. Given the benefits afforded to migrating cancer cells by collective dissemination, cells with a more hierarchical expression of collective dissemination-associated genes are likely to be over-represented in cancer metastases to distant organs as compared to metastases to local lymph nodes. Using the gene expression data from the study by Kimbung et al. (59), we calculated the CCC of collective disseminated-associated genes in samples from breast cancer metastases to different organs and observed a higher CCC for metastases to skin as compared to metastases to lymph nodes and liver, **Figure 7A**. A similar trend was observed on calculating the CCC of the IBC-associated gene network, **Figure 7B**.

We further explored whether the CCCs for the collective dissemination-associated gene network and the IBC-associated gene network were different in breast cancer patients with metastatic relapse within 5 years posttreatment and those with no metastasis during this follow-up period as computed from the data in the study by Wang et al. (53). Intriguingly, we observed that the CCCs of both networks were significantly higher, *p*-value = 0.03 in each case, for patients with no metastasis during the 5-year follow-up period as compared to the patients with metastatic relapse during the follow-up, **Figures 8A,B**. A similar trend was observed for breast tumor samples from The Cancer Genome Atlas (TCGA) for patients who exhibited relapse during the follow-up period and those who did not (60), **Figures 8C,D**. Given that healthy breast cells are inherently epithelial, a higher CCC for the patient

Figure 6 | (A) Normalized cophenetic correlation coefficient of the collective dissemination-associated gene network for breast cancer patients with metastatic relapse within a 30-month period posttreatment (*T* < 30; *n* = 56) or between 30 and 60 months posttreatment (*T* ≥ 30; *n* = 51). Gene expression data from the study by Wang et al. (53). (B) Normalized CCC of the inflammatory breast cancer-associated gene network for the same groups of breast cancer patients as in (A). (C) Normalized CCC of the collective dissemination-associated gene network for small cell lung cancer (SCLC) patients with less than 10 months of disease-free survival posttreatment (*T* < 10; *n* = 11) and SCLC patients with longer than 10 months of disease-free survival posttreatment but death during the follow-up period (*T* ≥ 10; *n* = 10). Gene expression data from the study by Rousseaux et al. (54). (D) Normalized CCC of the IBC-associated gene network for the same SCLC patient groups as in (C). Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \**p*-Value < 0.05.

group with no metastatic relapse during the follow-up period may be a consequence of the tumor being at initial stages of progression toward a metastatic phenotype at the time of diagnosis and sample collection in these patient groups. However, upon grouping the breast cancer patients by their estrogen-receptor status, no consistent trend was observed between patients with no relapse during the follow-up period and patients with metastatic relapse during the follow-up period for both gene sets, Figure S4 in Supplementary Material. These results indicate that the collective dissemination pathway in breast cancer patients with differing receptor statuses warrants further study.

### The CCC Provides Additional Information Regarding the Underlying Complexity of Collective Gene Expression

We next investigated if the insights described above can be obtained from an analysis of expression levels of collective

Figure 8 | (A) Normalized cophenetic correlation coefficient of the collective dissemination-associated gene network for tumor samples from 107 breast cancer patients who did not exhibit breast cancer relapse during the follow-up period and for tumor samples from 179 patients who exhibited metastatic relapse during the follow-up period. Gene expression data from the study by Wang et al. (53). (B) Normalized CCC of the Inflammatory breast cancer-associated gene network for the same patient groups as in (A). (C) Normalized CCC of the collective dissemination-associated gene network for tumor samples from 527 breast cancer patients with no metastasis during the follow-up period and for tumor samples from 13 patients with breast cancer metastasis during the follow-up period. Gene expression data from the cancer genome atlas (TCGA) (60). (D) Normalized CCC of the IBC-associated gene network for the same patient groups as in (C). Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \**p*-value < 0.05 and \*\**p*-value < 0.01.

dissemination-associated and IBC-associated genes. To determine how the CCCs of different gene networks correlate with the expression levels of these genes in different phenotypic groups, we carried out gene set enrichment analysis (GSEA) for different sets of genes on the epithelial and mesenchymal cell lines from the study by Grosse-Wilde et al. (39) and on the tumor samples from IBC patients and non-IBC breast cancer patients from the study by Iwamoto et al. (45). Using the GSEA software provided by the Broad Institute (61), we tested for enrichment in the expression of collective dissemination-associated genes, IBC-associated genes, and the canonical drivers of EMT in different phenotypic groups, i.e., epithelial versus mesenchymal cell lines in the data from Grosse-Wilde et al. (39) and IBC versus non-IBC patients in the data from Iwamoto et al. (45). The results are shown in **Figures 9A–F**. The expression of collective disseminationassociated genes was significantly enriched in epithelial cell lines as compared to mesenchymal cell lines (*p*-value < 0.001), **Figure 9A**, while IBC-associated genes and canonical EMT drivers did not show any such significant enrichment when compared across these two phenotypic groups. On the other hand, expression of IBC-associated genes was significantly enriched in tumor samples from IBC patients (*p*-value = 0.035), **Figure 9E**, while the collective dissemination-associated genes and canonical EMT drivers did not show significant enrichment on comparing IBC tumor samples with non-IBC breast tumor samples. We further divided the set of collective dissemination-associated genes into

two groups, genes with enriched expression levels in K14+ cells and genes with depleted expression levels in K14+ cells. Neither of these two subsets exhibited significant enrichment when carrying out IBC tumor samples versus non-IBC breast tumor samples GSEA, Figure S5 in Supplementary Material.

Previous studies have suggested a strong association between expression of the E-cadherin protein in tumor cells and IBC (62, 63). To investigate if the level of *CDH1* (E-cadherin) gene expression in tumor samples from breast cancer patients is also associated with IBC, we compared the levels of *CDH1* gene expression in tumor samples from IBC and non-IBC patients. There was no significant difference in the expression levels of *CDH1* gene between the two patient groups in any of the three breast cancer patient datasets, Iwamoto et al. (45), Boersma et al. (46), and Woodward et al. (24), **Figures 9G–I**.

Finally, to test the specificity of the collective disseminationassociated and IBC-associated gene sets in characterizing IBC behavior, we generated 100 random gene sets. Each gene set consisted of 83 genes, average of the sizes of the collective dissemination-associated and IBC-associated gene sets. We calculated the normalized CCC of these gene sets in tumor samples from IBC and non-IBC breast cancer patients from the study by Iwamoto et al. (45). Only for 2 of the 100 randomly generated gene sets, the CCC was significantly higher for the IBC group as compared to the non-IBC group (*p*-value < 0.05), Figure S6 in Supplementary Material. This indicates that our hypothesis of a

Figure 9 | (A–F) Gene set enrichment analysis using collective dissemination-associated genes (A,D), inflammatory breast cancer (IBC)-associated genes (B–E), and canonical epithelial-to-mesenchymal transition driving transcription factors (C–F) on: (A–C) gene expression data for epithelial and mesenchymal cell lines from the study by Grosse-Wilde et al. (39) and on (D–F) gene expression data for tumor samples from IBC and non-IBC breast cancer patients from the study by Iwamoto et al. (45). In (A–C), genes are ordered from left to right in decreasing order of correlation of expression with the epithelial phenotype. In (D–F), genes are ordered from left to right in decreasing order of correlation of expression with the IBC phenotype. Black bars along the top of each plot indicate the positions of hits to the gene set along the ordered list of genes. Nominal *p*-values of enrichment are indicated at the bottom of each plot. (G–I) Mean expression of *CDH1* (E-cadherin) gene in tumor samples from IBC and non-IBC breast cancer patients in studies by (G) Iwamoto et al. (45), (H) Boersma et al. (46), and (I) Woodward et al. (24).

more hierarchically organized gene expression in IBC samples as compared to non-IBC breast cancer samples is specific to collective dissemination-associated and IBC-associated gene sets and is not applicable to randomly chosen sets of genes.

### DISCUSSION

Cancer metastasis *via* migrating clusters of CTCs has emerged as a critical mechanism of seeding secondary tumors in recent studies (9–12). Although rare in comparison with individually disseminated cancer cells, CTC clusters can efficiently seed secondary tumors at distant organ sites (11, 12), and their presence in the bloodstream of cancer patients has been shown to be associated with poor disease prognosis, i.e., worse overall survival and worse disease-free survival (11). Understanding the molecular mechanisms underlying collective dissemination of tumor cells is, therefore, important for predicting metastasis, which remains the principal cause of cancer-associated mortalities. Determinants of single cell versus collective dissemination of tumor cells, however, remain elusive. Here, we have analyzed the topology of the network of genes implicated in the collective dissemination of tumor cells. We also investigated the topology of the network of genes reported to be differentially expressed in tumor samples from IBC patients as compared to tumor samples from non-IBC breast cancer patients. Taken together, our analysis suggests that maintenance of the epithelial phenotype in cancer cells disseminating from the primary tumor contributes toward metastasis *via* collective migration of tumor cells as CTC clusters.

Results suggest that the expression of genes differentially expressed in tumor cells migrating as clusters as compared to individually migrating tumor cells (12) exhibits a more hierarchical organization in epithelial cell lines as compared to mesenchymal cell lines among both, immortalized breast cancer cell lines (39) and the cancer cell lines in the NCI60 panel (41, 42). The importance of expression of such genes involved in cell migration, cell–extracellular matrix interaction, and cell–cell adhesion in the classification of NCI60 cell lines has been observed previously (64). Retention of some epithelial characteristics by cancer cells disseminating from the primary tumor has been reported to contribute toward collective invasion by tumor cells as CTC clusters (12, 65, 66). A more hierarchical organization in the expression of these genes may contribute toward a more robust epithelial phenotype in these cell lines (28–31). Higher hierarchical organization in the expression of these genes is also observed in tumor samples from IBC patients as compared to tumor samples from non-IBC breast cancer patients. This difference may contribute toward the strengthened presentation of epithelial characteristics such as cell–cell adhesion and juxtracrine signaling in tumor cells from IBC patients. The retention of these characteristics can foster the collective migration of tumor cells from the primary breast lesion (65). Further, our results reveal that hierarchical expression of collective dissemination-associated genes is of diagnostic relevance in IBC, thereby strengthening the case for IBC as a model system for the study of collective dissemination of tumor cells (22) and indicating the potential usefulness of mechanistic studies of tumor cell dissemination in determining the principles underlying IBC.

Next, we investigated the hierarchical organization in the expression of genes previously reported to be differentially expressed in tumor samples from IBC patients as compared to tumor samples from non-IBC breast cancer patients (38). The expression of these genes was more hierarchically organized in IBC samples as compared to non-IBC samples across multiple independent datasets. Further, epithelial cell lines exhibited a more hierarchical expression of these genes as compared to mesenchymal cell lines among immortalized breast cell lines (39) and among the cell lines in the NCI60 panel composed of nine different tumor types (41, 42). Thus, both collective dissemination-associated and IBC-associated genes exhibited a similar trend of higher CCC in immortalized breast cell lines or cancer cell lines as well as in tumor samples from IBC patients, adding to the existing evidence on collective dissemination *via* tumor emboli as the predominant mode of IBC metastasis and consequent aggressiveness. Intriguingly, the expression of canonical EMT-inducing transcription factors (52) was more hierarchically organized in non-IBC breast cancer samples as compared to IBC samples. Taken together, these results reinforce the notion that a complete EMT is not involved in IBC metastasis. Rather, it is the collective migration of tumor cells that are able to retain some epithelial characteristics that contributes toward the metastatic aggressiveness of IBC. The results presented here further strengthen the emerging notion that a complete EMT followed by MET is not necessarily as prevalent during cancer metastasis (7, 21) as posited earlier (67).

Both collective dissemination-associated and IBC-associated gene sets exhibited a higher CCC in breast cancer patients with faster posttreatment metastatic relapse as compared to patients with slower posttreatment relapse (53). A similar trend was observed in our calculations of the CCC for patients with SCLC (54), another metastatically aggressive cancer reported to metastasize *via* clusters of CTCs (55–57). These results indicate that a more hierarchical organization in the expression of genes involved in the collective dissemination of tumor cells may contribute toward a more aggressive behavior in metastatically aggressive tumors such as IBC and SCLC, which predominantly metastasize *via* clusters of CTCs. A mechanism-based investigation of the cross-talk between collective dissemination-associated and IBCassociated genes may, therefore, be a promising next step.

Further, samples from breast cancer metastases to lymph nodes and liver (59) exhibited a lower CCC as compared to breast cancer metastases to skin for collective dissemination-associated and IBC-associated gene sets (59). While metastasis of tumor cells to distant organs such as the skin is a complex, multi-step, and highly inefficient process, migration of tumor cells from the primary lesion to the local lymph nodes is likely to be a more facile process and can be brought about by the passive flow of the lymph. Metastasis to the liver is facilitated by the extravasation of migrating tumor cells into the liver *via* the fenestrated hepatic vascular epithelium (58). Correlation of the CCC for both gene sets, collective dissemination-associated and IBC-associated, with a higher rate of and propensity for metastasis to distant organs clearly speaks of the survival advantage afforded to migrating tumor cells by collective dissemination as clusters of CTCs. These advantages include enhanced ability to resist anoikis (cell death upon detachment from the substrate), evasion from immune system recognition, potential polyclonality, and enhanced ability to seed secondary tumors (68). In fact, CTC clusters can include non-tumor cells such as immune cells, platelets, and cancer-associated fibroblasts, thereby reproducing the primary tumor microenvironment conditions. Such an environment may contribute toward the survival of disseminating tumor cells in transit, promoting cancer metastasis (69).

A commonly used approach to determine if an *a priori* defined set of genes is associated with phenotypic differences between two groups is GSEA (70, 71). This method involves finding if the given set of genes is over-represented among genes that are differentially expressed in the two phenotypic groups. To determine if insights similar to those described above can be obtained *via* GSEA for the collective dissemination-associated gene set and for the IBC-associated gene set, we used the GSEA software provided by the Broad Institute (61) to calculate enrichment scores for the two gene sets in the data from Grosse-Wilde et al. (39), i.e., epithelial versus mesenchymal cell lines, and in the data from Iwamoto et al. (45), i.e., IBC versus non-IBC breast cancer patients. While we consistently obtained a higher CCC for collective dissemination-associated and IBC-associated gene sets in epithelial cell lines and in tumor samples from IBC patients, the expression of genes in these sets was not always enriched in epithelial versus mesenchymal cell lines or in IBC versus non-IBC patient samples. Together, these results indicate that the CCC need not correlate with GSEA. In fact, the CCC of a set of genes for two samples with a *k*-fold change in the expression of all genes in the set will be the same. The CCC can thus provide insights in addition to those that may be obtained from a direct analysis of gene expression data by using GSEA. The CCC of a gene network can be a robust metric of functional significance of a set of genes in different phenotypic groups, independent of the enrichment score calculated for the given gene set. It provides a prognostic measure based on the collective expression of genes in cells exhibiting different phenotypes beyond that provided by GSEA.

The classical view of cancer is that it involves de-differentiation of host cell pathways (36, 72). Since IBC is more metastatically aggressive as compared to non-IBC breast cancer, host cell pathways are likely to be more disrupted in tumor samples from IBC patients. This is indeed observed for breast tumor samples from the study by Iwamoto et al. (45). Of the 100 randomly generated gene sets, 41 exhibited a significantly higher CCC in the non-IBC breast cancer group as compared to the IBC group. This indicates that the host cell pathways are disrupted to a greater extent in IBC as compared to non-IBC breast cancer. However, structure in the pathways involving genes that promote cancer progression may be selected for as the disease advances. We previously showed that the expression of adult acute myeloid leukemia-associated genes is more hierarchically organized in samples from patients in whom the disease relapsed during the follow-up period as compared to patients that underwent complete remission upon treatment (37). Similarly, for breast cancer metastasis-associated genes, hierarchical organization was higher in patients who developed distant metastases during the follow-up period as compared to patients who did not (36). Here, we propose that due to the role of maintenance of an epithelial phenotype in collective dissemination of tumor cells and the subsequent metastatic efficiency of CTC clusters, a hierarchical organization in the expression of these genes may be selected for in metastatically aggressive cancers like IBC. A measure of hierarchical organization, here the CCC, can thus be a useful biomarker in cancer prognosis, particularly in the case of IBC.

### CONCLUSION

We have shown that a set of genes previously reported to be associated with the collective dissemination of tumor cells (12) is more hierarchically expressed in epithelial cell lines as compared to mesenchymal cell lines, thereby indicating a role for epithelial characteristics in the collective migration of tumor cells as clusters of CTCs. We further showed that IBC, an aggressive breast cancer subtype that metastasizes primarily *via* CTC clusters, exhibits a more hierarchical organization in the expression of these collective dissemination-associated genes as compared to non-IBC type breast cancer. Along similar lines, we showed that for genes differentially expressed in IBC as compared to non-IBC tumor samples, the expression is more hierarchical in tumor samples from IBC patients and in phenotypically epithelial cell lines, suggesting a a role for the retention of some epithelial traits in the metastatically aggressive nature of IBC. Taken together, our work indicates that at least some maintenance of the epithelial phenotype in disseminating tumor cells during disease progression plays a key role in successful metastasis of cancer to distant organs, and that IBC can be a suitable model system for studying mechanisms of collective migration of tumor cells as CTC clusters. Further, we have introduced the CCC as a quantitative metric for analyzing the collective migration of circulating tumor cell clusters, which can be useful in cancer prognosis, particularly in the case of IBC.

#### AUTHOR CONTRIBUTIONS

ST designed the study, carried out the analysis, and wrote the manuscript. MKJ designed the study, analyzed the results, and wrote the manuscript. WW and HL analyzed the results and edited the manuscript. MD supervised the study, analyzed the results, and edited the manuscript.

#### FUNDING

This work was supported by the Center for Theoretical Biological Physics, funded by the National Science Foundation (PHY-1427654). MKJ has a training fellowship from the Gulf Coast Consortia on the Computational Cancer Biology Training Program (CPRIT grant no. RP170593).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at https://www.frontiersin.org/articles/10.3389/fonc.2018.00244/ full#supplementary-material.

Figure S1 | Representative collective dissemination-associated gene networks for (A) epithelial cell lines and (B) mesenchymal cell lines from the study by Grosse-Wilde et al. (39), and for (C) tumor samples from IBC patients and (D) tumor samples from non-IBC breast cancer patients from the study by Iwamoto et al. (45). The nodes are collective dissemination-associated genes and the weights of the edges between different nodes were defined using Eq. 1.

Figure S2 | Normalized CCC of different gene sets calculated for different phenotypic groups using the resistance distance (2). Normalized CCC for 13 epithelial (E) and 11 mesenchymal cell lines from the study by Grosse-Wilde et al. (39), calculated using (A) collective dissemination-associated genes and (B) IBC-associated genes. Normalized CCC for 11 epithelial (E) and 47 epithelial–mesenchymal hybrid (E/M)+ mesenchymal (M) cell lines from the NCI60 dataset (41, 42), calculated using (C) collective disseminationassociated genes and (D) IBC-associated genes. Normalized CCC for tumor samples from the study by Iwamoto et al. (45) with 25 IBC and 57 non-IBC breast cancer patients, calculated using (E) collective disseminationassociated genes and (F) IBC-associated genes. Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \**p*-Value < 0.05. The trend in CCC values observed here is same as the trend when calculating CCC using the Euclidean commute time distance, Figures 1, 2A, 3, and 4A.

Figure S3 | Normalized CCC for estrogen-receptor-positive (ER+) breast cancer patients with metastatic relapse within a 30-month period posttreatment (*T* < 30; *n* = 36) or between 30 and 60 months posttreatment (*T* ≥ 30; *n* = 44): (A) normalized CCC of the collective dissemination-associated gene network and (B) normalized CCC of the IBC-associated gene network. Gene expression data from the study by Wang et al. (53). Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \*\**p*-value < 0.01. There were too few estrogen-receptor-negative (ER−) patients in the data set for similar analysis. The trend here is similar to the trend in Figure 6A.

Figure S4 | Normalized CCC for breast cancer patients with different estrogen-receptor statuses. Patients with estrogen-receptor-positive status: (A) normalized CCC of the collective dissemination-associated gene network and (B) normalized CCC of the IBC-associated gene network. In this group, there were 129 patients with no relapse during the 5-year follow-up period and 80 patients with metastatic relapse posttreatment during the follow-up period. Patients with estrogen-receptor-negative status: (C) normalized CCC of the collective dissemination-associated gene network and (D) normalized CCC of the IBC-associated gene network. In this group, there were 50 patients with no relapse during the 5-year follow-up period and 27 patients with metastatic relapse posttreatment during the follow-up period.

Gene expression data from the study by Wang et al. (53). Error bars indicate the SE in the estimate of CCCnorm calculated using the bootstrap method. \**p*-Value < 0.05 and \*\**p*-Value < 0.01. In (C), 10,000 bootstrap samples were drawn to calculate the normalized CCCs, obtain the error bars, and estimate the *p*-value.

Figure S5 | Gene set enrichment analysis on gene expression data for tumor samples from IBC and non-IBC breast cancer patients from the study by Iwamoto et al. (45) using (A) genes upregulated in cells in circulating tumor cell clusters, and (B) genes downregulated in cells in circulating tumor cell clusters. Genes are ordered from left to right in decreasing order of correlation of

#### REFERENCES


expression with the IBC phenotype. Black bars along the top of each plot indicate the positions of hits to the gene set along the ordered list of genes. Nominal *p*-values of enrichment are indicated at the bottom of each plot.

Figure S6 | Histogram of *p*-values calculated for the null hypothesis that the normalized CCC of a randomly generated gene set is higher in tumor samples from non-IBC breast cancer patients as compared to samples from IBC patients. Gene expression data from the study by Iwamoto et al. (45). Normalized CCC was calculated for 100 randomly generated gene sets consisting of 83 genes each. The red dotted line indicates *p*-value = 0.05 while the green dotted line indicates p-value = 0.95.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Tripathi, Jolly, Woodward, Levine and Deem. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Functional assay of cancer cell invasion Potential Based on Mechanotransduction of Focused Ultrasound

*Andrew C. Weitz1,2,3,4†, Nan Sook Lee1,3†, Chi Woo Yoon1 , Adrineh Bonyad3 , Kyo Suk Goo1 , Seaok Kim1 , Sunho Moon1 , Hayong Jung1 , Qifa Zhou1,4, Robert H. Chow3 \* and K. Kirk Shung1 \**

*1Ultrasonic Transducer Resource Center, University of Southern California, Los Angeles, CA, United States, 2 Institute for Biomedical Therapeutics, University of Southern California, Los Angeles, CA, United States, 3Zilkha Neurogenetic Institute, University of Southern California, Los Angeles, CA, United States, 4USC Roski Eye Institute, University of Southern California, Los Angeles, CA, United States*

#### *Edited by:*

*Triantafyllos Stylianopoulos, University of Cyprus, Cyprus*

#### *Reviewed by:*

*Fabio Grizzi, Humanitas Clinical and Research Center, Italy Balaji Krishnamachary, Johns Hopkins University, United States*

#### *\*Correspondence:*

*Robert H. Chow rchow@med.usc.edu; K. Kirk Shung kkshung@usc.edu*

*† These authors have contributed equally to this work.*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

> *Received: 04 May 2017 Accepted: 13 July 2017 Published: 07 August 2017*

#### *Citation:*

*Weitz AC, Lee NS, Yoon CW, Bonyad A, Goo KS, Kim S, Moon S, Jung H, Zhou Q, Chow RH and Shung KK (2017) Functional Assay of Cancer Cell Invasion Potential Based on Mechanotransduction of Focused Ultrasound. Front. Oncol. 7:161. doi: 10.3389/fonc.2017.00161*

Cancer cells undergo a number of biophysical changes as they transform from an indolent to an aggressive state. These changes, which include altered mechanical and electrical properties, can reveal important diagnostic information about disease status. Here, we introduce a high-throughput, functional technique for assessing cancer cell invasion potential, which works by probing for the mechanically excitable phenotype exhibited by invasive cancer cells. Cells are labeled with fluorescent calcium dye and imaged during stimulation with low-intensity focused ultrasound, a non-contact mechanical stimulus. We show that cells located at the focus of the stimulus exhibit calcium elevation for invasive prostate (PC-3 and DU-145) and bladder (T24/83) cancer cell lines, but not for non-invasive cell lines (BPH-1, PNT1A, and RT112/84). In invasive cells, ultrasound stimulation initiates a calcium wave that propagates from the cells at the transducer focus to other cells, over distances greater than 1 mm. We demonstrate that this wave is mediated by extracellular signaling molecules and can be abolished through inhibition of transient receptor potential channels and inositol trisphosphate receptors, implicating these proteins in the mechanotransduction process. If validated clinically, our technology could provide a means to assess tumor invasion potential in cytology specimens, which is not currently possible. It may therefore have applications in diseases such as bladder cancer, where cytologic diagnosis of tumor invasion could improve clinical decision-making.

Keywords: cancer invasion, focused ultrasound, calcium imaging, bladder cancer, prostate cancer

### INTRODUCTION

Cancer staging determines both patient prognosis and treatment protocol. To stage a biopsied tumor, the pathologist must determine the extent to which the tumor has invaded the surrounding tissue. In many instances, however, the intact tissue needed to assess invasion cannot be obtained from the patient. In such cases, fine-needle aspirations, washings, or brushings can be performed to collect cells from the tumor for cytologic diagnosis. These cells allow the cytopathologist to

**61**

determine whether the tumor is benign or malignant, but not whether it is invasive.

The inability to assess invasion can have devastating consequences, for example in the case of bladder cancer. Bladder tumors that have begun to invade the muscle wall must be treated promptly with cystectomy (surgical removal of the bladder), or they may become life-threatening (1). Carcinoma *in situ* (CIS) is an early form of bladder cancer that is considered high grade, as these tumors frequently recur as muscle-invasive disease (2). CIS is normally treated with bacillus Calmette–Guérin (BCG) immunotherapy upon initial diagnosis and recurrence (3). However, BCG inflames the bladder epithelium, making it difficult endoscopically to identify recurrent tumors for biopsy (4). In such cases, bladder wash cytology can be used to detect recurrence, but the invasion status of the recurrent malignancy cannot be determined. The inability to detect invasion precludes the use of preventative cystectomy, which can have fatal consequences if the cancer is indeed invasive. Thus, a method for assessing tumor invasion cytologically (e.g., in bladder washings) would enable appropriate and timely treatments that improve patient outcomes.

Classical cytology relies on examining cell morphology to identify the presence and appearance of malignant cells. Biophysical properties of tumor cells could reveal information about their malignant status that might escape detection in morphological studies. Recent work has revealed a number of biophysical changes that occur during cancer transformation and progression (5). For example, metastatic cells often express voltage-gated ion channels, including the Na<sup>+</sup> (6), K<sup>+</sup> (7), and Ca2<sup>+</sup> (8) types, rendering them electrically excitable. Similarly, mechanical properties of metastatic cancer cells differ from those of benign cells; metastatic cells are generally "softer" (9). Because these biophysical changes can serve as diagnostic markers, assays to measure biophysical properties of tumor cells have been proposed and are under development (9–11). These assays are typically intended to differentiate malignant from non-malignant cells in suspension, thus providing information similar to that obtained *via* standard cytological analysis.

Here, we present a new biophysical cancer assay that, to our knowledge, is the first that can assess the invasion potential of isolated tumor cells. The assay leverages the fact that cancer progression is accompanied by remodeling of calcium channels (8, 12, 13) and cellular mechanosensors (14). The assay applies a non-contact, mechanical stimulus—low-intensity focused ultrasound—to probe for presence of these proteins while monitoring their activity *via* calcium imaging. The stimulus elicits marked calcium elevations in invasive cancer cells, but not in non-invasive cells. We previously validated this assay in four breast cancer cell lines and demonstrated its effectiveness in quantifying invasion potential (15). However, analysis was limited to single cells, making the assay impractical for a clinical setting.

In the present study, we show that our assay can be used for high-throughput, functional analysis of cancer invasion in cell populations. A single ultrasound stimulus is applied while monitoring calcium activity in hundreds or thousands of cells simultaneously. We validate the technique in prostate and bladder cancer cell lines while exploring the role of different stimulus parameters. We also investigate the mechanism by which ultrasound elicits calcium elevations in invasive cells. If validated through further testing, this technology may lend itself to cytological assessment of tumor invasion, thus having important implications for diagnosis and management of diseases such as bladder cancer.

#### MATERIALS AND METHODS

#### Cell Lines

Six human cancer cell lines were used in this study: four prostate cancer lines (PC-3, DU-145, BPH-1, and PNT1A) and two bladder cancer lines (T24/83 and RT112/84). PC-3 cells were obtained from Frank Markland (University of Southern California), DU-145 and PNT1A from Mitchell Gross (University of Southern California), BPH-1 from Simon Hayward (NorthShore University HealthSystem), and T24/83 and RT112/84 from Sigma. PC-3 and DU-145 cells were cultured in DMEM, PNT1A, and BPH-1 cells in RPMI-1640, T24/83 cells in McCoy's 5A medium, and RT112/84 cells in EMEM. All medium was supplemented with 10% FBS and 2 mM l-glutamine, and RT112/84 medium was additionally supplemented with 1% non-essential amino acids. All cell lines were tested to be free of mycoplasma contamination using a mycoplasma PCR detection kit (Sigma). Cell authentication of all lines was performed with Promega's PowerPlex 16 System within 6 months of use.

PC-3, DU-145, and T24/83 are highly invasive cell lines, while BPH-1, PNT1A, and RT112/84 are weakly invasive (16–19). We confirmed the invasion status of each cell line with Matrigel Boyden chamber assays (20). Invasion potential was measured in BioCoat Matrigel Invasion Chambers (Corning), as described in our previous study (15). In brief, cells (1–1.5 × 105 ) were added to chambers and incubated for 1–2 days at 37°C. Non-invasive cells at the top of the chamber were removed by a cotton swab, and invasive cells that had passed through the Matrigel were stained with 0.2% crystal violet in 10% ethanol. Three independent fields of invasive cells per well were photographed under a microscope.

Preliminary results showed that BPH-1 cells exhibited variable levels of invasiveness over time, as measured by the Matrigel assay (Figure S1 in Supplementary Material) and our ultrasound assay (data not shown). In order to obtain a weakly invasive, homogeneous cell population, BPH-1 cells that passed through the Boyden chamber were removed after 24 h, and cells that did not pass through were selected and propagated. All BPH-1 experiments were performed with the selected cells, which were confirmed to be weakly invasive (see **Figures 2** and **3**). Experiments in the other five cell lines did not require selection, as these lines exhibited consistent levels of invasiveness over time.

#### Cell Preparation

Cells were plated on 35-mm polystyrene culture dishes to a density of 106 cells per dish. In one experiment, the substrate was coated with Cell-Tak Cell and Tissue Adhesive (Corning) to facilitate immediate cell adhesion. In another experiment, the culture dish substrate was replaced with an acoustically transparent, 50-μm-thick Mylar film (#48-2F-36; CS Hyde Company) to minimize reflection at the surface and eliminate ultrasonic surface waves. All cells were stained with cellmembrane permeant Fluo-4 AM (Thermo Fisher Scientific), a fluorescent reporter of intracellular calcium activity. Staining was performed by incubating dishes with 1 µM Fluo-4 AM for 30–60 min immediately prior to imaging. Following calcium dye loading, cells were washed with and maintained in external buffer solution consisting of 140 mM NaCl, 2.8 mM KCl, 1 mM MgCl2, 2 mM CaCl2, 10 mM HEPES, and 10 mM d-glucose, adjusted to pH 7.3 and 300 mOsm.

#### Ultrasound Transducers

Single-element, lithium niobate (LiNbO3), press-focused transducers were fabricated in house as described previously (21). A transducer with a center frequency of 38 MHz (f-number = 2, focal length = 8 mm) was used in most experiments (**Figure 1A**); a 3-MHz transducer (f-number = 2, focal length = 6 mm) was also tested. To drive the transducers, sinusoidal bursts from a signal generator (SG382; Stanford Research Systems) were fed to a 50-dB power amplifier (525LA; Electronics & Innovation) whose output was used to excite the transducer. Unless specified otherwise, amplitude was fixed at 16 Vp–p, pulse repetition frequency (PRF) at 1 kHz, and duty cycle at 5% (**Figure 1C**).

The acoustic output of the 38-MHz transducer was measured with a needle hydrophone (HGL-0085; Onda). **Figure 1B** shows the two-dimensional intensity beam profile, indicating that the diameter of the ultrasound focus is approximately 150 µm. Using the standard cell stimulation parameters provided above

FIGURE 1 | Cancer cells were stimulated with single-element, LiNbO3, press-focused transducers. (A) Photograph of the 38-MHz transducer used in most experiments. (B) Two-dimensional beam profile of *I*spta at the ultrasound focus, as measured by a hydrophone. (C) Typical voltage waveform used to drive the transducer. Carrier frequency is 38 MHz, amplitude is 16 Vp–p, pulse repetition frequency is 1 kHz, and duty cycle is 5%.

(16 Vp–p amplitude, 1 kHz PRF, and 5% duty cycle), the intensity and pressure at the focus were measured by the hydrophone to be 353 mW/cm2 spatial-peak temporal-average intensity (*I*spta), 7.0 W/cm2 spatial-peak pulse-average intensity (*I*sppa), and 394 kPa peak pressure. The mechanical index was measured to be 0.02. These values are below the FDA safety limit for diagnostic ultrasound (22).

#### Ultrasound Stimulation and Fluorescence Imaging

A custom microscope system was used to image cellular fluorescence while performing simultaneous ultrasonic stimulation (15). Petri dishes containing cells were placed on an inverted epifluorescence microscope (Olympus IX70), and the ultrasound transducer was lowered into the external buffer solution. A motorized three-axis micromanipulator was used to position the transducer in focus with the cell monolayer.

In each experiment, live-cell fluorescence imaging was performed for 300 s, with the ultrasound stimulus being delivered continuously between *t* = 50 and 200 s. Excitation light was provided by a mercury arc lamp and filtered through an excitation bandpass filter (488 ± 20 nm). Fluorescence emitted from the calcium dye was filtered through an emission bandpass filter (530 ± 20 nm) and recorded at 1 Hz (30% exposure duty cycle) with a digital CMOS camera (ORCA-Flash2.8; Hamamatsu). All imaging was performed at 4× magnification in order to capture activity from hundreds or thousands of cells simultaneously. For each cell line, simulation and imaging experiments were replicated in at least two different dishes of cells, and over least three independent fields of view per dish. (Experiments involving pharmacological blockers were limited to a single field of view per dish.) Figures show representative data obtained from one field of view.

#### Data Processing

Data were post-processed to determine the calcium response of every imaged cell. Cell locations were identified automatically with CellProfiler image analysis software (23) and used to extract the raw fluorescence intensities of each cell. These intensities were exported to MATLAB (MathWorks) in order to calculate each cell's normalized change in fluorescence (Δ*F*/*F*) during every imaging frame. Responding cells were defined as those that exhibited a Δ*F*/*F*max greater than 3.5 times the pre-stimulus root-mean-square noise level.

Two types of plots were generated for each 300-s experiment: a histogram showing the percentage of responding cells over time and a scatter plot indicating the time at which each cell first responded to the stimulus. Responding cells in these plots were arranged with respect to their distance from the transducer focus.

#### Pharmacology

To investigate the mechanism of ultrasound-induced calcium rise in invasive cancer cells, PC-3 and T24/83 cells were stimulated in the presence of various pharmacological agents. We tested five different blockers, each applied separately (**Table 1**). Blockers were dissolved in the external buffer solution 15–30 min before TABLE 1 | Pharmacological agents used to investigate the mechanism of ultrasound-induced calcium rise in invasive cancer cells.


performing imaging and ultrasound stimulation. Cellular responses were measured before adding the blockers, in the presence of blockers, and after washout.

#### RESULTS

#### Matrigel Boyden Chamber Assays

The Matrigel Boyden chamber assay (20) is the standard technique for assessing cancer cell invasion potential in the laboratory setting. We used this assay to measure the invasion potential of four prostate cancer cell lines (PC-3, DU-145, BPH-1, and PNT1A) and two bladder cancer cell lines (T24/83 and RT112/84) with varying levels of invasiveness. As expected, PC-3, DU-145, and T24/83 cells exhibited strong Matrigel invasion, while BPH-1, PNT1A, and RT112/84 cells did not (**Figure 2**; also see Figure S1 in Supplementary Material) (16–19).

Despite the simplicity of the Matrigel Boyden chamber assay, it is not suitable for clinical applications. The assay requires at least 10–25 thousand cells, necessitating time-consuming establishment of cell cultures from tumor biopsies (29, 30), which is not always possible. Furthermore, the Matrigel assay is relatively slow, taking at least 24 h to complete. Our mechanotransduction assay overcomes these limitations, as it can be completed in a matter of minutes and does not require large numbers of cells.

#### Calcium Responses to Ultrasound Stimulation Differ between Strongly and Weakly Invasive Prostate Cancer Cells

To investigate the use of mechanotransduction for determining the invasion potential of tumor cell populations, we imaged calcium activity of PC-3, DU-145, BPH-1, and PNT1A prostate cancer cells during stimulation with 38-MHz low-intensity focused ultrasound. The stimulus consistently evoked strong calcium responses in highly invasive PC-3 and DU-145 cells, but not in weakly invasive BPH-1 and PNT1A cells (**Figure 3**; Videos S1–S4 in Supplementary Material). In most cases, PC-3 and DU-145 stimulation evoked an intercellular calcium wave that emanated from the transducer focus and spread at a speed of ~50–100 μm/s. The wave began 5–25 s after stimulation onset and persisted for 1–2.5 min. Despite the acoustic focus being localized to an area approximately 150 µm in diameter (see **Figure 1B**), the calcium wave propagated over distances greater than 1 mm. Most cells that responded to ultrasound stimulation exhibited a single calcium transient, though double transients and calcium oscillations were also observed (**Figure 3B**).

FIGURE 2 | Assessment of prostate cancer (A) and bladder cancer (B) cell line invasion potential with Matrigel Boyden chamber assays. Cells that passed through the Matrigel barrier were stained with crystal violet. Results indicate that PC-3, DU-145, and T24/83 cells are highly invasive, while BPH-1, PNT1A, and RT112/84 cells are not. (A) Photomicrographs showing assay results for BPH-1, PNT1A, PC-3, and DU-145 prostate cancer cells. (B) Photomicrographs showing assay results for RT112/84 and T24/83 bladder cancer cells. One representative image for each cell line is shown.

#### Calcium Responses of Bladder Cancer Cell Lines Mirror Those of Prostate Cancer Cell Lines

As discussed above, a cytological assay of tumor cell invasion potential could have important implications for bladder cancer diagnosis and management. We therefore investigated whether the mechanotransduction assay could be applied to bladder cancer cells. We stimulated T24/83 (highly invasive) and RT112/84 (weakly invasive) bladder cancer cell lines with 38-MHz ultrasound while performing calcium imaging. As shown in **Figure 4A** and

(third column), and DU-145 (fourth column) cells were stimulated with 38-MHz low-intensity focused ultrasound while monitoring their activity with calcium imaging. Top row, baseline (pre-stimulation) fluorescence images of the cells. Red asterisks indicate the center position of the ultrasound focus. Scale bar is 500 µm. Second row, background-subtracted fluorescence images revealing all the cells that responded to the ultrasound stimulus (also see Videos S1–S4 in Supplementary Material). Only the PC-3 and DU-145 cells responded strongly, indicating that they are invasive. Third row, two-dimensional histograms showing the percentage of responding cells over time with respect to their distance from the transducer focus. Dotted green and red lines indicate stimulus onset and offset times, respectively. Bottom row, scatter plots showing the time at which each cell first responded to the stimulus (each dot represents a cell). The yellow shaded area indicates the times during which ultrasound was applied. Fitted blue lines indicate the speed of the calcium wave (77 µm/s in PC-3; 106 µm/s in DU-145). The histogram and scatter plots provide an intuitive means to visualize responses of the cell populations over time, making it easy to distinguish responding (invasive) cells from non-responding (non-invasive or weakly invasive) ones. (B) Fluorescence responses of three PC-3 cells showing the different types of calcium transients: single (left), double (middle), and oscillating (right).

Videos S5 and S6 in Supplementary Material, the stimulus evoked a calcium wave in T24/83 cells, but no responses were observed in RT112/84 cells. The pattern and timing of T24/83 responses were similar to those of PC-3 prostate cancer cells (see **Figure 3**).

#### Varying Stimulus Amplitude

We hypothesized that there was an acoustic activation threshold (i.e., intensity) below which invasive cancer cells would not exhibit calcium responses to ultrasound stimulation. To test for such a threshold, we stimulated T24/83 cells at amplitudes lower than 16 Vp–p (the voltage used in all prior experiments; **Figure 5A**). The smallest amplitude we tested, 2 Vp–p, did not evoke any detectable calcium activity. At 4 Vp–p, the stimulus activated just a few cells in the area where the ultrasound was focused. Stimulation at 8 Vp–p also elicited responses near the focus, but more cells responded at this amplitude. Increasing the

FIGURE 5 | Effect of ultrasound stimulus amplitude on bladder cancer cell calcium responses. Standard stimulus parameters (see Materials and Methods) were used while varying the transducer input voltage. All stated voltages represent peak-to-peak amplitude. Values in parentheses indicate the *I*spta at each voltage, as measured by a hydrophone. (A) Histogram and scatter plots from T24/83 cell stimulation. (B) Histogram and scatter plots from RT112/84 cell stimulation.

voltage to 16 Vp–p evoked hundreds of responses in the form of an intercellular calcium wave.

Given that RT112/84 bladder cancer cells are weakly invasive (19), we hypothesized that they might exhibit calcium responses at intensities greater than 16 Vp–p. As shown in **Figure 5B**, stimulation at 16 Vp–p did not evoke any calcium activity in these cells (the few data points on the histogram and scatter plots are false positives arising from movement of hyperfluorescent debris in the Petri dish; see Video S5 in Supplementary Material). At 32 Vp–p, three cells at the center of the transducer focus responded. At 48 Vp–p, 10 –15 cells at the focus responded, but there was still no calcium wave. Amplitudes higher than 48 Vp–p were not tested, as they likely would have damaged the ultrasound transducer.

These results reveal a dose–response relationship between stimulus amplitude and the strength of the calcium responses. For a given invasion potential, there appears to be an acoustic activation threshold (inflection point of the dose–response curve) below which no or a few cells respond to ultrasound stimulation. Given the inverse relationship between acoustic activation threshold and invasion potential, it may be possible to quantify the invasion potential of a tumor cell population by measuring its acoustic activation threshold.

#### Varying Stimulus Frequency

In this study, we have shown that 38-MHz focused ultrasound evokes calcium responses in invasive cancer cells, and we previously reported a similar effect for 200-MHz ultrasound (15). As shown in **Figure 6**, 3-MHz stimulation was also effective in eliciting calcium responses in invasive cells. Stimulating at 3 MHz evoked a calcium wave that propagated at 101 µm/s, similar to the speed of the calcium waves induced by 38-MHz stimulation (~50–100 μm/s; see **Figures 3A** and **4A**). These results suggest that the mechanism of ultrasound-induced calcium rise in invasive cancer cells is at least partly independent of stimulus frequency.

#### Mechanism of Ultrasound Stimulation

The fact that our assay can be applied to more than one cancer type (prostate and bladder) implies a conserved mechanism by which invasive cancer cells transduce ultrasound stimuli. To elucidate this mechanism, we applied pharmacological blockers of proteins we suspected were involved in the mechanotransduction process. We stimulated PC-3 and T24/83 cells in the presence of five different blockers, each applied separately (**Table 1**). We blocked voltage-gated Ca2<sup>+</sup> channels, which are known to be mechanosensitive (31), as well as stretch-activated Ca2<sup>+</sup> channels. We also blocked BKCa channels, which are stretch activatable (32) and expressed in PC-3 cells (33). Finally, we applied a drug that blocks both transient receptor potential (TRP) channels and inositol trisphosphate (IP3) receptors. TRP channels are mechanosensitive, permeable to calcium, and exhibit altered expression levels in several types of cancer (13, 34, 35). IP3 receptors are found on the endoplasmic reticulum and mediate Ca2<sup>+</sup> release from intracellular stores.

Of the five blockers we tested (**Table 1**), only 2-aminoethoxydiphenyl borate (2-APB) had an effect on the calcium responses. Application of 2-APB abolished all ultrasound-induced calcium activity, an effect that was partially reversed upon washout of the drug (**Figure 7**). This indicates that TRP channels and/or IP3 receptors are involved in mediating invasive cancer cell responses to ultrasound stimulation (see Discussion).

#### Mechanism of the Calcium Waves

Intercellular calcium waves are common biological phenomena that occur in many cell types. They play a role in a variety of cellular activities including migration and mechanotransduction (36). Calcium waves typically occur *via* transmission of an intracellular messenger such as IP3 through gap junctions (37), or by release of an extracellular messenger such as adenosine

baseline, before the drug was applied. Center column, cells were stimulated again 15 min after 2-APB (100 µM) application, and no responses were observed (the few data points on the histogram and scatter plots are false positives arising from movement of hyperfluorescent debris in the Petri dish). Right column, 30 min after 2-APB washout, ultrasound-stimulated calcium responses were partially restored. Results in T24/83 bladder cancer cells were similar.

triphosphate that diffuses to surrounding cells (38). In the case of our experiments, there was a third possibility—that ultrasound energy impinging upon the Petri dish surface was being reflected by the substrate and generating a surface wave that was activating distant cells.

To rule out the possibility of ultrasonic surface wave induced activation, we stimulated PC-3 cells seeded on an acoustically transparent Mylar film, thus minimizing ultrasound reflection by the substrate. As shown in **Figure 8** (left column), the ultrasound stimulus still evoked a calcium wave in these Mylar-seeded cells. This indicates that the wave was likely caused by cell-to-cell signaling, either through gap junctions or *via* release of an extracellular messenger.

To determine whether the calcium wave was propagating through gap junctions, we seeded PC-3 cells on a Petri dish coated with Cell-Tak Cell and Tissue Adhesive, in order to facilitate immediate cell adhesion. We then stimulated the cells within minutes, before they could form gap junctions. The calcium wave was not eliminated under these circumstances (**Figure 8**, center column), ruling out the possibility that it was propagating through gap junctions.

To test whether extracellular signaling molecules were generating the calcium wave, we used a pipette tip to scrape a channel of PC-3 cells off the substrate, leaving two regions of confluent cells separated by a 400-µm gap (**Figure 8**, top right). We then stimulated the cells on one side of the gap with ultrasound. The stimulus evoked a calcium wave that propagated across the gap, eliciting responses in cells on the other side (**Figure 8**, bottom right). These results support the notion that the calcium waves observed in invasive cancer cells are mediated by signaling molecules that diffuse through the extracellular solution.

#### DISCUSSION

We have demonstrated a new method for assessing the invasion potential of cancer cell populations and have validated the approach in prostate and bladder cancer cell lines. To our knowledge, this is the first technique other than the Matrigel Boyden chamber assay that can assess invasion of isolated tumor cells (as opposed to intact tissue). Unlike the Matrigel assay, our assay is rapid and does not require tens of thousands of cells. It can measure invasion potential of single cells as we have shown previously (15), or as shown in the present study, it can be applied to cancer cell populations.

We are not the first to investigate how low-intensity ultrasound interacts with invasive cancer cells. Tran et al. stimulated invasive MDA-MB-231 breast cancer cells with 1-MHz ultrasound (up to 500 kPa) while using the patch clamp technique to monitor membrane potential (32, 39). The cells immediately became hyperpolarized upon ultrasound application, but only when they were in direct contact with gas-filled microbubbles. By applying the BKCa channel blocker iberiotoxin, the authors determined that the microbubbles were activating mechanically sensitive BKCa channels, causing K<sup>+</sup> efflux and consequent hyperpolarization (32).

The BKCa-based mechanotransduction mechanism reported by the Tran study is fundamentally different from the mechanism

substrate. This result indicates that the wave is not caused by ultrasound reflection by the substrate. Center column, histogram, and scatter plots showing ultrasound responses of cells freshly seeded on a Cell-Tak-coated Petri dish. A calcium wave occurs even in the absence of gap junctions. Right column, baseline (top), and background-subtracted (bottom) fluorescence images of ultrasound responses occurring after scraping a channel of cells off the substrate. Red asterisks indicate the center position of the ultrasound focus. Calcium responses occurred on both sides of the gap, suggesting that calcium waves in invasive cancer cells are caused by release of signaling molecules that diffuse through the extracellular space. Scale bar is 500 µm.

observed in our study: We observed ultrasound responses in invasive cancer cells without using microbubbles. Furthermore, we found that PC-3 and T24/83 responses were not affected by iberiotoxin (**Table 1**), thus excluding involvement of BKCa channels in the mechanotransduction process. Instead, we found that the drug 2-APB blocked ultrasound-induced calcium activity (**Figure 7**), implicating a role of TRP channels and/or IP3 receptors. We also observed a delay of several seconds between stimulus onset and the calcium responses, suggesting a second messenger effect.

Several lines of evidence support the possibility that TRP channels are mediating the ultrasound responses observed in our study. TRP channels are a family of non-selective cation channels that can be activated by mechanical force, either directly or through a second messenger pathway (34). They have been implicated in several types of cancer including prostate (40) and bladder (41) cancer, regulating behaviors such as proliferation, differentiation, and migration [see Ref. (13, 35) for review]. Invasive and metastatic cancers are known to overexpress certain TRP channel isoforms, and silencing expression of these isoforms in cell lines significantly reduces migration and invasion (42, 43). It was recently reported that TRP-4, the *C. elegans* TRP channel homolog, transduces low-pressure (≤900 kPa) ultrasound stimuli in this model organism [though only in the presence of microbubbles; Ref. (44)].

Because 2-APB blocks both TRP channels and IP3 receptors, it is possible, however, that IP3 receptors are involved in transducing the ultrasound stimulus (instead of or in addition to TRP channels). IP3 and its receptors play a dominant role in transducing external stimuli into calcium signals by evoking Ca2<sup>+</sup> release from intracellular stores. They are also known to mediate oscillatory changes in cytosolic calcium concentration (45), such as those observed in this study (see **Figures 3B** and **4B**). In ongoing work, we aim to determine the precise mechanism of ultrasound-induced calcium rise *via* small interfering RNA-mediated downregulation of IP3 receptors versus TRP channels.

Although we found that ultrasound stimulation consistently evoked calcium waves in invasive cancer cells, the timing of those waves was somewhat variable. Calcium waves typically began 5–25 s after stimulation onset, traveled at a speed of ~50–100 μm/s, and persisted for 1–2.5 min. Though the reason for this variability is not clear, it is possible that the calcium wave timing correlates with invasion potential (for example, a faster wave could indicate greater invasion potential). In support of this theory, we did observe spontaneous variations in the degree of PC-3 and T24/83 Matrigel invasiveness that occurred during passaging in culture (data not shown). Likewise, others have reported that prostate and bladder cancer cells undergo cyclical, population-wide changes in tumorigenicity as they are passaged (46). Future studies will explore whether these changes in tumorigenicity and Matrigel invasiveness correlate with the pattern of ultrasound responses (such as calcium wave speed, percentage of responding cells, etc.).

If validated using clinical specimens, our mechanotransduction assay could provide pathologists with a means to detect tumor invasion in cytology specimens, a feat that is not currently possible. At present, pathologists can assess invasion only through histological analysis of biopsied tissue. However, there are many instances when intact tissue cannot be obtained from the patient. In such cases, diagnosis often relies exclusively on cytology, which cannot assess invasion. This can have devastating consequences, for example in the case of recurrent bladder cancer. As discussed above, recurrent invasive bladder malignancies sometimes go undetected by cystoscopy, which can be life-threatening (1, 4). An assay for identifying invasion in bladder wash cytology specimens could inform treatment decisions that improve patient outcomes. Esophageal carcinoma is another disease that could benefit from cytological assessment of tumor invasion. Detecting invasion in cells collected from esophageal brushings could limit the need for endoscopic tumor biopsies, which carry the serious risk of esophageal perforation (47).

In addition to assessing invasion potential, our technology could also be as a high-throughput optical screen for drugs that target the mechanotransduction pathway (48). Cells could be placed in multiwell plates, treated with different drugs, and then stimulated with ultrasound while imaging calcium activity. Drug effects would be indicated by any changes induced in the pattern of ultrasound responses.

A functional assay of cancer cell invasion potential, as demonstrated herein, could provide key advantages over other types of cancer diagnostic assays. Genomic tests, for example, are intended to predict tumor aggression or recurrence and have gained widespread use in recent years. These tests are limited, however, in that they provide information about nucleic acid expression, rather than protein expression/ translation or protein functional state, which can be altered post-translationally. Protein activity is what controls a tumor's behavior, including its ability to invade and metastasize (5). By probing cell function as a measure of invasion potential, our assay may provide a more precise measure of a tumor's propensity to spread.

It may eventually be possible to develop an *in vivo* version of the assay that could mitigate the need for tumor biopsies. In the case of bladder cancer, a confocal laser endoscope (49) with an integrated ultrasound transducer could deliver calcium dye to the tumor, stimulate it with ultrasound, and image cellular responses. For other types of solid cancers, the tumor could be stimulated percutaneously with ultrasound while using calcium-sensitive MRI contrast agents (50) to image the response. This approach would provide an entirely non-invasive means to assess tumor invasion potential.

#### AUTHOR CONTRIBUTIONS

AW, NL, CY, AB, KG, SM, and HJ collected and analyzed data. SK analyzed data. AW, NL, CY, QZ, RC, and KS designed the research. AW and NL wrote the manuscript. All authors read, contributed to, and approved the final manuscript.

#### ACKNOWLEDGMENTS

The authors thank Changyang Lee and Ruimin Chen for their assistance with ultrasound transducer characterization.

#### FUNDING

This work was supported by NIH grant numbers R01 EB012058, P41 EB002182, and R01 GM85791.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fonc.2017.00161/ full#supplementary-material.

FIGURE S1 | Matrigel Boyden chamber assays of BPH-1 cells at different time points. Cells that passed through the Matrigel barrier were stained with crystal violet. As cells were passaged in culture, their level of invasiveness changed spontaneously over time. To obtain a weakly invasive homogeneous population for use in ultrasound stimulation experiments, BPH-1 cells that did not pass through the Boyden chamber were selected and propagated (see Materials and Methods).

VIDEO S1 | Calcium responses of weakly invasive BPH-1 prostate cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 3, first column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× realtime speed.

VIDEO S2 | Calcium responses of weakly invasive PNT1A prostate cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 3, second column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× realtime speed.

VIDEO S3 | Calcium responses of strongly invasive PC-3 prostate cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 3, third column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× real-time speed.

VIDEO S4 | Calcium responses of strongly invasive DU-145 prostate cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 3, fourth column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× realtime speed.

VIDEO S5 | Calcium responses of weakly invasive RT112/84 bladder cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 4, left column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× real-time speed.

VIDEO S6 | Calcium responses of strongly invasive T24/83 bladder cancer cells to stimulation with 38-MHz low-intensity focused ultrasound (video corresponds to images in Figure 4, right column). The red asterisk indicates the center position of the ultrasound focus; its appearance and disappearance coincide with stimulation onset and offset, respectively. Video is played back at 60× realtime speed.

## REFERENCES


**Conflict of Interest Statement:** AW, NL, RC, and KS have applied for a patent (US 14/040,253) related to this work.

*Copyright © 2017 Weitz, Lee, Yoon, Bonyad, Goo, Kim, Moon, Jung, Zhou, Chow and Shung. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## The expression and Prognostic impact of immune cytolytic activity-related Markers in human Malignancies: a comprehensive Meta-analysis

*Constantinos Roufas1,2, Dimitrios Chasiotis1 , Anestis Makris1 , Christodoulos Efstathiades2 , Christos Dimopoulos2 and Apostolos Zaravinos1 \**

*1Department of Life Sciences, Biomedical Sciences Program, School of Sciences, European University Cyprus, Nicosia, Cyprus, 2 The Center for Risk and Decision Sciences (CERIDES), Department of Computer Sciences, School of Sciences, European University Cyprus, Nicosia, Cyprus*

#### *Edited by:*

*Triantafyllos Stylianopoulos, University of Cyprus, Cyprus*

#### *Reviewed by:*

*Abhishek D. Garg, KU Leuven, Belgium Gabriele Multhoff, Technische Universität München, Germany*

*\*Correspondence:*

*Apostolos Zaravinos a.zaravinos@euc.ac.cy*

#### *Specialty section:*

*This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology*

*Received: 21 November 2017 Accepted: 29 January 2018 Published: 21 February 2018*

#### *Citation:*

*Roufas C, Chasiotis D, Makris A, Efstathiades C, Dimopoulos C and Zaravinos A (2018) The Expression and Prognostic Impact of Immune Cytolytic Activity-Related Markers in Human Malignancies: A Comprehensive Meta-analysis. Front. Oncol. 8:27. doi: 10.3389/fonc.2018.00027*

Background: Recently, immune-checkpoint blockade has shown striking clinical results in different cancer patients. However, a significant inter-individual and inter-tumor variability exists among different cancers. The expression of the toxins granzyme A (GZMA) and perforin 1 (PRF1), secreted by effector cytotoxic T cells and natural killer (NK) cells, were recently used as a denominator of the intratumoral immune cytolytic activity (CYT). These levels are significantly elevated upon CD8+ T-cell activation as well as during a productive clinical response against immune-checkpoint blockade therapies. Still, it is not completely understood how different tumors induce and adapt to immune responses.

Methods: Here, we calculated the CYT across different cancer types and focused on differences between primary and metastatic tumors. Using data from 10,355, primary tumor resection samples and 2,787 normal samples that we extracted from The Cancer Genome Atlas and Genotype-Tissue Expression project databases, we screened the variation of CYT across 32 different cancer types and 28 different normal tissue types. We correlated the cytolytic levels in each cancer type with the corresponding patient group's overall survival, the expression of several immune-checkpoint molecules, as well as with the load of tumor-infiltrating lymphocytes (TILs), and tumor-associated neutrophils (TANs) in these tumors.

results: We found diverse levels of CYT across different cancer types, with highest levels in kidney, lung, and cervical cancers, and lowest levels in glioma, adrenocortical carcinoma (ACC), and uveal melanoma. GZMA protein was either lowly expressed or absent in at least half of these tumors; whereas PRF1 protein was not detected in almost any of the different tumor types, analyzing tissue microarrays from 20 different tumor types. CYT was significantly higher in metastatic skin melanoma and correlated significantly to the TIL load. In TCGA-ACC, skin melanoma, and bladder cancer, CYT was associated with an improved patient outcome and high levels of both GZMA and PRF1 synergistically affected patient survival in these cancers. In bladder, breast, colon, esophageal, kidney, ovarian, pancreatic, testicular, and thyroid cancers, high CYT was accompanied by upregulation of at least one immune-checkpoint molecule, indicating that similar to melanoma and prostate cancer, immune responses in cytolytic-high tumors elicit immune suppression in the tumor microenvironment.

conclusion: Overall, our data highlight the existence of diverse levels of CYT across different cancer types and suggest that along with the existence of complicated associations among various tumor-infiltrated immune cells, it is capable to promote or inhibit the establishment of a permissive tumor microenvironment, depending on the cancer type. High levels of immunosuppression seem to exist in several tumor types.

Keywords: granzyme A, perforin 1, immune cytolytic activity, metastasis, cancer immunotherapy, survival rate, tumor-infiltrating lymphocytes, tumor-associated neutrophils

#### INTRODUCTION

In normal cells, the role of cytotoxic T lymphocyte antigen-4 (CTLA-4 or CD152), programmed death-1 (PD-1 or CD279), or other similar immune-checkpoint molecules is to inhibit an autoimmune response and restrict an immune cell-mediated tissue damage. Cancer cells on the other hand, regularly use these immune-checkpoint molecules to escape from being detected and eliminated by the cells of the immune system (1–3).

Cytotoxic T cells (CTLs) and natural killer (NK) cells release perforin 1 (PRF1), granzymes, and granulysin, upon their expose to infected or dysfunctional somatic cells. The first cytotoxin polymerizes and creates a channel in the membrane of the target cell. Through these pores, granzymes will then enter the cytoplasm and trigger a caspase cascade, composed of cysteine proteases that will ultimately lead to apoptosis (4, 5). However, apoptosis can also be induced *via* cell–surface interaction between the CTL and the infected cell. Upon the activation of a CTL, the FAS ligand (FasL or CD95L) is expressed on its surface, and it binds to Fas (CD95) being expressed on the target cell (6). Furthermore, the TNF-related apoptosis-inducing ligand (TRAIL) and its receptors (TRAILR1/2) constitute another important axis of immune cytolytic activity (CYT) that leads to apoptosis (7).

Apart from tumor cells, the tumor microenvironment contains many different immune cell types, including neutrophils, macrophages, dendritic cells (DCs), NK cells, T and B cells (8–10). Spontaneous tumor immunity due to the infiltration of such immune cells to the tumor site (11) and immunotherapy can be used to predict the patient outcome in cancer (12–14). However, it is now known that these nonmalignant tumor-infiltrating immune cells can also contribute to cancer by taking part in the modulation of the tumor microenvironment together with other nonimmune stromal cells, including fibroblasts and endothelial cells (15–17).

Immunotherapies that depend on the blockade of such immunecheckpoint molecules can stimulate an anticancer response (18–21). Among them, PD-1 targeting drugs (Pembrolizumab and Nivolumab) or PD-L1 (Atezolizumab, Avelumab, and Durvalumab), and CTLA-4 inhibitors (Ipilimumab) can benefit treatment of several cancer types, comprising skin melanoma, non-small cell lung cancer, kidney cancer, bladder cancer, head and neck cancer, and Hodgkin lymphoma (22–24). Nevertheless, success rate varies from one tumor type to other and some cancers do not respond to therapy or they gradually develop resistance to it.

The interactions between cancer cells and cells of the immune system can be further understood using high-dimensional genomic and transcriptomic datasets stored in online repositories. One such publically available repository is The Cancer Genome Atlas (TCGA),1 which contains comprehensive, multidimensional maps of the key genomic changes in 33 different cancer types. Latest analysis of the TCGA datasets has linked the genomic landscape of tumors with tumor immunity, implicating neoantigen load in driving T-cell responses (25), and identifying somatic mutations associated with immune infiltrates (26). The Human Protein Atlas (HPA)2 (27–30) is another open access platform that provides a map to all the human proteins in cells, tissues, and organs, and integrates different "omics" technologies, such as antibody-based imaging, mass spectrometry-based proteomics, transcriptomics, and systems biology.

Here, we have used a large number of TCGA and HPA datasets containing thousands of solid tumor samples to understand how different cancers induce and adapt to immune responses. RNA-seq data for the genes of interest were extracted from different datasets in Fragments Per Kilobase Million (FPKM) and subsequently transformed to Transcripts Per Kilobase Million (TPM) values using the formula TPMi = FPKMi/sum(FPKMj) × 106 . We have further supported the RNA-level information using protein-level data across all cancer datasets. The CYT from each dataset has been further associated with the corresponding patient group's overall survival. To associate the CYT with patient survival both in primary and metastatic cancers, we have focused our attention on skin melanoma, breast, and thyroid cancers. We have also evaluated the density of tumor-infiltrating lymphocytes (TILs) and tumor-associated neutrophils (TANs) using hematoxylin and eosin (H&E)-stained sections of primary and metastatic tumors and made associations of their load with patient survival in each type of cancer.

<sup>1</sup>https://cancergenome.nih.gov/.

<sup>2</sup>https://www.proteinatlas.org/.

#### MATERIALS AND METHODS

#### Cancer Datasets

Using the Genomic Data Commons (GDC) Data Portal (The Cancer Genome Atlas, TCGA program3 ) and the GTEx web portal (Genotype-Tissue Expression project4 ), we extracted data from a total of 10,355 tumor resection samples and 2,935 normal samples and screened the variation of CYT across these 32 different cancer types and 28 different normal solid tissue types. TCGA-derived data represent mainly untreated primary tumors (*n* = 9,913). In addition, we extracted 47 recurrent and 395 metastatic cancer cases. The Skin Cutaneous Melanoma (SKCM) dataset included the majority of these metastatic cases (*n* = 368). Patients who received neoadjuvant therapy were excluded from the analysis. Where available, TCGA tumor samples were paired with their corresponding normal tissues, providing a germline reference.

In specific, the following tumor types were selected: diffuse large B-cell lymphoma (DLBCL, *n* = 48), kidney clear cell cancer (KIRC, *n* = 539), kidney papillary cancer (KIRP, *n* = 289), kidney chromophobe cancer (KIRCH, *n* = 65), testicular germ cell cancer (TGCT, *n* = 156), lung adenocarcinoma (LUAD, *n* = 535), lung squamous cell carcinoma (LUSC, *n* = 502), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, *n* = 306), thymoma (THYM, *n* = 119), (SKCM, *n* = 471), acute myeloid leukemia (LAML, *n* = 151), head and neck squamous cell carcinoma (HNSC, *n* = 502), pleural mesothelioma (MESO, *n* = 86), sarcoma (SARC, *n* = 263), stomach adenocarcinoma (STAD, *n* = 375), colorectal cancer (COAD, *n* = 480), and rectum adenocarcinoma (READ, *n* = 167), uterine corpus endometrial carcinoma (UCEC, *n* = 552), uterine carcinosarcoma (UCS, *n* = 56), bladder cancer (BLCA, *n* = 414), pancreatic cancer (*n* = 178), breast cancer (BRCA, *n* = 1109), bile duct cancer (*n* = 36), ovarian serous cystadenocarcinoma (OV, *n* = 379), liver hepatocellular carcinoma (LIHC, *n* = 374), thyroid carcinoma (THCA, *n* = 510), esophageal cancer (*n* = 162), prostate adenocarcinoma (PRAD, *n*= 499), glioblastoma (GBM, *n*= 169), brain lower grade glioma (LGG, *n* = 529), pheochromocytoma and paraganglioma (PCPG, *n* = 183), adrenocortical carcinoma (ACC, *n* = 79), and uveal melanoma (UVM, *n* = 80) (where each acronym denotes the corresponding project's code and "*n*" is the number of cancer tissue samples).

"Level 3" mRNA-Seq expression data of the genes of interest, along with the corresponding patient clinical information for each disease type (tumors and normals) were extracted from TCGA public access web portal [launch data portal3 ] and GTEx4 (for normal samples only). Gene expression data were additionally accessed from the Fantom5 Consortium5 and were used to evaluate gene expression markers.

We also retrieved protein expression data derived from antibody-based protein profiling using immunohistochemistry (IHC) from the Tissue Atlas of The Human Protein Atlas (HPA) (27–29). Information regarding the cellular distribution of each cytolytic protein (GZMA and PRF1) was also retrieved across all major cancers from the same repository. In total, we extracted IHC data from 19 different tumor types, among them BRCA (*n* = 12), cervical cancer (*n* = 11), colorectal cancer (*n* = 11), endometrial cancer (*n* = 12), glioma (*n* = 12), head and neck cancer (*n* = 3), liver cancer (*n* = 11), lung cancer (*n* = 12), lymphoma (*n* = 12), melanoma (*n* = 12), ovarian cancer (*n* = 12), pancreatic cancer (*n* = 10), prostate cancer (*n* = 10), renal cancer (*n* = 11), skin cancer (*n* = 11), stomach cancer (*n* = 11), testis cancer (*n* = 9), thyroid cancer (*n* = 4), and urothelial cancer (*n* = 11).

#### Calculation of CYT Followed by Downstream RNA-seq and Protein Profiling Analyses

We calculated the CYT (or "cytolytic index") as the geometric mean of GZMA and PRF1, as formerly defined (31). Briefly, we divided the total raw read counts per gene by the gene's maximum transcript length to signify a coverage depth estimate. Coverage estimates were then scaled to sum to a total depth of 1e6 per sample and inferred as Transcripts Per Kilobase Million (TPM). We compared the cytolytic index between metastatic and non-metastatic (primary) cancers, wherever a sufficient number of metastatic tumor cases were available (TCGA-BRCA, TCGA-SKCM, and TCGA-THCA datasets). We also calculated the expression of several other CTL/NK or non-CTL/NK expressing genes, including immunosuppressive factors, the C1Q complex, and interferon-stimulated chemokines, all of which were previously shown to associate with an increased CYT in cancer. We further correlated the cytolytic index with the expression of immune-checkpoint molecules, including CTLA-4, PD-1, CD274 (PD-L1), PDCD1LG2 (PD-L2), LAG3, IDO1, CD73 (NT5E), and ENTPD1 (CD39), across all TCGA datasets. The *p*-values from the comparisons of the CYT between tumor and normal samples or between metastatic and primary cancer samples were FDRadjusted. Loess regression was applied to diminish the noise of the variables during correlation analysis.

We also extracted GZMA and PRF1 protein expression data from the Tissue Atlas of HPA, and further analyzed them. GZMA was stained with an anti-GZMA antibody produced in rabbit (HPA054134, 1:200 dilution, Sigma-Aldrich) and PRF1 using two different antibodies produced in rabbit (either HPA037940, 1:29 dilution, or CAB002436, 1:10 dilution, Sigma-Aldrich) (27–29).

#### Overall Survival and Synergistic Target Analysis on the TCGA Datasets

We performed Kaplan–Meier curves analysis to calculate the overall survival of each TCGA-dataset's patient group, based on their cytolytic index, TIL, and TAN load, or specific tumor subtype (e.g., triple negative vs triple positive BRCA). In total, we assessed overall survival data of patients suffering from 25 different cancer types (37 TCGA-datasets). Analysis was performed using the log-rank (Mantel Cox) test with a statistical significance at the 95% level (*p* < 0.05). We further tested the synergistic effect of the genes PRF1 and GZMA on each dataset's patient survival outcome, using SynTarget (32).

<sup>3</sup>https://portal.gdc.cancer.gov/.

<sup>4</sup>https://www.gtexportal.org/home/.

<sup>5</sup>http://fantom.gsc.riken.jp/5/.

#### Detection and Quantification of Lymphocyte and Neutrophil Infiltration among Primary and Metastatic Cancers

We extracted digital slide images with H&E-stained histological slides of skin melanoma, breast, and thyroid cancer from The Cancer Digital Slide Archive (CDSA)6 and compared the load of TILs and TANs between metastatic and primary cancers. TILs were distinguished by the typical features of lymphocytes (33), including size, shape, and staining of the nucleus. The percentage (%) of lymphocyte and neutrophil infiltration was compared to the information extracted from the corresponding datasets at the GDC Data Portal. We further compared the percentage of necrosis between primary and metastatic cancers, as well as the percentage of tumor, normal, and stromal cells. We correlated the levels of immune CYT (TPM counts) with the load of TILs and TANs, as well as with the percentage of necrosis found among metastatic and non-metastatic breast, skin melanoma, and thyroid cancers, using Pearson's correlation test.

### RESULTS

### Immune CYT across Different Tumor Types

To assess the intratumoral immune cytolytic T-cell activity across various tumor types, we quantified the transcript levels of GZMA and PRF1, as previously done by Rooney et al. (31). GZMA is a tryptase leading to apoptosis through the caspase pathway, whereas PRF1 is a pore-forming enzyme facilitating the entrance of granzymes into the target cells. Both effector molecules are considerably overexpressed upon CD8+ T-cell activation (34) and during productive clinical responses to anti-CTLA-4 or anti-PD-L1 immunotherapy (12, 13). CTL/NK cells can kill cancer cells by overexpressing GZMA and PRF1. We now know that effector T cells at the tumor site are good predictors of a favorable outcome across various cancer types (35–40).

Although Rooney et al. previously measured the immune CYT of the local immune infiltrate across various tumor types (31), some datasets did not contain enough data at the time (e.g., there were only three normal cervix samples in the TCGA-CESC dataset). Given the increased number of tumor samples in the TCGA platform since 2014, we have now significantly enlarged the total number of different cancer types, from 18 to 32. We have also considerably increased the sample number in many datasets, thus providing an opportunity to better estimate the different cytolytic levels across diverse tumors.

Consistent with previous findings (31, 41), we found that the cytolytic index was highest in the kidney (in clear cell and papillary renal cell carcinoma, but not in chromophobe carcinoma), lung, and cervical cancers. Importantly, we show for the first time that DLBCL and testicular cancer also rank among the top cytolytic active tumors, with DLBCL exhibiting even higher cytolytic levels compared to KIRC (>100 TPM). In addition, melanoma and head and neck cancer exhibited significantly higher CYT compared to the corresponding normal tissues. Acute myeloid leukemia, pleural mesothelioma, sarcoma, and stomach cancer, also exhibited high tendency in CYT. On the contrary, ovarian, liver, thyroid, esophageal, and prostate cancers, as well as glioblastoma, glioma, pheochromocytoma and paraganglioma, adrenocortical carcinoma and uveal melanoma, exhibited the lowest cytolytic indexes (**Figure 1A**).

Although most normal tissues (11 tissues from TCGA or GTEx) showed significantly lower CYT compared to their corresponding tumors, some of them exhibited significantly higher activity. Specifically, lung cancer, thymoma, stomach, colorectal, uterine, bladder, breast, liver, and thyroid cancers, all exhibited lower CYT compared to their corresponding normal tissues. In the cases of lung adenocarcinoma, colorectal, uterine, liver, and thyroid cancers, the differences between cancer, and the normal tissues were statistically significant (**Figure 1A**). The vast range in CYT across different cancers and compared to their corresponding normal tissues reveals the existence of a combination of tissue- and tumor-specific mechanisms that control local immunity. In line with their synchronized roles, the expression of GZMA and PRF1 was strongly coordinated across the different cancer samples (Spearman rank correlation, rho = 0.87) (**Figure 1B**).

At the protein level, we analyzed tissue microarray (TMA) data from 20 different tumor types, and found that GZMA was either lowly expressed or absent in at least half of these tumors, whereas, PRF1 was not detected in almost any of the different tumor types (**Figure 1C**). GZMA exhibited medium protein expression in the majority of the pancreatic cancers (70%), in <35% of breast, cervical, liver, ovarian, prostate, renal, stomach, testis, and urothelial cancers, as well as in <10% of lymphomas and melanomas. These data are consistent with the low TPM values derived from our RNA-seq analysis (**Table 1**). Further information regarding anti-GZMA and anti-PRF1 antibody staining, intensity, quantity, and location are provided in Table S1 in Supplementary Material.

#### Immune CYT in Primary and Metastatic Cancers

Next, we focused our attention on whether the cytolytic index differs between primary and metastatic cancers. Among all TCGA datasets, the metastatic tumors were composed of 368 SKCM, seven BRCA, eight THCAs, one PRAD, two cervical cancers (CESC), one colorectal adenocarcinoma (COAD), one esophageal carcinoma (ESCA), two HNSCs, one pancreatic adenocarcinoma (PAAD), two PCPGs, one PRAD, and one sarcoma (SARC) sample. Therefore, since the majority of metastatic tumors were composed mainly of skin melanomas, breast and thyroid carcinomas we focused our downstream analysis on the corresponding datasets of these tumors.

To assess the cytolytic index in them, we obtained RNA-seq data from TCGA for 103 primary and 368 metastatic skin resection melanomas, 1,102 primary and seven metastatic BRCAs, as well as for 502 primary and eight metastatic THCAs. Although all metastatic tumors had higher CYT compared to their corresponding primary tumors, the difference was statistically significant only for the skin melanoma dataset. This is obviously due to the significantly higher number of metastatic melanoma cases (*n* = 368) (**Figure 2**).

<sup>6</sup>http://cancer.digitalslidearchive.net/.

all cancers, a Spearman rank correlation (*r*) of 0.87 was observed. (C) Low levels GZMA and PRF1 protein expression detected in tissue microarrays of 20 different

tumor types. All representative immunohistochemistry images of each tumor type derived from the Human Protein Atlas.

Similarly, we investigated the expression of various suppressive factors previously shown to be associated with CYT, and compared their expression levels between metastatic and primary tumors. These genes included the immune-checkpoint molecules, CTLA-4, PD-1 (PDCD1), PD-L1 (CD274), PD-L2 (PDCD1LG2), LAG3, CD73 (NT5E)/CD39 (ENTPD1), IDO1/2, DOK3, the GMCSF receptors (CSF2RA, CSF2RB) (42), CD70, UBD, DOC3, NKG7, PLA2G2D, and the C1Q complex. We also included interferon-stimulated chemokines that attract T cells (CXCL9, CLCL10, and CXCL11) (11). We further investigated the expression of alternative genes through which T cells can induce cytolysis of cancer cells, including CD95-CD95L (FAS-FASLG) and TRAIL-TRAILR (TNFSF10, TNFRSF10A/B). Among the investigated genes, CD247, GZMK, GZMH, NKG7, PRF1, GZMA, GZMB, GZMH, GZMK, CD3E, and CD2 are expressed in CTL/NK cells; whereas CSF2RB, LTA, DOK3, PDCD1LG2, IDO1, PLA2G2D, CXCL9, CXCL10, CXCL11, CXCL13, UBD, C1QA, C1QB, C1QC, BATF2, and CSF2RA are expressed in non-CTL/NK cells (31).

In TCGA-SKCM, all genes (apart from CD70) exhibited significantly higher levels in metastatic skin melanomas compared to primary tumors. We also noticed a similar, but non-significant trend in datasets TCGA-BRCA and TCGA-THCA, presumably due to the small sample number of metastatic cases (**Figures 3**–**5**).

#### Kaplan–Meier and Synergistic Survival Analysis of GZMA and PRF1 across TCGA-Datasets

We next performed Kaplan–Meier survival analysis on 37 TCGA-datasets deriving from 25 different cancer types in order to estimate the risk of individual and/or simultaneous high (or low) PRF1 and GZMA expression on patient overall survival.

In TCGA-ACC, non-metastatic cutaneous melanoma ("m0" TCGA-SKCM), and bladder urothelial carcinoma (TCGA-BLCA but not the GSE32894 dataset), both individual and simultaneous high levels of PRF1 and GZMA were significantly associated with better prognosis. On the reverse, simultaneous low expression of both genes led to a significant shift toward negative effect vs all other ACC (or SKCM) patients. As expected, metastatic melanoma sufferers succumbed much earlier than non-metastatic skin melanoma patients did. These data provide significant evidence that high expression of both cytolytic genes in these cancer types, synergistically affects patient survival (**Figure 6A**).

Table 1 | Protein expression profiles of granzyme A (GZMA) and perforin 1 (PRF1) across 19 different cancer types, using antibody-based protein profiling data from immunohistochemistry (the Human Protein Atlas).


*Numbers indicate the tumor patients expressing each gene out of the total number (percentage, %).*

In TCGA-LIHC, only the individual high levels of PRF1 and GZMA were significantly associated with a positive effect on patient survival. A similar non-significant association of (individual or simultaneous) high GZMA and PRF1 expression with better effect on patient survival could also be observed in TCGA-MESO, ovarian cancer (GSE13876 and GSE49997), TCGA-STAD, TCGA-THCA, and TCGA-UCEC (Figure S1 in Supplementary Material). These data suggest that high CYT is widely associated with an improved prognosis among the above-mentioned cancer types.

On the contrary, across TCGA-LGG, BRCAs (GSE25066), and TCGA-THYM, both individual and simultaneous high levels of GZMA and PRF1 were significantly associated with a worse prognosis, whereas the simultaneous low levels of both genes led to a significant shift toward positive effect (**Figure 6B**). Regarding BRCA, though, we could not confirm these results using the independent datasets METABRIC and TCGA-BRCA, which revealed a tendency for the opposite effects of both cytolytic genes on patient survival. Regarding the METABRIC dataset, we separated BRCA patients who were subjected to hormonal therapy plus radiotherapy (HT/ RT) (*n* = 605) from the untreated patients; however, an association of high levels of GZMA and PRF1 with a worse prognosis could not be confirmed (Figure S6 in Supplementary Material).

Analogous non-significant associations of (individual or simultaneous) high cytolytic levels with worse effect on patient survival were also observed in lung cancer (GSE30219, TCGA-LUAD, and TCGA-LUSC), TCGA-PAAD, TCGA-PRAD and GSE16560, and TCGA-READ (Figure S2 in Supplementary Material).

In colon cancer, neither the individual nor the simultaneous high levels of the two genes were associated with a better prognosis, although simultaneous low levels of GZMA and PRF1 tended to shift toward a negative effect. Depending on the probe used, it seemed that a combination of high PRF1 and low GZMA levels yields a better patient outcome (GSE39582, TCGA-COAD, TCGA-COADREAD). Among metastatic colon cancer patients ("M1" patients in the TCGA-COAD dataset), simultaneous high levels of both genes were marginally significantly associated with worse prognosis, but simultaneous low levels of both genes could not provide the reverse trend (Figure S3 in Supplementary Material). We did not notice the same trend in the TCGA-COADREAD colorectal cancer patient cohort, though, implying that the aforementioned results are specific for colon (but not rectum) cancers.

Among clear-cell (TCGA-KIRC) and papillary renal cell carcinomas (TCGA-KIRP), we could not deduce any similar association among metastatic or non-metastatic tumors. In chromophobe renal carcinoma (TCGA-KICH) though, individual and simultaneous high levels of both genes tended to associate with better patient survival. On the other hand, concurrent low levels of both cytolytic genes, tended to associate with a worse prognosis. Interestingly, in the TCGA-KIPAN dataset, both the individual and synchronized high levels of GZMA and PRF1 significantly connected with worse patient survival. The simultaneous low expression of both genes exhibited reverse outcome (Figure S4 in Supplementary Material).

In DLBCL (GSE10846, and GSE32918), using various combinations of distinct molecular probes for the two cytolytic genes (PRF1, 214617\_AT, 1553681\_A\_AT, or ILMN\_1740633; GZMA, 205488\_AT, or ILMN\_1779324), we could not provide any

significant association with patient survival. A similar absence of significant associations was also detected in glioblastoma (GSE4271, GSE13041, and TCGA-GBM) and non-metastatic HNSCs. We could not deduce any further association or trend between the expression of both cytolytic genes and the survival of TCGA-TGCT and uterine carcinosarcoma patients (Figure S5 in Supplementary Material).

#### Infiltration of Lymphocytes and Neutrophils in Primary and Metastatic TCGA-Datasets

We further evaluated the infiltration of lymphocytes (TILs) and neutrophils (TANs) to the tumor site of primary and metastatic cancer samples across the TCGA-SKCM, TCGA-BRCA, and TCGA-THCA datasets, using the Cancer Digital Slide Archive (see text footnote 6). TILs contained both stromal- and intratumoralcompartment lymphocytes, as previously defined (43). Both of them were mainly composed of T cells and a smaller number of B cells, NK cells, and macrophages (44, 45).

In the TCGA-BRCA dataset, the number of TILs appeared enriched in the stroma of the primary tumors compared to the corresponding areas on the slide of the metastatic BRCAs. However, this might probably be due to the higher number of stroma cells detected in the primary breast tumors (percentage of stromal cells in primary vs metastatic BRCA, 21.15 ± 0.520 vs 7.143 ± 3.595; *p* = 0.032).

Although the number of TILs and neutrophils was higher in several cases of primary BRCA, the overall difference was not statistically significant (mean% of TILs ± SD in primary vs metastatic BRCAs, 6.102 ± 0.403 vs 4.714 ± 2.179; *p* = 0.78 and mean% of neutrophil infiltration ± SD in primary vs metastatic BRCAs, 1.625 ± 0.167 vs 0 ± 0; *p*= 0.44). Among primary tumors, comparing between triple negative (ER−, PR−, Her2/neu−, or TNBC), and triple positive (ER+, PR+, Her2/neu+, or TPBC) BRCAs, the load of TILs (and TANs) was not significantly different and was not significantly associated with a worse outcome, in argument with previous observations (46–48). In addition, the percentage of necrosis did not differ between metastatic and primary skin melanoma (<2%) (**Figure 7**).

In the TCGA-SKCM dataset, although in several cases the number of TILs was more enriched in the stroma of primary melanomas (as opposed to metastatic cancers), the overall load of TIL and TAN did not differ significantly between them. It is also worth noticing that the number of stroma cells counted in metastatic melanomas was higher compared to primary skin tumors (percentage of stromal cells in primary vs metastatic melanoma, 5.835 ± 1.083 vs 9.043 ± 0.571; *p* = 0.009). In addition, the rate of necrosis was marginally higher in metastatic skin melanoma

compared to primary tumors (*p* = 0.042). The overall survival did not differ between high TIL load (>1% TILs) or low TIL load (<1% TILs) in primary skin melanoma patients. However, among metastatic patients, a high percentage of lymphocytic infiltration shifted toward a better prognosis (**Figure 8**). According to recent data, the number of TILs in stage III metastatic melanoma associates with the response to Ipilimumab once these patients progress to stage IV disease (49).

In the TCGA-THCA dataset, the infiltration of lymphocytes was significantly higher in metastatic thyroid tumors and the high TIL load (>2% TILs) was associated with a better prognosis within the primary tumor group (mean% of TILs ± SD in primary vs metastatic cancers, 1.597 ± 0.160 vs 8.375 ± 3.59, *p* < 0.0001). The infiltration of neutrophils was minor (<0.1%) and did not differ between primary and metastatic THCAs. The necrotic rate was equally low between the two groups (**Figure 9**).

#### Correlation of the Cytolytic Index with Immune-Checkpoint Molecules and TILs in Primary and Metastatic TCGA-Datasets

In order to understand the context of PRF1/GZMA deregulation relative to the expression of various immune-checkpoint molecules, we correlated the cytolytic index with the expression of CTLA-4, PD-1, CD274 (PD-L1), PDCD1LG2 (PD-L2), LAG3, IDO1, CD73 (NT5E), and CD39 (ENTPD1) across all TCGA datasets (Figure S7 in Supplementary Material). In the majority of the cancers, a high cytolytic index was accompanied by upregulation of at least one immune-checkpoint molecule, indicating that similar to melanoma (42) and prostate cancer (41), immune response in CYT-high tumors elicits multiple host and tumor mechanisms of immune suppression in the tumor microenvironment (Figure S7 in Supplementary Material). For example, in TCGA-BRCA, CTLA-4, and PD-1 expression was significantly associated with a high cytolytic index (CTLA-4, *p* = 8.75e−199, Pearson's rho = 0.75; PD-1, 4.09e−309, Pearson's rho = 0.85). As expected, this correlation was 20 times stronger compared to the normal breast, due to absence of immunosuppression in the latter (CTLA-4, *p* = 2.44e−12, Pearson's rho = 0.60; PD-1, 1.12e−14, Pearson's rho = 0.65). Importantly, this association was even stronger among metastatic melanomas, suggesting the existence of a more intense immunosuppression in these tumors (e.g., in primary melanoma, PD-1, *p* = 1.14e−35, Pearson's rho = 0.887; LAG3, *p* = 1.03e−45, Pearson's rho = 0.930; IDO1, *p* = 3.05e−08, Pearson's rho = 0.513. In metastatic melanoma, PD-1, *p* = 1.48e−148, Pearson's rho = 0.917; LAG3, *p* = 4.28e−163, Pearson's rho = 0.931; IDO1, *p*= 2.38e−53, Pearson's rho = 0.690) (Figure S8 in Supplementary Material).

The cytolytic index was significantly correlated with lymphocyte infiltration in BRCA, thyroid cancer, and skin melanoma. The association between a high TIL load and CYT was stronger among primary breast and THCAs, but not in melanomas.

Consistent with the fact that apoptosis is a hallmark of CYT, we scored no further correlation between CYT and necrosis or between CYT and infiltration of neutrophils (**Figure 2B**).

#### Number of Tumor and Normal Cells across Metastatic and Non-Metastatic TCGA-Datasets

Expression analysis can be hampered due to a different number of cells within each tumor, thus reducing the ability to confidently measure the cytolytic index and correlate it with the expression of immune-checkpoint molecules in each dataset, as well as to compare gene expression between primary and metastatic cancers. To address this, we calculated the number of tumor and normal cells within each tumor. Overall, the primary and metastatic cancer samples across the three datasets contained an equal number of tumor cells (70–90% tumor cells, *p* > 0.05) (**Figures 7**–**9**). Thus, the detected differences should not be the result of enrichment in tumor cells in one group or the other. Similarly, the percentage in normal cells did not differ between primary and metastatic BRCAs (3.48 ± 0.26, in primary cancers vs 7.86 ± 3.56, in metastatic cancers; *p*= 0.188). On the other hand, primary melanomas had higher percentage of normal cells compared to metastatic tumors (7.485 ± 1.263 vs 1.155 ± 0.313, *p* < 0.001), and metastatic THCAs had a higher percentage of normal cells compared to their primary counterparts (8.125 ± 8.125 vs 2.126 ± 0.305, *p* = 0.021).

### DISCUSSION

In this study, we quantified the cytolytic index based on the expression of GZMA and PRF1, both of which mediate cytolysis. This index is strongly associated with CTLs, plasmacytoid dendritic cells, counter-regulatory Tregs, and known T-cell co-inhibitory receptors (31).

In agreement with Rooney et al. (31), we found a great variation in the immune CYT across different types of cancer, which possibly reflects the existence of merged tissue- and tumorspecific mechanisms orchestrating the local immunity. Cancers of the ovaries, liver, thyroid, esophagus, and prostate, as well as glioblastoma, glioma, pheochromocytoma and paraganglioma, adrenocortical carcinoma and uveal melanoma all exhibited minimal levels of CYT. On the contrary, DLBCL, clear-cell renal cell carcinoma, testicular cancer, cervical cancer, skin melanoma, and head and neck carcinoma exhibited increased levels of CYT.

The tumor-intrinsic resistance to CYT has been suggested to be due to different mechanisms. Among them, recurrent mutations in immune-related genes have been proposed, such as B2M, HLA-A, -B, and -C, and CASP8, as well as copy number aberrations in loci containing immunosuppressive factors, including

the receptors PD-L1/2 and CTLA-4 (31). PD-1 is transiently induced in activated T cells (50) and its expression is preserved in TILs (51–53). PD-L1 expression is high in several human malignancies, such as skin melanoma, lung, head and neck, and ovarian cancers (54, 55). PD-L1 expression is also correlated with a bad prognosis among patients with esophageal, colon, ovarian, or kidney cancer (56–60). The PD-1/PD-L1 axis is significant in tumor-induced immune evasion and both molecules are hopeful target candidates for immunotherapy. Actually, recent clinical trials have demonstrated that blockage of this signaling can benefit patients with advanced melanoma, kidney, or non-small cell lung cancer (2, 61–63). In metastatic melanoma, PD-L1 expression on peripheral T cells was recently shown to be prognostic on overall and progression-free survival (64).

CTLA-4 expression levels are low on resting T cells, but increase upon T-cell activation. In acute infection, CTLA-4 is transiently induced and binds to B7-1/2, thus competing with CD28 and weakening the T-cell response (65). On the other hand, CTLA-4 is constitutively expressed in T cells during chronic infection and cancer due to chronic antigen exposure. CTLA-4 is also constitutively expressed on antigen-experienced memory CD4+ and CD8+ T cells, as well as Tregs (65). Similarly, B7-1 is not expressed on resting antigen-presenting cells (APCs) (as opposed to B7-2) and is induced after APC activation. Anti-CTLA-4 therapy (Ipilimumab) was shown to induce cancer regression in metastatic kidney cancer (22, 66) and melanoma (67–70). Importantly, CTLA-4 blockade was reported to associate with bowel inflammation in melanoma patients (71), signifying that its signaling is crucial for the preservation of immune homeostasis in the gut.

Another example of immune-inhibitory molecule is indoleamine-pyrrole 2,3-dioxygenase (IDO). This molecule is constitutively expressed in the tumor microenvironment either by tumor cells or by host immune cells and is stimulated by inflammatory cytokines as IFN-γ, leading to host immune inhibition through increased Treg and effector T-cell proliferation blockade. A combination of IDO inhibition and immunecheckpoint blockade are currently under clinical investigation, with promising results (72).

Arginase is also an immune-inhibitory metabolic enzyme being expressed by both tumor cells as well as infiltrating myeloid cells (73). Both IDO and arginase inhibit immune responses by locally depleting the essential amino acids for anabolic functions in T cells or synthesizing specific natural ligands for cytosolic receptors, which can change the functions of lymphocytes. Inhibition of both IDO and arginase can enhance intratumoral inflammation (74, 75).

We found that high levels of several immune-checkpoint molecules, including CTLA-4, PD-1, PD-L1/2, LAG3, IDO1, CD73, and CD39 are associated with an increased cytolytic

index, across many cancers; and we expect that a combinatorial targeting of such immune-checkpoint molecules can provide a synergistic effect in cancer immunotherapy. Garg et al. found that predictive biomarkers of responsiveness to immune-checkpoint inhibitors in glioblastoma (GBM) exhibited inconsistent patterns among patients, predicting either resistance or susceptibility to therapeutic targeting of CTLA-4 or IDO1 (76).

Furthermore, different levels of tumor-intrinsic resistance to CYT can be attributed to the diverse levels of neoepitopes in these tumor types. Neoepitopes are tumor-specific antigens produced from DNA mutations occurring in cancer cells. Such mutations can be missense mutations, indels (insertions/deletions), and/or gene fusions. Increasing evidence shows that neoepitope-specific antitumor immune responses occur naturally in cancer cells and have great potentials as immunotherapeutic agents (77). Theoretically, immune responses to neoepitopes are not diminished by host central tolerance in the thymus and cannot trigger an autoimmune reaction (77, 78). These neoepitopes were lately shown to facilitate recognition of a tumor as foreign (78, 79), and an increased load of them is associated with effective immune responses to immune-checkpoint therapy (80). Currently, strategies to selectively enhance T-cell reactivity against genetically defined neoepitopes are under development (78, 81–84). Furthermore, recent findings identified target neoepitopes which can be helpful in the design of a vaccine against murine melanoma (85). Importantly, the immunogenicity and specificity of these neoepitopes was validated *in vivo*, after administering mice either mutated or wild type synthetic peptides. Further advance in the field was made by Verdegaal et al. who analyzed the stability of neoantigen-specific T-cell responses and the antigens they recognize in melanoma patients treated by adoptive T-cell transfer. This study demonstrated that T cells mediate neoantigen immunoediting, indicating that the therapeutic induction of broad neoantigen-specific T-cell responses should be used to avoid tumor resistance (86).

In comparison to melanoma, the immune CYT in breast cancer, the burden of nonsynonymous mutations, and the predicted load of neoepitopes were previously found to be relatively modest, suggesting that a combination of immune agents with nonredundant mechanisms of action should be of high-priority (87). Recently, Vonderheide et al. highlight the critical steps that need to be followed for a more successful immunotherapy

part); overall patient survival with respect to the percentage of tumor-infiltrating lymphocytes (TILs) (>5%, high TILs, <5%, low TILs); and representative hematoxylin and eosin slides of ER+, PR+, Her2/neu+ (TPBC) and triple negative breast cancer (lower part).

in breast cancer, including immune suppression in the tumor microenvironment and failed or suboptimal T-cell priming (87).

Also Chen et al. (88) categorized the tumor microenvironment into four types, depending on the expression of PD-L1, as well as the ratio CD8A/CYT, and proposed that this classification can serve the design of more suitable immunotherapeutic strategies.

A very interesting improvement in the field was further made by Riaz et al. (89), who showed that the mutation burden in melanoma patients decreases with successful anti-PD-1 blockade therapy, suggesting that the selection against mutant neoepitopes is a critical mechanism of action of this immunotherapy.

All these advances, show that neoepitopes can be used as biomarkers to predict the clinical response to immunotherapy and the outcome, as well as to serve as immunotherapy targets (25, 90). Besides epitope selection, the reduction of gene expression heterogeneity within tumor cells, the definition of the optimum number of simultaneously targeted neoantigens, of the patient profile that can benefit from neoantigen-based immunotherapy and escape the risk of adverse effects, and a synergistic combination of immune-checkpoint blockade and/or adoptive T-cell therapy, are all issues that need to be successfully addressed in order to select potent neoantigens for cancer immunotherapy (77).

Adding to the variability in cytolytic levels that we detected among different cancer types, CYT has also been previously shown to correlate with oncogenic viruses in certain tumor types. For example, CYT is associated with HPV infection in cervical cancer, and head and neck cancer, with EBV infection in stomach cancer, and with HBV and HCV infection in liver cancer (31). Overall, it seems that CYT is part of an inflammatory environment in a premalignant state of certain tumor types, whereas, in others, oncogenic mutations, copy number aberrations, or viral infection can induce a tumor-promoting inflammatory microenvironment, within which complex interactions between different cell types regulate cancer development and metastasis (91, 92).

In the context of metastasis, we observed that CYT was significantly higher in metastatic skin melanoma compared to primary skin tumors. The increased cytolytic levels could be

further observed in metastatic breast and thyroid cancers suggesting that although initially regarded as an indicator of a failed immune response, CTLs/NK cells (among other inflammatory cells) also support tumor development (93, 94). This observation is in agreement with previous reports supporting that regardless of the tumor's origin, an inflamed tumor microenvironment has many tumor-promoting effects (91, 92). In line with this, we found significantly elevated expression of various suppressive factors, correlating with a high cytolytic index, in metastatic skin melanomas, breast, and thyroid cancers. For example, high levels of CTLA-4, PD-1, PD-L1/2, LAG3, and IDO1 that we detected in metastatic melanoma, were many fold times more significantly associated with high cytolytic levels, pointing towards the existence of immunosuppression in these metastatic tumors (Figure S8 in Supplementary Material).

The above-mentioned vast range in the CYT and the different levels of infiltration of inflammatory cells (T cells and neutrophils) is best reflected by the different survival curves produced among different types of cancer. In some tumor types (ACC, SKCM, BLCA, LIHC, MESO, OV, STAD, THCA, and UCEC), high CYT was associated with an improved outcome; whereas in others (LGG, BRCA, THYM, LUAD/LUSC, PAAD, PRAD, and READ) it is correlated with a worse outcome. Among LGG, THYM, and BRCA, we showed that both individual and simultaneous high levels of GZMA and PRF1 were significantly associated with a worse prognosis, whereas the simultaneous low levels of both cytolytic genes led to a significant shift toward a positive effect. Nevertheless, we could not observe this across different breast cancer datasets. Furthermore, contrasting results mentioning a worse effect of PRF1 on survival of BRCA patients were also recently reported in another large-scale meta-analysis (95). The difference between the two studies might be due to cohort-specific bias or power-related discrepancies. Of interest, among certain tumor types including ACC, SKCM and BRCA, the simultaneous expression of both cytolytic genes synergistically affected patient survival.

Tumor-infiltrating lymphocytes are mononuclear cells of the immune system that intrude the tumor tissue, and their presence

carcinomas (lower part).

has been reported in solid tumors, such as breast, colon, lung, and cervical cancers, as well as in melanoma (43, 96–98). Low levels of CD8+ TILs are related with the likelihood of response and may escalate during therapy in responding tumors (2, 99). Further, the location of CD8+ TILs at the invasive margin of tumors may indicate an effective immune response (42, 99, 100). The tumor microenvironment may limit extravasation of effector T cells into the tumor, diminish their expansion, or reduce their viability (101). In BRCA, an increased TIL load in the stroma of the tumor was reported to associate with a higher prospect of therapy in early stage TNBC and Her2+ patients (98). Assessing light microscopy data of tissue slides, we found a higher TIL load in primary BRCA compared to the metastatic counterparts, but the differences were not significant (*p* > 0.05). These lymphocytic infiltrates mirror favorable host antitumor immune responses within these samples. Although the presence of high TIL levels has been previously linked with a more favorable prognosis in patients with Her2+ and early stage TNBC (46–48), we found no significant difference in the outcome of TNBC or TPBC between high and low TIL load.

Tumor-associated neutrophils also compose a significant part of the inflammatory cell infiltrate in several tumor types (102–105), but the mechanisms by which they affect tumor progression are only now being investigated. Recent studies point toward the tumor-promoting effects of neutrophils. Histologic studies performed on a variety of tumor types have shown that the increased TAN load correlates with unfavorable recurrence-free, cancer-specific and overall patient survival in kidney cancer, skin melanoma, colorectal cancer, and head and neck cancer (106). It has also been suggested that TANs can drive the metastasis of breast cancer cells to the liver or the lung (107, 108), activating angiogenesis (109, 110). In contrast, older reports suggested that neutrophils have antitumoral effects, by inducing direct cytotoxicity of target cells, and decreasing the size and the number of lung metastatic foci (111–113). Interestingly, the anticancer activity of TANs was reported to mostly culminate into anticancer activity *via* oxidative burst (114, 115). Despite the heavily debated role in favor or against cancer, the latest research shows that TANs do play a key role in various aspects of tumor development, from malignant transformation to tumor progression, modification of the extracellular matrix, angiogenesis, cell migration, and immunosuppression (116–121). Due to their contradictory roles in cancer, neutrophils are now classified into two subpopulations, antitumor and pro-tumor TANs (117). We detected very low percentage of TANs in primary breast and thyroid carcinomas and almost null levels in their metastatic counterparts. We noticed higher neutrophilic infiltration in TNBC compared to TPBC, but without reaching statistical significance (*p* > 005). In skin melanoma, we noticed even less neutrophilic infiltration, being slightly higher in the metastatic tumors.

Overall, the multiple crosstalks among different tumor-infiltrating immune cells, including TILs and TANs, was suggested to promote or inhibit the establishment of a permissive tumor microenvironment (17). A better understanding of the role of these cells will provide opportunities for the immunomodulation and the improvement of the existing antitumor therapies.

To conclude, we have measured the CYT in terms of RNA and protein levels in a large number of TCGA datasets, in order to understand how different cancers induce and adapt to immune responses. We associated each cancer's CYT with patient survival both in primary and metastatic cases and evaluated the tumorinfiltration of lymphocytes and neutrophils in H&E-stained sections of the same tumors. Our data suggest that the cytolytic index along with the existence of complicated associations among various tumor-infiltrated immune cells is capable to promote evasion from immunosurveillance in certain cancers.

#### REFERENCES


#### AUTHOR CONTRIBUTIONS

CR extracted data, analyzed them, and performed statistical computations. DC, AM, and CE extracted data. CD critically read the manuscript and approved its final version. AZ supervised the study, extracted, and analyzed data, interpreted the results, and wrote the manuscript.

#### ACKNOWLEDGMENTS

All authors would like to thank the staff members of the TCGA, GTEx and HPA platforms and all the corresponding patients for the data of whom they extracted and analyzed.

#### FUNDING

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fonc.2018.00027/ full#supplementary-material.


translational blockade of immune checkpoints. *Int J Mol Sci* (2016) 17:E1151. doi:10.3390/ijms17071151


inflammatory activity, modulation by cancer cells and expansion in advanced disease. *Int J Cancer* (2011) 129:2183–93. doi:10.1002/ijc.25892


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Roufas, Chasiotis, Makris, Efstathiades, Dimopoulos and Zaravinos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*