Immunogenomic Landscape and Immune-Related Gene-Based Prognostic Signature in Asian Gastric Cancer

Background Asians have the highest incidence of gastric cancer (GC), and the prognosis of Asian GC is poor. Furthermore, the therapeutics for Asian GC is limited because of genetic heterogeneity and screening difficulty at the early stage. This study aimed to develop an immune-related gene (IRG)-based prognostic signature and to explore prognosis-related regulatory mechanism and therapeutic target for Asian GC. Methods To elucidate the prognostic value of IRGs in Asian GC, a comprehensive analysis of IRG expression profiles and overall survival times in 364 Asian GC patients from the Asian Cancer Research Group (ACRG) and The Cancer Genome Atlas (TCGA) databases was performed, and a novel prognostic index was established. To further explore regulatory prognosis mechanisms and therapeutic targets, a tumor immunogenomic landscape analysis, including stromal and immune subcomponents, cell types, panimmune gene sets, and immunomodulatory genes, was performed. Result Our analysis allowed the creation of an optimal risk assessment model, the Asian-specific IRG-based prognostic index (ASIRGPI), which showed a high accuracy in predicting survival in Asian GC. We also developed an ASIRGPI-based nomogram to predict the 3- and 5-year overall survival (OS) of Asian GC patients. The impact of the ASIRGPI on the worse prognosis of Asian GC was possibly related to the stromal component remodeling. Specifically, TGFβ gene sets were significantly associated with the ASIRGPI and worse prognosis. Immunomodulatory gene analysis further revealed that TGFβ1 and EDNRB may be the novel potential therapeutic targets for Asian GC. Conclusions As a tumor microenvironment-relevant gene set-based prognostic signature, the ASIRGPI model provides an effective approach for evaluating the prognosis of Asian GC and may even prolong OS by enabling the selection of individualized therapy with the novel targets.


INTRODUCTION
Gastric cancer (GC) was the fourth most common malignant tumor and the third leading cause of cancer mortality in 2019 according to global cancer statistics (1,2). The morbidity of GC varies greatly across regions, with nearly 60% of GC occurring in East Asia (3), which has the highest cancer-related mortality for GC worldwide. Additionally, studies have also shown that Asian Americans have a higher GC incidence than non-Hispanic whites (4)(5)(6), which indicates that Asian GC (AGC) is associated not only with lifestyle and culture but also with genetics. Previous research (5)(6)(7) reported that the prognosis of AGC is better as it is often diagnosed at earlier tumor stages, at a more distal anatomic site, and at a younger age and receives more aggressive treatment. Considering that the prognosis of GCs depends not only on tumor stage but also on heterogeneous and epigenetic molecular features (8,9), the differences in GC survival patterns and the possible causes of survival disparities among different ethnic groups have not yet been clarified. Thus, elucidating the mechanisms of AGC will offer new insights for prognosis prediction and treatment of GC.
Over the past decade, immunotherapy, especially immune check inhibitors (ICIs), has become a promising treatment strategy (10,11). Chen et al. demonstrated that anti-PD1, anti-PD-L1, and anti-CTLA-4 ICIs could improve some survival endpoints in advanced GCs (12), indicating that immunotherapeutic approaches have promising prospects for long-term and durable remission of GCs. However, the effect of immunotherapy for GC is limited, and there are great individual differences in GC immunotherapy. Additionally, the tumor microenvironment (TME), including the extracellular matrix (ECM), stromal cells, and immunoinfiltrating cells, was found to play significant roles in the progression, metastasis, and therapeutic response of a variety of tumors (13,14). Li et al. developed a TME-based risk score as an independent prognostic factor for GC (15). Studies have also shown that a high M2 macrophage level is related to the status of peritoneal dissemination, angiogenesis, immune evasion, and poor prognosis in GCs (16,17). Thus, TME-related gene set-based prognostic signatures might be immunotherapic signs in GCs.
In this study, we systematically characterized the immunogenomic landscape and immune-related gene (IRG) signatures of Asian and white GC by investigating the TCGA-STAD and ACRG transcriptional profiles. We further determined the clinical role of immune genes as tools for classifying the prognoses of AGC patients and developed and validated an individualized Asian-specific IRG-based prognostic index (ASIRGPI). Furthermore, the relationship between ASIRGPI and the TME was analyzed to explore ASIRGPIrelated survival mechanisms and therapeutic targets.

