Integrated analysis of the underlying immunomodulator and survival signature of STON1 gene in kidney renal clear cell carcinoma

Axiu Zheng Guangdong Medical University Jianrong Bai Guangdong Medical University Yanping Ha Guangdong Medical University Bingshu Wang Hainan Medical University Yuan Zou Guangdong Medical University Caixia Wu Guangdong Medical University Jingci Xing Guangdong Medical University Shaojiang Zheng Hainan Medical University Zhihua Shen Guangdong Medical University Botao Luo Guangdong Medical University Wei Jie (  wei_jie@hainmc.edu.cn ) Hainan Medical University https://orcid.org/0000-0002-1869-1458

including colitis/diarrhea, hepatitis, pneumonitis, myocarditis, myositis, encephalitis, and myasthenia gravis [8]. Thus, there is an urgent need to identify novel biomarkers to guide prognosis prediction and clinical management.
Stonins are an evolutionally conserved family that mediates the recovery and circulation of vesicles at neuromuscular junctions, and consists of Stonin1 (STON1, also known as germ line-speci c transcription factor) and Stonin2 (STON2) [9]. The human STON1 protein consists of 735 amino acids and has a predicted molecular mass of 83 kDa [10]. Additionally, the level of STON1 gene methylation was veri ed to be a prognostic marker for the progression of obesity and in the development of personalized treatment [11]. High expression of STON1 as a potential mRNA biomarker in sperm was identi ed in cigarette smoke-exposed individuals [12]. Furthermore, STON1 has an irreplaceable role in regulating the endocytosis of neuron-glial antigen 2 (NG2). NG2 appears to be oncogene of glioblastoma and melanoma [13], but the potential functions of STON1 in mediating speci c molecular mechanisms in various carcinomas remain entirely unexplored.
Over the past few years, various prognostic models have been constructed for the underlying mechanism of KIRC based on metabolism-related genes, long non-coding RNAs, alternative splicing variants, immunerelated genes, RNA binding proteins, microRNA related genes, and autophagy-related genes [14][15][16][17][18][19][20]. In our current study, the independent prognostic value of STON1 was validated using online bioinformatics databases. Subsequently, we identi ed several STON1-related immunomodulators with prognostic indication and used them to build a risk score model. Using this model, risk score was estimated in KIRC patients. Finally, a nomogram integrating clinicopathological parameters and risk score was established to maximize the convenience and e ciency for physicians in predicting the overall survival (OS) time of patients with KIRC.
Analysis of the relationships between STON1, clinical phenotype, and prognosis in KIRC RNA sequencing and related clinicopathological data were downloaded from TCGA using UCSC Xena (https://xena.ucsc.edu/), an online explorer that allows users to explore functional genomic data sets for correlations between genomic and phenotypic variables. A total of 491 patients were enrolled in the cohort after exclusion criteria were applied. The exclusion criteria were as follows: 1) patients without follow-up records; 2) patients without a diagnosis of KIRC; 3) primary tumor could not be evaluated; and 4) stage and grade were not reported. Finally, 491 patients were divided into high and low STON1 expression groups according to the median expression value of STON1. The OncoLnc database (http://www.oncolnc.org/) is a platform for immediate survival analysis with TCGA data [24]. The UALCAN database (http://ualcan.path.uab.edu/) can provide graphs and plots depicting expression pro les and patient survival information for protein-coding, miRNA-coding, and lincRNA-coding genes [23]. Kaplan-Meier Plotter (https://kmplot.com/analysis/) is an accessible online website with the purpose of identifying survival biomarkers [25]. We also analyzed the prognostic value of STON1 expression in KIRC patients using the Gene Expression Pro ling Interactive Analysis (GEPIA) browser (http://gepia2.cancer- Correlation of tumor mutational burden, microsatellite instability, and STON1 The raw data of mRNA pro les in KIRC were downloaded from TCGA dataset. "R.utils", "rjson", and "edgeR" R packages were used to process the raw expression matrix. A marginal scatter plot was acquired by the "ggpubr" R package for the KIRC mRNA pro le. The relationship between STON1 and tumor mutational burden (TMB) was displayed by the Assistant for Clinical Bioinformatics, an open integrated database (https://www.aclbi.com/).

Analysis of STON1 mRNA levels, immune cells, and immunomodulators in KIRC
The TIMER database was employed to display the correlation between gene copy numbers of STON1 and immune cells [27]. The Tumor and Immune System Interaction Database (TISIDB) (http://cis.hku.hk/TISIDB/) was used to evaluate the relationship between STON1 and 28 immune cell types, as well as 69 types of immunomodulators [28]. Heatmaps and bubble graphs produced by the "ggplot2" R package intuitively presented P-values and Rho values for immune cells and STON1 expression, while a violin diagram showed the relevance of STON1 and immune subtype in KIRC.

GO, KEGG analysis, and Integrated Network in KIRC
A total of 31 immunoregulators with statistical signi cance were uploaded into cBioPortal (http://www.cbioportal.org/) to identify genes that are co-expressed with immunomodulators [29]. Using the screened 57 co-expressed genes and 31 immunoregulators, gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis were performed on the DAVID online tool [30] (https://david.ncifcrf.gov/). A P-value of less than 0.05 was considered statistically signi cant. Results are presented as a bubble plot generated using the R package "ggplot2". A protein interaction network was constructed by Cytoscape software, an open source for visualizing complex networks [31].
Univariate and multivariate Cox regression, Kaplan-Meier survival analysis, and establishment of a STON1-releated immunoregulator prognostic signature in KIRC The processed RNA-sequencing data of KIRC, which contained 504 samples, was randomly divided into a training set (n=328) and an internal validation set (n=176) at a proportion of 65% via the "caret" Rpackage. Using the "survminer" and "survival" R packages, univariate Cox analysis was applied to evaluate the correlation between immunomodulators and OS in the training set. Next, we entered 15 immunomodulators with statistical signi cance into a stepwise multivariate Cox regression model with the purpose of identifying independent prognostic biomarkers.
Risk score was calculated using the following formula: Risk score= h 0 (t) × exp (∑coe cient of gene (i) × expression of gene (i)) Subsequently, the patients in the training set were strati ed into a high-risk group and a low-risk group according to the median risk score. The Kaplan-Meier survival curve was applied to compare the survival outcomes among the two groups, employing R packages of survival, as well as the log rank test.
Furthermore, we used a time-dependent ROC curve to evaluate the model accuracy with the "timeROC" R package. Using the "ggrisk" package of R software, the risk graph presented a landscape of the STON1releated immunoregulator prognostic signature.

Validation of STON1-related immunoregulator prognostic signature
To validate the stability of our signature, the risk score was calculated using the coe cient described above. Next, the validation set was also divided into high-risk and low-risk groups based on the median risk score. The Kaplan-Meier curve and log rank test presented the survival difference between the two groups, as well as a risk graph. A time-dependent receiver operating characteristic (ROC) curve further validated the accuracy of the STON1-releated immunoregulator prognostic signature.
Construction and estimation of a nomogram based on risk score and clinicopathological parameters Univariate Cox regression analysis was performed on the training set to screen variables with statistical signi cance using the "survival" and "survminer" packages of R software. Then, stepwise multivariate Cox regression analysis was executed with the screened variables that were derived from the univariate Cox regression analysis. The results were displayed in a forest plot using the "forest plot" of R package. An integrated nomogram was established with the training set based on stepwise multivariate Cox regression analysis to predict the survival probability at 1, 3, and 5 years, employing the "rms" package of R software. Subsequently, time-dependent ROC curves and calibration curves were used to measure the accuracy and stability of the nomogram in the training set and validation set, as well as the c-index. Although the ROC curve and calibration curve could thoroughly predict the accuracy and stability of the nomogram model, it could not reveal the bene t rate of patients in the model. Decision Curve Analysis (DCA) is a statistical method for evaluating the degree of bene t that KIRC patients can obtain. It is therefore an indispensable validation tool in addition to time-dependent ROC curves and calibration curves [32]. Using the "ggDCA" package in R software, DCA curves were plotted to validate our nomogram better than other clinical parameters for clinical bene t in the training and validation sets.

Statistical analysis
All mRNA expression was normalized by log2 transition. A columnar scatter plot was created by GraphPad Prism (version 8, GraphPad Software, CA, USA). Chi-squared test was performed using IBM SPSS Statistics (version 25) to analyze the association between STON1 mRNA level and clinicopathological parameters. The correlation analysis between STON1, mismatch system, and TMB was performed using Pearson's test in RStudio software (version 4.0.3). For the Kaplan-Meier curve, logrank test and Cox proportional hazard regression model analysis were performed using RStudio software (version 4.0.3). P < 0.05 was considered statistically signi cant.

Correlation between STON1 expression and clinicopathological characteristics of KIRC patients
The KIRC GDC Expression Matrix was downloaded from the UCSC Xena, which is one of the most comprehensive clinical oncology databases [33]. A total of 491 KIRC patients were included in the cohort after excluding non-KIRC diagnoses and missing clinical information, and were classi ed into a high expression group (n = 246) and a low expression group (n = 245) according to the median value of STON1. Notably, STON1 mRNA levels were related to grade (P = 0.022), TNM stage (P = 0.039), distant metastasis (P = 0.019), and vital status (P = 5.46 e−6 ). However, no correlation was found between STON1 mRNA levels and age (P = 0.136), sex (P = 0.949), invasion depth (P = 0.074), and lymph node metastasis (P = 0.385) ( Table 1).

STON1 was a protective prognostic factor in KIRC
To further analyze the prognostic value of STON1, we performed a survival analysis using four oncological databases. The results indicated that STON1 expression levels were strongly associated with OS in KIRC. Patients with high expression of STON1 had a better outcome compared with those with low levels of STON1 in the OncoLnc database (log rank P = 9.78 e−8 ) (Figure 2a). Similar results were acquired from other oncology databases, including the UALCAN database, the Kaplan-Meier database, and GEPIA2 with TCGA data (Figure 2b, 2c, 2d).
Correlations of STON1 levels with tumor mutation burden and mismatch repair genes in KIRC Given the sensitive connection between TMB, microsatellite instability (MSI), and immunotherapeutic response, we focused on their relationship with STON1 expression levels. Marginal scatter plots revealed that alterations in STON1 expression were always accompanied by changes in mismatch repair proteins (MLH1, PMS2, MSH2, MSH6, and EpCAM). Moreover, the four mismatch proteins seem to be positively associated with STON1 ( Figure 3a, P < 2.2 e−16 ). Additionally, STON1 mRNA level was negatively associated with the TMB score (P = 0.012) (Figure 3b).

Association between STON1 expression and immune cells in KIRC
STON1 gene copy numbers were accompanied by changes in immune cell in ltration levels, including that of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells (Figure 4a). Subsequently, we examined the relationship between STON1 expression levels and the in ltration levels of 26 immune-related cell types. Our data illustrated that the abundance of different subtypes of immune cells was signi cantly positively (effector memory CD8 T cells, effector memory CD4 T cells, type 1 T helper cells, type 2 T helper cells, regulatory T cells, memory B cells, natural killer cells, immature dendritic cells, eosinophils, mast cells, and neutrophils) or negatively (activated CD8 T cells, CD56 bright natural killer cells, CD56 dim natural killer cells, myeloid-derived suppressor cells, activated dendritic cells, and plasmacytoid dendritic cells) related to STON1 expression (Figure 4b, 4C). These results demonstrated that STON1 was positively linked with most immune cells. Furthermore, STON1 mRNA expression was obviously highest in the C3 subtype (in ammatory) (Figure 4d), which was reported to have the best prognosis in KIRC [34]. This result indicated that STON1 may have a protective role in KIRC by exhibiting upregulated expression in the C3 immune subtype.

STON1 may mediate immune response through immunomodulators
Thirty-one immunomodulators were obviously related to STON1 in KIRC, including 21 immunostimulators (C10orf54, CD27, CD80, CXCL12, CXCR4, ENTPD1, HHLA2, IL6, IL6R, LTA, NT5E, PVR, RAET1E, TMIGD2,  TNFRSF14, TNFRSF17, TNFRSF18, TNFRSF25, TNFSF14, TNFSF15, and TNFSF9) (Figure 5a) and 10 immunoinhibitors (ADORA2A, CD274, CTLA4, IDO1, KDR, LAG3, PDCD1, PDCD1LG2, PVRL2, and TGFBR1) (Figure 5b). We subsequently explored the possible signaling pathways through which STON1 may modulate the immune response in KIRC using GO and KEGG enrichment analyses. Thirty-one immunomodulators and 57 closely related genes were uploaded onto the DAVID. Biological process analysis revealed that these genes were enriched in the tumor necrosis factor (TNF)-mediated signaling pathway, T cell co-stimulation, and immune response. Molecular functional analysis revealed that they were enriched in TNF receptor binding, TNF-activated receptor activity, and cytokine activity. Cellular component analysis indicated that STON1 was mainly enriched in the plasma membrane (Figure 6a). KEGG enrichment analysis further con rmed that these genes are tightly related to cytokine-cytokine receptor interactions (Figure 6b). The associations between 31 immunomodulators and 57 closely related genes are presented in a protein-to-protein interaction network (Figure 6c) Prognostic value of STON1-related immunoregulators in KIRC In total, 504 patients were randomly divided into a training set (n = 328) and a validation set (n = 176) at a ratio of 65%-35%, which was downloaded from GDC. To explore the prognostic value of STON1-related immunomodulators, univariate Cox analysis was rst performed on the training set, which identi ed 15 of the abovementioned immunoregulators with statistical signi cance ( Table 2). We then input these 15 immunomodulators into a stepwise multivariate Cox regression model using the training set. CTLA4 (P = 0.0177), KDR (P = 0.0135), PDCD1 (P = 0.024), and HHLA2 (P < 0.001) were independent prognostic factors, as shown in Figure 7a. Here, risk scores were obtained using the following formula:  Figure 9a). Meanwhile, multivariate Cox proportional hazards regression analysis was implemented to assist with the identi cation of independent prognostic factors. The independent prognostic factors were risk score (HR = 1.197, 95% CI = 1.127-1.271, P = 6.71 e−9 ), age (HR = 1.038, 95% CI = 1.020-1.057, P = 3.52 e−5 ), and stage (HR = 1.699, 95% CI = 1.420-2.033, P = 6.76 e−9 ) (Figure 9b). Therefore, a comprehensive nomogram prognostic model was constructed according to all independent prognostic factors screened by multivariate Cox proportional hazards regression analysis of the training set to predict the probability of 1-, 3-, and 5-year OS (Figure 9c). In the training set, AUC of 1-, 3-, and 5-year OS were 0.889, 0.804, 0.790, respectively (Figure 9d). In the validation set, AUC of 1-, 3-, and 5-year OS were 0.882, 0.847, 0.809, respectively (Figure 9e). The AUC of the training set and validation set illustrated our model with good predictive ability for OS. Furthermore, the calibration curves, which re ected the consistency between the predicted risk and the actual risk, showed that the nomogram-predicted probability was well-matched with the ideal line for 1-, 3-,5-year OS in the training set and validation set (Figure 10a). The same methods were applied to the validation set to further validate the authenticity of the above results ( Figure  10b). The concordance probability (C-index) can be used to measure and compare the discriminative power of a risk prediction model. The C-indexes of our nomogram model were 0.783 and 0.795 in the training set and validation set, respectively. Moreover, we also discovered that our nomogram prognostic model showed a better net bene t for predicting OS compared with age, stage, and risk score in the training set and validation set (Figure 10c, 10d).

Discussion
Although early-stage KIRC patients can be cured by partial or complete nephrectomy with a better outcome, approximately 90% of KIRC patients still suffer from recurrence and metastasis [1,5]. Recently, immune checkpoint blockers, such as nivolumab plus ipilimumab combination immunotherapy, have shown promise in KIRC patients with rapid progression in clinical individualized treatment based on targets of the tumor microenvironment [35], but a considerable number of patients displayed a low or absent immune response [36]. Therefore, further identi cation of novel immunotherapeutic agents and biomarkers to guide clinical therapy and predict survival of KIRC patients remains challenging. In our current study, we focused on STON1, a protein-coding gene involved in vesicle transport. We demonstrated that STON1 was decreased in most human cancer types including KIRC, and was strongly associated with grade, stage, distant metastasis, and vital status in KIRC. Furthermore, STON1 overexpression bodes well for prognosis. The above results indicated that STON1 may as a tumor suppressor and inhibit tumorigenesis and progression. STON1 was identi ed as target gene at the rs13405728 locus in polycystic ovary syndrome [37]. However, to date, no reports have published the underlying functions of STON1 in cancer. A previous study demonstrated that alterations in focal adhesion dynamics, cellular motility, and signaling were induced by the absence of STON1. The above processes are intimately correlated with tumor migration. Moreover, as a speci c adaptor for NG2 endocytosis, STON1 de ciency resulted in NG2 accumulation on the cell surface, thus promoting tumor deterioration [13]. Our results indicate that STON1 may inhibit tumor growth and progression in KIRC.
Mismatch repair genes are related to the human mismatch repair response, which is responsible for the repair of base mismatches that occur during DNA replication [38]. Tumor mutational burden is a novel indicator of mutation quantity [39]. A previous review demonstrated that tumors with mismatch repair possess the capacity to increase the mutational burden, resulting in increased immune checkpoint protein expression including PD-1 and PD-L1. Furthermore, the immunogenic neoepitopes generated by the imbalance of the mismatch system signi cantly improved the immune response rate [40]. Li et al.
revealed that several microsatellite site alterations indicate poor survival of KIRC patients [41], while Chalmers et al. stated that the median mutation per MB was 2.7 in renal clear cell carcinoma, and microsatellite instability could increase TMB to some extent [42]. Interestingly, as our study revealed, the mRNA level of STON1 was positively associated with mismatch repair (MMR) proteins (MLH1, PMS2, MSH2, MSH6, and EpCAM). Therefore, we speculated that TMB aggravation is partly due to increased immunogenic neoepitopes created by the mismatch repair system in the absence of STON1 in KIRC. Our subsequent results corresponded with previous studies and reinforced our assumption. As the STON1 expression decreased, the TMB score increased. Hence, STON1 may be a candidate prognostic biomarker for evaluating the e cacy of immunotherapy responses and further experiments should be conducted to con rm this theory.
KIRC is one of the most immune-in ltrated cancers among all human solid cancer types [43]. The emergence of immune checkpoint targets provides a new direction for tumor immunotherapy. In recent years, CTLA4 inhibitors, PD-l inhibitors, PD-L1 inhibitors, and PD-L1 inhibitors plus tyrosine kinase inhibitor (TKI) in progressed KIRC patients were approved by the FDA and achieved great curative effect [44]. We conducted a comprehensive bioinformatics analysis to reveal the details of the STON1associated immune response in KIRC. The results of the TIMER database analysis indicated that armlevel deletion of STON1 was strongly associated with the in ltration level of six immune cell types: B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and dendritic cells. Furthermore, our study revealed that the high abundance of most immune cells was positively associated with STON1 mRNA expression in KIRC, including regulatory T cells (Tregs), type 1 T helper cells (Th1), type 2 T helper cells (Th2), effector memory CD4 T cells (Tem_CD4), effector memory CD8 T cells (Tem_CD8), natural killer cells (NK), neutrophil cells, memory B cells (Mem_B), mast cells, and immature dendritic cells (iDC).
Tumor-in ltrating Tem_CD4 was associated with downregulation of CD127, a classic marker of exhausted T-cells [45]. Interferon-γ producing Th1 cells secrete TNF, resulting in tumor destruction. Th1 and Th2 cells can not only recruit NK cells and macrophages near to the tumor site to exert antitumor effects, but also act as enhancers to promote tumor-speci c cytotoxic T-lymphocyte (CTL) responses [46]. It is reported that high mast cell density is associated with a longer OS and progression-free survival in metastatic renal cell carcinoma [47]. However, myeloid-derived suppressor cells contribute to immune suppression in the tumor microenvironment[48]. The above results indicated that STON1 may indirectly inhibit tumor progression by regulating most immune cells in the immune microenvironment.
Our violin plot showed that the mRNA level of STON1 in KIRC was highest in the C3 subtype but the lowest in C5 subtype. C1-C6 immune subtypes were previously identi ed from ve representative signatures across multiple tumors. Furthermore, the C3 subtype seems to possess a better outcome than the other immune types [34]. We speculated that higher STON1 expression was remarkably associated with longer OS in C3-dominated tumors. Considering the crucial role of immunomodulators in KIRC, the results obtained from the TISDB database identi ed 21 STON1-associated immune-stimulators and 10 STON1-associated immunoinhibitors. Subsequently, GO terms analysis indicated that TNF receptor binding was highly enriched. Taken together, these results further validated that STON1 is strongly associated with the activation of the tumor immune microenvironment.
Immune prognostic signatures have been reported previously. Xu et al. identi ed seven immune-related genes and constructed a gene signature, of which the AUC value was 0.751. The signature showed moderate prognostic accuracy in predicting survival time [49]. Consistently, we conducted stepwise multivariate Cox regression with STON1-related immunomodulators screened by univariate Cox regression analysis in the training test. A total of seven STON1-related immune genes were screened, and CTLA4, KDR, PDCD1, and HHLA2 displayed a dependent prognostic value in KIRC. All of them had been proven to be tightly associated with human cancer biology. CTLA4 inhibits T cell activation to some extent [50]; thus, anti-CTLA4 antibody immunotherapy can enhance the anti-tumor immune effects of T cells [51]. Wu et al. reported that autophagy induces KDR phosphorylation, which leads to angiogenic mimicry that promotes tumor growth [52]. An anti-PD1 (PDCD1) agent series was approved by the FDA and used in various carcinomas with better immune response [53]. Xu et al. revealed that HHLA2 was downregulated and had protective effects in epithelial ovarian cancer [54]. On the basis of this seven-gene immunoregulator signature, KIRC patients were divided into low-risk and high-risk groups according to the median value of the risk score derived from the signature. The Kaplan-Meier survival curve of the training set and validation set revealed that the high-risk score group had poor OS in KIRC. In our survival predictive system, the AUC values for 1-, 3-, and 5-year OS for both the training set and validation set were greater than 0.70, indicating the relative stability of our signature. As for clinicopathological characteristics, we integrated a nomogram consisting of age, stage, and risk score according to the stepwise regression model in the training set. With our nomogram, clinicians could acquire great

Conclusions
Collectively, lower STON1 expression in KIRC indicates tumor progression and worse survival outcome. STON1 may act as a stimulator to induce the activation of the anti-tumor immune system. Our current study provides new insights into the function of STON1 in the tumor microenvironment. Importantly, our integrated nomogram derived from the risk score of STON1-associated immunoregulators and clinical parameters could be used to predict survival in KIRC patients. A limitation of our study is that it lacked external validation because of missing survival data, and cytological functional studies are required to further verify the speci c role and molecular mechanism of STON1 in KIRC.

Declarations
Ethics approval and consent to participate Not applicable. The data of KIRC patients were obtained from public databases.

Consent for publication
Not applicable.

Availability of data and materials
The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.           Supplementary Files