Associations of individual and joint expressions of ERCC6 and ERCC8 with clinicopathological parameters and prognosis of gastric cancer

Human tissue specimens

Paired gastric cancer and para-cancerous normal tissue samples were collected from 145 individuals with gastric cancer, who were diagnosed at the anorectal department of the First Hospital of China Medical University. Each selected patient had not received neoadjuvant chemoradiotherapy or any other treatment from 2012 to 2015. Histological diagnoses were completed following the updated Sydney System Classification (Stolte & Meining, 2001) and the World Health Organization criteria (Hamilton & Aaltonen, 2000) for gastritis and GC, respectively. The 2010 7th edition of the TNM staging system of the International Union Against Cancer/American Joint Committee on Cancer was selected to stage the tumor (Santiago, Sasako & Osorio, 2009), based on postoperative pathologic examination. This study was approved by the ethics committee of the First Hospital of China Medical University. Written informed consent was obtained from all participants.

Immunohistochemistry staining and evaluation

Every single paraffin-embedded sample was sectioned into a 4-µm-thick slide. All the deparaffinized and rehydrated slides were boiled at 95 °C for 30 min in citrate antigen retrieval solution (pH 6.0). Immunohistochemistry (IHC) staining was performed with the avidin-biotin complex method as our previous study described (Jing et al., 2017). CSB (ERCC6) rabbit polyclonal antibody (1:300 dilution, Origene, TA313375, USA) and AB1(anti-ERCC8) antibody produced in rabbit (1:500 dilution, Sigma, AV31542, USA) were used.

Two experienced pathologists who were blinded to patient-related information evaluated and scored the staining results independently. A semi-quantitative scoring criterion was applied to analyze the staining extent of each slide (Detre, Jotti & Dowsett, 1995). The staining score was categorized on the basis of intensity: 0 - no staining, 1 - light staining, 2 - moderate staining and 3 - strong staining and coloring ratio: 0 (≤5%), 1 (5%–25%), 2 (25%–50%), 3 (50%–75%), 4 (≥75%) of the IHC results. An immunoreactivity score (IS) was generated by multiplying the two scores for each sample. All the scores were applied to indicate certain extent of positive staining except the score of 0, which suggested a negative protein expressed level.

RNA isolation and qRT-PCR

Total RNA were extracted from 36 pairs of tissue specimens with RNAiso Plus reagent based on the experimental guidelines (TaKaRa, Kusatsu, Japan). The cDNA was synthesized by PrimeScript RT Master Mix following the protocol (Perfect Real Time; TaKaRa, Kusatsu, Japan). Relative expressions of ERCC6 and ERCC8 were detected with TB Green Premic EX Taq II (TliRNaseH Plu; TaKaRa, Kusatsu, Japan). β-actin was selected as the endogenous reference control. Each melting curve was with single peak in the study. A 2−ΔCt method was applied to assess the relative expression levels of ERCC6 and ERCC8. Primers for ERCC6, ERCC8 and β-actin were tabulated in Table S1.

Clinical information and RNA-seq data collection

Clinical information of 109 IHC cases concerning age, gender, smoking, family history and alcohol consumption were obtained via questionnaire. Clinical characteristics were collected from medical records, including TNM stage, Lauren’s classification, Borrmann classification, tumor size, phase of progression, lymph node metastasis, perineural invasion, and vascular invasion. The final follow-up was completed on July 2016. Full data concerning prognosis were obtained from 97 participants.

