Mining The Cancer Genome Atlas database for tumor mutation burden and its clinical implications in gastric cancer

Abstract

BACKGROUND

Tumor mutational burden (TMB) is an important independent biomarker for the response to immunotherapy in multiple cancers. However, the clinical implications of TMB in gastric cancer (GC) have not been fully elucidated.

AIM

To explore the landscape of mutation profiles and determine the correlation between TMB and microRNA (miRNA) expression in GC.

METHODS

Genomic, transcriptomic, and clinical data from The Cancer Genome Atlas were used to obtain mutational profiles and investigate the statistical correlation between mutational burden and the overall survival of GC patients. The difference in immune infiltration between high- and low-TMB subgroups was evaluated by Wilcoxon rank-sum test. Furthermore, miRNAs differentially expressed between the high- and low-TMB subgroups were identified and the least absolute shrinkage and selection operator method was employed to construct a miRNA-based signature for TMB prediction. The biological functions of the predictive miRNAs were identified with DIANA-miRPath v3.0.

RESULTS

C>T single nucleotide mutations exhibited the highest mutation incidence, and the top three mutated genes were TTN, TP53, and MUC16 in GC. High TMB values (top 20%) were markedly correlated with better survival outcome, and multivariable regression analysis indicated that TMB remained prognostic independent of TNM stage, histological grade, age, and gender. Different TMB levels exhibited different immune infiltration patterns. Significant differences between the high- and low-TMB subgroups were observed in the infiltration of CD8+ T cells, M1 macrophages, regulatory T cells, and CD4+ T cells. In addition, we developed a miRNA-based signature using 23 differentially expressed miRNAs to predict TMB values of GC patients. The predictive performance of the signature was confirmed in the testing and the whole set. Receiver operating characteristic curve analysis demonstrated the optimal performance of the signature. Finally, enrichment analysis demonstrated that the set of miRNAs was significantly enriched in many key cancer and immune-related pathways.

CONCLUSION

TMB

Key Words: Tumor mutational burden, Gastric cancer, Prognosis, Immune infiltration, microRNA, Immunotherapy

Core Tip: Whether tumor mutation burden (TMB) is associated with a favorable prognosis remains controversial in various cancers. Accumulating evidence highlights that it is necessary to explore clinical impact of TMB in gastric cancer (GC). We defined the highest mutation load quintile (top 20%) in GC as the high-TMB group and found that high TMB values were associated with improved clinical outcomes, which might be attributed to the induction of antitumor immune responses in the microenvironment. We developed a microRNA-based signature to predict TMB values, which might serve as a surrogate biomarker for TMB in GC and aid physicians in clinical medical decision-making.

INTRODUCTION

Gastric cancer (GC) represents the fifth most frequent malignant disease around the globe. GC alone accounts for 8.2% of cancer-related deaths worldwide and represents a heavy economic burden and a serious public health concern, especially in China[1,2]. Despite the considerable benefits of current therapies, such as targeted biological agents and combination therapies, prognosis for advanced GC is still poor with a median survival of less than 12 mo[3,4]. Thus, it is urgently important to develop new therapeutic approaches to prolong patient life. More recently, immunotherapy with immune checkpoint inhibitors (ICIs) has emerged as one of the most promising therapeutic approaches for various solid tumors[5-7]. ICIs, specifically programmed death receptor-1/ligand 1 (PD-1/L1) antibodies, have been approved for the treatment of advanced and refractory GC[8,9]. However, only a small subset of these patients have shown a response to ICIs owing to the complexity of immunosuppressive mechanisms and genetic heterogeneity among tumors[10-13]. Therefore, a reliable biomarker is required to determine which patients can respond to ICIs and guide the selection of GC patients for immunotherapy.

Currently, tumor mutational burden (TMB), referring to the amount of nonsynonymous mutations per one million bases, is in the spotlight as a novel biomarker and a rational target for predicting response to ICIs. High TMB may be a response biomarker for immunotherapy, based on the established notion that high mutational burden could facilitate neoantigen accumulation on tumor cells, enhancing immune cell activities in the microenvironment, subsequently eliciting T-cell-dependent immune responses, and thereby inhibiting tumor development[14,15]. ICIs can restore neoantigen-mediated antitumor immune responses, thus patients with high TMB are more likely to respond to immunotherapy and exhibit improved clinical outcomes. For example, patients with high mutational load, including bladder cancer, melanoma, and lung adenocarcinoma, have appeared to benefit from ICIs[16-18]. A positive correlation between high incidence of TMB and overall survival benefit has been found in a small cohort of patients suffering from refractory GC treated with PD-1 antibody[19]. However, studies of TMB in GC patients treated with ICIs are limited in number to date[19,20]. The association of mutational load with clinical characteristics/ outcomes and immune infiltration in the microenvironment also remains lacking. Continued research is needed to delineate the somatic mutation profile of GC and explore the correlations of TMB with immune cell fractions.

