Nuclear Factor of Activated T-cell Activity Is Associated with Metastatic Capacity in Colon Cancer

Colorectal carcinoma is a leading cause of cancer-related deaths in the United States (1, 2). Conventional staging is insufficient to adequately predict metastasis or recurrence in colorectal cancer, particularly for patients with stage II or stage III disease (3). Over the past several years, gene expression profiling has been used to predict patient outcomes independently of conventional staging; however, the biologic mechanisms driving poor clinical outcomes remain incompletely understood. Furthermore, because published prognostic colorectal cancer gene expression signatures are nonoverlapping, inference of regulatory molecular targets driving clinical outcomes poses a significant challenge (4, 5). Nevertheless, it is likely that the combined data underlying these studies provide enough biologic information for inference of potential transcriptional regulators of colorectal cancer progression and metastasis.

In the present study, we reasoned that the gene expression signatures of metastasis in mouse models and poor prognosis in humans reflect a functional consequence of the deregulation of coregulated gene networks. Accordingly, we integrated eleven human microarray datasets, representing data from 1,295 colorectal cancer tumor specimens, to produce a set of coexpressed gene modules. These coexpression modules were filtered using gene expression data from an immunocompetent mouse model of metastatic colon cancer to identify a metastasis-associated coexpression module. We further predicted the nuclear factor of activated T cell (NFAT) family of transcription factors as a primary regulator of 21 genes in the module.

NFAT is a cell signaling molecule involved in complex adaptive systems particular to vertebrate biology related to the Rel family of cell signal proteins (6, 7). NFATs are regulated by receptor-mediated calcium/calmodulin signaling driving nuclear transport where they control transcriptional events governing functions as diverse as cell proliferation, survival and differentiation (8, 9), epithelial-to-mesenchymal transition (EMT; refs. 6, 10, 11), and stem cell quiescence (12). NFATs also modulate inflammatory cell function, angiogenesis, cell migration and invasion (11, 13), and regulation of tumor stromal cell function in both a cell-type–specific and context-dependent manner (14).

Here, we used a systems biology approach to identify NFAT as a potential regulator of colon cancer progression and metastasis. A schematic of our overall workflow appears in Supplementary Fig. S1. We experimentally determined that expression of NFATc1, specifically, promotes increased colon cancer cell invasion and metastasis, and that expression of the NFATc1-driven transcriptional program correlates with poor prognosis in stage II and III colorectal cancer patients.

