RNA sequencing (RNA-seq) data of 371 STAD samples and their clinical information were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), while the RNA-seq data of 868 STAD samples (GSE84437, n = 433; GSE62254, n = 300; GSE183136, n = 135) and the corresponding survival information were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Their clinical information was collected in Supplementary Table S1.

First, a univariate Cox regression analysis was performed to investigate the correlation between the gene expression value of SLC family members and the overall survival (OS) of patients with STAD. Subsequently, a stepwise multivariate Cox regression analysis was performed with the genes that significantly affect OS to develop a formula for the risk score. Patients were classified into high- and low-risk groups based on the risk score. Moreover, the prognostic power of the risk score was evaluated using Kaplan-Meier (K-M) survival curves with both TCGA and GEO cohort data. Furthermore, univariate and multivariate Cox regression analyses were performed to validate the independent prognostic value of the risk score.

To explore the functional signaling pathway, GSEA was performed. Gene sets with |normalized enrichment score (NES)|>1 more than 1, p < 0.05, and false discovery rate (FDR) <0.25 were considered to be significantly enriched [23].

Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm was used to estimate the relative proportion of 22 immune cell types in the groups. Single sample GSEA (ssGSEA) was performed to further define the immune and molecular functions and compare the scores of the groups.

The TIDE score, T cell dysfunction score, and T cell exclusion score of patients with STAD from the TCGA dataset were computed online (http://tide.dfci.harvard.edu) after uploading the transcriptome profiles. T cell dysfunction and T cell exclusion are the two primary mechanisms to model tumor immune evasion. TIDE integrates the expression signatures of T cell dysfunction and T cell exclusion, two primary mechanisms of immune evasion, to predict ICI therapy response [24].

All analyses were performed using R software (version 4.1.1, www.r-project.org) and codes are uploaded as Supplementary Datasheet S1. The independent t-test was performed to compare continuous values between the two groups, while the χ² test was performed to examine categorical data. Spearman’s analysis was utilized to check the correlation between the groups. K-M survival analysis was utilized for univariate survival analysis, while the Cox regression model was utilized for multivariate survival analysis. Results with two-sides of p < 0.05 were considered statistically significant, providing credibility for the data analysis.

A total of 58 amino acid transporter-related SLC family genes were enrolled in this study. Univariate Cox analysis was used to evaluate the correlation between gene expression and OS (Table 1). Nine genes were found to be significantly associated with OS. SLC1A7, SLC7A2, SLC7A3, and SLC17A7 contributed to the poorer OS with a hazard ratio (HR) of more than 1. On the contrary, HRs of SLC1A4, SLC1A5, SLC6A9, SLC7A1, and SLC25A15 were less than 1 (Figure 1A). Subsequently, a stepwise multivariate Cox regression analysis was performed on nine genes, and four genes were ultimately confirmed to be significantly associated with prognosis (SLC6A9, SLC7A2, SLC7A3, and SLC25A15) (Supplementary Table S2). We also performed K-M analysis to identify the prediction potential of these genes and utilized the “surv_cutpoint” function to automatically determine optimal cutoff values (SLC6A9: 1.2744, SLC7A2: 0.6353, SLC7A3: 1.3351, SLC25A15: 0.9835). Results showed that SLC6A9 and SLC25A15 were positively correlated with OS, while SLC7A2 and SLC7A3 showed a negative correlation (Figure 1B). Consequently, a risk score formula was constructed, “risk score = 0.2384 * SLC7A2 expression + 0.4544 * SLC7A3 expression − 0.3645 * SLC6A9 expression − 0.3365 * SLC25A15 expression.” Tumor samples had significantly higher risk scores than normal samples (Supplementary Figure S1).

The median risk score was selected as a cut-off to divide all patients into the high- and low-risk groups. Univariate Cox regression analysis showed that age (≥65), TNM stage, T stage, lymphatic metastasis, and risk score were significantly correlated with the prognosis of STAD. Multivariate Cox regression analysis confirmed the independent prognostic value of the risk score (HR = 1.5623, 95% CI: 1.3210–1.8478, p < 0.0001) after adjustment of other clinicopathologic factors (Table 2). K-M analysis demonstrated that the high-risk group had a worse prognosis (Figures 1C,D). Results from the validation cohorts GSE62254 and GSE183136 were consistent with those from the training cohort (Figures 2A,B,E,F). Univariate Cox regression analysis showed that age (≥65), TNM stage, and risk score were significantly correlated with the prognosis of STAD in GSE62254. Multivariate Cox regression analysis further confirmed the independent prognostic value of the risk score (HR = 2.212, 95%CI: 1.58–3.095, p < 0.0001) (Supplementary Table S3). Although the OS difference was not significant at GSE84437 (Figures 2C,D, p = 0.095), a significant difference could be obtained if the upper and lower quartiles were compared (Supplementary Figure S2A, p = 0.008). Moreover, ROC analysis demonstrated that the AUC of TCGA cohort, GSE62254, GSE84437, and GSE183136 cohorts were 0.62, 0.67, 0.57, and 0.61 (Supplementary Figure S2B, p < 0.05).

The baseline clinical factors were similar for the high- and low-risk groups, as shown in the TCGA clinical heatmap (Figure 3A). However, Lauren’s age and Lauren’s classification were significantly different in the two risk groups in the GSE62254 cohort (Supplementary Figure S3A). Patients with intestinal-type gastric cancer in the low-risk group tended to be older (Supplementary Figures S3B,C). The status of the top 20 commonly mutated genes in the TCGA cohort was checked. The most frequent mutation type was missense mutation, followed by multi-hit mutation. The mutation frequency of those genes in the low-risk group was equal to that of the high-risk group (Figure 3B). GSEA was utilized to clarify the pathway enrichment of two risk groups. Gene sets of the low-risk group were enriched in base excision repair, mismatch repair, and apoptosis pathways. Contrarily, the high-risk groups were enriched for cell adhesion molecules, extracellular matrix (ECM) receptor interaction, focal adhesion, and Mitogen-Activated Protein Kinase (MAPK) signaling pathways, which gave the samples more active proliferation and significant aggressiveness (Figure 3C).

A total of 447 differential genes were identified between the two risk groups with log₂FC > 1, and the high-risk group had an enriched distribution of highly expressed genes. KEGG analysis was performed to compare the difference in signaling pathways. Focal adhesion, cell adhesion molecules, ECM receptor interaction, and Wnt signaling pathways were significantly expressed genes in high-risk groups, suggesting the aggressiveness of tumors (Figure 4).

The correlations between the risk score and the proportion of immune cells were examined by the Wilcoxon test. The high-risk group had higher levels of naïve B cells, resting CD4⁺T memory cells, monocytes, resting DCs, and resting mast cells, while the low-risk group had higher levels of activated CD4⁺ T memory cells, follicular T helper cells, resting NK cells, M0 macrophages, M1 macrophages, and neutrophils (Figure 5A). Combined with the results of three GEO cohorts, activated CD4⁺ T memory cells, M0 macrophages, and neutrophils were significantly higher in the low-risk group than in the high-risk group, while resting CD4⁺T memory cells and naïve B cells showed an opposite trend (Figures 5B–D). The immune functions of two risk groups were assessed based on the TCGA cohort (Figure 5E). The high-risk group had higher functional scores of B cells and mast cells and type II interferon (IFN) response. The low-risk group had more significant APC co-inhibition and stronger class I major histocompatibility complex (MHC) molecules, Th1 cells, and Th2 cells. Similarly, patients with the top 25% of risk scores in three GEO cohorts had significantly higher functional scores of B cells and type II IFN response than patients with the lowest 25% of risk scores (Supplementary Figure S4).

Three widely used biomarkers related to ICI therapy response; microsatellite instability (MSI), CD274 expression, and tumor mutation burden (TMB), were tested in two risk groups in the TCGA cohort. We discovered that patients in the low-risk group had higher MSI, CD274 expression, and TMB, suggesting that ICI therapy might be more beneficial for these patients (Figures 6A–C). Additionally, the TIDE score was introduced into our analysis. The low-risk group, as predicted, had lower TIDE scores, which indicated that patients in this group might acquire more benefits from ICI therapy (Figure 6D). T cell dysfunction and exclusion scores, which together made up the TIDE score, were both lower in the low-risk group (Figures 6E,F). The low-risk group had higher levels of other co-inhibitory molecules such as CTLA4, CEACAM1, IDO1, LGALS9, PDCD1, TNFRSF9, TNFRSF14, and TNFRSF18 (Figure 6G). We then used a cohort (NCT#02589496) comprising 44 GC patients who received anti-PD1 therapy with available clinical information to examine the relationship between ICI therapy response and risk score [26]. Unfortunately, due to the small sample size, we only managed to establish a numerical difference rather than a statistical one (Supplementary Figure S5).