Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma

1The incidence (accounting for 5.6% of all cancer cases) and mortality (accounting for 7.7% of all cancer deaths) of gastric cancer (GC) ranks fifth and fourth worldwide among all malignancies, which critically threatens human health. It was estimated that over one million new cases of GC were reported worldwide in 2020. Gastric adenocarcinomas derived from gastric glandular epithelial cells accounts for more than 90% of all GCs. GC is a multifactorial disease, which is contributed by both environmental and genetic factors, such as smoking, family history, Epstein–Barr virus (EBV) infection, alcohol consumption, and diet. Most GCs are diagnosed at late stage of progression due to limited premalignant indications and symptoms. Nowadays, the application of the endoscopic examination largely improved the survival rate of GC patients, and a 30% reduction in GC mortality using endoscopic screening. However, stomach cancer is still one of the most lethal malignant tumors, with a 5 year survival rate of around 20%. Therefore, it is necessary to identify more reliable biomarkers for predicting the prognosis and exploring more potential therapeutic targets in GC.

2Accumulating evidences indicated that mitochondria plays essential roles in regulation of cell growth, cell death, and cell metabolism during the whole process of tumor progression. Mitochondria are involved in bioenergetics metabolism, such as ATP production, reactive oxygen species (ROS) production, apoptosis, and calcium homeostasis. Moreover, mitochondrial dysfunction may contribute to the chemoresistance. Therefore, mitochondrial-targeting therapies may be applied for the treatment of GC, including ROS production and elimination, mitochondrial fission and fusion, ATP production, and apoptosis. For instance, nanohybrid-induced oxidative stress triggered mitochondria-mediated autophagy, which inhibited cell growth in cancer cell.

3Considering that mitochondrial dysfunction was a risk factor for the tumorigenesis of GC, identifying effective mitochondrial-related biomarkers for the prognosis of GC patients should be an encouraging direction of research. Several studies have constructed GC prognosis-related models to predict patient survival. However, rare studies have been applied for the establishment of prognostic models for GC associated with mitochondria.

4Tumor microenvironment (TME) was mainly composed of the stromal cells, immune cells and cytokines. The components of TME affected the immune cell evasion or inhibition, and drug resistance in malignancies. For example, immune progenitors in the complex microenvironment of the TME were more likely to differentiate into M2 macrophages and Treg cells, but not to play their tumor-inhibiting functions as fully mature immune cells. In addition, the response to immune checkpoint blockade (ICB) was closely related to the constitution of the TME. ICB revived an effective anti-tumor immune response. It was reported that PD-1/PD-L1 inhibitors immunotherapy has an impact on the therapy of patients with advanced gastrointestinal malignancy. Moreover, a study reported that tumors with a higher tumor mutation burden (TMB), had a better immunotherapy response, especially with PD-1/PD-L1 blockade. Thus, to figure out the correlation between the risk score and TME, we explored the TME signatures, the immune cell infiltration in TME, the expression level of immune checkpoints, and the response to immunotherapy.

5In summary, a novel mitochondrial-related risk model was constructed using NOX4, FKBP10, ALDH3A2, and MAOA gene set, which could effectively predict the prognosis and immunotherapy responsiveness for patients with STAD. In addition, we estimated the drug sensitivity of STAD patients to 138 drugs, including chemotherapy drugs, immunotherapy drugs, and targeted drugs, et al., and found that patients in high-risk group was more sensitive to rapamycin, PD-0325901 and dasatinib, whereas patients in low-risk group was more sensitive to AZD7762, CEP-701 and methotrexate. Taken together, our mitochondrial-related risk model could be a reliable prognostic biomarker for personalized treatment of STAD patients.

