Multiscale network analysis reveals molecular mechanisms and key regulators of the tumor microenvironment in gastric cancer

Gastric cancer (GC) is the third leading cause of cancer deaths and the fourth most prevalent malignancy worldwide. The high incidence and mortality rates of gastric cancer result from multiple factors such as ineffective screening, diagnosis, and limited treatment options. In our study, we sought to systematically identify predictive molecular networks and key regulators to elucidate complex interacting signaling pathways in GC. We performed an integrative network analysis of the transcriptomic data in The Cancer Genome Atlas (TCGA) gastric cancer cohort and then comprehensively characterized the predictive subnetworks and key regulators by the matched genetic and epigenetic data. We identified 221 gene subnetworks (modules) in GC. The most prognostic subnetworks captured multiple aspects of the tumor microenvironment in GC involving interactions among stromal, epithelial and immune cells. We revealed the genetic and epigenetic underpinnings of those subnetworks and their key transcriptional regulators. We computationally predicted and experimentally validated specific mechanisms of anticancer effects of GKN2 in gastric cancer proliferation and invasion in vitro. The network models and the key regulators of the tumor microenvironment in GC identified here pave a way for developing novel therapeutic strategies for GC.


Introduction
Gastric cancer (GC) is the third leading cause of cancer death and the fourth most prevalent malignancy worldwide, accounting for 8% of cancer incidence and 10% of cancer deaths, and approximately 21,320 cases of GC (13,020 men and 8,300 women) were diagnosed and 10,540 patients died from GC in 2012 in the United States. 1 The molecular mechanisms driving tumorigenesis of GC include a number of biological and cellular processes activated in tumor pathogenesis such as proliferation, angiogenesis, the bypass of senescence and cell death pathways, evasion of tumor suppressing mechanisms, immortality, invasion and so on. 2 leading inconsistent treatment responses and marginal improvements. 3 High-throughput molecular profiling data makes it possible to dissect the heterogeneity of GC in a comprehensive and unbiased manner. 3 A number of gene signatures have been identified for diagnosis and classification of GC as well as prediction of therapeutic response. 3 However, the reproducibility of such gene signature-based models are usually poor due to multiple factors including limited cross validation of predictive gene lists per tumor type and outcome. 4 Recently, The Cancer Genome Atlas (TCGA), a comprehensive multi-Omics cohort for studying multiple cancers was developed. TCGA includes genomic, transcriptomic and epigenomic molecular data of primary gastric adenocarcinomas. Analysis of the TCGA gastric cancer data stratifies five distinct molecular subtypes GC, specifically, Epstein-Barr virus positive (EBV), microsatellite instability (MSI), genomically stable (GS) and chromosomal instability (CIN) and these subtypes complement histopathological classifications by their distinctive patterns of DNA methylation, somatic genomic alterations and gene/protein expressions. 5 Network biology has been successfully established to systematically model molecular interactions underlying complex human diseases. 6 Here, we employed an integrative multiscale gene network analysis framework to the gastric cancer data in The Cancer Genome Atlas (TCGA-GCC) 5 to reveal key molecular mechanisms underlying GC prognosis in a data-driven manner. We postulate that these key pathways emerge as coherent modules and associate to respective key genetic and/or epigenetic alterations. Combination of these key molecular events will further generate hypothesis on multiple molecular "niche" exploited in GC etiology. Note that all the p-values reported in the manuscript were corrected for multiple testing unless stated otherwise.