There are two traditional approaches used to assess TMB in formalin-fixed, paraffin-embedded tissue of patients in most studies to date, including whole genome sequencing and whole exome sequencing (WES). WES is generally regarded as the definitive standard for mutation load assessment, but it is currently still clinically impractical because of high cost, turnaround time, and tissue heterogeneity[21,22]. Furthermore, WES requires the analysis of the matched normal tissue to remove germline mutations, and it is difficult for clinical doctors to use complex bioinform-atics algorithms to quantify TMB[23]. Several targeted gene panels focusing on cancer-related regions have recently developed as new methods to determine TMB; however, the demand for larger amounts of tumor DNA limits their use in clinical practice[24]. Thus, to find surrogate biomarkers that can predict TMB status accurately in GC is highly desirable. Lv et al[25] have identified a classifier based on microRNA (miRNA) expression patterns to predict mutational load in lung adenocarcinoma. MiRNAs are endogenous non-coding RNAs with the capacity to modulate many biological processes through gene regulation and have the potential to be biomarkers in GC immunotherapy[26,27]. Therefore, we hypothesized that miRNAs can also be surrogate biomarkers that highly correlate with the TMB status in GC.

In the current study, we determined the tumor mutational profiles of patients with GC by using the Cancer Genome Atlas (TCGA) data (https://portal.gdc.cancer.gov/repository). Specifically, we attempted to answer the following questions: Is the TMB an independent predictive biomarker for GC patients? Is the abundance of immune cell fractions in the microenvironment of high mutational load subgroup different from that of low subgroup? Do miRNAs have the potential to predict TMB values in GC?

MATERIALS AND METHODS

Data acquisition

Data on somatic mutations, RNA-seq, and miRNA expression profiles for GC were obtained from the TCGA database. For mutation data, we chose the “Masked Somatic Mutation” data which were based on VarScan software and subsequently applied the Maftools package for mutational analysis and comprehensive visual presentation[28]. The TCGA database is freely available and open to the public; therefore, there is no requirement for additional ethical approval.

Estimation of TMB values

The TMB values for each sample were determined by measuring the total amount of nonsynonymous mutations, including somatic substitutions, coding deletions and insertions, and coding errors of genes via Perl scripts and represented as the amount of mutations per mega-base (Mb) of the genomic region being sequenced. Based on a previous study, we used a cut-off of the top 20% of the TMB (9 mutations/Mb) in this study as the cut-off value. Samples with TMB ≥ 9 mutations/Mb were defined as high mutational burden (TMB-H), whereas samples with TMB < 9 mutations/Mb were defined as low mutational burden (TMB-L)[19].

Prognosis analysis

Based on the requirement for the prognosis analysis, we excluded the GC patients having a survival time of less than 1 mo, potentially implying death caused by other disease. The Kaplan-Meier curve was utilized to explore the association of somatic mutation count with overall survival, and the survival difference between the TMB-H and TMB-L group was assessed by log-rank test. We also performed Wilcoxon rank-sum tests to evaluate the relationships of TMB levels with clinicopathological parameters, including TNM status, histological grade, age, gender, and American Joint Committee on Cancer (AJCC) stage. In addition, in order to determine whether the prognostic value of TMB was independent of other clinicopathological parameters, Cox proportional hazards regression analyses (univariate and multivariable) were utilized to determine statistical significance, representing results as hazard ratios (HR) and 95% confidence intervals (95%CI).

Immune infiltration analysis

Based on RNA-seq expression data, the R package “CIBERSORT” was utilized to quantify the levels of 22 immune cells in GC patients with a threshold P value < 0.05. CIBERSORT is a deconvolution algorithm that requires an input matrix of a known reference set and accurately determines the relative levels of leukocyte subtypes from their gene expression profiles[29]. Next, we used the R package “pheatmap” to visualize the distributions of immune cell fractions in the high- and low-TMB subgroups. The difference in immune cell abundance between the high- and low-TMB groups was compared by Wilcoxon rank-sum test and visualized with the R package “vioplot”.

Construction of a miRNA-based signature

The patients diagnosed as GC were randomly sorted into a training dataset (60%) and a testing dataset (40%). The training set was used to construct a miRNA-based signature for TMB prediction. First, the miRNAs differentially expressed between the high mutation load and low mutation load groups were identified by analyzing the training set using the R package ‘limma’, and only the miRNAs with a P value < 0.01 and fold change (FC) > 1.5 were selected for subsequent analysis. We illustrated the differential expression patterns of miRNAs in the TMB-L and TMB-H subgroups by performing bidirectional hierarchical clustering and generating a heatmap plot. To construct the miRNA-based signature for predicting TMB values, we used the package “glmnet” in R software to conduct least absolute shrinkage and selection operator (LASSO) regression analysis. LASSO analysis is a powerful method that can improve prediction accuracy by constructing a penalty function, shrinking some coefficients to zero, reducing the number of variables, and finally selecting only a subset of variables into the model[30]. Then, ten-fold cross-validation was performed with type.measure = “auc” to identify the optimal tuning parameters (minimum value of lambda) in the LASSO model. The final model retained all predictors with coefficients not equal to zero. Finally, the individual index score was constructed as a linear combination of the expression level of miRNA multiplied with a regression coefficient (β) for the corresponding miRNA obtained from the LASSO regression model: The index = (β_miRNA1 * expression level of miRNA1) + (β_miRNA2 * expression level of miRNA2) + (β_miRNA3 * expression level of miRNA3) + (β_miRNAn * expression level of miRNAn).

Performance of the miRNA-based signature

The testing set and the whole set were utilized to validate the robustness and predictive performance of this signature. The ability of the signature to predict TMB value was assessed by using receiver operating characteristic (ROC) curve methodology and calculating the area under the curve (AUC) with the R package “survival ROC”. Sensitivity and specificity of this signature, as well as positive predictive value (PPV) and negative predictive value (NPV), were determined in all the sets. The relationship between the signature index of each sample and TMB value was determined by the Spearman method. With the aim to investigate the association between this signature and immune checkpoint molecules, we explored the correlation of the signature index with gene expression levels of immune checkpoints including cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), PD-L1, and PD-1.

