Comprehensive analysis of COMMD10 as a novel prognostic biomarker for gastric cancer

Data download and analysis

mRNA expression data and clinical information were downloaded from the TCGA database (https://cancergenome.nih.gov/) and Genotype-Tissue Expression (GTEx) database. The gene amplification and mutation status of COMMD10 was obtained from the cBioPortal (http://www.cbioportal.org/) for Cancer Genomics. The Human Protein Atlas (HPA) (http://www.proteinatlas.org/) database was obtained for subcellular localization and immunofluorescence images of COMMD10 expression in GC cells.

Survival analysis

The GEPIA database (http://gepia.cancer-pku.cn/) and the TCGA database were used to estimate the correlation between COMMD10 expression and survival in STAD patients with different clinical characteristics. The R package “survival” (version 3.6) was used to obtain an OS survival map for COMMD10. A threshold value of 50% was chosen as the division threshold to divide the cohort into high and low expression groups. Visualization was performed using ggplot2.

Construction and evaluation of prognostic nomogram

All independent clinicopathological prognostic factors from a Cox regression analysis were selected to construct a line plot to assess the 1-, 3-, and 5-year OS probabilities for patients with STAD. The accuracy of the nomogram was validated by comparing the predicted probabilities of the line graph with the observed actual probabilities through calibration curves. The overlapping reference lines indicate that the model is accurate. The RMS R package (https://cran.r-project.org/web/packages/rms/index.html) and survival R package were used to generate a nomogram. C-index was used to determine the discrimination of the nomogram, which was calculated by a bootstrap approach with 1,000 resamples.

Correlation and gene set enrichment analysis

To understand the biological processes and pathways in which COMMD10 may be involved, correlation analyses between COMMD10 and other mRNAs in GC were performed using data collected from TCGA. Genomic ontology (GO) terminology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway studies were performed for significantly co-expressed genes using the clusterProfiler package in R language.

Association of COMMD10 expression levels with m6A modifications in STAD

The association of COMMD10 expression levels with m6A-related genes, including YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, WTAP, RBM15, RBM15B, ZC3H13, HNRNPC, METTL14, METTL3, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, HNRNPA2B1, VIRMA, FTO and ALKBH5 was analyzed utilizing the R statistical computing language. Correlations between genes were analyzed using the corrlot software package.

Cell culture

Human normal gastric epithelial cells (Ges-1) were obtained from Guangdong Hybribio Biotech Co., Ltd (Guangdong, China). And GC cell lines AGS, MNK-45, HGC-27, SGC-7901 were obtained from Hunan Fenghui Biotechnology Co., Ltd (Hunan, China). Ges-1, HGC-27, MNK-45 were cultured in RPMI-1640 medium (Gibco, USA) plus 10% fetal bovine serum (Procell, Wuhan, Hubei, China) and 1% penicillin and streptomycin (ABT920; G-clone, Beijing, China). The basal medium for AGS and SGC-7901 were Ham’s F12 (Procell, Las Vegas, NV, USA) and Dulbecco’s Modified Eagle Medium (DMEM) (Gibco, Waltham, MA, USA), respectively, and the rest of the culture conditions were the same. All were incubated at 37 °C in a humidified incubator with 5% CO₂.

RNA extraction, reverse-transcription RNA, and quantitative real-time polymerase chain reaction

Total cellular RNA was extracted using TRIzol reagent (YEASEN, Shanghai, China), according to the instructions. Reverse transcription was performed using HiScript II Q RT SuperMix for qPCR kit (R223-01; Vazyme, Nanjing, Jiangsu, China) and cDNA was synthesized according to the instructions. cDNA was synthesized using LightCycler 480 II system (Roche, Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal reference. COMMD10 and GAPDH primers were synthesized by Tsingke Biotechnology Co., Ltd (Beijing, China). The primer sequences were as follows:

COMMD10-forward: 5′-ATGGCGGTCCCCGCGGCGCT-3′

COMMD10-reverse: 5′-TCATGTAAGGGAATCCAGCTG-3′

GAPDH-forward: 5′-GGTCACCAGGGCTGCTTTA-3′

GAPDH-reverse: 5′-GGATCTCGCTCCTGGAAGATG-3′

Western blot

A total of 10 paired cancer and paracancerous tissues were collected from December 2021 to March 2022 from gastric cancer patients who underwent surgical resection in the Department of Gastrointestinal Surgery at the Third Xiangya Hospital of Central South University, and written informed consent was obtained from all patients prior to the study. The research was approved by Ethics Committee of Xiangya Third Hospital, Central South University (22101). The methods for tissue protein extraction were as follows: tissue clipping, tissue homogenizer homogenization; lysis on ice with RIPA lysis solution; centrifugation; protein quantification using a BCA kit (KGPBCA; KeyGEN BioTECH, Jiangsu, China); appropriate amount of protein up-sampling; electrophoresis; gel cutting; membrane transfer; milk closure; addition of 1:800 dilution of rabbit anti-COMMD10 (123867; ZEN BIO, Chengdu, Sichuan, China); 1:1,000 dilution of rabbit anti-GAPDH monoclonal antibody (GB11002; Servicebio, Wuhan, Hubei, China); incubated overnight at 4 °C; washed the membrane; and added goat anti-rabbit IgG secondary antibody (HRP; Proteintech, Rosemont, IL, USA); incubated at 37 °C for 90 min; membrane washed; developed by ECL luminescence kit (BL520A; Biosharp, Beijing, China); developed and fixed; and the expression level of COMMD10 in cells and tissues was analyzed with GAPDH as the internal reference protein.

Immunohistochemical staining

Ten paired GC cancer and paracancerous tissues were used for IHC staining. Briefly, paraffin sections were dewaxed, antigenically repaired, endogenous peroxidase blocked by hydrogen peroxide, and serum closed, then incubated overnight at 4 °C with 1:100 dilution of anti-COMMD10 primary antibody (123867, ZEN BIO, Chengdu, Sichuan, China). This was followed by the addition of a 1:200 dilution of HRP-labeled goat anti-rabbit secondary antibody (GB23303; Servicebio, Ghent, Belgium). Then the sample was incubated at room temperature for 50 min, and the DAB kit (Servicebio, Ghent, Belgium) developed color. The positive expression of DAB is brownish yellow.

m6A dot blot analysis

Total gastric cancer RNA was extracted and RNA concentration was detected using Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). The sample underwent heat denaturation at 65 °C for 10 min in a thermal cycler and was quickly cooled on ice. A total of 1 ug RNA was loaded on nylon membranes. Next, the membrane was cross-linked by UV, blocked with 5% non-fat milk, and incubated with M6A antibody (1:600, HA601049) at 4 °C overnight. Subsequently, the membranes were incubated with HRP-conjugated goat anti-mouse IgG dilution (1:2,000, Proteintech, Rosemont, IL, USA). Finally, the membrane was visualized using a chemiluminescence imaging analysis system. To determine the consistency of loading, membranes were stained with 0.02% methylene blue for 30 min, then rinsed twice with ultrapure water and photographed.

Statistical analysis

Box line plot and scatter plot methods were used to detect the expression level of COMMD10 gene in GC patients. The cut-off value of COMMD10 expression was selected as the median method of gene expression. Spearman correlation analysis was used to explore the degree of correlation between COMMD10-related genes and COMMD10. Univariate and multivariate Cox analyses were used to screen for potential prognostic factors. In all analyses, *, ** and *** indicate p < 0.05, p < 0.01 and p < 0.001, respectively.

COMMD10 expression in GC

Data downloaded from TCGA and GTEx were used to analyze the expression of COMMD10 in pan-cancer. The results revealed significant differences in the expression of COMMD10 in a variety of tumors (Fig. 1A). Based on the role of COMMD10 in GI-related tumors such as liver and colon cancers, we sought to explore whether it also plays an important role in GC. We further evaluated the mRNA expression levels of COMMD10 in cancer and paracancerous tissues of STAD patients, and the results showed that the expression levels of COMMD10 were significantly higher in tumor tissues than in paracancerous tissues (p < 0.01) (Fig. 1B). In addition, the expression of COMMD10 mRNA in normal cell lines and GC cell lines were detected by qRT-PCR. The expression levels in MNK-45 and AGS cells were found to be significantly higher than those in the normal gastric epithelial cell line GSE-1 (Fig. 1C). Western blotting was also used to detect the expression of COMMD10 in 10 pairs of cancer and paracancerous tissues. We observed that COMMD10 protein levels were higher in most cancer tissues than in paired paracancerous tissues (Figs. 1D, 1E). Finally, IHC staining on clinical specimens were performed to further validate the expression level of COMMD10 in GC. The results showed that the expression of COMMD10 was higher in cancer than in paracancerous tissues (Figs. 1F, 1G). In addition, to understand the mutation level of COMMD10 in STAD, we analyzed its genome and copy number. OncoPrint mapping of the COMMD10 gene in STAD patients in the TCGA dataset was analyzed using cBioPortal plots (Fig. 1H), which showed less than 3% gene amplification, missense and deep deletion mutations in COMMD10. These data confirmed that COMMD10 was highly expressed in GC.

Correlation between COMMD10 expression levels and clinicopathological characteristics of GC patients

TCGA database was used to analyze the correlation between different clinical characteristics and COMMD10 expression levels in patients with STAD. Specifically, the relationship between COMMD10 expression and ten clinicopathological characteristics of GC patients were investigated, including age, gender, H. pylori infection, T-stage, N-stage, M-stage, pathological stage, residual tumor and OS events (Fig. 2). The results showed that high expression of COMMD10 was significantly associated with T stage (p < 0.05) (Fig. 2D), residual tumor (p < 0.001) (Fig. 2H), and OS events (p < 0.05) (Fig. 2I) in these patients. These data suggested that COMMD10 was significantly upregulated in STAD and correlated with different clinical features.

High COMMD10 affects the prognosis of patients with STAD in different clinicopathological states

To determine whether COMMD10 expression affects patient survival, GC patients were divided into a high COMMD10 expression group (the top 50% of the highest expression samples) and a low COMMD10 expression group (the remaining 50% of samples) using the GEPIA database. Survival analysis was then performed based on the mean COMMD10 expression value. Kaplan–Meier survival analysis showed that high COMMD10 expression was associated with poor prognosis in STAD patients in terms of OS (HR = 1.7, p = 0.00093) and progression-free disease survival (HR = 1.8, p = 0.0025) (Figs. 3A, 3B). Subsequently another Kaplan–Meier analysis was applied based on the TCGA database to analyze the relationship between COMMD10 expression and prognosis in STAD patients. As seen Fig. 3C, patients with high COMMD10 expression had a worse prognosis than those with low COMMD10 expression (HR = 1.56, p = 0.007). However, it showed no significant value in progress free interval (HR = 1.38, p = 0.071) (Fig. 3D). Then subgroup analyses were performed of patients with different clinicopathological status STAD and showed that high COMMD10 expression was significantly associated with poor prognosis in STAD in: patients over 65 years of female patients (HR = 1.85, p = 0.039), T3 (HR = 1.65, p = 0.037), tissue grade G3 (HR = 1.81, p = 0.005), tumor anatomical site: cardia/proximal (HR = 2.46, p = 0.037) and age (HR = 1.73, p = 0.01). These data are shown in (Figs. 3E–3I).

High COMMD10 expression is an independent risk factor for overall survival

Univariate and multivariate Cox risk regression analyses were used to identify the independent prognostic factors of STAD in terms of age, sex, T stage, N stage, M stage and COMMD10. Univariate and multivariate analyses showed that high COMMD10 expression was an independent prognostic factor for OS in STAD patients (HR = 1.492, 95% CI [1.053–2.113], p = 0.025). In addition, age (HR = 1.792, 95% CI [1.247–2.575], p = 0.002), N stage (HR = 1.782, 95% CI [1.134–2.801], p = 0.012), M stage (HR = 2.181, 95% CI [1.220–3.899], p = 0.009) were also independent prognostic factors (Table 1).

Construction of prognostic line graphs based on independent prognostic factors

Then these independent prognostic factors including age, N-stage, M-stage, and COMMD10 were used to construct a prognostic nomogram and a calibration curve was plotted to test the validity of the nomogram. The relationship between the four clinicopathological variables (age, N stage, M stage and COMMD10) and the 1-year, 3-year and 5-year survival probabilities of OS was visually displayed (Fig. 4A). As shown in (Fig. 4B), the calibration curves for 1, 3 and 5 years showed the consistency of our results and predicted values, indicating that this COMMD10-based nomogram performed satisfactorily. The c index of 0.645 showed moderate predictive accuracy. In conclusion, this nomogram may be a better model for predicting survival of GC patients compared to individual prognostic factors.

Functional enrichment of COMMD10-related genes

We perform differential expression analysis between the low and high COMMD10 in STAD using the Limma package. Only protein-coding genes were retained. A total of 900 differential genes, including 204 up-regulated genes and 696 down-regulated genes, were selected according to the |logFC| >1 and p.adj <0.05 screening conditions (Fig. 5A). Next, we analyse the co-expressed genes associated with COMMD10 expression in the TCGA STAD dataset. A total of 2,728 genes were obtained under |cor spearman| >0.3 and p < 0.05, with 2,724 positively associated genes and four negatively associated genes. The top 50 positively correlated genes were shown in the (Fig. 5B). The 2,728 co-expressed genes were studied for GO terms and KEGG pathway using the clusterProfiler package in R language. COMMD10 co-expressed genes were involved in 709 biological processes, 231 cellular components, 157 molecular functions and 40 KEGGs (p.adj < 0.05 and q value < 0.05). The bubble plots showed the top five information of biological processes, cellular components, molecular functions and KEGGs, respectively. GO term annotations indicated that these genes were mainly involved in RNA localization, RNA splicing, ncRNA processing, tRNA metabolic process; ubiquitin ligase complex, chromosomal region, trophoblast region, mitochondrial matrix; ubiquitin-like protein transferase activity, histone binding (Figs. 5C–5E). The KEGG pathway analysis indicated that these genes were mainly involved in signaling pathways such as RNA transport, ubiquitin-mediated protein hydrolysis, autophagy, mRNA surveillance pathway, and spliceosome. (Fig. 5F). Table S1 summarized the GO terms and KEGG pathway details for COMMD10 co-expression enrichment analyses.

Association of COMMD10 expression levels with m6A modifications in STAD

The functional enrichment of COMMD10-related genes were mainly involved in biological processes such as RNA localization, RNA splicing, RNA transport, mRNA monitoring pathway, and spliceosome. And m6A RNA methylation plays a crucial role in RNA splicing, export, processing, translation and decay. This suggested that the expression level of COMMD10 may correlate with m6A modification. The level of m6A modification in GC tissues with high COMMD10 expression were assessed using dot blot analysis. The results showed that the expression of m6A was higher in GC tissues with high COMMD10 expression compared with adjacent tissues (Fig. 6A). Next the association between COMMD10 expression levels in STAD and 20 m6A-related genes was explored. As shown in Fig. 6B, COMMD10 expression was significantly and positively correlated with YTHDF3, HNRNPC, METTL14, YTHDC2, RBMX, WTAP, METTL3, YTHDF2, FTO, VIRMA, YTHDC1, RBM15, ALKBH5, HNRNPA2B1, YTHDF1, and RBM15B (p < 0.001). Figure 6C showed the scatter plot of COMMD10 expression with 11 m6A-related genes (Spearman r > 0.3, p < 0.001). As shown, COMMD10 expression were significantly positively correlated with YTHDF3 (r = 0.441, p < 0.001), HNRNPC (r = 0.434, p < 0.001), METTL14 (r = 0.431, p < 0.001), YTHDC2 (r = 0.415, p < 0.001), RBMX (r = 0.395, p < 0.001), WTAP (r = 0.377, p < 0.001), METTL3 (r = 0.365, p < 0.001), YTHDF2 (r = 0.360, p < 0.001), FTO (r = 0.355, p < 0.001), VIRMA (r = 0.345, p < 0.001), and YTHDC1 (r = 0.326, p < 0.001). We also observed that YTHDF3, METTL3, FTO, YTHDC1 and YTHDF1 had a poor prognosis in GC (Fig. 6D). This suggests that there may be a regulatory relationship between YTHDF3, METTL3, FTO, YTHDC1 or YTHDF1 and COMMD10, thus affecting the prognosis of patients.

Correlation of COMMD10 expression with immune characteristics

To explore the correlation between COMMD10 expression levels and tumor immune response, TCGA database was used to investigate immune infiltration in STAD with different COMMD10 expression levels. The results showed that COMMD10 expression in STAD patients was positively correlated with macrophages, dendritic cells (DCs), immature dendritic cells (iDC), eosinophils, central memory T (Tcm) cells, effective memory T (Tem) cells, helper T (Th) cells, Th1 cells, Th2 cells, and negatively correlated with NK CD56 leukocytes (Fig. 7A). Further studies showed that COMMD10 expression was significantly positively correlated with macrophage (r = 0.238, p < 0.001) (Fig. 7B) and iDC infiltration levels (r = 0.197, p < 0.001) (Fig. 7C). This prompted us to explore the relationship between the level of COMMD10 expression and immune infiltration. We found significant differences (p < 0.05) in the levels of infiltration of helper T cells, iDCs and macrophages into immune cells when COMMD10 expression was divided into high and low expression groups (Figs. 7D–7G). It is noteworthy that COMMD10 was closely related to M2 macrophages based on the CIBERSOFT algorithm (Fig. S1).