Bioinformatics analysis
Gene expression data processing. Ilumina HiSeq RNA Sequencing data, processed by Reads per kilo base per million (RPKM) method from TCGA, have been downloaded and comprehensive data quality control (QC) has been performed: log2(RPKM +1) transform, quantile-normalization, correction for batch effects by center, platform and tissue source site (TSS) ids from TCGA sample barcodes, and correction for confounding factors including ethnic group, age and gender, resulting to 218 annotated primary tumor tissue samples across 26,539 genes.
Integrative network analyses of TCGA GC cohort. In order to handle multifaceted, large-scale -omics data for TCGA GC cohort (TCGA-GCC; see Supporting Information Table S1 for cohort description), we developed an integrative network analysis framework to prioritize altered pathways in GC etiology (illustrated in Supporting Information Fig. S1). The backbone of the framework work lies in constructing multiscale coexpression network (Supporting Information Fig. S1C) by Multiscale Embedded Gene Co-expression Network Analysis (MEGENA). 7 Briefly, MEGENA first selects gene pairs with significant correlations (FDR < 0.05) and then embeds them onto a three-dimensional topological sphere, leading to planar filtered network (PFN). Multiscale clustering is performed on PFN to unveil gene modules varying compactness, leading to a hierarchy of parent and child modules. The gene modules are annotated by enrichments of known pathway and functions from MSigDB 8 (Supporting Information Fig. S1d).
Then, we generated a compendium of gene signatures reflecting genomic and epigenomic alterations (Supporting Information Fig. S1b; see Supporting Materials and Methods in Supporting Material 1). These gene signatures were intersected with gene modules to test enrichments. For instance, mutation signature is identified as differentially expressed genes in somatic mutants of a gene through limma 9 with FDR < 0.05 and fold change >1.2. Gene modules enriched for the mutation signatures are then associated to this somatic gene mutation.
The gene modules were prioritized by associations to overall and recurrence-free survival across the whole TCGA-GCC as well as within clinical subtypes (Supporting Information Fig. S2). For each group of patients, the association between a module and an outcome was tested by two approaches. First, the module eigen-gene 10 (i.e., first principal component of module gene expressions) was modeled by univariate Cox proportional hazard model. 11 Second, median expression of the module eigen-gene was used to stratify patients into low-and high-expression groups, and these were tested for significant difference in outcomes by logrank test (see Supporting Materials and Methods in Supporting Material 1). Overall, the importance of each gene module is determined by a score MS m that summarizes the module's correlations to GC survival by where m = m th module, T = number of tests, p t = Cox or logrank p-values.
To validate the network interactions of the nominated drivers, differentially expressed genes by respective gene perturbation from in vivo/in vitro experiments were projected to the network neighbor genes to test for enrichments (Supporting Information Fig. S1f; see Supporting Materials 1).

Experimental procedures
Patients and specimens. One hundred and four patients who underwent gastric carcinoma resection between July 2012 and April 2015 were recruited in our study. The study protocol was approved by the Research Ethics Committee of the Fujian Provincial Cancer Hospital, and informed consent was obtained from all participants (namely, Fujian cohort; Approval number: SQ2015-068-01). There were 82 men and 22 women with a median age of 58.2 years (interquartile range, 22.0-82.0 years) (see Supporting Information Table S2). Histologic type was determined according to Lauren. 12 None of the patients had received any chemotherapy prior to surgery. Fresh tumor tissues and adjacent nontumorous stomach tissues were obtained immediately after tumor resection. One part of the tissues was immediately snap-frozen in liquid nitrogen and stored at −80 C, and the other part was fixed in 10% buffered formalin and embedded in paraffin.
Cell lines. Human gastric cancer cell line (AGS -RRID: CVCL_0139) was purchased from the Cell Banks, the Chinese Academy of Sciences (Shanghai, China). AGS cell line was authenticated by DNA finger printing analysis. AGS was grown in RPMI1640 (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA), supplemented with 10% fetal bovine serum (FBS) (Gibco; Thermo Fisher Scientific, Inc.). AGS was cultured in a 5% CO 2 incubator at 37 C. When the cells reached the logarithmic growth phase, succeeding experiments were performed. All experiments were performed with mycoplasma-free cells.

Ethics approval and consent to participate
The study protocol for Fujian cohort was approved by the Research Ethics Committee of the Fujian Provincial Cancer Hospital, and informed consent was obtained from all participants (approval number: SQ2015-068-01).

