A three‐gene signature from protein–protein interaction network of LOXL2‐ and actin‐related proteins for esophageal squamous cell carcinoma prognosis

Human ESCC cell lines KYSE150 14 and SHEEC 15 were used in the study. KYSE150 cells were cultured in RPMI‐1640 medium with 10% fetal bovine serum, whereas SHEEC cells were maintained in high‐glucose DMEM medium with 4.0 mmol/L l‐Glutamine and Ham's Nutrient Mixture F12 with 10% newborn bovine serum.

Human ESCC tumors and adjacent nontumorous esophageal epithelial tissues were collected directly after surgical resection between January 2000 and January 2013, at Shantou Central Hospital (Sun Yat‐sen University, Shantou, China). Fifteen paired tissues between 2012 and 2013 were used for genome‐wide transcriptome sequencing (RNA‐seq, SPA accession number: SRP064894↗, Shanghai Biotechnology Corporation, Shanghai, China) 16. We validated the three‐gene signature using a validation tissue microarray cohort of 124 randomly selected patients with clear pathologic diagnosis and follow‐up information (Table S1). All tumors were confirmed as ESCC by pathologists in the Clinical Pathology Department of the hospital, and the tissue specimens were classified according to the tumor‐node metastasis (TNM) classification of the International Union against Cancer. Written informed consent was obtained from all patients or their appropriate surrogates. The study was approved by the ethical committee of Shantou Central Hospital and Shantou University Medical College.

Published gene expression datasets, comprising 316 patients of ESCC in total, were from different centers and platforms in the study. RNA‐seq data and clinical information of 84 patients with ESCC were downloaded from The Cancer Genome Atlas (TCGA) database 17. The gene microarray expression profiles and corresponding clinical information of patients with ESCC were obtained from the publicly available Gene Expression Omnibus (GEO) database, including 53 patients from GSE23400↗ (Affymetrix HG‐U133A array) 18, 119 patients from GSE53624↗ (Agilent human lncRNA+mRNA array V.2.0), and 60 patients from GSE53622↗ (Agilent human lncRNA+mRNA array V.2.0) 19. For prognostic gene signature analysis, the 119 specimens in GSE53624↗ were randomly divided into two groups by the stratified randomization function from SAS (r) Proprietary Software 9.2 (Table S2). One group was used as a training cohort (n = 60) to develop a multigene signature for ESCC specimens, whereas the other one was a testing cohort (n = 59). GSE53622↗ (n = 60) was provided as a validation cohort for assessing the prediction performance of the signature. The analyses were done blinded from clinical and biologic data and sample identification.

KYSE150 and SHEEC cells were cultured for 24–48 h. When cells reached 80–90% confluence, they were fixed in 4% paraformaldehyde in phosphate‐buffered saline (PBS) for 10 min, rinsed three times with PBS, and permeabilized in 0.1% Triton X‐100 for 10 min or less. Cells were washed with PBS and blocked in blocking serum (Life Technology, Shanghai, China) for 60 min at room temperature. Cells were stained overnight with anti‐LOXL2 antibody (1:100; Novus, Littleton, CO), anti‐LOXL2 C‐terminal antibody (1:100; Abcam, Cambridge, UK), or anti‐FlagM2 (1:200; Sigma‐Aldrich Guangzhou, China) in blocking buffer, and then probed with Alexa Fluor 488‐conjugated Affinipure donkey antirabbit secondary antibody (1:200; Jackson ImmunoResearch, West Grove, PA) and Actin‐stain^™ 555 phalloidin (1:200; Cytoskeleton, Denver, CO), followed by counterstaining with 0.1 μg/mL 4′,6‐diamidino‐2‐phenylindole (Sigma Aldrich, Guangzhou, China). Confocal imaging was performed on a laser‐scanning confocal microscope (LSM 880; Carl Zeiss, Oberkochen, Germany).

Genes involved in regulation of the actin cytoskeleton were obtained from KEGG data (Release 76.0, KEGG database/map04810, gene numbers = 181). Gene set enrichment analysis (GSEA2.2.0, Cambridge, MA), a powerful analytical method, was performed using Pearson's correlation to determine the degree of linear relationship between LOXL2 and actin cytoskeleton genes.