6RNA-seq data and microsatellite status information for 407 STAD samples were downloaded from the TCGA database (https://www.cancer.gov/tcga). The clinical information was extracted from the UCSC Xena (http://xena.ucsc.edu). 61 samples were excluded due to incomplete clinical information or survival less than 30 days. In total, 346 samples, comprising 317 tumor samples and 29 healthy samples, were analyzed in the present study. Two validation cohorts, GSE66229 and GSE15459, were applied in the present study (https://www.ncbi.nlm.nih.gov/geo/). In GSE15459, 10 samples were excluded due to survival less than 30 days. Altogether, 300 and 182 samples were analyzed in GSE66229 and GSE15459, respectively. The list of mitochondrial-related genes was collected from MitoCarta 3.0 database (https://www.broadinstitute.org/mitocarta/mitocarta30-inventory-mammalian-mitochondrial-proteins-and-pathways) and the Gene set enrichment analyses (GSEA, http://www.gsea-msigdb.org/gsea/index.jsp) (Additional file 10: Table S1).

7The “limma” package of R (version 3.5.1) was applied to produce DEGs between normal and tumor samples, or between high-risk and low-risk groups from the training set. |Log (2) fold change|> 2 and adjusted P < 0.01 were the criteria for defining DEGs. “GdcVolcanoPlot” packages in R were employed to generate volcano map to visualize the DEGs, and a Venn plot was exploited to display the common DEGs in both DEGs groups and mitochondrial-related genes.

8The mitochondrial-related genes were screened by univariate Cox regression, Lasso regression analysis and multivariable Cox regression analysis to construct a novel prognostic gene signature. Each sample’s risk score was calculated using the following formula:where expgene, i, and βi represent the expression level of gene, the number of signature genes, and the coefficient index, respectively. In all participated cohorts, the samples were divided into low-risk and high-risk groups based on the risk score (median cut-off value). To analyze the survival conditions for the prognosis signature, the optimized cutoff and the Kaplan–Meier (K–M) survival curve were conducted by R package “survival” and “survminer”. The predictive performance was presented by ROC curve, risk plot and concordance index (C-index). Detailed information for prognostic genes was obtained from The Human Protein Atlas (HPA, https://www.proteinatlas.org/) and National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/).

9Risk score and clinical factors including age, gender, T stage, N stage, M stage, tumor stage, family history, H pylori infection, grade, reflux history, and disease types were analyzed using univariate Cox regression analysis to screen the factors significantly related to survival (P < 0.1). Then, multivariate COX regression analysis was applied to identify the candidate predictors significantly related to survival (P < 0.05). Based on this, nomograms were constructed using these predictors, and scores in nomogram model were assigned for these variables. By adding the scores of the predictors enrolled in nomogram model, the total score of each patient was obtained. Finally, the patient's survival outcome in 1, 3 and 5 years can be calculated using the total score and the probability of survival outcome. ROC curve, calibration curves and decision curve analysis (DCA) were applied to estimate the discrimination and accuracy of the nomogram model.

10In the present study, R “clusterProfiler”, “org.Hs.eg.db”, “enrichplot” and “ggplot2” package (R version: 3.5.1) were employed to analyze the function of mitochondrial-related DEGs, or the DEGs between high-risk and low-risk groups. Furthermore, adjusted P < 0.05 was used to filter the functional candidates.

11Curated sets v7.4 collections were obtained from the Molecular Signatures Database as the target sets with which GSEA was performed by using GSEA 4.2.1 software. The total transcriptome of tumor samples was used for the GSEA, and only gene sets with P < 0.001 and FDR, q < 0.001 were regarded to be statistically significant.

12Stromal scores were calculated using the ESTIMATE algorithm by R (version 3.5.1) package “estimate”. The list of TME-related biomarkers was extracted from the Gene set enrichment analyses (GSEA, http://www.gsea-msigdb.org/gsea/index.jsp) (Additional file 10: Table S2–S5).

13The abundance of 22 tumor-infiltrating immune cells (TIICs) in STAD samples was calculated using the CIBERSORT algorithm by R package. CIBERSORT is a deconvolution algorithm that can infer 22 kinds of TIICs and harnesses the ability to predict the relative abundance of each immune cell population by calculating the expression of specific marker. The relative abundance of the TIICs between high-risk and low-risk groups was compared using the Wilcox text. The list of immune check points was referenced from a published study. Immune score and tumor purity were also calculated using the ESTIMATE algorithm by R (version 3.5.1) package “estimate”. The list of immune cell signatures was downloaded from TISIDB (http://cis.hku.hk/TISIDB/download.php).

14The capability of risk score in predicting the response to immunotherapy or 138 drugs for chemotherapies/targeted therapies was explored in the present study. The 50% inhibiting concentration (IC50) values of the 138 drugs were calculated using the “pRRophetic” package of R (version 3.5.1) and the value was normally transformed. The detailed information of 138 drugs was acquired from Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/). The potential response to immunotherapy was inferred by the tumor immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu) score.

15The somatic mutation data were downloaded from cBioportal database (https://www.cbioportal.org/). The R (version 3.5.1) package “maftools” was then used to draw a waterfall plot to illustrate the mutation landscape in STAD patients with the high- and low-risk group and calculate the TMB score for each sample.

16The normal gastric epithelial cell line GES-1 and human gastric cancer cell lines SGC-7901 and HGC-27 were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). Cells were cultured in DMEM medium (Gibco, Thermo Fisher Scientific, Inc., Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS, Hyclone). Cells were routinely cultured in a humidified atmosphere containg 5% CO2 at 37 ℃. A total of 10 fresh tumor and paired adjacent normal tissues from patients with STAD were collected in the First Affiliated Hospital of Shanxi Medical University (Taiyuan, Shanxi, China). All patients provided written informed consent, and this study was approved by the ethics committee of Shanxi Medical University.

17RNAs in tissues and cell lines of STAD were extracted with a RNAiso Plus (Takara,Tokyo, Japan) and were reversely transcribed into cDNA usin PrimeScriptTM RT Master Mix (Takara). Quantitative real-time PCR (qRT-PCR) was performed by TB Green ®Premix Ex TaqTM II (Takara). β-actin was used as reference genes. The primers were listed in Additional file 10: Table S6.

18R (version 3.5.1) and the GraphPad Prism 8 software were applied for statistical analysis. A Student’s t-test was used to analyze the expression and the distribution of risk score, stromal score, immune score, tumor purity and TMB in different groups. Chi-square test is applied to evaluate the difference in immunotherapy response, status of top 5 mutant genes and clinical factors in different groups. The correlation was evaluated using the Spearman method. C-index was used to estimate the predictive power of age and risk score to OS. P < 0.05 was defined as statistically significant.

19The general workflow of our current study was presented in Fig. 1. As shown in Additional file 10: Table S7, 2381 DEGs, including 2145 protein-coding genes, were screened and visualized via volcano maps between normal and tumor groups (Fig. 2A, Additional file 10: Table S8). Next, combined analysis for selected mitochondrial-related genes from the GSEA and 2145 DEGs from our study were performed to filter out 183 candidate mitochondrial-related DEGs in STAD (Fig. 2B, Additional file 10: Table S9).

20GO enrichment analysis were then carried out to uncover important roles of mitochondrial-related DEGs in STAD. These DEGs potentially participated in small molecule catabolic process, regulation of mitochondrial organization, et al. Regarding the cellular component, they were mainly related to mitochondrial matrix, mitochondrial outer membrane, et al. In terms of molecular function, these DEGs were involved in tubulin binding, and ubiquitin-like protein ligase binding, et al. (Additional file 1: Fig. S1A, Additional file 10: Table S10). Moreover, KEGG pathway analysis was also applied to demonstrate important pathways being involved by these DEGs, such as lipid and atherosclerosis, Hepatitis B infection, Diabetic cardiomyopathy, and apoptosis, et al. (Additional file 1: Fig. S1B, Additional file 10: Table S11).

21Based on above 183 mitochondrial-related DEGs, 19 genes were further selected as potential risk factors for the prognosis of patients with STAD through univariate Cox regression analysis (P < 0.05, Fig. 2C). The gene number was further narrowed down to 9 according to LASSO regression analysis and to 4 by multivariable Cox regression analysis (Fig. 2D, E). Finally, 4 mitochondrial-related DEGs, including NOX4, FKBP10, ALDH3A2 and MAOA, were utilized to establish a prognostic model for patients with STAD (Table 1).

22As shown in Fig. 2F, higher expressions of NOX4 and FKBP10, and lower expressions of ALDH3A2 and MAOA were observed in tumor samples compared with the normal tissues, respectively. Moreover, the immunohistochemistry results from HPA database showed that FKBP10 was upregulated in gastric cancer tissues, while ALDH3A2 and MAOA were downregulated in gastric cancer tissues, when compared with corresponding non-cancerous tissues (Additional file 1: Fig. S1C). Further K-M analysis demonstrated that patients with higher expression of NOX4 (P = 0.030), FKBP10 (P = 0.040), and MAOA (P = 0.018) had a shorter OS than those with lower expression, respectively (Additional file 1: Fig. S1D). However, the patients with higher level of ALDH3A2 had a better OS than those with lower expression, even though it is a bit beyond statistically significant difference (P = 0.052, Additional file 1: Fig. S1D).

23Then, the risk score for each patient with STAD in both training and validation cohorts was computed based on the following formula:Patients were divided into high-risk and low-risk subgroups based on the median risk score. K–M curves showed that patients in high-risk group had worse OS (P = 0.0009, Fig. 3A). To assess the accuracy of prognostic risk models in predicting 1-, 3-, and 5-year OS, ROC curves were plotted with AUC values of 0.635, 0.640, and 0.793, respectively (Fig. 3B). The relationship between the risk score and the survival time, survival status, and risk ranking, and a heatmap of the expressions of the 4 genes were shown in Fig. 3C. Taken together, these results demonstrated the robustness of our risk model in predicting the prognosis of patients with STAD.

24The robustness of the prognostic risk model was further validated in GSE66229 and GSE15459 datasets. In line with that of the training cohort (TCGA-STAD), patients in high-risk group also had worse prognosis in the validation cohorts (Fig. 3D, Additional file 2: Fig. S2A). In GSE66229 dataset, the AUC values of the ROC curve for 1-year, 3-year and 5-year survival were 0.620, 0.625, and 0.601, respectively (Fig. 3E). Corresponding AUC values of 0.620, 0.647, and 0.657 were observed in GSE15459 dataset (Additional file 2: Fig. S2B). The higher the risk score, the worse the survival (Fig. 3F, Additional file 2: Fig. S2C). The heatmaps of the expressions of the 4 genes were shown in Fig. 3F, Additional file 2: Fig. S2C. Consistent with the TCGA-STAD training cohort, the expressions of the NOX4 and FKBP10 were significantly up-regulated, while the expressions of the MAOA and ALDH3A2 were significantly down-regulated in STAD in GSE66229 validation cohort (Additional file 2: Fig. S2D). Next, we systematically analyzed the relationship between the risk score and clinical characteristics in STAD. The risk scores were remarkably higher in patients with H pylori infection, and cystic, mucinous and serous neoplasms (Additional file 3: Fig. S3). Nevertheless, no differences were observed in the mean of risk score among the groups of age, gender, T stage, N stage, M stage, tumor stage, family history of GC, grade, and reflux history (Additional file 3: Fig. S3). The clinical characteristics of the low-risk and high-risk subgroups were then compared, and the difference of Gender (P = 0.006), N stage (P = 0.043), H pylori infection (P = 0.035), disease type (P = 0.046) and survival status (P = 0.008) among the two risk subgroups reached statistical significance (Table 2).

25The nomogram integrated the risk score and all important clinical features, which can be used to quantitatively predict the prognosis of patients and provide a reference for clinical decision making. In our study, risk score (P = 0.0005) and age (P = 0.020) were finally identified as prognostic indicators by using univariate and multivariate Cox regression analysis to construct nomogram (Table 3). As a result, a predictive nomogram integrating risk score (a score of 100) and age (a score of 67.5) for prognosis was constructed (Additional file 4: Fig. S4A). ROC curves showed that the AUC values of the nomogram were 0.651, 0.664, and 0.749 for 1-, 3-, and 5-years OS, respectively (Additional file 4: Fig. S4B). The calibration curve showed that the actual survival probabilities at 1-, 3- and 5-year were almost in accordance with the survival probabilities predicted by the nomogram model (Additional file 4: Fig. S4C). The decision curves showed that the nomogram model was better than other factors in predicting the prognosis in STAD (Additional file 4: Fig. S4D).

26We further conducted functional enrichment analyses of the 298 DEGs in high-risk and low-risk groups (Additional file 10: Table S12). GO enrichment analysis indicated that the differential genes annotated to biological processes were involved in extracellular matrix (ECM) organization and extracellular structure organization. Differential genes annotated to cellular component categories were mainly enriched in collagen-containing ECM and collagen trimer. Differential genes annotated to molecular function categories were mainly enriched in ECM structural constituent and collagen binding (Fig. 4A, Additional file 10: Table S13). The top 10 pathways obtained by KEGG analysis were: protein digestion and absorption, proteoglycans in cancer, focal adhesion, human papillomavirus infection, PI3K-Akt signaling pathway, ECM-receptor interaction, cell adhesion molecules, axon guidance, cAMP signaling pathway, and vascular smooth muscle contraction (Fig. 4B, Additional file 10: Table S14). GSEA results showed that risk score was significantly associated with ECM glycoproteins, core matrisome, ECM organization in high-risk group (Fig. 4C). The detailed GSEA results for the high-risk and low-risk groups were presented in Additional file 5: Fig. S5.

27Given the TME-associated signal pathway was enriched through functional enrichment analyses, we explored the relationship between the risk score and the TME signatures. As shown in Fig. 5A, risk score was highly positively correlated with stromal score in STAD, and stromal score was higher in high-risk group compared to that in the low-risk group. We further investigated the relationship between risk score and TME components. Our results indicated that risk score had a significant and positive correlation with the expressions of the majority of carcinoma associated fibroblast (CAF) signatures (Fig. 5B), as well as ECM-collagen and matrisome signatures (Fig. 5C, D). Taken together, these results suggested a close relationship between mitochondrial-related risk score and TME signatures in STAD.

28The tumor immune microenvironment was closely related with the therapeutic effects and prognosis of patients with malignant tumor. It is reasonable to check the relationship between the risk score and the immune cell infiltration in STAD. The contents of naive B cells, regulatory T cell (Tregs), M0 macrophage, and M2 macrophage were remarkably higher in the high-risk group. In contrast, CD8+ T cells and resting CD4+ T cells were higher in low-risk group (Fig. 6A, Additional file 6: Fig. S6A, B). Consistent with the above results, risk score was positively correlated with the expressions of the majority signatures of M2 macrophage, while negatively correlated with the expressions of the majority signatures of activated CD8+ T cell (Fig. 6B). The correlations between the risk score and the signatures of other immune cells were presented in Additional file 8: Fig. S8. Consistent with the previous studies, the ESTIMATE results showed that the patients in high-risk had higher immune score, and significantly lower tumor purity (P = 0.370), than those in the low-risk groups (Fig. 6C, D). Due to the positive correlation between the risk score, and matrisome and CAF signatures, as well as the negative correlation between the risk score and activated CD8+ T cells signatures, we speculated activated CD8+ T cells signatures was negatively correlated with the matrisome and CAF signatures. Interestingly, our results showed that the expression of activated CD8+ T cell signatures were negatively correlated with both matrisome and CAF signatures (Additional file 6: Fig. S6C). Taken together, these results strongly suggested the tumor immunosuppressive microenvironment might contributed to the worse prognosis of the patients with STAD in high-risk group, which needs to be validated in further study.

29Nowadays, immune checkpoint inhibitors were studied and well applied in cancer immunotherapy. In the present study. Our results indicated that 43 immune checkpoints were considerably modulated in high-risk group (Fig. 6E). In addition, risk score was significantly positively correlated with the expression level of 7 immune checkpoints, including CD200, NRP1, TNFSF4, B7-H3, TNFSF18, LAIR1 and OX40 (r > 0.2, Additional file 6: Fig. S6D). Currently, the inhibitors for PD-1 and CTLA-4 are research hotspots in the treatment of advanced STAD. As shown in Fig. 6C, the expressions of PD-1, PD-L1 and CTLA-4 were significantly down-regulated in the high-risk group. Consistently, risk score was significantly negatively correlated with the expressions of PD1, PD-L1 and CTLA-4, respectively (Additional file 6: Fig. S6D).

30Given the above results, we further used the TIDE algorithm to evaluate the ability of risk score in predicting the responses to immunotherapy in STAD. Our results showed that risk score had a significantly positive correlation with TIDE score (Additional file 7: Fig. S7), indicating that patients in low-risk group received better response to immunotherapy. The immunotherapy response rate in high-risk group (32.28%) was significantly lower than that in low-risk group (66.67%) (Fig. 7A).

31Next, PRJEB25780 cohort (PD-L1 inhibitor treatment for 45 patients with advanced gastric cancer) was applied to validate whether the mitochondrial-related risk signature could accurately predict the responses to immunotherapy for patients with STAD. Consistent with the prediction results by TIDE, the immunotherapy respond rate in high-risk group (13.04%) was significantly lower than that in low-risk group (40.91%) in PRJEB25780 validation cohort (Fig. 7B).

32As shown in Fig. 7C, the immunotherapy response rate in the low-immune subgroup (34.59%) was remarkably lower than that in the high-immune subgroup (64.56%). Interestingly, the immunotherapy response rate in low-risk group (45.56%) was remarkably higher than that in high-risk group (20.29%) in the subgroup with low-immune score. Moreover, the immunotherapy response rate in low-risk group (94.20%) was significantly higher than that in high-risk group (41.57%) in the subgroup with high-immune score, strongly suggesting that combined risk score and immune score was a robust indicator to predict the responses to immunotherapy in STAD (Fig. 7D). In addition, the immunotherapy response rate in the high stromal group (37.97%) was significantly lower than that in the low stromal group (61.01%) (Fig. 7E). In the low stromal subgroup, the immunotherapy response rates were 65.35% and 53.45% in low risk and high risk subgroup, respectively, which were similar to the low-stromal group (61.01%), indicating combined risk score and stromal score was not better than stromal score alone in predicting the response to immunotherapy in STAD patients with low-stromal score. However, the immunotherapy response rate in the low-risk + high-stromal group (68.97%) was significantly higher than that in high-stromal group (37.97%), whereas the immunotherapy response rate in the high-risk + high-stromal group (20.00%) was significantly lower than that in the high-stromal group (37.97%), strongly suggesting that combined risk score and stromal score can more accurately predict response to immunotherapy in STAD patients with high-stromal score (Fig. 7F).

33The phenotype for microsatellite instability–high (MSI-H) is a distinct tumor subclass that is highly susceptible to immunotherapy. Consistent with the previous studies, the immunotherapy response rate MSI-H subgroup (72.73%) was remarkably higher than that in the MSS subgroup (46.30%) and MSI-L subgroup (36.96%) (all P = 0.0005, Fig. 7G). Furthermore, our results also showed that the immunotherapy response rate in low-risk group (65.35%) was remarkably higher than that in high-risk group (29.57%) in the subgroup with MSS. The immunotherapy response rate in low-risk group (57.89%) was also significantly higher than that in high-risk group (22.22%) in the subgroup with MSI-L. The above results strongly suggested that combination of risk-score and MSS/MSI-L can be used as a robust indicator to predict the response to immunotherapy in STAD (Fig. 7H).

34Progressive accumulation of mutations throughout life can lead to cancer. Genome sequencing has revolutionized our understanding of somatic mutation in cancer, providing a detailed view of the mutational processes and genes that drive cancer. Therefore, we mapped the mutation landscape of STAD in both the high risk and low risk groups, and analyzed the relationship between risk score and mutation profile. As shown in Fig. 8A, the top 20 high-frequency mutated genes in high-risk and low-risk group were presented. TP53, TTN, MUC16, LRP1B, SYNE1, CSMD3, FAT4, OBSCN, ARID1A, FLG, CSMD1, DNAH5, SPTA1, PCLO and RYR2 were the common high-frequency mutation genes in both groups. In top 5 mutant genes, the mutation rates of TTN and MUC16 were significantly decreased in high-risk group (Fig. 8B).

35Besides, accumulating evidences supported that TMB functioned as a potentially predictive biomarker for multiple applications, including the biomarker for response to immunotherapy in malignancies. Our data showed that TMB in high-risk group was significantly lower than that in low-risk group, and the risk score had a significant negative correlation with TMB in STAD (Fig. 8C, D). Interestingly, the higher TMB tended to have a better OS compared with the lower TMB, but without a statistically significant difference (P = 0.153, Fig. 8E). Moreover, low-risk group tended to have a better OS compared with high-risk group in high TMB subgroup. The similar results were obtained in the low TMB subgroup, but without a statistically significant difference (P = 0.099, Fig. 8E). Taken together, these results strongly suggested that combination of risk score and TMB might be a valuable biomarker for predicting the prognosis for STAD patients (Fig. 8E).

36To find the potency of risk score as an index for predicting the response to drugs (including chemotherapy, targeted therapy, and immunotherapy) in STAD, we inferred the IC50 value of the 138 drugs in TCGA-STAD patients. We found that patients in high-risk group might be more sensitive to rapamycin, PD-0325901, dasatinib, et al., whereas patients in low-risk group might be more sensitive to AZD7762, CEP-701, methotrexate, et al., which could provide a reliable reference for clinical treatment (Fig. 9). Detailed information on the top 10 sensitive drugs for the high-risk and low-risk subgroup was illustrated in Tables 4 and 5.

37To verify the expressions of hub genes in STAD samples, qRT-PCR was conducted on 10 pairs of STAD tumor and paired adjacent normal tissues. As shown in Fig. 10A, the expressions of NOX4 and FKBP10 were significantly higher in STAD tumor tissues than those in the paired adjacent normal tissues, respectively (all P < 0.05). The expressions of ALDH3A2 and MAOA were significantly lower in STAD tumor tissues than those in the paired adjacent normal tissues, respectively (all P < 0.05, Fig. 10A). We also verified the expressions of the four hub genes in a human normal gastric epithelial cell line GES-1 and human gastric cancer cell lines through qRT-PCR. Our results showed that the expressions of NOX4 and FKBP10 were significantly higher than in the normal gastric epithelial cell line GSE-1 (all P < 0.05, Fig. 10B), while the expressions of ALDH3A2 and MAOA were significantly lower than in the normal gastric epithelial cell line GSE-1 (all P < 0.05, Fig. 10B). These results supported our hypothesis and provide solid evidence for the rationality of choosing these four genes for prognostic model construction (see Fig. 11).

38GC remains a frequent cancer worldwide with high incidence and mortality globally. Effective biomarkers are still missing even though 3 biomarkers (HER2, MSI-H and PD-L1) have been proven to predict the responses of targeted therapy in GC. Therefore, identifying more effective biomarkers for targeted therapy and prognosis prediction is highly demand. Mitochondria were important pharmacological targets due to their critical roles in cell proliferation and death. The mitochondrial energy metabolisms are now known to be reprogrammed to meet the challenges of high energy demand, with better use of available fuels for malignant cell growth and migration. Thus, mitochondria play a vital and multifunctional role in tumor occurrence and development, and targeting mitochondria provides therapeutic opportunities. A growing body of research showed that mitochondrial-related genes can be used as biomarkers for the diagnosis and treatment of malignancies.

39NOX4 has been identified as a biomarker and therapeutic target for a variety of human cancers. NOX4 was upregulated in pancreatic cancer and was involved in the development of pancreatic cancer by promoting cell proliferation, regulating cell metabolism, and mediating angiogenesis, suggesting NOX4 was a potential therapeutic target for pancreatic cancer. Up-regulation of NOX4 predicted worse prognosis and accelerated tumor growth in colorectal carcinoma. Moreover, it has been shown that NOX4 recruited M2-macrophages via ROS/PI3K signaling pathway-dependent cytokines production, thus contributing to the cell division in NSCLC. Consistent with the previous studies, NOX4 in our model was significantly up-regulated in STAD, and high level of NOX4 was associated with worse prognosis in patients with STAD.

40It has been reported that a cancer-specific molecular mechanism for NSCLC was related with FKBP10-dependent protein translation. The expression of FKBP10 was positive in cancer lesions. Li et al. reported that FKBP10 silencing decreased the expression of integrin αV and integrin α6, and P-AKT, suggesting that FKBP10 might promote metastasis. FKBP10 was up-regulated in GC tissues and might be a reliable therapeutic target in GC. Consistent with the above studies, our results showed that the high expression of FKBP10 was related with high risk, and predicated poor prognosis of STAD patients. Consistent with the previous studies, FKBP10 in our model was significantly up-regulated in STAD, and high level of FKBP10 was associated with worse prognosis in patients with STAD.

41It has been shown that ALDH3A2 was overexpressed in low-grade GC compared with high-grade GC, and patients with low expression of ALDH3A2 had worse OS than those with high ALDH3A2 expression. ALDH3A2 was reported as a reliable biomarker for the immunotherapy, as well as an independent predictor for the prognosis of GC. In renal clear cell carcinoma, low level of ALDH3A2 was related with shorter survival. Consistently, ALDH3A2 was significantly down-regulated in STAD. However, ALDH3A2 couldn’t effectively predict the prognosis for patients with STAD, and this inconsistent effects of ALDH3A2 on prognosis may be caused by the different inclusion and exclusion criteria across different study, which needs more comprehensive investigation.

42MAOA exerted different biological effects in different tumors. MAOA was found to be involved in mitochondrial dysfunction, and promoted malignant growth and metastasis in gastric cancer. MAOA promoted prostate cancer progression by increasing cell growth and cancer stem cells, which suggested that MAOA might be a potential therapeutic target for the treatment of prostate cancer. In the present study, MAOA was significantly down-regulated in STAD. However, the high level of MAOA was associated with worse prognosis in patients with STAD, which need to be validated in more samples in further study.

43Currently, many biomarkers were applied for prognostic prediction of GC, such as NOX4, FKBP10 and ALDH3A2, but most of them are studied for a single biomarker. Increasing evidences indicated that prognostic model constructed by multi-genes as a prognostic index was more comprehensive and effective than single gene in kinds of malignancies. For instance, Nie et al. constructed a GC prognosis model based on metabolic signature, which were mainly related to the dysregulation of the metabolic microenvironment. Wu et al. constructed a immune-related prognostic model. As the dysfunction and dysregulation of mitochondria have been associated with cancer, we constructed a STAD prognostic model based on mitochondrial-related genes, which could effectively predict the prognosis for patients with STAD.

44The DEGs between the high-risk and low-risk groups were mainly enriched to the extracellular matrix (ECM) and focal adhesion pathway. ECM accumulation was a classical characteristic feature of tumors, and a higher ECM content predicted a poorer prognosis in a broad range of cancer types. The TME is a composition of cancer cells, non-cancerous stromal cells, soluble growth factors, cytokines, proteases, and ECM, which provides essential signals for tumor survival, growth, and acquisition of invasiveness, while hindering antitumor immunity. Fibroblasts constitute one of the most vital cells in the stroma and turn into cancer-associated fibroblasts (CAFs) in TME. CAFs not only play active roles in tumorigenesis and progression both by soluble factors and direct cell-to-cell contact, but also sculpt TME by suppressing anti-tumor immune responses or by recruiting immunosuppressive cells. Matrisome defined as the compendium of genes encoding core ECM proteins, or structural component of the ECM. Consistently, our results indicated that risk score was positively correlated with the CAF signature, ECM signature, and Matrisome signatures. Moreover, risk score had a positive correlation with stromal score and a negative correlation with tumor purity, which could be on behalf of higher infiltration degrees of stromal cells in the TME of the high-risk group.

45The immune cells play essential roles in TME. The success of cancer immunotherapy relies on the comprehensive understanding of the tumor microenvironment and immune evasion mechanisms in which the tumor, stroma, and infiltrating immune cells coordinated in a complex network. The main benefit of immunotherapy is to generate memory CD8+ T cells for sustained protection against metastasis and preventing recurrence of the disease. Active immune cells could enter into the tumor parenchyma and perform their anti-tumor function. Therefore, the ultimate goal of immunotherapy is to convert an immunodorminant TME into an immunostimulatory TME, which allows the immune system to clear tumor lesions. Treg cells are involved in tumor progression by inhibiting antitumor immunity. High Treg cell infiltration in the TME was involved in unfavorable prognosis in patients with various types of cancer. M2-polarized macrophages, commonly deemed tumor-associated macrophages (TAMs), were contributors to many pro-tumorigenic outcomes in cancer. Macrophage type 2 (M2) cells, and Tregs cells could make immunologic barriers against CD8+ T cell‐mediated antitumor immune responses. Our results indicated that naive B cells, regulatory T cell (Tregs), M0 macrophage and M2 macrophage were significantly enriched in high-risk group, whereas CD8+ T cells and resting CD4+ T cells were remarkably enriched in low-risk group. Taken together, these results suggested that high-risk group was suffused with immunosuppressive cells such as Tregs, M2 macrophages, producing the immunosuppressive microenvironment to hamper CD8+ T cells-mediated eradication for tumor cells.

46Monoclonal antibodies against immunological checkpoint molecules provided a vast breakthrough in cancer therapeutics. For instance, PD-1/PD-L1 and CTLA-4 inhibitors showed promising therapeutic outcomes. In our study, high-risk group had a considerably lower rate of immunotherapy response than that in the low-risk group, which was consistent with the expression levels of PD-1/PD-L1 and CTLA-4 in the high-risk and low-risk groups. Interestingly, combination of risk score and immune score, or stromal score can more accurately predict the responsiveness to immunotherapy of patients with STAD. These findings further demonstrated the effectiveness of the risk score as a biomarker in predicating the response to immunotherapy. The TMB was associated with the formation of neoantigens which activated antitumor immunity, which was a reliable biomarker to predict the response to PD-L1 therapy. In our study, the TMB in high-risk group was lower than that in low-risk group, which strongly suggested that the lower response rate to immunotherapy in high-risk group might be due to the lower TMB.

47In addition, risk score might be helpful in screening the therapeutics drugs for patients with STAD. For instance, an independent study showed that high-risk STAD patients showed higher sensitivity to the chemotherapy agents, including rapamycin. Another study found that methotrexate is suitable to inhibit the function of Early B-cell factors (EBFs) in gastric cancer. In the present study, rapamycin, PD-0325901 and dasatinib were found to be more effective for patients in the high-risk group, whereas AZD7762, CEP-701 and methotrexate were predicted to be more effective for patients in the low-risk group. However, the toxicities of the screened drugs was uncertain. For instance, due to the unpredictable cardiac toxicity, the development of AZD7762 was not going forward in patients with advanced solid tumors. 5-fluorouracil, doxorubicin, high-dose methotrexate (FAMTX) schedule was reported to be active in advanced gastric cancer, and the main toxicity was myelosuppression. Oxidative stress is a component of many diseases, including cancer. Although numerous small molecule drugs evaluated as antioxidants have exhibited potential therapeutic ability in preclinical studies, results from clinical trial was disappointed. A greater understanding of the pharmacological mechanisms through which anti oxidative drug act might provide a rational usage would lead to greater therapeutic success in malignancies.

48Our research has some unique advantages. First, the combination of multigene has robust predictive capability for cancer prognosis than single genes. An integrated mitochondria-related gene prognostic risk model would play more vital roles in the diagnosis and prognosis of STAD patients. Second, the results of the study provide us with a more accessible method to determine whether patients belong to the high- or low-risk group, which is simple and feasible. In addition, we evaluated its predictive value, chemotherapy efficacy, immunotherapy efficacy and immune cell infiltration for patients with STAD, which could provide individualized management and treatment for STAD patients.

49Although our findings in this study have important clinical consequences, there are still some limitations. Firstly, this is a retrospective study, and an independent prospective cohort is needed to verify the risk model constructed in this study. Secondly, this study heavily relied on datasets and computational predictions while validation component is poor. More experimental studies are needed to validate in further studies. Finally, the carcinogenic effects of the prognostic genes in the model and the mechanisms of interaction between prognostic genes and mitochondrial dysfunction in STAD are mainly unknown and need to be further explored. Based on above information, our future direction will focus on three aspects: (1) Applying mouse model to verify our current hypothesis; (2) Collecting gastric cancer cases and clinical information to validate the risk score model, and compare it with the gold standard of clinical diagnosis, and ensure that the risk model constructed in the present study can be applicable for clinical practice; (3) To screen more novel genes associated with mitochondria in more datasets.

50We established a STAD patient risk score model including NOX4, FKBP10, ALDH3A2, and MAOA. Functionally, the risk score was highly correlated to the TME and immune cell infiltration of STAD patients. Combined analysis for risk score and stromal score, or immune score, or MSS/MSI can predict the response to immunotherapy more accurately than single index in STAD. Regarding drug sensitivity, patients in high-risk group was more sensitive to rapamycin, PD-0325901 and dasatinib, whereas patients in low-risk group was more sensitive to AZD7762, CEP-701 and methotrexate. Taken together, our mitochondrial-related risk model could be a reliable prognostic biomarker for personalized treatment of STAD patients.