Data availability
The original gene expression data, analyzed as "TCGA-GCC" and "vTCGA-GCC," that support findings of our study are available from TCGA data portal (https://portal.gdc.cancer.gov/). Respective patient barcodes for TCGA-GCC and vTCGA-GCC are provided in Supporting Material 2C within this manuscript. Another independent validation data, GSE84437, that support findings our study, are available from Gene Expression Omnibus (GEO) with the accession code "GSE84437." The primary software MEGENA (version 1.3.7) for network analysis is publicly available as R package in The Comprehensive R Archive Network (CRAN) (https://cran.r-project.org/web/ packages/MEGENA/index.html). The development version of MEGENA is available from GitHub repository (https://github. com/songw01/MEGENA).

Multiscale gene coexpression network of gastric cancer
We developed an integrative network analysis framework to analyze multi-Omics data in the TCGA gastric cancer cohort (denoted TCGA-GCC; see Materials and Methods), using Multiscale Embedded Gene co-Expression Network Analysis (MEGENA) 7 (Fig. 1a). Two hundred and twenty-one modules pertaining a parent-child hierarchy were identified and were associated with known pathways (e.g., cell cycles, extracellular matrix, focal adhesions, immune system process, DNA metabolic process, proteasome, biological oxidations) or reflected some unknown biological processes (see Supporting Material 2A). These modules were prioritized by their prognostic power of survival in the whole TCGA-GCC population as well as various subtypes via Cox proportional hazard model (  Taking account of module overlaps by the hierarchy, the biology of top five modules was effectively captured within M121/ M666 (associated with intestinal brush border), M434 (associated with digestion) and M28 (associated with epithelial cell differentiation) and respective module hub genes were nominated as key drivers (Fig. 1c).
Reproducibility of the TCGA-GCC modules was assessed using two additional independent cohorts. As TCGA-GCC has expanded with more samples during our study, we collected the RNA-seq data from an additional 135 primary gastric tumor samples (termed vTCGA-GCC) to verify the findings from the original TCGA-GCC study. Also, we collected a microarray based gene expression study of 433 gastric cancer patients (Illumina HumanHT-12 V3.0 expression beadchip) from GEO with an accession number GSE84437 (denoted as GSE84437). The module preservation analysis 13 showed that 76.0% (168) and 94.6% (209) of the 221 TCGA-GCC modules were significantly preserved in GSE84437 and vTCGA-GCC, respectively ( Fig. 1b and Supporting Material 2B).

The most prognostic modules capture multiple facets of GC tumor microenvironment
To characterize the GC network modules, we identified genetic, epigenetic and cell type gene signatures associated with GC. We systematically identified gene signatures as potential functional manifestation of nonsynonymous somatic mutations as differentially expressed genes between the mutants and wildtypes (Somatic Mutation associated Genes (SMGs); see the section Identification of Somatic Mutation Associated Gene Signatures in Supporting Material 1). Also, gene signatures correlated to cis−/transmethylation sites were identified (expression associated methylation site gene (eMSGs); see Extraction of cis−/trans-methylation significantly correlated gene signatures in Supporting Material 1). These signatures were then used to characterize genetic and/or epigenetic alterations in the previously identified gene modules and key regulators using Fisher's Exact Test (FET) ( Fig. 1b; SMGs: Supporting Material 2B, eMSGs: Supporting Material 7). We further inferred abundances of stromal and immune cell components in the primary tumors by using ESTIMATE 14 and CIBERSORT, 15 respectively. Gene signatures significantly correlated with cell type compositions were identified and then projected onto gene modules to determine their cell type specificity ( Fig. 1b; see Supporting Material 4).
Interestingly, the top modules of TCGA-GCC network suggest strong cross-talks among different cell types in GC tumor microenvironment (Fig. 5). As summarized in Figure 5a, our data reveal these modules are regulated by somatic mutations in NEB and PIK3CA, epigenetic alterations (methylation changes by EBV infection) and loss of GKN1/2 and TFF1/2 expressions. The complex signaling network across stromal, epithelial and immune cells demonstrates activated or suppressed key oncogenic/tumor suppressive pathways in GC such as epithelial-mesenchymal (EMT), gastrointestinal mucosal barrier and cytotoxic CD8+ Tcell/Natural Killer (NK) cell activation (Fig. 5a, Table 1).
Furthermore, eigen-genes of these modules (i.e., first principal components of gene modules) 10 are significantly correlated with each other, suggesting cross-talks (Fig. 5b).
In the rest of the article, we will focus on the modules shown in Figure 1c and comprehensively investigate their network structures and biological implications as well as their interactions in the context of GC tumor microenvironment.

Deactivation of gastric mucosal barrier drives GC proliferation and invasion
The second ranked module, M434, represents downregulation of tumor suppressors such as TFF2, GKN1 and GKN2 in gastric mucosal barrier (GMB) homeostasis in GC (Fig. 1c). M434 was highly enriched for the genes downregulated in tumor compared to adjacent normal tissue (Bonferroni corrected FET p = 1.04E-25, 9.25 fold enrichment (FE)) and highly reproduced in vTCGA-GCC and GSE84437 with high overlaps (% overlap with vTCGA-GCC = 43%, GSE84437 = 48%). Suppression of the genes in M434 is associated with poor overall survival in GC and increased tumor sizes. For instance, lower expression levels of the key drivers GKN1 and GKN2 are associated with poor survival in node positive cases in TCGA-GCC (Fig. 2a). Within node positive patients, tumor weights in the group with median-low GKN1 or GKN2 expressions were significantly higher than those in the group with median-high GKN1 expression (Supporting Information Fig. S3b), but are independent of pathological T stage (Supporting Information Fig. S3c). We confirmed prognostic significance of GKN1/2 in an independent km-plotter cohort 17 (Supporting Information Fig. S3e).
Novel downstream functions of overexpressing GKN2 in GC. M434 captures a tumor suppressive axis TFF1-TFF2-GKN1-GKN2 (Fig. 1c), where TFF2, GKN1/2 are found in M434, and TFF1 is found at a neighboring module directly connected to M434 (Supporting Information Fig. S3f). Within the TFF1-TFF2-GKN1-GKN2 axis, we investigated functional impacts of overexpressing GKN2 in GC, as detailed downstream mechanisms of GKN2 are largely unknown, compared to other genes in the axis. We first comprehensively validated clinicpathological features of GKN2 expression in an independent cohort of 104 patients who underwent gastric carcinoma resection in the Fujian Provincial Cancer Hospital (Fujian cohort; see Patients and Specimens in Materials and Methods, Supporting Material 1). The correlations between GKN2 gene expression in cancer tissues and pathological features are shown in Supporting Information Table S2. GKN2 mRNA expression was significantly correlated with tumor size and location with p = 0.023 and 0.006, respectively. The average mRNA level of GKN2 in the group with large tumors was 5.2-fold lower than the group smaller tumors while the average expression in the patients with non-EGJ (esophagogastric junction) was 3.8-fold higher than EGJ.
Next The altered genes by GKN2-LV were captured in GKN2's neighborhood in the TCGA-GCC network. The downregulated genes were significantly enriched in GKN2's network neighborhood (Fig. 2b). Interestingly, TFF1 was upregulated by GKN2 overexpression (Fig. 2a,b), known synergistic heterodimeric interaction partner of GKN2 conferring antiproliferative and proapoptotic effects on gastric cancer cells. 19 This suggests GKN2 alone is capable of conferring antitumor effect.
We systematically identified de novo downstream functions of GKN2 as enriched gene modules by the GKN2 overexpression signatures (Fig. 2c). M121 and M126 were most significantly downregulated whereas M293 was most significantly upregulated. M121 represents intestinal niche/fibrosis in GC as discussed in later sections, and M293 is enriched for cholesterol biosynthesis pathway.
M28 contains several oncogenic pathways. For instance, M28 significantly intersects with cell viability signatures from KRAS mutant epithelial cells 20 (Bonferroni corrected FET p-value = 5.421E-2, 21.8 FE). Comparing by KRAS mutation status, KRAS mutant viability signatures genes in M28, MST1R and SH2D3A, were consistently upregulated in nonsynonymous KRAS mutants in TCGA-GCC and vTCGA-GCC. Furthermore, M28 is epithelial tumor specific. Several key drivers of M28 such as EPCAM 14 and ESRP1 21 are epithelial tumor cell specific markers (Fig. 1c). Expression of the genes in M28 is negatively correlated with ESTIMATE inferred stromal score (Fig. 3a).
To some extent, the association between M28 and stromal score can be explained by nonsynonymous somatic mutations in an actin-binding protein, Nebulin (NEB). NEB somatic nonsynonymous mutations were present in 26 (~11.6%) of the 224 samples, and associated to key stromal phenotypes such as lower stromal scores (Fig. 3b) and less occurrence of N3 patients (Supporting Information Fig. S4c) in NEB mutants. On the contrary, the downregulated NEB SMGs are enriched in an epithelial-mesenchymal transition (EMT) module, M109 (Supporting Information Fig. S4a,B), a key pathway mediating epithelial-stromal interaction within tumor microenvironment. 22 The epithelial-stromal interaction may involve key tumorigenic pathways in gastrointestinal stromal tumor (GIST). KIT/PDGFRA mutations serve as characteristic oncogenic mutations, present iñ 98% of GISTs. 23 The genes downregulated in active KIT/PDGFRA mutations from GISTs, curated from the CREED database (initial source data: GSE17743), 24 are significantly enriched with the stromal score correlated gene signature (Fig. 3c).
Taken together, these results suggest NEB may serve as a gateway for epithelial-stromal interactions in advanced GC microenvironment (Fig. 3d), mediating the cross-talk between tumorigenic pathways of epithelial and stromal tumor cells, and activating (NEB wildtype) or suppressing (NEB mutant) EMT, a critical pathway for cancerous cell dissemination.

GC epithelium leverages intestinal mucosal barrier niche against EBV infection
M28 strongly interacts with the intestinal epithelium specific/goblet cell-like modules, M666 and its parent module, M121 (Fig. 1c). Eigen-gene of M666 strongly correlates to M28 eigen-gene (Fig. 5b), and M121/M666 contains genes in several key intestinal mucosal barrier (IMB) functions including the members of the intermicrovillar adhesion complex (IMAC) such as MYO7B, ANKS4B, CDHR2 and CDHR5, which are critical for enterocyte brush border formation 25 ESPN involved in actin bundling 26 mucins (TFF3 and MUC2) secreted by intestinal goblet cells 27 (Fig. 2e,f ). M121/M666 are also conserved in vTCGA-GCC and GSE84437 ( Fig. 1b; Supporting Material 2B, C). M121/M666 is specifically silenced in the EBV subtype via hypermethylation. By comparing the EBV subtype and the rest of the samples, we identified the genes downregulated in the EBV subtype (using the threshold of FDR < 0.05 and fold change >1.2) and the cis-hyper-methylated genes in the EBV subtype (FDR < 0.05) and then intersected them to derive an EBV-specific hypermethylation gene signature (EBV-HEMG; see Supporting Material 6A, B). These EBV-HEMGs are significantly enriched in M121 (Bonferroni corrected FET p < 2.63E-22, 6.49 FE). Projection of M121 and M28 onto the protein-protein interaction network from the String database 16 (Fig. 2e) revealed that many key features of IMB maintenance involve interactions with GC epithelium specific pathways (M28), and are suppressed in the EBV subtype. Furthermore, suppression of key genes in M121/M666 including MYO7B and PPP1R14D, were associated with poor recurrence-free survival (see Supporting Information Fig. S5a,b).
Enrichment of PIK3CA mutations in the EBV subtype has been well established within TCGA-GCC. 5 Not surprisingly, the upregulated PIK3CA SMGs in the TCGA-GCC tumors are highly enriched in M102 (Bonferroni FET p = 1.89E-110, 14.3 FE), and M102 key drivers including TIGIT (co-inhibitory Ig receptor) show higher expressions in EBV subtype, compared to other subtypes (Fig. 4b). These results indicate PIK3CA mutations and EBV infection synergistically upregulate NK−/T-cell cytotoxic pathways including cosignaling T-cell activation receptors (Fig. 4c).

Discussion
Our integrative network approach not only revealed robust network structures and key regulators of GC, but also presented detailed signaling circuits underlying survival-associated GC microenvironments involving cross-talks among epithelial, stromal and immune cells (Table 1; Fig. 5a,b). We validated GKN2, an emerging driver of gastric mucosa, through in silico and in vitro perturbation experiments as well as pathological examination of a large number of human gastric cancer specimens.

Gastrointestinal mucosal barrier: different modes of protections for different contexts of GC progression
Several modules involved gastric-specific mucosa (M434) and intestinal-specific mucosa (M121 and M666) emerged as key tumor suppressive pathways in GC. Our analysis reveals different modes of operations for these mucosal protective functions. The gastric mucosal barrier (GMB) module is specifically associated with GC proliferation and invasion while the intestinal mucosal barrier (IMB) niche is specifically associated with the EBV infection.
GKN2 is an independent regulator of GC proliferation and invasion in GMB niche. The GMB module (M434), ranked as the second module most predictive of survival time, reveals the antitumor effect of the TFF1-TFF2-GKN1-GKN2 axis. These genes are robustly downregulated in GC compared to normal gastric tissue 29 and their downregulation in tumor is associated with increased tumor size and poor overall survival in GC. These findings are confirmed in the independent Fujian cohort.
The TFF1-TFF2-GKN1-GKN2 axis exerts pressures on GC tumorigenesis via gastric mucosal defense 30 silencing oncogenic PTEN/PI3K/AKT/mTOR 31 JAK/STAT pathway 32 ERK1/2 axes 33 and eventually inhibits GC proliferation and invasion. 19,32 GKN1 is known to play an important role in the progression of gastric cancers via inhibition of EMT and cancer cell migration 34 and TFF2 forms heteromer with MUC6, which assembles and stabilizes the laminated structure of gastric mucus and has antibiotic activity against Helicobacter pylori. 30 Heterodimeric interaction between GKN2 and TFF1 has synergistic antiproliferative and proapoptotic effects in GC. 19 In our study, in vitro overexpression of GKN2 in GC cells silences PTEN/PI3K/AKT/mTOR and JAK/STAT pathway, and inhibits GC proliferation and invasion as previously reported. Interestingly, overexpression of GKN2 alone recovered expression of its heterodimeric interaction partner, TFF1 to confer antitumor effects.
On the contrary, GKN2 overexpression upregulated cholesterol synthesis module, M293, includes HMGCR, which promotes GC growth and migration. 35 Although we observed a clear antitumor effect of GKN2 overexpression in our clinical and in vitro studies, the unexpected activation of a GC oncogenic pathway requires further in-depth investigation to clarify its roles conditioned on recovered GKN2 functions in vivo.
EBV infection silences the IMB niche. Two top ranked modules, M121 and M666, include many of intestinal epithelia specific genes in the IMB niche and are specifically downregulated by hypermethylation in the EBV subtype.
The IMB niche in GC acquires several intestinal features such as brush border formation and intestinal cell specific goblet cells. For instance, IMAC is a critical adhesion complex for brush border formation, and consists of USH1C, MYO7B, ANKS4B and intestinal specific protocadherins CDHR2/5 as essential components for the complex formation. 25 USH1C belongs to GC epithelial specific module M28 while MYO7B, ANKS4B and CDHR2/5 fall into M121. These data suggest that IMAC formation requires contributions from both actively expressed intestinal components and GC specific epithelium to become fully functional.
On the other hand, the IMB niche involving M121 and M666 is not limited to the intestinal subtype in TCGA-GCC. For instance, an intestinal brush border specific gene, ESPN, which falls into both M121 and M666, is not differentially expressed between intestinal and diffuse subtypes in TCGA-GCC (t-test p = 0.93). Many other genes in the modules show similar results. Although intestinal phenotypes such as gland formation and fibrosis are observed in early stages of GC and atropic gastritis 36 and constitute a major subtype of GC as intestinal subtype by Lauren classification, 12 the IMB niche defined by M121 and M666 is generally applicable to a broader spectrum of GC beyond the intestinal subtype.
In summary, the IMB niche represents an emergent protective mechanism in advanced GC. With the known protective roles of IMB in host defense in intestine against gut microbiome 37 the observed suppression of gastric IMB may indicate a prerequisite for EBV infection to function in GC development and progression.

Nebulin mutations: blocking stromal-epithelial cell interactions in GC
Stromal cells constitute large portions of GC tumors, and cancer-stromal cell interactions promote tumor growth and metastasis. 38 Coexistence of often benign gastro-stromal intestinal tumors (GISTs) with gastric adenocarcinoma has been frequently observed. 23 In our TCGA-GCC analysis, a KRAS dependent GC tumorigenic module (M28) associates to varying degree of "stromal-ness." This stromal-ness corresponded to GIST specific oncogenic mutations, KIT and PDGFRA, 23 pertaining stromal-epithelial interactions in GC.
In our study, somatic nonsynonymous mutations in an actin-binding protein, Nebulin (NEB) is associated to this stromal-epithelial interactions. Nebulin is a giant 600-to 900-kDa filamentous protein constituting cytoskeletal matrix that coexists with the thick and thin filaments within the sarcomeres of skeletal muscle 39 and is an overexpressed protein in gastric cancer tissues compared to normal. 40 However, its roles in GC progression have remained largely unknown.
We present NEB-centered molecular mechanisms of GC invasion by leveraging differentially expressed genes in NEB somatic nonsynonymous mutation. Specifically, upregulated genes in NEB mutants coincided with GC epithelial specific M28, and the respective downregulated genes were associated to EMT module (M109), which bears ZEB1 as a top key driver, a key transcription factor in EMT and metastasis 41 (Supporting Information Fig. S4a,b). Another interesting key driver of M109 is MYLK (myosin light chain kinase), which phosphorylates myosin regulatory light chains to facilitate myosin interaction with actin filaments to produce contractile activity, 42 and involved in cell motility and morphology.
Together with NEB functions in maintaining cytoskeletal matrix, the downregulation of M109 in NEB mutants strongly suggests disruption of dynamic changes in actin cytoskeleton via loss-offunction mutation in NEB, a necessary process of EMT for cancerous cells to acquire cell motility and morphological changes. 22 This also explains the underrepresentation of N3 patients in NEB mutant samples in TCGA-GCC (Supporting Information Fig. S4c).
Overall, NEB-centered GC mechanisms illustrate the complex cross-talks between GC stroma and epithelium via gastric actin cytoskeleton, and propose NEB mutation as a protective mechanism by disrupting EMT.
PIK3CA and EBV driven methylation: recipes for GC immune evasion M102 is NK−/T-cell cytotoxic pathway module, including key receptors for T-cell and NK-cell activation for immune surveillance of GC tumor. Interestingly, M102 is subject to combination of genetic alteration (PIK3CA mutation) and epigenetic alteration (hypomethylation in EBV subtype), which synergistically upregulated M102 compared to other GC subtypes in TCGA-GCC.
Among the key driver genes of M102, TIGIT is a poliovirus receptor (PVR)-like protein, an immunoreceptor expressed in T-cells that acts as inhibitory checkpoint on both of T-cells and NK-cells. 43 Recently, TIGIT emerged as a promising immune checkpoint blockade target in GC. TIGIT+ CD8+ T-cell population increases in GC, and these cells show functional exhaustions impairing antitumoral activities. 44 Blocking TIGIT showed synergistic effect in recovering antitumoral CD8+ T-cell functions with anti PD-L1 treatment. 44 As M102 also includes CD8+ T-cell markers such as CD8A and NKG7 as key drivers, these strongly suggest M102 depicts exhausted CD8+ T-cells by exploiting coinhibitory TIGIT signaling.
The top key driver of M102, CD96 (also known as T Cell-Activated Increased Late Expression Protein), is expressed in T-and NK-cells with adhesive functions to modulate their interactions and enhance cytotoxicity. 45 In HIV-1 infected adults, CD96 is associated with different cell effector functions of CD8+ T cells. 45 It has been shown that targeting NK cell activity via CD96 has promising anticancer potential, complementary to the existing PD-1 and PD-L1 targeted therapeutics 46 and suppresses metastasis in melanoma lung metastasis mouse model (B16F10). 47 To our knowledge, the effectiveness of CD96 blockade, thereby recovering NK-cell antitumoral activity in GC has not been explored.
These results suggest targeted immunotherapy for PIK3CAmutant and/or EBV-infected GC patients, via targeting adaptive and innate immune response systems. Especially, CD96 and TIGIT are attractive blockade targets in such GC patients though future investigation is needed to confirm the finding.
In summary, the most prognostic coexpression network modules reveal a series of pathways representing complex cross-talks among epithelial, stromal and immune cells in GC (Fig. 5a), which are potentially driven by environmental factor (i.e., EBV infection), genetic alterations (i.e., NEB and PIK3CA mutations) and epigenetic alterations (i.e. hypo-/hyper-methylations). We further systematically investigate key regulators and proposed interactions of the top modules. In particular, we validate antitumor effects of one key network driver GKN2 overexpression and identified downstream pathways of GKN2 in GC. These findings enable generation of novel hypotheses regarding complex interplay among multifaceted axes (e.g., different pathways and different cell types) in GC progression. For instance, we propose the IMB niche in M121 and M666 as a suppressor of EBV driven immune evasion in M102, where host defense function of the IMB niche may be capable of inhibiting EBV activation.
Our integrative network analysis of TCGA-GCC reveals some fundamental patterns of molecular interactions and specific mechanisms in GC progression. The network models and the key regulators identified here pave a way for defining novel therapeutic strategies. Future work will include validation of prioritized key drivers and subnetworks in GC and develop network-based biomarkers for stratifying GC patients.