A PPI human reference set of 224,522 PPIs between 18,130 proteins was obtained from the Biological General Repository for Interaction Datasets (BioGRID, Release 3.4), and was mapped into Cytoscape software (Version 3.2.1, http://www.cytoscape.org↗) as a visualizing parent PPI network. LOXL2 PPIs from variable resources and PPIs of ACTB or ACTG1 were extracted from the parent PPIN to generate a LOXL2‐ACTB/ACTG1 PPIN. The network then created a LOXL2‐ACTB/ACTG1 PPIN in ESCC with a filtering method according to reads per kilobase of exon per million mapped reads (RPKM) of the RNA‐seq of ESCC specimens. When the RPKM of one gene in 15 tumor samples is higher than 10, the gene will be regarded as an effective detecting candidate for the following analysis. Self‐loops and duplicated edges of proteins in the networks were removed. Gene Ontology (GO) and KEGG pathway assignments were performed using Gene Functional Classification Tool from the DAVID Bioinformatics Resources database (Release 6.7). The GO terms in the Biological Process (GOTERM‐BP‐FAT) and KEGG pathways with P < 0.05 and enrichment score > 2 were considered as significantly enriched function annotations.

Tissue microarray (TMA) construction of esophageal carcinoma tissue has been described in our previous study 20. The SuperPicTure Polymer Detection Kit and the Liquid DAB Substrate Kit (Zymed/Invitrogen, Carlsbad, CA) were used to conduct immunohistochemistry according to the manufacturer's instructions. Primary antibodies used for immunohistochemistry staining in this study are listed below: anti‐LOXL2 antibody (1:100; Abnova, Taipei, China), anti‐ Fibronectin antibody (1:200; Santa Cruz Biotechnology, Dallas, Texas), and anti‐E‐cadherin (1:100; BD Biosciences, Franklin Lakes, NJ). The TMAs were intelligently scanned by a multimodal and high‐throughput imaging system, Vectra 2.0.8 (PerkinElmer Inc., Hopkinton, MA). Following the manufacturer's protocols, the expressions of three proteins were quantified using the Nuance v3.0 incorporated with proprietary image analysis software, inForm v1.2 (PerkinElmer Inc., Hopkinton, MA). First, we built a tissue segmentation algorithm to separate cancer from its surrounding stroma in multispectral image of stained TMAs sections. Then, the tumor cells were circled individually and the amount of staining intensity was measured and digitized. Through a series of analysis above, the sum scores of each protein expression in the three‐gene signature of individual tissue specimens were obtained (range from 0 to 300) and subjected to statistical analysis.

The freely available gene expression data of ESCC specimens in GSE53624↗ and GSE53622↗ were already log2 scale transformed 19. High expression or low expression of a given gene was defined on the basis of the level of mRNA log2 expressions by the automatic optimum selection method in X‐tile software (Release 3.6.1), and the same was applied to immunohistochemistry scores of the protein. In all cohorts in the study, each gene of each patient was assigned a numerical value in which high expression equaled 1 and low expression equaled 0. To estimate overall survival of a signature formed by LOXL2, CDH1, and FN1 genes, we adopted a developed strategy using a univariate Cox proportional hazards regression analysis coefficient for the three genes 21. A given patient's prognosis risk score was then calculated as the sum of each gene's score, which was derived by multiplying the gene expression level by its corresponding regression coefficient (prognostic risk score = ΣCox coefficient of gene Gi × expression value of gene Gi). Finally, the patients were separated into two groups (high‐risk signature and low‐risk signature) using the 50th percentile cut‐off of the risk score of the three‐gene signature as a threshold. To reduce the influence of heterogeneity among different patients or cohorts, the risk scores and the cut‐off value derived from the training cohort were reestimated and then applied to the other three cohorts.

Unless otherwise stated, statistical analyses were performed with SPSS 22.0 software (IBM, Chicago, IL) or GraphPad Prism 5 software (La Jolla, CA). The differences between groups were analyzed by the Student t test or paired‐samples t test. The chi‐squared test or Fisher's exact test was used to assess the correlations between category variables. Overall survival (OS) was estimated by the Kaplan–Meier method. Multivariate Cox proportional hazards regression (HR) analysis was applied to evaluate independent prognostic factors associated with overall survival, and the three‐gene signature, age, gender, and tumor stage, was used as covariates. To estimate the sensitivity and specificity of the signature and area under the curve (AUC) value, receiver operating characteristic (ROC) analysis was adopted in three cohorts. The differences between the AUC of the three‐gene signature, TNM stage, and a combination of both were detected by MedCalc software (Vision 15.8, Ostend, Belgium). A P < 0.05 was considered to indicate statistical significance, and all statistical tests were two‐sided.

To explore the subcellular location of LOXL2 in ESCC, we studied the human ESCC cell lines KYSE150 and SHEEC. LOXL2 in SHEEC cells was mainly localized at the median and basal cytoplasm and nucleus, but weak staining was also observed at the apical region (Fig. 1A upper panel). Similar staining of LOXL2 was detected in KYSE150 cells, except for a strong signal in the nucleus (Fig. 1A lower panel). To our great surprise, LOXL2 and filamentous actin (F‐actin) were colocalized in ESCC cells (Fig. 1A and S1). Although F‐actin staining was abundant throughout the cells, colocalization of LOXL2 and F‐actin was primarily restricted to the cytoplasm near the cell membrane and a few cell–cell contacts. The gene sets enrichment analysis (GSEA) plots indicated that LOXL2 expression was positively correlated with regulation of actin cytoskeleton genes in patients with ESCC (Fig. 1B). Altogether, these observations point to a new relationship in which LOXL2 is related to the actin cytoskeleton in ESCC.

Based on the work above, we assumed that the relevance of LOXL2 and actin cytoskeleton in ESCC is through a number of their interacting proteins. We gathered certain and potential LOXL2 protein–protein interactors (PPIs) from variable protein databases and published reports (Fig. 2A). All these and acknowledged human PPI data from the BioGRID database were combined to provide the original data for subsequent analyses. A LOXL2‐ACTB/ACTG1 PPIN, from all human tissues, containing 504 nodes and 3568 edges was generated by mapping total LOXL2, ACTB, and ACTG1 PPIs to the parental PPI network (Fig. 2C left panel). To create an ESCC‐specific subnetwork, we examined 21,157 gene expressions in RNA‐seq of 15 patients with ESCC. The low‐expressional ones out of 504 proteins in ESCC tissues were excluded through an RPKM‐filtering method, and the remaining 2580 PPIs between 362 proteins were simultaneously selected for the following survey (Fig. 2B, C right panel).

We found that the genes in the ESCC‐specific network were intimately connected to regulation and reorganization of the cytoskeleton, as well as multiple typical carcinoma‐related cell behaviors, such as cell adhesion, cell invasion, and cell migration (Fig. 2D and Table S3). To better know the network composition of the ESCC‐specific network, we grouped the proteins in three subcategories upon the network topological properties: central proteins consisting of LOXL2, ACTB, and ACTG1, core proteins comprising 14 interacting proteins, and other 345 common proteins (Fig. 2C).

The 14 LOXL2‐ACTB/ACTG1 PPIN core proteins are listed as below: A2M, ANGPLT4, CAP1, CDH1, CTGF, CTNND1, DAB2, EZH2, FN1, HSP90AA1, PFN1, TAF11, TMSB4X, and YWHAZ. Due to its absence in microarray expression profiles in our study (GSE53624↗ and GSE53622↗), the core protein PFN1 was left out at first. Central gene LOXL2 and the other 13 core genes were considered as potential ESCC prognostic biomarkers and were further studied.

For prognostic gene signature analysis, we obtained microarray data from 119 patients with ESCC (GEO accession number: GSE53624↗) and randomly assigned it to either the training cohort or the testing cohort. Among the 60 cases in the training cohort, potentially crucial genes, excluding CTGF, DAB2, and TMSB4X, were significantly differentially expressed between ESCC tumors and paired normal specimens (Fig. 3A and Fig. S2). Five from the differentially expressed genes were upregulated in tumor tissues, whereas the others were downregulated in contrast. Moreover, high‐level expression of five genes was found to be significantly associated (P <0.05 by the log‐rank test) with outcome of patients with ESCC in the training cohort: four genes (FN1, LOXL2, TAF11, and TMSB4X) were associated with shorter OS and high‐level expression of one gene, CDH1, was associated with a better clinical outcome (Fig. 3B and S3).

Summarizing the analytic results in the training cohort, LOXL2, CDH1, FN1, and TAF11 were in accordance with the inclusive criteria. On account of the biological functions of these genes, we ultimately selected LOXL2 (lysyl oxidase‐like 2), CDH1 (E‐cadherin 1, type 1), and FN1 (fibronectin 1) to integrate a three‐gene signature in ESCC. Afterward, the three‐gene signature in the patients among the training cohort was evaluated with a risk score algorithm provided by univariate Cox regression. The Cox coefficients for LOXL2, CDH1, and FN1 were 1.089, 0.980, and −0.827, respectively. The 50th percentile of the prognostic risk scores was 0 (range from −0.827 to 2.069), therefore, all individuals were classified into two groups: high‐risk signature (risk score > 0) and low‐risk signature (risk score ≤ 0). Patients with a low‐risk gene signature (n = 45) had a longer median overall survival than patients with a high‐risk gene signature (n = 15) (42.0 months vs. 23.3 months, P =0.003) (Fig. 3B). Multivariate Cox proportional hazard regression analyses revealed that mortality risk assessed by the three‐gene signature was, along with tumor stage, an independent predictor of survival (P <0.001, HR = 4.48, 95% CI = 2.03–9.90) among ESCC patients in the training cohort (Table 1). In conclusion, high risk, predicted by the three‐gene signature, predicts poor survival in ESCC patients.

To test the three‐gene signature for its prognostic value, we applied it to the microarray data of the testing cohort. The same risk score algorithm and criteria, as which were derived from the training cohort, classified 19 and 40 ESCC patients of the testing cohort into the high‐risk and low‐risk groups, respectively. The high‐risk gene signature was associated with a median OS of 23.8 months, whereas the low‐risk signature was associated with a median OS of 43.0 months (P =0.001 by the log‐rank test) (Fig. 4A and S4). The three‐gene signature was then validated in another cohort comprising 60 patients with ESCC (GEO accession number: GSE53622↗). Similar to the above, patients in the low‐risk group had significantly longer survival time than those in the high‐risk group (median survival months: 37.7 vs. 26.2, P =0.021) (Fig. 4B and S5A). The high‐risk three‐gene signature remained significantly associated with shorter overall survival in the testing cohort (P =0.002, HR = 2.87, 95% CI = 1.47–5.61) and the validation cohort (P =0.032, HR = 2.21, 95% CI = 1.07–4.58) (Table 1).

Not only in the mRNA level, but the three‐gene signature was further validated in the protein level with a tissue microarray cohort of 124 patients with ESCC. Each protein displayed two immunostaining phenotypes; that is, low staining and high diffuse staining (Fig. 4C). The OS in the high‐risk signature group was significantly shorter than that in the low‐risk signature group (median survival months: 46.5 vs. 73.6, P =0.007) (Fig. 4D and Fig. S5B). Most importantly, the high‐risk group predicted by three‐gene signature also represented a higher mortality risk than the low‐risk group (P =0.016, HR = 1.88, 95% CI = 1.13‐3.14) (Table 1). Therefore, the three‐gene signature is a robust high‐risk factor for ESCC prognosis.

In current clinical analysis, TNM stage classification is deemed to be the optimal prognostic indicator for cancer. As demonstrated in Table S4 and S5, there was no statistically evident relationship between expression level of the three‐gene signature and any ESCC clinical characteristic in not only the testing cohort but also in the three validation cohorts. A prognostic model combining TNM stage and the three‐gene signature was conducted and receiver operating characteristic (ROC) analysis was performed. The combined model had a higher area under curve (AUC) than both the signature and the TNM stage alone in the four cohorts (Fig. 5). In applying this criteria to ESCC patients, we found that the predictive ability of the combined model was significantly better than the signature or the TNM stage alone. As an example, the predictive ability between the combined model and TNM stage alone was 0.696 versus 0.616 (P <0.05, Fig. 5B).