Determination of miRNA functions

To identify the biological functions of the list of miRNAs generated to construct the predictive signature, pathway enrichment analysis including KEGG pathways and GO was conducted using DIANA-miRPath v3.0, an online software suite capable of deciphering miRNA function with experimental support (http://www.microrna.gr/miRPathv3)[31]. The results of GO and KEGG pathway analyses were considered to indicate significance at a threshold of P value < 0.01 and were visualized with the R package “ggplot2”.

Statistical analysis

R software (version 3.6.1) was employed to implement the statistical analyses in the study. P values < 0.05 were considered significant unless otherwise specified.

RESULTS

Mutational genomic landscape in GC

We included 433 GC patients for mutational analysis in this study. As depicted in Figure 1A-F missense mutation was the most common mutational category, and single nucleotide polymorphism accounted for the most frequent variant type. We classified single nucleotide variants into six classes and the results revealed that C>T mutations exhibited the highest incidence (91939) in GC. The number of variants per sample ranged from 0 to 5612, and the median number was 89. The waterfall map presented the top 30 mutated genes and their status with respect to mutational categories (Figure 1G). The top 20 mutated genes were as follows: TTN, TP53, MUC16, ARID1A, LRP1B, SYNE1, FLG, FAT4, CSMD3, PCLO, DNAH5, KMT2D, FAT3, OBSCN, HMCN1, RYR2, ZEFX4, SPTA1, CSMD1, and PIK3CA. We used the interaction plot to display the co-occurrence and exclusive correlations among the top 20 mutated genes (Supplementary Figure 1). Blackish green represents the coincident associations across mutated genes, whereas yellow represents the exclusive associations. A gene-cloud plot is used to present mutation information for genes in Supplementary Figure 2.

Figure 1 Mutational genomic landscape in The Cancer Genome Atlas gastric cancer cohort. A: Variant categories; B: Variant types; C: Single nucleotide variant types; D: Number of variants per sample; E: Summary of variant categories; F: Top 10 mutated genes. G: Waterfall map of the top 30 mutated genes and their status of variant categories. Various colors with annotations at the bottom represent the different variant categories while the barplot above the legend exhibits the value of tumor mutational burden. SNV: Single nucleotide variants; SNP: Single nucleotide polymorphism; INS: Insertion; DEL: Deletion; TCGA: The Cancer Genome Atlas.

Prognosis value of TMB in GC

First, the value of TMB was calculated for each sample in the whole set and patients were classified into the TMB-H (n = 87) and TMB-L group (n = 346) on the basis of the cut-off of the top 20% of TMB value. Then, the prognosis capacity of TMB was determined by Kaplan-Meier analysis and log-rank test, which indicated that patients in the high mutational burden subgroup showed a significantly better overall survival than those in the low subgroup (P = 0.020, Figure 2A). In addition, we included age as a continuous variable and gender, histological grade, and TNM stage as categorical variables for univariate and multivariable Cox regression analyses to further investigate the clinical value of TMB (Figure 3). Results of the Cox regression indicated that TMB value was an independent and favorable prognostic biomarker for overall survival in GC (HR = 0.982, 95%CI: 0.967−0.997, P = 0.021), as were gender (HR = 1.693, 95%CI: 1.119−2.561, P = 0.013) and lymph node metastasis (HR = 1.280, 95%CI: 1.022−1.604, P = 0.032) (Figure 3B). High TMB was considered a favorable prognostic factor, while other parameters were confirmed to be unfavorable prognostic factors for GC. We further compared the differences of TMB among different subgroups. High levels of TMB were observed in the GC patients with the following characteristics (Figure 2): Over the age of 65 years (P < 0.001, Figure 2B), female gender (P = 0.033, Figure 2C), early AJCC stage (P = 0.044, Figure 2E), AJCC T1-2 stage (P = 0.039, Figure 2F), and lack of lymph node metastasis (P = 0.018, Figure 2G). However, there were no significant differences observed in the correlations of TMB value with AJCC M stage (P = 0.104, Figure 2H) or histological grade (P = 0.051, Figure 2D).

Figure 2 Clinical significance of tumor mutation burden in gastric cancer patients. A: Survival analysis to explore the overall survival of gastric cancer patients between the high and low tumor mutation burden groups; B-H: Correlation between tumor mutation burden values and clinical characteristics in gastric cancer. TMB: Tumor mutation burden.

Figure 3 Identifying independent prognostic parameters in gastric cancer. A: Forrest plot of the univariate Cox regression analysis in gastric cancer; B: Forrest plot of the multivariate Cox regression analysis in gastric cancer. Solid squares indicate the hazard ratios of death, and close-ended horizontal lines represent the 95% confidence intervals. TMB: Tumor mutation burden.

TMB and immune cell infiltration

To further explore the potential relationships between TMB value and immune infiltration in the tumor microenvironment, we used a deconvolution algorithm to calculate the 22 immune cell fractions and present the immune infiltration landscape for each sample in Figure 4A. We further used Wilcoxon rank-sum tests to compare the difference in immune infiltration between the high mutational load and low mutational load subgroups, and the results demonstrated that patients with high mutational load exhibited significantly increased abundance of CD8+ T cells (P = 0.023), T follicular helper cells (Tfh, P < 0.001), M1 macrophages (P < 0.001), and activated CD4+ T memory cells (P < 0.001). However, resting CD4+ memory T cells (P = 0.026), regulatory T cells (Tregs, P < 0.001), and naïve B cells (P < 0.001) displayed notably decreased proportions in the TMB-H group (Figure 4B).

Figure 4 Quantitative analysis of immune infiltration in two groups based on tumor mutation burden status. A: The box plot of 22 tumor-infiltrating immune cells in each sample; B: Expression comparison of 22 tumor-infiltrating immune cells between the high tumor mutation burden (TMB) and low TMB groups. TMB: Tumor mutation burden.

Identification of a miRNA-based signature for TMB prediction from the training dataset

Patients with complete miRNA expression information were enrolled for subsequent study (n = 425) and randomly sorted into a training dataset (n = 255) and a testing dataset (n = 170). No clinicopathological characteristics, including age, gender, histological grade, or AJCC stage, were significantly different between the training and testing cohorts, as shown in Supplementary Table 1. Next, we conducted differential expression analysis in the training set, and 70 miRNAs were identified based on a cut-off point (|log₂FC| > 0.585 and P < 0.01) to be differentially expressed between the TMB-L and TMB-H groups. Among these miRNAs, 22 differentially expressed miRNAs were downregulated and 48 miRNAs were upregulated. A heatmap of the top 50 differentially expressed miRNAs is shown in Figure 5. Next, the differentially expressed miRNAs were input into LASSO analysis. Ten-fold cross-validation was used for selecting parameters in the LASSO model by minimum criteria. At the optimal values log (λ), dotted vertical lines were set via the minimum criteria, where 23 features were selected to establish the prognostic model (Figure 6A). Subsequently, the index of each patient was measured to predict the value of TMB as follows: Index = -4.801 + miR-452-5p*(-0.773) + miR-203b-3p*(0.081) + miR-582-3p*(0.194) + miR-582-5p*(0.140) + miR-27a-5p*(0.032) + miR-651-5p*(0.057) + miR-508-3p*(0.273) + miR-410-3p*(0.139) + miR-181d-5p*(-0.760) + miR-96-5p*(0.467) + miR-30a-3p*(-0.238) + miR-155-5p*(0.773) + miR-4662a-5p*(-0.130) + miR-196b-5p*(0.264) + miR-3913-5p*(-0.205) + let-7g-3p*(-0.338) + miR-210-3p*(0.352) + miR-497-5p*(-0.727) + miR-9-5p*(-0.037) + miR-625-5p*(-0.061) + miR-181b-5p*(-0.869) + miR-100-5p*(0.527) + miR-338-5p*(0.453).

Figure 5 Heatmap of top 50 differentially expressed miRNAs between the high and low tumor mutation burden groups.

Figure 6 Least absolute shrinkage and selection operator analysis for prognostic features screening and receiver operating characteristic curves for the miRNA-based signature. A: Least absolute shrinkage and selection operator regression with ten-fold cross-validation obtained 23 prognostic parameters using minimum lambda value; B: Receiver operating characteristic analysis showed that the areas under curves for the training set, the testing set, and the whole set were 0.982, 0.887, and 0.947, respectively. λ: Lambda; AUC: Area under curve.

Performance of the miRNA signature for TMB prediction

In the study, we used AUC, sensitivity, specificity, PPV, and NPV to describe the performance of the miRNA signature. The performance of the signature for predicting TMB in each set is displayed in Figure 6B and Table 1. The AUCs reached 0.982, 0.887, and 0.947 for the training set, the testing set, and the whole set, respectively, demonstrating the competitive power of the signature for predicting the TMB values of GC patients. We obtained an accuracy of 0.934 with a sensitivity of 0.802, specificity of 0.968, PPV of 0.863, and NPV of 0.951 for the whole set. The accuracy in the training set was 0.953 and 0.906 in the testing set, respectively, with a sensitivity of 0.750 and specificity of 0.937. These results revealed that the miRNA-based signature predicted TMB values with high efficiency and could be used as a predictor. Previous studies have shown PD-L1 expression and TMB to be independent of each other in most cancer subtypes[32,33]. Hence, it would be meaningful to explore the association between the miRNA signature index and immune checkpoint molecules. Clearly, the signature index exhibited a significant strong correlation with the level of TMB (r = 0.57, P < 0.001). However, the index showed a low correlation with the gene expression of PD-1 (r = 0.13, P = 0.013) and CTLA4 (r = 0.13, P = 0.011) and a moderate correlation with that of PD-L1 (r = 0.33, P < 0.001) (Figure 7).

Figure 7 Correlation of the miRNA-based signature with tumor mutation burden values and immune checkpoint molecules (programmed death-1, programmed death ligand-1, and cytotoxic T lymphocyte associated antigen 4). TMB: Tumor mutation burden; PD-1: Programmed death-1; PD-L1: Programmed death ligand-1; CTLA4: Cytotoxic T lymphocyte associated antigen 4.

Table 1 Performance of the miRNA-based signature to predict tumor mutation burden in gastric cancer.

Type	Training set	Testing set	Whole set
sensitivity	0.828	0.750	0.802
specificity	0.990	0.937	0.968
PPV	0.960	0.700	0.863
NPV	0.951	0.950	0.951
Accuracy	0.953	0.906	0.934
AUC	0.982	0.887	0.947

PPV: Positive predictive value; NPV: Negative predictive value; AUC: Area under curve.

Determination of miRNA function

Functional enrichment analysis showed that 48 KEGG pathways (Supplementary Table 1) and 104 GO terms (Supplementary Table 2) were enriched for the 23 miRNAs (P < 0.01). It was noted that several KEGG pathways are involved in the development and progression of GC, including the HIPPO, PI3K-Akt, WNT, ERBB, and transforming growth factor-beta (TGF-β) signaling pathways (Figure 8A). GO enrichment analysis demonstrated significant enrichment of immune related pathways, such as Toll-like receptor signaling pathway, Fc-epsilon receptor signaling pathway, immune system process, innate immune response, and Fc-gamma receptor signaling pathway involved in phagocytosis (Figure 8B). Taken together, the functional enrichment analysis revealed the potential roles of the miRNAs in cancer-related immune processes.

Figure 8 Functional enrichment analysis of the 23 miRNAs. A: Significant enriched cancer-related KEGG pathways; B: Significant enriched immune-related GO terms. X axis refers to the number of targeted genes. Y axis refers to GO or KEGG entry name. Bubble color refers to the enrichment P value. Bubble size refers to the number of targeting miRNAs.

DISCUSSION

Biological behaviors of cancers, including tumor initiation, angiogenesis, tumor invasion, and metastasis, are driven by genome instability, expression-level modulation, and immune cells present in the tumor microenvironment[34]. With the widespread use of high-throughput molecular technologies, TMB, one of the manifestations of genetic instability, has attracted the attention of researchers. Currently, TMB is considered an innovative biomarker of immunotherapy response for multiple cancers, including GC. Further understanding of TMB and its relationship with immune cells becomes more important in the field of personalized medicine. Therefore, the present study focused on the clinical implications of TMB in GC.

In the current study, we summarize the mutational genomic landscape in GC patients based on the TCGA dataset. The top 3 mutated genes were TTN, TP53, and MUC16. TP53 is one of the most extensively studied tumor suppressor genes, and its mutation not only inhibits the suppression of tumor development, but also produces certain cancer-promoting proteins[35]. TTN, the longest known gene, originally known for its effects on the development and regulation of cardiac and skeletal muscles, has been confirmed to be highly correlated with TMB levels and the responsiveness to ICIs in solid tumors[36]. Despite the length of the gene contributing to a greater number of somatic mutations, the molecular mechanisms of TTN in the production of tumor mutations require further study. MUC16 is one of the most frequently mutated genes in various cancers and contributes to tumor proliferation and metastasis by regulating immune response to cancer[37]. In a recent study, it has been suggested that MUC16 was associated with high TMB and favorable prognosis in GC[38].

Our study analyzed the correlation between TMB values and overall survival rates in patients with GC. The Kaplan-Meier analysis in the GC cohort from the TCGA dataset was performed by defining the highest TMB quintile (top 20%) as the high TMB group. Based on this approach, we found that GC patients with high mutational burden had a favorable prognosis. In addition, higher mutational load in GC was also significantly associated with older ages, female gender, earlier tumor stage, and lack of lymph node metastasis. However, the prognostic impact of mutational burden in different types of cancer remains controversial. Yuan et al[39] demonstrated that patients with higher mutational load had a worse prognosis in esophageal cancer[39], consistent with previous studies in clear cell renal cell carcinoma and prostate cancer[40,41], whereas bladder cancer with higher mutational load was associated with better outcomes[42]. In these previous studies, they used the median value as the cut-off to distinguish TMB-H and TMB-L groups.

Currently, a universal cut-off value for defining high TMB is lacking. Previous studies have revealed that there are several factors that influence the threshold for high TMB, such as cancer types, sample types, pre-analytic variables, and detection methods[43]. Some studies have been exploratory, distributing TMB into three groups. For example, a study on a total of 908 resected lung cancer species measured TMB levels by using targeted gene panels and divided TMB into terciles: High (> 8 mutations/Mb), low (4 mutations/Mb), and intermediate subgroup (> 4 and ≤ 8 mutations/Mb)[44]. In a large cohort focusing on diverse cancers, TMB levels was separated into three subgroups: Low-(50% of patients), intermediate-(40% of patients), and high-TMB subgroup (10% of patients)[23]. Most published clinical studies were likely to distribute TMB into two groups. In a study with small sample size, researchers defined the higher mutation load quintile (top 20%, 14.31 mutations/Mb) in advanced GC treated with ICIs as the high TMB group and found that TMB was correlated with clinical outcomes[20]. Wang et al[19] have studied the correlation of TMB with survival in chemo-refractory GC under treatment with toripalimab and demonstrated that high TMB (top 20%, 12 mutations/Mb) was significantly correlated with improved survival[19]. A large-scale study also used an upper 20^th percentile cut-off to define high TMB and analyzed the association between tumor mutational load and the clinical responses to ICIs across multiple cancer types[45]. Therefore, we took the same approach of separating TMB-H and TMB-L groups (top 20%, 9 mutations/Mb) in order to dichotomize the data, but this was not a universal number of high TMB in GC and might not have any clinical significance. When we set the cut-off at the median of mutational load (2.671 mutations/Mb) in the current study, the difference of overall survival rates between the high mutational load and low mutational load subgroups was still statistically significant (P = 0.016, Supplementary Figure 3). Thus, it is critical and urgent to develop a satisfactory and reproducible definition of the predictive TMB cut-off value before implementing TMB as a biomarker of immunotherapy response in individualized treatment. To further clarify whether TMB is a favorable dependent factor in GC, Cox proportional hazards regression analysis was implemented and demonstrated that TMB was significantly correlated with overall survival as a continuous variable, indicating that it was able to predict the survival of GC patients without consideration of other conventional clinicopathological variables.

Immune cells, a large proportion of infiltrating cells in the tumor microenvironment, interact with tumor cells by releasing inflammatory cytokines and chemokines, which drive biological behaviors of cancers and influence their therapeutic responses to immunotherapy. In the present study, we explored the potential relationships between TMB values and immune infiltration in GC patients. The abundance of CD8+ T cells, CD4+ T cells, and Tfh cells were markedly higher in the TMB-H group than in the TMB-L group. Accumulating evidence has shown that high mutational load tends to cause neoantigens to accumulate in cancer and results in the activation of CD8+ cytotoxic T cells and subsequent initiation of tumor cell lysis, consistent with our results[14,15]. CD4+ T cells were considered a component of anticancer immunomech-anisms and an independent indicator of favorable prognosis in multiple cancers[46]. It is to be noted that not all mutations would generate high neoantigens on the surfaces of tumor cells, and further research is warranted to determine which gene mutations are responsible for the induction of immune response. Different from the conventional view on immune cells, TAMs play a dual role in tumor development depending on their polarization status. The M1 macrophages are involved in antitumor response whereas M2 macrophages have pro-tumorigenic properties[47]. In our study, M1 macrophages, but not M2, were significantly elevated in the tumor microenvironment with high TMB levels. Furthermore, we found that patients with high mutational load showed significantly decreased infiltrating naïve B cells and Tregs. Tregs have been considered to have a central role in inhibiting effective antitumor immunity and correlate with an unfavorable prognosis in many cancers[48]. The mechanisms underlying the decreased level of naïve B cells in the TMB-H group remain unclear. We hypothesize that B-cell differentiation factors secreted by tumor cells may be responsible for decreased infiltration of naïve B cells. In summary, these results indicate that tumor microenvironment exhibited a significant antitumor response in patients with high TMB.

Growing evidence has revealed that aberrant expression of miRNAs can be found in various cancer types and plays a crucial role in the carcinogenesis, migration, and invasion of tumor cells by regulating adaptive and innate immune responses in the tumor microenvironment[26,27]. It has been suggested that miRNA-targeted immunotherapeutics have great potential in clinical practice. However, the relationship between miRNA and mutational load in GC has not previously been explored. In the present study, we screened the differentially expressed miRNAs between the TMB-H and TMB-L groups in the GC cohort and found that different mutational load values were correlated with different miRNA profiles. Next, we conducted LASSO analysis to select parameters from differentially expressed miRNAs and established a 23-miRNA classifier to predict TMB values based on the training dataset, which was further validated in the whole dataset and testing dataset. The accuracies of this predictive model for the training set, the testing set, and the whole set were 0.951, 0.906, and 0.934, respectively, and specificities were 0.990, 0.937, and 0.968, respectively, which implied that this signature was very effective with a high accuracy and specificity in predicting TMB values. The efficiency of this signature was also evaluated by ROC analysis and the results revealed that this signature was credible in predictive performance throughout the testing set and whole set. Moreover, the strong correlation between the signature index and TMB values in GC patients further confirmed the robustness of this signature. However, PPV in the testing set was lower than that in the training and whole set, indicating that the ability of the signature to recognize high mutational load needs to be improved. We also explored the association of the signature with immune checkpoint molecules and demonstrated that the signature showed a low correlation with PD-1, PD-L1, and CTLA4. These results were consistent with previous studies, which demonstrated that TMB was independent of the expression of immune checkpoint molecules[32,33].

With respect to the biological functions of the miRNAs in the predictive signature, functional annotation was conducted. KEGG pathway analysis indicated that the functions of the 23 miRNAs were potentially associated with HIPPO, PI3K-Akt, WNT, ERBB, and TGF-β signaling pathways, and those were supposed to play a critical role in the tumorigenesis of GC[49]. Moreover, the miRNA sets were found to be involved in immune-related pathways, including immune system process and innate immune response. Increasing evidence has demonstrated that high mutational burden might lead to the activation of antitumor immune responses[14,15]. The enrichment analysis indicated that the 23 miRNAs contributed to vital cancer and immune pathways, which might provide strong biological evidence for the feasibility of the miRNA-based signature in predicting TMB values.

This study has several limitations. First, all of our samples and clinical data were based on the TCGA dataset, and most of patients were Westerners. Cohorts with larger sample sizes from other regions are needed to confirm our results, and external validation of the miRNA-based signature is necessary in the future. Second, the cut-off definition used in the present study to distinguish the TMB-H and TMB-L groups was not uniform, and multi-center randomized controlled studies focusing on immunotherapy in GC to identify TMB cut-off values are proposed. Third, all the results in the present study were based on a bioinformatics analysis and description. Mechanistic investigation should be performed to clarify the underlying mechanism of high mutation load in the activation of antitumor immune responses in GC patients. In addition, the functions of 23 miRNAs in immune responses were not characterized using in vitro or in vivo experimentation.

CONCLUSION

Taken together, mutational burden is considered an independent and favorable prognostic biomarker in patients suffering from GC. High TMB is notably correlated with a good survival and might lead to the activation of antitumor immune cells in the tumor microenvironment. Moreover, different mutational load is associated with different miRNA expression patterns, and a miRNA-based signature was established to predict TMB values in GC, which might aid physicians in clinical medical decision-making.

ARTICLE HIGHLIGHTS

Research background

Tumor mutational burden (TMB) is in the spotlight as a novel biomarker and a rational target for predicting response to immunotherapy in multiple cancers. Gastric cancer (GC) is one of the most common gastrointestinal malignant tumors worldwide. Accumulating evidence highlights that it is necessary to further explore clinical impact of TMB in GC.

Research motivation

The association of TMB with clinical outcomes and immune infiltration in the tumor microenvironment in GC patients has not yet been elucidated. MicroRNAs (miRNAs) have a crucial role in the carcinogenesis, migration, and invasion of tumor cells by regulating adaptive and innate immune responses, but the relationship between miRNA expression patterns and mutational load is not clear in GC.

Research objectives

This study aimed to explore the clinical impact of TMB and establish a miRNA-based signature for TMB prediction in GC patients.

Research methods

The Kaplan-Meier analysis in the GC cohort from The Cancer Genome Atlas dataset was performed by defining the highest TMB quintile (top 20%) as the high-TMB group. The difference in immune infiltration between the high- and low-TMB subgroups was evaluated by Wilcoxon rank-sum test. The least absolute shrinkage and selection operator analysis was conducted to select parameters from differentially expressed miRNAs between the high- and low-TMB subgroups and construct a miRNA-based signature classifier for TMB prediction.

Research results

Higher mutational load in GC was significantly associated with better prognosis, older ages, female gender, earlier tumor stage, and lack of lymph node metastasis. Different mutational load levels exhibited different immune infiltration patterns and different miRNA expression patterns. In addition, we developed a miRNA-based signature using 23 differentially expressed miRNAs to predict TMB values of GC patients.

Research conclusions

High TMB is notably correlated with good survival and might lead to the activation of antitumor immune cells in the tumor microenvironment in GC. The miRNA-based signature might be developed as a surrogate biomarker for TMB in GC.

Research perspectives

The miRNA-based signature for TMB prediction might help develop treatment strategies for GC patients and have an impact on the clinical practice in the course of GC.

ACKNOWLEDGEMENTS

The authors thank The Cancer Genome Atlas (TCGA) project for providing invaluable datasets for statistical analyses.