Data of 415 stomach cancer patients, which included RNA-seq data and survival information, were downloaded from TCGA (https://cancergenome.nih.gov/) (Weinstein et al., 2013; Akbani et al., 2014). Detailed information of clinicopathological parameters including age, gender, grade, TNM stage and histological type were also downloaded for further analyzing. GEPIA (Gene Expression Profiling Interactive Analysis) (http://gepia.cancer-pku.cn/) was applied to obtain the mRNA expression levels of ERCC6 and ERCC8 and the data of disease-free survival (DFS) (Tang et al., 2017).

Interaction and functional analysis of ERCC6 and ERCC8

Given the clinical significance of ERCC6/ERCC8 in GC, we further investigated their biological functions. First we used ERCC6 and ERCC8 as core genes to construct protein-protein interaction (PPI) networks by Search Tool for the Retrieval of Interacting Genes (STRING v.11.0; https://string-db.org/; accessed on August 27, 2020), to mine proteins that have functional interactions with ERCC6/ERCC8 (Szklarczyk et al., 2019). Analytic information including nodes degrees and biological networks was visualized with Cytoscape platform (v.3.7.2) (Su et al., 2014). And ten most associated proteins were showed in the diagrams.

Then we conducted enrichment analyses of Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) to explore the biological functions of ERCC6 and ERCC8 with the Database for Annotation, Visualization and Integrated Discovery (DAVID; v.6.8; https://david.ncifcrf.gov/home.jsp; accessed on August 31, 2020), a user friendly database providing comprehensive analysis of gene annotation (Huang da, Sherman & Lempicki, 2009). R language (Version 3.6.3) and the ggplot2 package were applied to visualize the analytic results. Terms with a P < 0.05 were deemed significant and for GO, only top ten terms of each group were selected to be visualized.

Identification of regulation networks of ERCC6 and ERCC8 by GSEA /KEGG

Gene set enrichment analysis was conducted on the GSEA platform (version 4.1.0; https://www.broadinstitute.org/gsea/) coupled with MSigDB database (Subramanian et al., 2005a; Subramanian et al., 2005b). Oncogenic signature gene sets (c6.all.v7.1.symbols.gmt) and TCGA expression data were included in this analysis. Through C6 Oncogenic Signatures of GSEA, we could clarify the gene sets associated regulation networks that were involved in ERCC6 and ERCC8. First TCGA expression data was grouped into ERCC6/ERCC8-high and ERCC6/ERCC8-low groups according to their expression levels. And then differentially expressed genes of the two groups were figured out. Finally significantly changed oncogenic regulation networks of these genes were identified through a thousand times of phenotype permutation test and the metric for the analysis was set as pearson. A normalized p value <0.01 and a false discovery rate (FDR) less than 0.25 were selected as criteria for significant enrichment results. Normalized enriched score (NES) was applied to rank the obtained results. KEGG pathway analysis was further conducted to identify pathways these gene sets mainly involved in.

Identifying gene-gene interactions between ERCC6 and ERCC8

To identify interactions between ERCC6 and ERCC8, we performed gene-gene interaction analysis with the Gene Multiple Association Network Integration Algorithm (GeneMANIA; https://www.genemania.org/; accessed on November 13, 2020), which is a user-friendly interface providing analysis with available genomics and proteomics data (Warde-Farley et al., 2010). To display interactions, nodes and links represented genes and networks respectively in the visualized results.

Validation of protein expression of ERCC6/ERCC8 related genes

CCLE (The Broad Institute Cancer Cell Line Encyclopedia; http://www.broadinstitute.org/ccle) is an online database which provides genomic data and protein expression of more than 900 cancer cell lines (Barretina et al., 2012). To validate the protein expression of the genes obtained from GSEA, we extracted mRNA expression data of ERCC6, ERCC8 and protein expression data of EIF4E, ERBB2, JAK2, Src, and Cyclin_D1, which were mainly involved in PI3K/AKT/mTOR pathway, from 37 gastric cancer cell lines. Pearson’s correlation was further calculated to evaluate the correlation between their expression levels. A P < 0.05 was deemed significant.

Statistical analysis

IBM SPSS Statistics for Windows, version 23.0 (IBM Corp., Armonk, N.Y., USA), GraphPad Prism V5.0 software (GraphPad software, USA) and R platform (Version 3.6.3) were used for statistical analyses. For IHC, Pearson χ² test was used when analyzing the correlations between ERCC6 and ERCC8 expressed levels and clinicopathological parameters; univariate and multivariate Cox regression analyses have been selected for the determination of their impact on overall survival (OS), and variables including age, TNM stage, perineural invasion, vascular invasion and lymph node metastasis were further adjusted in the multivariate model to evaluate the independent prognostic value. For qRT-PCR data, Wilcoxon matched-pairs signed rank test was used to assess the differential expression of ERCC6 and ERCC8 among different groups. For RNA-seq data, Wilcoxon test and Kruskal-Wallis H test have been employed when calculating the interrelationships between ERCC6/8 expressions and clinical characteristics; two-sided Log-rank test and multivariate Cox proportional-hazards model adjusted by gender, age, grade, stage, T, N, and M were selected to clarify the prognostic value of ERCC6 and ERCC8. A P < 0.05 suggested statistical difference.

Expression of ERCC6 and ERCC8 in GC and adjacent normal mucosa

In this study we compared the expressed levels of ERCC6 and ERCC8 between GC and adjacent normal mucosa. Representative ERCC6 and ERCC8 staining were present in Fig. 1. Our results suggested that individual and joint expressions of ERCC6 and ERCC8 were obviously higher in adjacent normal mucosa than in GC tissues (all P < 0.001) (Table 1). Specifically, the ERCC6-ERCC8 double positive rate dropped to 16.5% in GC and the double negative rate was only 1.8% in adjacent normal mucosa. As shown in Fig. 2, no significant difference was observed with ERCC6 mRNA expression in GC (P = 0.300) and ERCC8 mRNA expression was lower in GC tissues when compared with normal tissues (P < 0.0001).

Correlations of ERCC6 and ERCC8 expression with clinicopathological characteristics in GC

We explored the associations of ERCC6/ERCC8 expressed levels with clinicopathological parameters in GC patients and the results were summarized in Table S2. ERCC6 expression was significantly related to Borrmman classification (P = 0.017), Lauren’s classification (P = 0.004) , TNM stage (P = 0.005) (P = 0.012 for T stage) and perineural invasion (P = 0.001). High ERCC6 expression was observed in gastric cancers of Borrmman class I–II, TNM stage I–II, intestinal-type and without perineural invasion. Expression of ERCC8 was statistically higher in patients with TNM stage I-II in comparison to stage III-IV (P < 0.001) (P < 0.001 for T stage; P = 0.002 for lymph node metastasis), and was higher in those with early-stage and small size (P = 0.031 and 0.007, respectively). Higher expression of ERCC8 was also observed in intestinal type GC (P = 0.008). As for the joint expression of ERCC6/ ERCC8, double positivity was related to small tumor size (P = 0.005), Borrmman I-II stage (P < 0.001), TNM I-II stage (P = 0.001) (P < 0.001 for T stage), Lauren intestinal type (P = 0.014) of GC and negative perineural invasion (P = 0.015). Double negativity was associated with TNM III–IV stage (P < 0.001) (P = 0.002 for T stage), positive perineural invasion (P = 0.002), advanced stage (P = 0.034) and diffuse type (P < 0.001) of GC.

Analysis results of RNA-seq data obtained from TCGA was shown in Fig. 3, higher ERCC6 expression was related with better T stage (P = 0.027), which is consistent with the IHC analysis results, while no statistical results were found between ERCC6 expression and age (P = 0.570), gender (P = 0.646), histological type (P = 0.425), grade (P = 0.072), stage (P = 0.091), N stage (P = 0.572) and M stage (P = 0.242). As for ERCC8, Fig. 4 revealed that overexpressed ERCC8 was closely related to worse grade (P = 0.018), advanced stage (P = 0.022), and worse N stage (P = 0.037); the associations between ERCC8 expression and age, gender, histological type, T stage and M stage were not significant (all P > 0.05).

The relationship between ERCC6 and ERCC8 expression and prognosis of GC

As shown in Table 2, univariate survival analysis showed no significant correlation between protein expressed levels of ERCC8 and GC prognosis (P = 0.211), while a significant correlation was observed between ERCC6 protein expressed levels and GC prognosis (P = 0.047, HR =3.416, 95% CI [1.017–11.475]). In addition, double negative expression of ERCC6 and ERCC8 was significantly associated with poorer prognosis (P = 0.022, HR = 2.603, 95% CI [1.148–5.905]). However, no statistical association was found between the expression levels of ERCC8/ ERCC6-ERCC8 and GC prognosis (both P > 0.05). Because TNM stage (P < 0.001), age (P = 0.039), perineural invasion (P = 0.043), vascular invasion (P = 0.007), and lymph node metastasis (P < 0.001) are all statistically associated with gastric cancer prognosis (Table S3), Cox’s proportional hazards model adjusted by perineural invasion, TNM stage, age and vascular invasion was further applied to evaluate the prognostic value. However, the multivariate analysis suggested that ERCC6 or ERCC8 or ERCC6-ERCC8 expressed level was not an independent factor for GC prognosis (all P > 0.05) (Table 2).

Survival analysis with RNA-seq data suggested that higher ERCC8 mRNA expression was related to better OS (P = 0.025; Fig. 4I) while ERCC6 expression made no sense (P = 0.15; Fig. 3I). And multivariate Cox proportional hazard models showed no significant results with ERCC6 (P = 0.969; HR = 0.995; 95% CI [0.761 –1.300]) and ERCC8 (P = 0.078; HR = 1.393; 95% CI [0.964–2.014]). And the results suggested that double lower ERCC6 and ERCC8 mRNA expression was not significantly associated with GC prognosis neither in univariate model (P = 0.242; HR = 1.364; 95% CI [0.810–2.297]) nor in multivariate model (P = 0.197; HR = 1.426; 95% CI [0.832 –2.446]). There existed no significant correlation between ERCC6 mRNA expression (P = 0.600; File S4A), ERCC8 mRNA expression (P = 0.780; Supplementary file 4B) and DFS, respectively.

Function analysis for ERCC6 and ERCC8

To explore the function of ERCC6 and ERCC8, we first constructed PPI networks using ERCC6 and ERCC8 as core genes. With a confidence score of more than 0.4, Fig. 5A showed that ten most associated proteins of ERCC6 were ERCC2, ERCC3, ERCC4, ERCC5, ERCC8, TP53, XPC, XPA, POLR2A and POLR2I; and for ERCC8, the related proteins were XAB2, DDB1, ERCC3, ERCC4, ERCC2, CUL4A, UVSSA, ERCC5, XPA and GTF2H2 (Fig. 5B).

Then we conducted analyses of GO and KEGG according to the network results we obtained from String. As revealed in Fig. 6, ERCC6 network genes showed enrichment in molecular functions of protein N-terminus binding, damaged DNA binding and DNA-dependent ATPase activity (Fig. 6A). They were mainly involved in nucleoplasm, transcription factor TFIID complex, and holo TFIIH complex according to cellular components analysis result (Fig. 6B). Figure 6C showed that ERCC6-interactive genes were significantly enriched in biological processes of nucleotide-excision repair and UV protection. Analysis results of KEGG suggested that these genes were closely related to nucleotide excision repair, Huntington’s disease, RNA polymerase and basal transcription factors (Fig. 6D). As for ERCC8 network genes, results of GO enrichment analysis showed that these genes were related to the composition of transcriptional initiation complexes and ubiquitin ligase complexes, and not surprisingly, their main molecular functions and participated biological processes bore a remarkable resemblance to ERCC6 network genes (Figs. 7A–7C). Similar KEGG results, nucleotide excision repair and basal transcription factors were also observed in ERCC8 network genes. And a special part of this analysis results was ubiquitin mediated proteolysis (Fig. 7D).

Identification of gene sets associated regulatory networks of ERCC6 and ERCC8

Further, we identified the most positive and negative related gene sets with ERCC6 and ERCC8, to figure out the cellular regulatory networks in GC that ERCC6 and ERCC8 were involved in. Oncogenic signatures analysis indicated that in GC there existed six and ten most significant gene sets for ERCC6 and ERCC8, respectively. Among the results, both ERCC6 and ERCC8 were associated with TBK1 and BCAT associated cellular regulatory networks; ERCC6 was also associated with EIF4E, MTOR, JAK2 and CSR related regulatory networks; ERCC8 was also associated with PIGF, RB, ERBB2, GCNP, SRC and CYCLIN D1 related regulatory networks. Detailed information were shown in Fig. 8 for ERCC6 and Fig. 9 for ERCC8. KEGG pathway further revealed that these gene sets were mainly involved in the PI3K/AKT/mTOR pathway (File S5).

Gene-gene interaction network between ERCC6 and ERCC8

Gene-gene interaction network accessed from GeneMANIA clarified the correlations of ERCC6 and ERCC8 among pathway, predicted, shared protein domains, physical interactions, co-localization, and co-expression. As shown in Fig. 10, there existed direct interactions including physical interactions and pathway and indirect interactions including prediction, co-expression, colocalization and shared protein domains between ERCC6 and ERCC8.

Validation of protein expression of ERCC6/ERCC8 related genes in gastric cancer cell lines

Pearson’s correlation analysis using data from 37 gastric cancer cell lines showed that ERCC6 mRNA expression was significantly correlated with JAK2 protein expression (r = −0.345; P = 0.037), and ERCC8 mRNA expression was associated with Src protein expression (r = −0.417; P = 0.010). No significant relationship was found between ERCC6 mRNA expression and protein expression of EIF4E (r = −0.069; P = 0.686), ERBB2 (r = 0.295; P = 0.076), or between ERCC8 mRNA expression and Cyclin_D1 protein expression (r = 0.171; P = 0.311). Detailed information were shown in Table S4.