Ethical Information and Study Cohorts
This study was approved by the local Research Ethics Board of The Wenzhou Medical University Second Affiliated Hospital, and informed patient consent was waived. All procedures performed in studies involving human participants were in accordance with the 1964 Helsinki declaration and its later amendments.
Gene sequencing data and the corresponding clinical data for stomach adenocarcinoma (STAD) and adjacent normal tissue samples were downloaded from The Cancer Genome Atlas database (TCGA; https://portal.gdc.cancer.gov). In addition, the raw GSE62254/ACRG (Asian Cancer Research Group) data, containing 300 AGC patients, including survival information, were downloaded from the Gene Expression Omnibus database (GEO: https://www.ncbi.nlm.nih.gov/geo/). Z-score normalization was performed for all gene sequencing data.

Prognostic Model Establishment
Differentially expressed genes (DEGs) between tumor and adjacent normal tissues were determined with the filtering condition of log2 | fold change | > 1 and a false discovery rate (FDR) < 0.05, using the limma package in R software 3.6.3. The IRG list, which has been validated to be involved in immunity, was downloaded from ImmPort, a platform that provides accurate and timely immunological data. Differentially expressed IRGs were thus extracted from the DEGs.
Patients who have an overall survival (OS) of <90 days as well as missing survival information were excluded from further study. Thus, a total of 364 patients were finally enrolled. Univariate Cox and further LASSO regression analysis was used to establish the prognostic model. The survival receiver operating characteristic (ROC) curve was performed to evaluate the performance of the ASIRGPI and to determine the cutofffor classifying AGC patients as low-or high-risk. Survival analysis associated with the prognostic model was carried out via Kaplan-Meier analysis.
Furthermore, multivariate Cox regression analyses were used to evaluate whether the prognostic model could independently predict the prognosis of AGCs. A nomogram was thus formulated using the coefficients of the multivariable Cox regression model via the rms package in R. Calibration curves were assessed graphically by plotting the observed rates against the nomogram-predicted probabilities.

TME Characterization
ESTIMATE, which uses gene expression signatures to infer stromal and immune cell fractions to determine stromal and immune scores, was performed to analyze the TME subcomponent, while xCELL, which also uses gene expression to infer the proportions of 64 tumor-infiltrating immune cell (TIIC) and stromal cell types, was further used to analyze the TME cell type. Additionally, gene set variation analysis (GSVA) was performed to estimate the enrichment scores of 110 immunoregulation-related pathways in AGC while the immunomodulatory gene was estimated via analyzing 78 immunomodulatory genes summarized by the TCGA immune response working group. A detailed flowchart is shown in Figure 1.

Immunohistochemistry Detection
Tissue microarrays from 154 local AGC patients (35 tissue pairs and 119 GC tissue) were constructed according to our previous research (18). IHC was performed in the Pathology Laboratory of  The Second Affiliated Hospital of Wenzhou Medical University, using rabbit anti-TGFB1 monoclonal antibody (Abcam, cat. no. ab27969) at a dilution of 1:250. The results were initially determined by two pathology experts and were accepted if a third expert also confirmed the result. Otherwise, the data were reviewed by all three experts and discussed until a consensus was reached. Finally, a computer-automated method was conducted (Image-Pro Plus 6.0; Media Cybernetics, Inc.) and the expression level was represented as numbers of positive cells per square millimeter (positive cells number/total area).

Statistical Analysis
All statistical analyses are carried out using R (version 3.6.3) and SPSS 22.0. Pearson correlation analysis was used to determine the correlation, and survival analysis was performed using the log-rank test. p < 0.05 was considered statistically different.

Landscape of IRGs and TME in Asian GC
To identify Asian-specific IRGs, we first selected all Asian tissues (74 tumor and 7 adjacent normal tissues) in the TCGA-STAD cohort and the differential expression patterns of 2,498 IRGs in AGC tissues were further analyzed to form an expression profile to analyze differentially expressed IRGs. We determined 685 DEGs, including 515 upregulated and 130 downregulated genes between GC and adjacent normal samples, using the limma package with cutoffs |log fold change| > 1 and FDR < 0.05. After removing genes not detected in the ACRG dataset, 183 upregulated and 67 downregulated IRGs were finally identified (Supplementary Figures S1A, B). The expression profiles of the 250 IRGs in the TCGA-STAD Asian cohort and TCGA-STAD white cohort are shown in Figure 2A. Further, gene set enrichment analysis (GSEA) revealed that TGFb signaling, coagulation, myogenesis, TNFA signaling, and epithelialmesenchymal transition (EMT) were enriched in the TCGA-STAD white cohort, while MYC targets V2, DNA repair, G2M checkpoint, MYC targets V1, and E2F targets were enriched in the TCGA-STAD Asian cohort ( Figure 2B and Supplementary Figure S2).
To further explore the survival differences between the different ethnicities, a Kaplan-Meier survival curve was performed. After excluding patients with an OS of <90 days and those without survival information, we found that white GC patients (TCGA-STAD white cohort) correlated with poor prognosis compared with AGC patients (ACRG cohort and TCGA-STAD Asian cohort; p = 0.004, Figure 2C).
Meanwhile, statistical differences of infiltrating TME cells between the Asian and white GC patients were investigated; we found that the fractions of B cells naïve, plasma cells, T cells CD8, T cells CD4 memory resting, T cells regulatory, NK cells resting, NK cells activated, monocytes, macrophages M0, and macrophages M2 were significantly increased in AGC, while T cells follicular helper, T cells gamma delta, macrophages M1, dendritic cells resting, dendritic cells activated, mast cells activated, eosinophils, and neutrophils were significantly decreased ( Figure 2D).

Generation and Validation of ASIRGPI in TCGA and ACRG Cohorts
To further explore the clinical significance of differentially expressed IRGs, univariate Cox regression analysis was performed and 35 IRGs were found to be survival-associated using the 64 AGC samples and clinical information for the TCGA-STAD cohort ( Figure 3A). LASSO regression analysis further reduced 35 potential IRGs to 7, which had non-zero coefficients in the regression model ( Figures 3B, C). An Asianspecific IRG-based prognostic index was calculated for each patient using the following formula: ASIRGPI = (0.0146 * expression of CRABP1) + (0.2543 * expression of F2R) + (0.0413 * expression of LTB) + (0.1680 * expression of PLSCR1) + (0.2501 * expression of S100B) + (0.0979 * expression of SEMG1) + (0.1692 * expression of TYROBP).
ROC analysis was carried out to assess the ASIRGPI. The high area under the curve (AUC of 0.903) confirmed the high prognostic performance of the ASIRGPI in survival surveillance. Additionally, the maximal Youden index value for an ROC curve was determined as the optimal cutoff point (0.3086; Supplementary Figure S3). AGC patients in the TCGA were thus assigned into low-and high-risk cohorts. High-risk patients (16,25.0%) had shorter OS (HR = 56.89, p < 0.001) than low-risk patients (48, 75.0%) among the 64 Asian patients with GC ( Figure 3E). The expression of the seven prognostic genes and the relationship between ASIRGPI distribution and the survival status are shown in Figures 3D, E.

Clinical Correlation Analysis and Construction of ASIRGPI-Based Nomogram
The clinical characteristics of the patients in the TCGA and ACRG cohorts are depicted in Table 1. To explore the prognostic value of ASIRGPI, multivariate Cox regression analyses were conducted for all 364 AGC patients. We found that the ASIRGPI, age, and TNM stage serve as independent predictors of AGC patient survival outcomes ( Figure 5A). Similarly, high-risk patients (88, 29.7%) had shorter OS (HR = 1.79, p < 0.001) than low-risk patients (256, 70.3%; Figure 5B).
To further elaborate the clinical significance of the ASIRGPI, the relationship between the ASIRGPI and clinical and demographic characteristics, including age and TNM stage (T stage, lymphatic invasion, and distant metastasis), according to the International Union against Cancer was analyzed. We found that the ASIRGPI was positively correlated with both age and TNM stage ( Figures 5C, D).
Finally, to provide a more intuitive clinical application tool, a nomogram was integrated with the ASIRGPI; age and TNM stage were constructed based on multivariate Cox analysis ( Figure 5E). Calibration plots showed that the nomogram could predict the probability of 3-and 5-year OS well ( Figure 5F).

Association of TME Subcomponents and Cell Types With ASIRGPI and AGC Patient Outcomes
To explore potential ASIRGPI-related survival mechanisms, ESTIMATE was performed and stromal/immune scores were inferred via ssGSEA for the entire cohort. Pearson's correlation analysis showed that both the stromal and immune scores were positively correlated with ASIRGPI [r = 0.74, p < 0.001 ( Figure 6A) and r = 0.70, p < 0.001 ( Figure 6B), respectively]. Additionally, using the median score as the cutoff values, survival analysis showed that high-stromalscored AGC patients had a worse OS than low-stromalscored patients (HR = 1.59, p = 0.003; Figure 6A), but no OS difference was observed between low-and high-immune-scored patients (HR = 1.16, p = 0.330; Figure 6B). Furthermore, TME cell-type analysis was performed and 64 TME cell types were inferred. Among these, nine cell types ( Figure 6C) were both significantly related to overall survival (log-rank test, p < 0.05)

TGFB1 Validation as the Therapeutic Target in AGC Patients
Next, to explore potential ASIRGPI-associated therapeutic targets, GSVA was performed to estimate the enrichment scores of the 110 immunoregulation related pathways. We found 13 panimmune gene sets to be significantly associated with overall survival (P < 0.05) and with ASIRGPI (|r| ≥ 0.40, P < 0.05); 11 were positively associated with worse outcomes and ASIRGPI risk scores ( Figure 6D). The remaining two genes were positively associated with worse outcomes and negatively associated with ASIRGPI risk scores, while the Module11_Prolif_score was the opposite. Interestingly, TGFb-related gene sets (TGFB_PCA_17349583 and TGFB_score_21050467), which are closely associated with the remodeling of stromal components in the TME, were highlighted during this screening.
To further analyze ASIRGPI-related molecular targets, 29 immunomodulatory genes that could be detected in the entire cohort of 364 patients were explored ; only two (EDNRB and TGFB1) were found to be positively correlated with the ASIRGPI (r ≥ 0.40, P < 0.05) and significantly associated with poor outcomes ( Figure 7A). OS analysis (Figures 7B, C) further proved that expression of the therapeutic target TGFb1 was only significantly associated with poor outcome of the AGC patients (TCGA-STAD Asian cohort and ACRG cohorts), while no such phenomenon was found in the TCGA-STAD White cohort ( Figure 7D).
To verify the effect of TGFB1 on Asian gastric cancer, IHC detection of TGFB1 was performed in the TMA. The TGFB1 expression in 35 pairs of matched tissues was analyzed (Supplementary Figure S4A). The TGFB1 expression was significantly higher in the tumor tissue (p = 0.031, Supplementary Figure S4B). Among all samples, 83.8% (129 of 154) were TGFB1-positive overall. A cutoff point of 8.68, which was optimized using a ROC, was chosen to categorize patients into TGFB1 high and low expressed subgroups. Further survival analysis found that the high-TGFB1-expressed patients had a shorter OS than the low-expressed patients (p < 0.001; Figures 7E, F).

DISCUSSION
GC patients often display heterogeneous clinical outcomes, with OS ranging from months to decades (2,8). Even among patients at the same TNM stage and receiving the same treatment, survival outcomes vary widely (15). In this regard, treatments are limited, as all GC patients receive a similar therapeutic regimen, lacking individual differences. Considering the striking differences in both the incidence rate and OS of the disease between Asian and Western countries (19), we first paid special attention to this difference and identified gastric cancer patients by ethnicity. Interestingly, we found that the activation   (24). Thus, we first attempted to analyze the prognostic value of IRGs and TIICs in Asian patients with GC. In this study, 35 differentially expressed IRGs in AGC that are significantly associated with OS were identified using public AGC cohorts. LASSO Cox regression analysis was performed, and seven IRGs were identified to be significantly associated with AGC progression. We developed a robust ASIRGPI model based on the seven IRGs and proved its efficacy in the ACRG cohort. Additionally, the unavailability of our ASIRGPI model for Caucasians further confirmed its Asian specificity. The correlations of overall survival with age, gender, TNM stage, and ASIRGPI were analyzed, and the ASIRGPI was determined as an independent predictor for outcomes, which could provide potential practical guidance for individual therapeutic regimens and improved antitumor immune responses in AGC.
Previous studies (24)(25)(26) have shown that TIICs are highly associated with tumorigenesis, invasion, and metastasis. The interactions between TIICs and tumor cells are considered to be directly associated with physical tumor cell destruction, tumor A C D B FIGURE 6 | Association of TME subcomponents and TME cell types with ASIRGPI and outcome in AGC patients. (A, B) Association of TME subcomponents with risk score and outcome in the entire 364 AGC cohort. (A) Scatter plots depicting the positive correlation between stromal score and ASIRGPI in AGC patients. Pearson's correlation coefficient is shown in the graphs (p < 0.001). Kaplan-Meier curves for overall survival of 364 AGC patients according to stromal score. Logrank test, p = 0.003. (B) Scatter plots depicting the positive correlation between immune score and ASIRGPI. Pearson's correlation coefficient is shown in the graphs (p < 0.001). Kaplan-Meier curves for overall survival of 364 AGC patients according to immune score. Log-rank test, p = 0.330. (C) TME cell type analysis, xCell, a method that uses gene expression signatures to infer the proportions of 64 immune and stromal cell types in samples, was utilized to determine the enrichment score of each cell type via ssGSEA. Among the 64 TME cell types, those significantly related to overall survival (log-rank test, p < 0.05) and risk score (Pearson's correlation test, |r| ≥ 0.40, p < 0.05) are listed. The circular data markers indicate estimated hazard ratios. The error bars represent 95% CIs. Pearson's correlation coefficients between nine TME cell types and stromal scores are also shown (p < 0.05). (D) Panimmune gene set analysis, gene set variation analysis (GSVA) was used to estimate the enrichment scores of 110 immunoregulation-related pathways in AGC samples. Among the 110 panimmune gene sets, those significantly related to overall survival (log-rank test, p < 0.05) and risk score (Pearson's correlation test, |r| ≥ 0.40, p < 0.05) are listed. The circular data markers indicate estimated hazard ratios. The error bars represent 95% CIs. Pearson's correlation coefficients between 13 panimmune gene sets and stromal scores are also shown (p < 0.05).
burden reduction, and clinical prognosis improvement. In this study, we also examined the relationship between TME subcomponents and ASIRGPI and outcomes in AGC patients. Considering that our ASIRGPI was based on IRGs, we found that the ASIRGPI was strongly and positively related to the stromal and immune scores. However, in spite of this correlation, only stromal score was relevant to the prognosis of AGC patients, which is consistent with the results of a previous study (27). This shows that the abundance of stromal components is independently related to the prognosis of GC. Hence, we infer that the effect of ASIRGPI on the worse survival of AGC patients may be related to the remodeling of stromal components; moreover, further TME cell type analysis also revealed that stromal cells accounted for the largest proportion of the nine cell types that were positively correlated with both prognosis and ASIRGPI. In summary, our results provide new insight into the mechanism by which ASIRGPI regulates the prognosis of AGC patients. To further explore potential ASIRGPI-based therapeutic targets for AGC patients with poor prognosis, a panimmune gene set was performed and the TGFb-related gene set was revealed to be significantly correlated with ASIRGPI and poor survival. Similar to previous studies, which demonstrated that high expression levels of TGFb1 were associated with poor prognosis in several cancers as well as GC (27,28), further immunomodulatory gene analysis also indicated that TGFb1 are significantly correlated with ASIRGPI and poor survival. Our experimental verification in GC tissue microarrays also confirmed this conclusion at the protein level. Additionally, since TGF-b signaling pathway activation is regarded as a symbol of extracellular matrix dysregulation and EMT (29), we also found that TGFB1 was highly expressed in extracellular matrices in the GC tissue microarrays. On this basis, we infer that the poor prognosis in both white GC patients and AGC patients with high ASIRGPI can be attributed to the remodeling of stromal components via TGF-b signaling pathway activation; however, the exact mechanism requires further study. It is worth mentioning that TGFb1 was only found to be significantly associated with poor outcome of the AGC patients. Galunisertib and M7824, molecules targeted to block the TGFb signaling pathway, have already been used in the clinical treatment of a variety of cancers (30-32); they may now become a new treatment for AGC rather than for all GC patients.
Despite our significant findings, this study had certain limitations. First, although a well-validated prognostic model was established, we only focused on AGC patients; the sample size was relatively small, and further validation of local data is necessary. Second, we gathered TCGA-STAD and ACRG/GSE62254 cohorts; sampling bias in sequencing methods on account of the differences in platforms used is inevitable. Finally, our study provides new insight into the AGC stromal microenvironment and related potential targets for individual therapy. However, this was a retrospective study; prospective studies and in-depth investigations into their functions are urgently needed to confirm these findings.
In summary, this study is the first to systematically demonstrate that AGC patients with good prognosis have lower TGFb signaling enrichment. We further analyzed the role of IRGs in monitoring the survival of AGC and developed and validated a survival-associated IRG-based ASIRGPI model. This may have important clinical implications for the survival outcomes of AGC. Finally, we comprehensively analyzed the TME characterization to explore related survival mechanisms based on ASIRGPI and found that TGFB1 may be a suitable novel target for individual AGC therapy.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the local Research Ethics Board of The Wenzhou Medical University Second Affiliated Hospital. The patients/ participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
XY and XS conceived and designed the study. CM, LM, YH, and XY were responsible for data acquisition. HH and WC have verified the underlying data. CM, LM, YH, XY, and AS were responsible for analysis of data. CM, LM, and RG were responsible for the interpretation of data. CM, XX, and XS prepared the manuscript with assistance from all authors. All authors were responsible for the revision of the manuscript, approval of the final version for publication, and accuracy and integrity of the work.