Human gene expression datasets were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) and the ArrayExpress Archive (http://www.ebi.ac.uk/microarray-as/ae/; Supplementary Table S1). Each gene expression dataset was processed separately to generate a gene by sample expression matrix with normalized, log-transformed, and standardized expression values (see Supplementary Methods). The mouse gene expression dataset was downloaded from the GEO GSE19073 (15) and the differential expression analysis was described in detail in Supplementary Methods.

For each human CRC dataset selected for integration, Pearsons correlation coefficients and the log likelihood ratio (LLR) for functional relevance were calculated for all gene pairs. Next, an integrated LLR score was calculated for each gene pair, and a score cutoff was selected to achieve a balance between consistency with existing knowledge and network coverage. An integrated gene coexpression network was then constructed in which each node is a gene while two nodes are connected by an edge if corresponding LLR score is above the selected cutoff (for additional detail, see Supplementary Methods). For the identification of coexpression modules, we used a modified version of the Iterative Clique Enumeration (ICE) algorithm (16). GO biologic process enrichment analysis of the coexpression modules and transcription factor target enrichment analysis were performed using WebGestalt (17). The ranked list of genes according to their differential expression levels in the mouse dataset was used to identify metastasis-related modules using the GseaPreranked analysis available in the GSEA software downloaded from http://www.broadinstitute.org/gsea/. Default parameters were used for the analysis.

Survival analysis was performed in R using the survival package. Specifically, survival curves were estimated using the Kaplan–Meier method, and survival comparisons among groups were made by the log-rank test. Two datasets were used for the survival analysis. The dataset GSE17536 was from the Moffitt Cancer Center (MCC; Tampa, FL) and has been previously published. The dataset GSE38832 was from the Vanderbilt University Medical Center (VUMC; Nashville, TN), which included 54 samples from the previously published dataset GSE17537 and 68 newly analyzed samples. The COMBAT approach based on empirical Bayes frameworks was applied to remove batch effects across the two batches of microarray experiments (18).

Human tissues used for microarray analysis were collected following written informed consent and clinically annotated according to established protocols and approved by the appropriate Institutional Review Boards at the Moffitt Cancer Center and Vanderbilt University as previously described (GSE17536 and GSE17537; ref. 15).

Total RNA from cells or tissues was isolated using QIAGEN kits (QIAGEN) and DNase-I treated, quantified by Nanodrop-1000 (Thermo Scientific) and assessed for quality on an Agilent Bioanalyzer as previously described (15). Chromosome immunoprecipitation (ChIP) studies were conducted using mouse anti-NFATc1–specific antibody from Santa Cruz Biotechnology and a kit from Millipore, according to the manufacturer's instructions. qRT-PCR (qPCR) was performed as described elsewhere (19), gene-specific primers are listed in Supplementary Tables S2 and S3.

The MC-38 mouse adenocarcinoma cell line and its derivatives were provided by Dr. Robert Coffey (Department of Internal Medicine, Vanderbilt University Medical Center, Nashville, TN) are described elsewhere (15). HCT116 and HT29 colon cancer cell lines were obtained from ATCC. All cell lines were maintained at low passage as monolayers in RPMI-1640 medium (Gibco Life Technologies), supplemented with 10% FBS (Atlanta Biologicals), 500 U/mL penicillin G, 500 μg/mL streptomycin (Gibco Life Technologies Inc.), and l-glutamine (Gibco Life Technologies Inc.). FK506 (Sigma) was used at 20 ng/mL. Integrity of human cell lines used in this study was tested by RNAseq analysis in May 2013. Cytoplasmic and nuclear extracts were prepared using Nuclear Extract kit (Active Motif), according to manufacturer's instructions.

Cells were harvested in RIPA lysis buffer (50 mmol/L Tris, pH7.5, 150 mmol/L NaCl, 1% NP-40, 0.5% Na-deoxy Cholate, 0.1% SDS) containing a cocktail of protease inhibitors (Roche Diagnostics), with a brief sonication. Samples were mixed with LDS buffer containing DTT (Invitrogen), and fractionated on 4% to 12% NuPAGE gels in MOPS-SDS buffer (Invitrogen). Antibodies against NFATc1-c4 and PARP1/2 were obtained from Santa Cruz Biotechnology, β-actin from Sigma Chemical, and α-tubulin from Abcam Scientific. Ramos cell (Burkitt lymphoma, B lymphocytes) lysate (Santa Cruz Biotechnology) was used as positive control.

Invasion assays were conducted using both Boyden chambers as described elsewhere (15) as well as the xCELLigence system from Roche Diagnostics (Supplementary Methods).

RNAi studies were performed as previously described using NFATc1-specific ON-target plus SMART pool siRNA or ON-target plus Non-targeting Pool (Thermo Scientific). Specific si-RNA sequences are given in Supplementary Table S4.

For preparing stable overexpressing cells, mouse wild-type NFATc1 (Addgene plasmid 11101; ref. 20) was cloned into vector pGP-Lenti3 (GenBank Accession no. JX861384) between unique XbaI and BamHI sites. Puromycin and GFP expressions were used as selection markers for creation of stable cell lines.

For gene knockdown experiments, NFATc1-specific GIPZ lentiviral shRNAmir (V3LMM_418820; Thermo Scientific) was selected for experimental use. Briefly, MC38 cells were electroporated using NEON (Invitrogen). Puromycin and GFP expressions were used as selection markers for creation of stable cell lines. For HT29 studies, cells were transiently cotransfected with pEF6-NFATc1 and pMaxGFP vectors, flow sorted to enrich for highly expressing GFP-positive cells, and confirmed NFATc1 overexpression by Western blot analysis.

All animal studies were approved by the Vanderbilt Institutional Animal Care and Use Committee and performed in accordance with the standards of the Association of Assessment and Accreditation of Laboratory Care. Procedures were performed using both parental and stable transfected cell lines as previously described (15).

Statistical analyses for all experimental procedures were conducted in GraphPad Prism 5 (GraphPad Software). P < 0.05 was considered statistically significant.

We first identified 14 publicly available human colorectal cancer gene expression datasets (Supplementary Table S1) and eliminated three sets, for which the association between pairwise gene coexpression and functional relevance was not clear (Supplementary Results and Supplementary Fig. S2). We also found that the 11 remaining datasets contain complementary functional information (Supplementary Results, Supplementary Fig. S3) and were therefore suitable for integration and development of a comprehensive coexpression network.

Following Ramani and colleagues (21), we constructed three coexpression networks based on three LLR cutoffs and evaluated the functional relevance of the networks using random gene pairs and interacting protein pairs (Supplementary Methods) as negative and positive references, respectively. As shown in Fig. 1A, all three networks (cyan, pink, and green curves) clearly outperformed the negative reference (black), and a higher LLR cutoff corresponded to higher functional relevance of gene pairs in the coexpression network. The functional relevance of the coexpression network derived from an LLR cutoff of 1.25 (green) approached that for the curated protein interaction network (red), suggesting a high biologic validity of the coexpression network. To capture novel functional relationships offered by the coexpression network analysis that may not be supported by existing knowledge, we used a slightly relaxed LLR cutoff of 1 and constructed the integrated colorectal cancer coexpression network with 2,285 genes and 13,083 edges (Supplementary File S1). A coexpression network derived from the GSE17536 dataset with an LLR cutoff of 0.74 showed a similar level of functional relevance as that of the integrated network (Fig. 1A, dashed pink curve). However, it included only 1,226 genes and 7,516 edges. Therefore, data integration helped achieve a significantly better network coverage with 86% more genes and 74% more edges. The integrated network represents one of the most comprehensive coexpression networks for a specific cancer type and can serve as a resource for other investigators for the study of cancer-associated gene function and identification of candidate cancer driver gene programs.

The integrated network had a high clustering coefficient of 0.58, indicating a highly modular network structure (22). Using a modified version of the ICE algorithm (16), we identified 441 coexpression modules (Supplementary Fig. S4 and Supplementary File S2) from the network, among which 38 included 10 or more genes.

Figure 1B depicts the five largest modules with 187, 107, 105, 85, and 72 genes, respectively. Four of the top five modules showed significant enrichment (FDR <0.01, hypergeometric test) in distinct biologic processes including developmental process, cell cycle, translational elongation, and immune system process, respectively. The fifth module was not enriched in any of the defined GO biologic processes (labeled as unknown function). Interestingly, among the 72 genes in this module, only 27 (37%) had a GO biologic process annotation, in contrast with an average of 85% annotated genes in the other four modules. It is possible that this module may represent a transcriptional program underlying an understudied biologic function.

We ranked all human genes based on the differential expression level of their mouse orthologs in the poorly metastatic parental MC38 cell line (designated as MC38^Par) and the highly metastatic MC38 cell-derived colon cancer cell line (MC38^Met; ref. 15). For each coexpression module, we used the GSEA to test whether genes in the module were significantly enriched in the top or bottom of the ranked list, indicating coordinate up- or downregulation of the module in the metastatic cell line.

The module consisting of 187 genes functionally linked to developmental process (Fig. 1B) was enriched with genes upregulated in the metastatic cells (FDR < 0.001), whereas the module consisting of 107 genes functionally linked to cell cycle was enriched with genes that are downregulated in the metastatic cells (FDR = 0.01). This pattern is consistent with reduced cell proliferation and increased motility in cells undergoing EMT, a developmental process implicated in cancer metastasis (23). Owing to its significant association to the metastatic cell gene expression profile, we focused on the 187-gene module for further study. We identified 63 core genes in this module that accounted for the enrichment signal in the metastatic cells (Fig. 1C; Supplementary File S3). As shown in Fig. 1D, mouse orthologs of the 63 core genes showed clear upregulation in MC38^Met compared with MC38^Par, whereas those of non-core genes were either not differentially expressed or downregulated in MC38^Met. This result indicates that a substantial portion, but not all, of the coexpression relationship was conserved between human tumors and the mouse model. The 63 core genes are positively correlated in human and mouse, and thus, serve as candidate effectors of a transcriptional program associated with a metastatic phenotype.

We performed enrichment analysis for the 63 core genes against the gene sets of transcription factor targets and found that one third of the genes in the module were putative targets of the NFAT family transcription factors, representing significant enrichment (FDR = 0.0008, Fisher exact test, Fig. 1E and Supplementary File S4). As a set, these genes are functionally enriched for skeletal and muscular system development, connective tissue development, and for cell movement, growth, and proliferation (Supplementary Table S5).

We next assessed the mRNA expression of NFAT family members and 21 putative target genes in 22 stage II and III fresh-frozen primary human colorectal cancer tumor specimens. We found that mRNA for 19 out of the 21 putative NFAT target genes were significantly upregulated in the tumor specimens as compared with normal adjacent tissue (Supplementary Fig. S5). We next assessed the expression of individual NFAT members and the 21 target genes using microarray data from the isogenic mouse cell lines and from primary tumor data. Among the four NFAT members examined, NFATc1 mRNA was the most differentially expressed (FDR = 0.0001) between MC38^Met and MC38^Par cells, with almost 3-fold (P < 3E-06, moderated t test) increase in metastatic cells (Supplementary Table S6). Moreover, two NFATc1 probe sets (21105_s_at and 210162_s_at) correlated best with the major cluster of the 21 metastasis-related NFAT transcriptional target genes in three independent patient cohorts (Supplementary Fig. S6B–S6D). These data suggested that NFATc1, specifically, is associated with a metastasis-associated transcriptional program in human colorectal cancer.

We found that pharmacologic inhibition of NFATs using FK506 reduced both cell invasiveness and expression of metastasis-associated target genes (Supplementary Fig. S7). We next wanted to assess the role of NFATc1 expression, specifically, in regulating invasiveness. We found that experimental inhibition of NFATc1 expression in MC38^Met cells by RNAi resulted in a specific decrease in NFATc1 mRNA and protein expression. In Fig. 2, four members of the NFAT family are evaluated for off-target effects of the RNAi at the mRNA and protein levels. First, the effect of siNFATc1 on NFATc1-specific mRNA levels is evaluated in MC38^Met cells relative to siSCR control and parental lines (Fig. 2A). Next, the effect of the RNAis are evaluated on NFATc2 (Fig. 2B) and NFATc3 (Fig. 2C)-specific mRNA expression, demonstrating no off-target effects relative to parental siSCR controls. In Fig. 2D, the specificity of the NFATc1-specific knockdown can be appreciated at the protein level, as well, with no off-target effects shown for NFAT family members c2-c4. Importantly, treatment of MC38^Met with NFATc1-specific RNAi caused a significant reduction in cell invasion of MC38^Met cells (Fig. 2E). Thus, the siNFATc1-specific RNAi (a pool of four distinct RNAis was used in this experiment, sequences are given in Supplementary Table S4) has no detectable off-target effects on NFAT family member-specific mRNA species or proteins, demonstrating the specificity of the observed effects on cell invasion.

To determine the relevance of NFATc1 expression in human cancer cells, we manipulated expression levels in HT29 (low to undetectable NFATc1 expression) and in HCT116 cells (high endogenous NFATc1 expression). We found that transient overexpression of NFATc1 in HT29 cells (Fig. 3A) had no effect on cell viability or proliferation (Fig. 3B) but significantly increased cell invasion (Fig. 3C). In contrast, when RNAi was used to specifically inhibit NFATc1 expression in HCT116 cells (Fig. 3D and E), we found that targeted RNAi treatment did not significantly affect cell proliferation (Fig. 3F) and was associated with a significant decrease in the rate of invasion through Matrigel as compared with the siSCR control (Fig. 3G).

NFATc1 activity is primarily regulated at the posttranslational level by phosphorylation, so variations in NFATc1 expression levels across human colon cancer specimens, particularly in the epithelial compartment, are low (data not shown). However, the range of expression of the putative targets of NFATc1 activity identified with our filtered integrated data network is diverse and so may be used as surrogate biomarkers for NFATc1 activity in human specimens. Combining the results of experiments conducted in MC38 cells, we found that eight out of 21 genes from the metastatic module are consistently altered by NFAT activity and NFATc1 expression in MC38^Met cells (Fig. 4A and B, Supplementary Fig. S7 and Supplementary Table S5). These eight genes include Angiopoeitin (ANGPTL2), Asporin (ASPN), Collagen 3A1 (COL3A1), Matrix Gla Protein (MGP), Mannose receptor (MRC2), Fibroblast activating protein (FAP), Polymerase I and Transfer Release Factor (PTRF), and TWIST1. ChIP assays confirmed enrichment of NFATc1 on promoters on seven out of eight of these genes in metastatic MC38 cells (Fig. 4C, Supplementary File S5). This 8-gene set is functionally enriched in cell movement (ANGPTL2, COL3A1, FAP, MGP, MRC2, TWIST1) and cancer (TWIST1, COL3A1, ANGPTL2, ASPN, FAP, MGP, MRC2, PTRF). Importantly, we also observed significant overlap between these eight genes and previously published colon cancer-specific prognostic signatures (15, 24–27), indicating that our integrative approach largely supports and builds upon a significant group of published studies in this area by providing a biologically based underpinning for the signature (Supplementary Table S7). Although developmental pathways and activation of cell movement processes have been strongly linked with cancer progression in other experimental models (15, 25–29), these findings are the first to link them with cell autonomous NFATc1 expression and invasive activity in tumor cells.

To test whether these eight genes can be used to identify colorectal cancers with poor prognosis, we examined their expression levels in two clinically annotated cohorts of stage II and III colorectal cancer tumors. As shown in Fig. 5A and B, increased expression of these genes is associated with increased number of recurrence and death events for stage II and III patients. To statistically evaluate the correlation of the eight genes with colorectal cancer outcomes, we performed a survival analysis based on a combined VUMC and MCC dataset. On the basis of the average of standardized gene expression levels of the eight genes, patients were classified into two groups (above median expression and below median expression). We found that the group of patients with above-median expression had significantly shorter overall survival (OS: P = 3.2e-05, HR = 3.56; Fig. 5C), worse disease-specific survival (DSS; P = 6.5E-06, HR = 7.89, Fig. 5D), and shorter disease-free survival (DFS; P = 1.1E-04, HR = 3.76, Fig. 5E) as compared with the below-median group. These data show that increased expression of the NFATc1-driven transcriptional program correlates with worse outcomes for stage II and III patients and suggest a novel role for NFATc1 in promoting tumor progression.

To determine whether the NFATc1-associated changes in cell migration and invasion observed in cell culture are also associated with changes in metastatic potential in vivo, we manipulated NFATc1 expression in MC38 cells and tested their metastatic potential in isogenic and immunocompetent mice as described previously (15). For these studies, we performed stable transfection of MC38^Parcells with constitutively expressed NFATc1 (MC38^Par+NFATc1; Fig. 6A and B) and showed that these cells exhibited a higher rate of trans-endothelial invasion as compared with the MC38^Par+Vec control cells (Fig. 6C). For the in vivo experiments, mice injected intrasplenically with MC38^Par+NFATc1 cells exhibited a greater metastatic liver tumor incidence and overall tumor burden as compared with those injected with MC38^Par+Vec control cells. Strikingly, liver metastasis from MC38^Par+NFATc1 cells were mostly diffuse and replaced most of the liver parenchyma in contrast with very little metastatic colonization in mice injected with MC38^Par+Vec control cells (n = 12-14 mice per cell line, Fig. 6D–G and Supplementary Fig. S8C and S8D). Moreover, target gene expression was predictably upregulated in these cells, consistent with a regulatory role for this network in metastasis in vivo (Supplementary Fig. S8E). Thus, overexpression of NFATc1 and the target gene network in MC38^Par cells increased their metastatic potential in vivo.

In complementary studies, stable NFATc1-specific shRNAmir expressing MC38^Met cells (MC38^Met+shNFATc1) or control scrambled shRNAmir (MC38^Met+shCtrl) cells were generated (Fig. 7A and B). We found that MC38^Met+shNFATc1 cells exhibited a lower rate of invasion as compared with the MC38^Met+shCtrl cells (Fig. 7C). Intrasplenic injection of the MC38^Met+shNFATc1 cells resulted in reduced metastatic liver tumor burden as measured by bioluminescence (Fig. 7D and E; P < 0.04), lower tumor incidence (Fig. 7F; P < 0.03), and reduced liver weight to body weight ratio (Fig. 7G; P < 0.0004) as compared with mice injected with MC38^Met+shCtrl cells (Supplementary Fig. S8A–S8C). Moreover, target gene expression was predictably downregulated in these cells and consistent with a role for this NFATc1 network in regulating metastasis in vivo. Supplementary Table S8 summarizes the statistics of the in vivo experiments in table format.

Thus, the above-described in vivo studies confirm a role for tumor cell expression of NFATc1 in regulating the metastatic behavior of colon cancer cells in an immunocompetent mouse model, strengthening the evidence that such a mechanism can also play a role in human colorectal cancer. However, our results neither rule out a cooperative role of other NFAT members in disease progression, as has been shown by others in pancreatic cancer (30–32), leukemia (33, 34), breast cancer (35, 36), and in melanoma (37), nor they rule out cooperative effects of other signaling pathways in contributing to these effects. Further studies are needed to determine the precise regulation of genes from this module and to determine whether or not non-cell autonomous signals may be required for full expression of the prognostic genes.

In summary, we combined multiple publicly available microarray datasets to identify a novel cell autonomous role for NFATc1 in regulating cell invasion and metastasis in colon cancer and to further link that activity to expression of eight target genes whose expression is associated with clinical outcomes in patient samples. A unique strength of this study is the quality evaluation of publicly available human microarray datasets and systematic integration of the data into a coexpression network for module analysis. Although the coexpression modules were identified in an unsupervised manner, filtering the modules using expression data from a mouse model allowed us to identify a group of coregulated genes whose expression is associated with poor outcomes in human colorectal cancer. Moreover, NFATc1, the regulator of a group of coexpressed genes associated with poor colorectal cancer patient outcomes, was demonstrated to promote metastasis in mice. Thus, the integrative approach addresses the sample size limitations of earlier prognostic gene signature studies and the use of the immunocompetent mouse model represents important aspects of a disease-appropriate microenvironment in the metastatic process.

The authors disclose no potential conflicts of interest.

Conception and design: N. Deane, H. An, S. Mima, K. Ciombor, B. Zhang, R.D. Beauchamp

Development of methodology: M.K. Tripathi, N. Deane, H. An, K. Ciombor, B. Zhang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.K. Tripathi, N. Deane, H. An, C. Padmanabhan, K. Ciombor, N. Prodduturi, M.K. Washington

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.K. Tripathi, N. Deane, J. Zhu, H. An, X. Wang, Z. Shi, X. Chen, M.K. Washington, B. Zhang, R.D. Beauchamp

Writing, review, and/or revision of the manuscript: M.K. Tripathi, N. Deane, J. Zhu, H. An, K. Ciombor, M.K. Washington, B. Zhang, R.D. Beauchamp

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.K. Tripathi, N. Deane, H. An, K. Ciombor, B. Zhang

Study supervision: N. Deane, B. Zhang, R.D. Beauchamp

The authors thank Tanner J. Freeman for critical input and Connie Weaver and Keeli Lewis for technical expertise.

This work was supported by grants GM088822 (B. Zhang, R.D. Beauchamp), CA095103 (R.D. Beauchamp), CA158472 (X. Chen and R.D. Beauchamp), K12 CA9060625 (X. Chen), the Vanderbilt-Ingram Cancer Center Support Grant (CA0684845), the Clinical and Translational Science Award to Vanderbilt (UL1RR024975), Vanderbilt Digestive Disease Research Center (DK058404), and the Vanderbilt Cooperative Human Tissue Network U01CA094664 and NIH.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.