Genome Instability-related MiRNAs Predict Survival, Immune Landscape, and Immunotherapy Responses In Gastric Cancer

doi:10.21203/rs.3.rs-553563/v1

Download PDF

Research

Genome Instability-related MiRNAs Predict Survival, Immune Landscape, and Immunotherapy Responses In Gastric Cancer

https://doi.org/10.21203/rs.3.rs-553563/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Backgroud Increasing evidence suggests that microRNAs (miRNAs) are involved in genome instability (GI) and drive the occurrence of tumors. However, the role of GI-related miRNAs in gastric cancer (GC) remains largely unknown. Herein, we developed a novel GI-related miRNA signature (GIMiSig) and further investigated its role in prognosis, the immune landscape, and immunotherapy responses in GC patients.

Methods An analysis of somatic mutation data on 434 gastric cancer cases from The Cancer Genome Atlas (TCGA) database was performed, thereby generating genome stability (GS) and GI groups. By detecting differentially expressed miRNAs between the GS and GI groups that were associated with overall survival, 8 miRNAs were identified and used to construct the GIMiSig.

Results The GIMiSig showed high accuracy in detecting GC patients. Using GIMiSig to stratify the patients into high- and low-risk subgroups to predict survival outperformed the use of regular clinical features such as age, gender, or disease stage. Patients with low risk had a more favorable survival time than those with high risk. More importantly, the high-risk patients were associated with decreased UBQLN4 expression, higher accumulation of immune cells, lower Titin (TTN) mutation frequency, worse immunotherapy efficacy, and cancer-associated pathways. Conversely, the low-risk patients were characterized by UBQLN4 overexpression, lower fraction of immune cells, higher TTN mutation frequency, better response to immunotherapy, and GI-related pathways.

Conclusion In summary, we constructed a novel GIMiSig that could stratify GC patients into distinct risk groups that have different survival outcomes and immunotherapy efficacy. The results may provide new clues for improving GC outcomes.

Cancer Biology

gastric cancer

microRNAs

genome instability

immune infiltrates

immunotherapy

Gastric cancer (GC) is the second most common cause of cancer-related mortality worldwide, causing severe economic and social burdens (1). Although multiple strategies such as surgery, chemotherapy, and radiotherapy have been used, the outcomes of GC patients remain poor, partly due to lack of early detection and the high frequency of recurrence (2). Recently, immune checkpoint inhibitors (ICIs) targeting immune checkpoints, especially programmed cell death-1 (PD-1) and its ligand (PD-L1), have shown encouraging results in several cancer types (3). However, immunotherapy efficacy in GC remains poor because of the heterogeneity among GC subtypes and the complex interactions within the immune system (4). Therefore, there is a need to identify novel biomarkers for risk stratification, immune landscape assessment, and prediction of therapeutic effects.

Genomic instability (GI) is regarded as a hallmark of malignant tumors because it results in genomic alteration of cells, thereby facilitating cell proliferation without restraint (5). Unsurprisingly, GI is also a crucial step in driving tumorigenesis in GC patients (6). The mechanisms of GI in cancer remain mysterious, but it is widely accepted that miRNAs are involved in the regulation of GI, and they can be used to predict cancer outcomes. For instance, Tianzhen et al built a risk model involving 10 GI-related miRNAs that predicted prognosis and chemotherapy responses in ovarian cancer (7). In malignant glioma, miR-100 led to GI by impairing cell responses to DNA damage (8). In hepatocellular carcinoma, miR-24, miR-103, and miR-107 are considered to be carcinogenic factors because they enhance GI in multiple ways, such as influencing the cell cycle and reducing DNA damage repair (9, 10). Additionally, in breast cancer, miR-155 overexpression provoked GI by disturbing telomere homeostasis, contributing to an unfavorable prognosis (11). Although emerging research has shown that miRNAs participate in GI in various tumor types, the roles of GI-related miRNAs in GC remain unclear.

Additionally, the tumor microenvironment (TME) has become a research hotspot in cancer. The TME is a complex system that consists of multiple cells, proteases, and cytokines, and it impacts tumorigenesis and immunotherapy responses. More and more evidence has confirmed the crucial roles of miRNAs in the TME. For example, Chou et al found that miR-29b modulated the TME to prevent the metastasis of breast cancer (12). Furthermore, miR-9 modulates the TME by stimulating angiogenesis (13). However, research on the function of GI-related miRNAs in the TME is very scarce. Thus, we aimed to investigate the interaction between GI-related miRNAs and the TME in GC.

In this study, we combined somatic mutation data and miRNA sequencing data to develop a risk signature associated with GI. Our investigation of how the signature influences diagnosis, prognosis, immune landscape, and immunotherapy responsiveness may offer insights into the relationships between miRNAs and GC. More importantly, our study indicates that these miRNAs may be exploited as predictive hallmarks and even drug targets in GC.

Data collection

The Cancer Genome Atlas (TCGA)-GC profiles in this study were retrieved from TCGA portal (https://portal.gdc.cancer.gov/), including miRNA sequencing data (446 GC samples), mRNA sequencing data (375 GC samples), clinical data (397 GC samples), and somatic mutation data (434 GC samples). The whole TCGA cohort was randomly split into two subsets. The subset used to build the risk model was designated the training set and the other subset, which was used to validate the performance of the risk model, was designated the testing set. The characteristics of the three sets are presented in Supplementary Table 1. The GSE112264 (normal samples = 41, and GC samples = 50) and GSE54129 datasets (normal samples = 21, and GC samples = 111) used in the validation phase were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Identification of GI-related miRNAs

Integrating miRNA sequencing data and somatic mutation data from TCGA, GI-related miRNAs were identified using the following steps: (i) the cumulative number of somatic mutations in each patient was defined as the somatic mutation value (SMV); (ii) patients were ranked by descending SMV and clustered into two groups: the first group, with the top 25% SMVs, was characterized by genome instability (GI) and the second group, with the bottom 25% SMVs, was characterized by with genome stability (GS); (iii) miRNA expression was compared between the GI and GS groups; and (v) miRNAs that met the selection criteria (log[fold change] = 1 and adjusted p < 0.05) were identified as the GI-related miRNAs, which were presented in a heatmap. The R packages "limma", "pheatmap", and "sparcl" were utilized to process and visualize the data. The top 40 differential miRNAs between GS and GI groups were selected according to the absolute value of fold change, which were shown in a heatmap. Based on FunRich (14), gene ontology (GO) functions were applied to identify the possible biological functions of differentially expressed miRNAs. GO analysis were made up of three parts: cellular component (CC), molecular function (MF), and biological process (BP).

Construction of GI-related miRNA signature (GIMiSig)

Univariate and multivariate Cox regression analyses were performed to identify GI-related miRNAs that influenced survival. Eight miRNAs were used to establish the GIMiSig, which was used to calculate the risk score. The method to generate the risk score was as follows: risk score = (ExpmiRNA₁×Coef₁ + ExpmiRNA₂×Coef₂ + ExpmiRNA₃×Coef₃ + …ExpmiRNA_n×Coef_n.). Patients with higher risk score (≥ mean value) and lower risk score (< mean value) were allocated to the high- and low-risk groups, respectively. Subsequently, a Kaplan–Meier analysis was used to determine the survival distributions in the high- and low-risk GC patients, and a log-rank test was used to compare the two groups, with p < 0.05 as the threshold. The prediction accuracy of the GIMiSig regarding survival time was assessed based on the area under the receiver operating characteristic (ROC) curve (AUC). Furthermore, we performed univariate and multivariate Cox regression analyses to identify whether the GIMiSig was an independent predictor. The hazard ratios (HRs) and p values were presented in forest plots. Additionally, ROC curve analysis was applied to determine the predictive performance of the GIMiSig and clinical variables regarding survival time.

Minimally invasive diagnostic value of the GIMiSig

The expression levels of the eight miRNAs were compared between the normal and tumor groups from TCGA. Next, ROC curve analysis was performed based on the miRNA expression data of TCGA and GSE112264, which examined the performance of the GIMiSig for the diagnosis of GC.

Relationship between the GIMiSig and Titin (TTN) mutation status

Using somatic mutation data in Mutation Annotation Format, we employed the R package "maftools" to visualize the gene mutation landscape in the whole TCGA set and in the high- and low-risk groups. The ratio of wildtype to mutated TTN was compared in the high- and low-risk groups. Furthermore, GIMiSig and TTN mutation status were used to classify patients into four subgroups and the survival distributions of the four groups were analyzed.

Immune landscape analysis and prediction of immunotherapy response

First, the relationship between the GIMiSig and immune cells was determined by a Spearman’s correlation analysis, which is presented in Supplementary Table 2. We merged the results of multiple methods, such as TIMER, CIBERSORT, and XCELL, into a bubble diagram. The R packages "ggtext", "ggplot2", and "scales" were used for visualization.

Next, ESTIMATE (15) and single-sample gene set enrichment analysis (ssGSEA) (16) methods were used to evaluate the immune landscape within each risk group. ESTIMATE is an effective tool for evaluating the composition of the TME, including stromal and immune cells. Four parameters are generated by ESTIMATE: stromal score (proportion of stromal cells), immune score (proportion of immune cells), ESTIMATE Score (sum of stromal and immune scores), and tumor purity (proportion of tumor cells in the tumor tissue). The analyses were completed using the R packages "estimate", "ggplot2", "reshape2", and "violin". However, ESTIMATE cannot identify specific cell subtypes. Therefore, ssGSEA was utilized to assess the differential abundance of intra-tumor immune cells using the R package “GSVA”. The following subtypes were included: dendritic cells (DCs), activated DCs (aDCs), immature DCs (iDCs), plasmacytoid DCs (pDCs), B cells, CD8 + T cells, macrophages, mast cells, neutrophils, natural killer (NK) cells, T helper cells, follicular helper T (Tfh) cells, type-1 T helper (Th1) cells, type-2 T helper (Th2) cells, tumor infiltrating lymphocytes (TILs), and regulatory T cells (Tregs). We used the Wilcoxon signed-rank test to compare the immune landscape between the high- and low-risk groups.

The immune checkpoint-related gene expression levels (PDL1, CTLA4, PD1, LAG3, and TIM3) were extracted from TCGA. The microsatellite instability (MSI) of each patient was determined based on somatic mutations. The aneuploidy score was retrieved from a previous study (17). The tumor mutation burden (TMB) value was generated for each patient (by using the Perl language to analyze the TCGA somatic mutation data), which was log2-transformed. The Tumor Immune Dysfunction and Exclusion (TIDE) score was calculated online (http://tide.dfci.harvard.edu//) (18). Thereafter, the immune checkpoint-related gene levels, MSI score, aneuploidy score, TMB, and TIDE score were compared between the high- and low-risk groups.

Biological significance of the GIMiSig

We carried out GSEA to elucidate the GIMiSig-related functions in GC, using the Kyoto Encyclopedia of Genes and Genomes (KEGG) as the reference. The R packages "plyr", "ggplot2", "grid", and "gridExtra" were used to visualize the pathways enriched in the high- and low-risk groups.

Target genes of the GIMiSig

The possible targets of the 8 miRNAs in the GIMiSig were predicted using three databases: (i) TargetScan (http://www.targetscan.org/), (ii) miRanda (http://mirdb.org/miRDB/), and (iii) miRTarBase(http://miRTarBase.mbc.nctu.edu.tw/). The overlapping results of the three platforms were regarded as putative target genes. The R package "VennDiagram" was used to generate Venn diagrams, which show the number of target genes of each miRNA according to the three databases. GO and KEGG pathway of target genes were performed using clusterProfiler tool (19). Based on a confidence score ≥ 0.9, target genes were selected to build a PPI network using the STRING database (http://string-db.org/). Subsequently, the network was inputted into Cytoscape and the hub target genes were selected using CytoHubba. Finally, expression levels of the top 5 hub genes were compared between the normal and tumor groups from TCGA and GSE54129.

Statistical analysis

All statistical analyses were conducted using R software (version 3.6.1). Data from different groups were compared using the Mann–Whitney–Wilcoxon test. The log-rank test was performed to assess differences in survival. P < 0.05 was the significance threshold.

Different genomic alterations in the GS and GI groups

Based on the genomic profiles in the whole TCGA set, we tried to differentiate the molecular subtypes of GC. SMV (sum of somatic mutations), in descending order, was used to rank the patients. Patients with the top 25% and the bottom 25% SMVs were identified as GI patients (n = 108) and GS patients (n = 110), respectively. To obtain insights into special characteristics in patients with GI, we analyzed miRNAs, SMV, and UBQLN4 expression profiles in the two subgroups. As shown in Fig. 1A, the clustering heatmap depicts the differentially expressed miRNAs between the GI and GS patients. SMV was significantly higher in GI patients than GS patients (p = 5.7e-16) (Fig. 1B). UBQLN4 is a novel hallmark, contributing to the occurrence of GI (20). Thus, we compared UBQLN4 expression between the two groups. The boxplot shows that UBQLN4 was overexpressed in the GI group compared to the GS group (p = 0.0089) (Fig. 1C). These results suggested that it was reasonable to distinguish the GS and GI groups based on somatic mutations. Subsequently, expression level of the top 40 differential miRNAs between the GS and GI groups were shown in the heatmap (Fig. 1D). Furthermore, GO analysis was performed to explore biological function of differentially expressed miRNAs. The results showed that these miRNAs mainly participated in nucleus (CC), transcription factor activity (MF), signal transduction ( BP), regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism ( BP), and cell communication ( BP) (Fig. 1E).

Establishment of a GI-related miRNA signature (GIMiSig)

To investigate the clinical significance of GI-related miRNAs in GC, a series of bioinformatics analyses were conducted to determine whether these miRNAs correlate with the prognosis of patients. A total of 389 GC patients with clinical information were divided into the training and test set, and the former was used to build the risk model. Using univariate and multivariate Cox regression analyses(Table 1), 8 GI-related miRNAs related to overall survival (OS) were selected to build the risk model in the training set. The expression of the 8 miRNAs in each patient and the corresponding coefficients were used to calculate a risk score for each patient as follows: Risk score = (-0.446)×miR-125b-5p + 0.170×miR-99a-3p + 0.0815×miR-548v + 0.338×miR-100-5p + (-0.074)×miR-196b-3p + (-0.061)×miR-1275 + 0.043×miR-380-3p + 0.109×miR-363-3p. After calculating the risk score, patients with higher and lower risk scores compared to the mean risk score within the three sets were assigned to the high- and low-risk groups, respectively.

Table 1

Multivariate Cox regression analyses in the training set.
id	coef	HR	HR.95L	HR.95H	p-value
miR-125b-5p	-0.446	0.640	0.409	1.000	0.050
miR-99a-3p	0.110	1.116	1.000	1.246	0.049
miR-548v	0.082	1.085	1.008	1.168	0.030
miR-100-5p	0.338	1.402	0.934	2.105	0.103
miR-196b-3p	-0.074	0.929	0.888	0.971	0.001
miR-1275	-0.061	0.940	0.900	0.982	0.006
miR-380-3p	0.043	1.044	0.995	1.096	0.080
miR-363-3p	0.109	1.116	0.969	1.285	0.129

When the GIMiSig was applied to the training set, the survival data suggested that high-risk patients exhibited shorter overall survival than low-risk patients (p < 0.01) (Fig. 2A). AUC at 1, 3, and 5 years was 0.814, 0.746, and 0.796, respectively (Fig. 2B). Expression levels of the 8 miRNAs in the high- and low-risk patients were showed in the heatmap (Fig. 2C). Moreover, low-risk patients had higher SMV (p = 2.4e-06) (Fig. 2D). No significant difference was observed in UBQLN4 expression between high- and low- risk groups in the training set (p = 0.066) (Fig. 2E).

To confirm the stability of the GIMiSig, we examined its performance in the testing and whole TCGA sets. Similarly, high-risk patients had poorer prognosis in the testing (p = 0.007) and whole TCGA sets (p < 0.001) (Fig. 3A,F). The ROC curves for the GIMiSig also indicated good accuracy in the testing and whole TCGA sets (Fig. 3B,H). Distribution of miRNA expression levels changed with increasing risk score in the test and whole sets (Fig. 3C,I). As shown in the boxplot, higher SMV (p < 0.01) and UBQLN4 overexpression (p < 0.05) were observed in low-risk patients in the testing and whole sets (Fig. 3D,E,J,K ). These results indicated that the GIMiSig could stratify GC patients into two groups differing in prognosis, SMV, and UBQLN4 expression.

Evaluation of independent prognostic factors

The risk score generated by the GIMiSig was identified as a prognostic factor by univariate analysis (p < 0.001, HR = 1.351, 95% CI 1.205–1.514) (Fig. 4A), and as a factor that was independent of the other clinical features by multivariate analysis (p < 0.001, HR = 1.436, 95% CI 1.259–1.636) (Fig. 4B). More importantly, the risk score yielded the highest AUC value (0.730) compared to other parameters, such as age (0.529) and grade (0.549), which indicated that the GIMiSig is an effective tool for risk assessment of GC patients (Fig. 4C).

Survival analysis of GC subgroups

We then examined the survival time of patients in the high- and low-risk subgroups stratified by age, gender, stage, and grade. Low-risk patients had better survival time than high-risk patients in the young (age < 65 years, p < 0.001) and old (age ≥ 65 years, p = 0.004) subgroups (Fig. 5A), female (p < 0.001) and male (p = 0.011) subgroups (Fig. 5B), stage I–II (p = 0.003) and stage III–IV (p < 0.001) (Fig. 5C) subgroups, and grade 1–2 (p < 0.001) and grade 3 (p < 0.001) subgroups (Fig. 5D).

Minimally invasive diagnostic value of the GIMiSig

In addition, we assessed the minimally invasive diagnostic value of the GIMiSig. First, expression levels of the eight miRNAs were compared between the normal and tumor groups from TCGA. We found that miR-125b-5p (p < 0.001), miR-99a-3p (p < 0.001), miR-100-5p (p < 0.001), miR-548v (p < 0.001), and miR-363-3p (p < 0.001) were downregulated in tumors compared with that in normal control (Fig. 6A). On the contrary, miR-380-3p (p < 0.05), miR-196b-3p (p < 0.001), and miR-1275 (p < 0.05) was highly expressed in tumors than that in normal control (Fig. 6A). The differential expression of these miRNAs indicated that the GIMiSig may have the potential for early diagnosis of GC. Next, we evaluated whether the GIMiSig could differentiate GC from normal control in TCGA and GSE112264. The results revealed that the AUC of TCGA and GSE112264 cohort were 0.700 (Fig. 6B) and 0.877 (Fig. 6C), indicating the GIMiSig was able to detect GC patients accurately.

Relationship between the GIMiSig and TTN mutation status

As shown in the waterfall plot (Fig. 7A), TTN mutation accounted for 48% of all gene variants in the whole TCGA set. Additionally, the proportion of mutated TTN was higher than the proportions of other gene mutations in both the high-risk (38%) and low-risk (58%) groups (Supplementary Fig. 1). Thus, we next examined the distribution of wildtype and mutated TTN in the high- and low-risk groups. The results suggested that low-risk patients had more frequent TTN mutation than high-risk patients in the three sets (Fig. 7B). In the training set, TTN mutation was detected in 30% of high-risk patients and 64% of low-risk patients (p < 0.001). Although there were no significant association of TTN mutation status with risk in the test set, there was a trend toward more frequent TTN mutation in the low-risk patients (p = 0.051). In the whole TCGA set, TTN mutation was more common in the low-risk group (65%) than the high-risk group (42%) (p < 0.001). We further assessed whether the GIMiSig could predict survival better than TTN mutation status. Using TTN mutation status, patients were classified into TTN mutation and TTN wildtype group. There were no significant differences in OS between TTN mutation and wildtype TTN among the high-risk patients or between TTN mutation and wildtype TTN among the low-risk patients (p > 0.05) (Fig. 7C). When the GIMiSig was used to divide patients into high- and low-risk group, we observed significant differences in OS between high- and low-risk group among patients with TTN mutation or with wildtype TTN (p < 0.001, Fig. 7C).

Relationship between the GIMiSig and immune landscape

Accumulating evidence has shown that the TME greatly impacts tumor occurrence and responses to immune-targeting strategies. Multiple methods were utilized to evaluate the association of the GIMiSig with the TME. First, correlation analysis was conducted between the GIMiSig and various intra-tumor cell subpopulations. The results revealed that high-risk patients exhibited positive associations with CD8 + T cells, fibroblasts, B cells, and macrophages while they had negative associations with activated CD4 + T cells, T helper cells, and Treg cells (Fig. 8A). Next, the ESTIMATE and ssGSEA algorithms were used to evaluate immune characteristics in the high- and low-risk groups. ESTIMATE was applied to calculate the stromal, immune, and ESTIMATE scores for each GC patient. Interestingly, compared to the low-risk patients, the high-risk patients had higher stromal, immune, and ESTIMATE scores but lower tumor purity (p < 0.001) (Fig. 8B). ssGSEA was performed to investigate the relationships between GIMiSig and 14 intra-tumoral immune cells. The high-risk patients had higher proportions of B cells (p < 0.001), DCs (p < 0.01), iDCs (p < 0.05), mast cells (p < 0.01), neutrophils (p < 0.001), pDCs (p < 0.01), T helper cells (p < 0.01), follicular helper T (Tfh) cells (p < 0.01), TILs (p < 0.001), and Treg cells (p < 0.05) (Fig. 8C). In aggregate, these data suggested that the GIMiSig may influence the immune microenvironment.

Prediction of immunotherapy responsiveness

Immune checkpoint-related gene expression levels, MSI score, aneuploidy score, TMB, and TIDE score were used to assess the potential clinical efficacy of immunotherapy in the high- and low-risk groups from TCGA. We found that PD-L2 (p < 0.01) and TIM3 (p < 0.05) expression were significantly higher in the high-risk group than the low-risk group (Fig. 9A). The GIMiSig had no significant influence on other immune checkpoint-related genes, although a trend towards elevated CTLA-4 (p > 0.05) expression was observed in high-risk patients. Additionally, the MSI score was lower in low-risk patients than high-risk patients (p < 0.001) (Fig. 9B). As for aneuploidy score, it was higher in low-risk patients than high-risk patients (p < 0.01) (Fig. 9C). Patients with high TMB are more likely to trigger immune responses by presenting more neoantigens, so they benefit more from ICI therapy. Tumor cells in patients with high TIDE scores are more likely to escape from the surveillance of the immune system, resulting in insensitivity to ICI therapy. In this study, high-risk patients had lower TMB ( p = 1.62e-08) (Fig. 9D) and higher TIDE score (p = 7.4e-05) (Fig. 9E) compared to low-risk patients, suggesting that high-risk patients may have unfavorable responses to immunotherapy.

GSEA

To elucidate the link between the GIMiSig and its functional properties in GC, GSEA was used to detect the GIMiSig-related pathways. The results showed that several cancer-associated and immune response pathways were discovered in the high-risk patients, such as cell adhesion molecules, cytokine–cytokine receptor interaction, and the MAPK signaling pathway (Fig. 10A). Additionally, several GI-related pathways were enriched in the low-risk patients, including cell cycle, nucleotide excision repair, and mismatch repair (Fig. 10B). This might explain why higher SMV and UBQLN4 overexpression were detected in low-risk patients. Taken together, these results indicated potential roles of the 8 aberrantly expressed miRNAs in the occurrence and progression of GC: (i) contributing to tumor metastasis via cell adhesion molecules; (ii) altering the immune microenvironment based on cytokine release; and (iii) disrupting gene repair.

Target genes of the GIMiSig

To ascertain the downstream targets of the 8 miRNAs, we predicted the target genes using three public databases (TargetScan, miRanda, and miRTarBase) and presented the overlapping results (i.e., target genes shared across all three databases) in Venn diagrams. There were 31, 12, 24, 135, 106, 2, and 59 target genes of miR-548v, miR-100-5p, miR-380-3p, miR-363-3p, miR-125b-5p, miR-196b-3p, and miR-1275, respectively (Supplementary Fig. 2). As there were no overlapping target genes for miR-99a-3p, only the remaining 7 miRNAs were retained to construct the miRNA–target regulatory network. After obtaining the target gene list, GO and KEGG analysis were conducted to identify enriched biological processes of target genes. We found that target genes may play roles in several GI-related biological pathways, such as DNA − binding transcription activator activity, RNA polymerase II − specific, chromatin DNA binding, and cell cycle (Fig. 11A, B). Next, the interactions among target genes were investigated using a PPI network analysis. The network was inputted into Cytoscape and hub target genes were selected by CytoHubba. Finally, based on the connectivity degree, the top 5 hub target genes were as follows: COL1A2 (degree = 6), COL5A1 (degree = 5), COL11A1 (degree = 5), COL27A1 (degree = 5), COL12A1 (degree = 4) (Fig. 11C). Thus, our analyses showed that the 8 miRNAs in the GIMiSig may target these hub genes to affect biological processes in GC. Besides, all the top 5 hub genes had significantly higher expression in GC compared with that in normal samples from TCGA (p < 0.001) (Fig. 11D) and GSE54129 (p < 0.001) (Fig. 11E).

To the best of our knowledge, this is the first study to establish a GIMiSig and explore its predictive value and biological function in GC. The GIMiSig consisted of 8 GI-related miRNAs (miR-99a-3p, miR-548v, miR-100-5p, miR-380-3p, miR-363-3p, miR-125b-5p, miR-196b-3p, and miR-1275). The risk score generated by the GIMiSig was utilized to distinguish the high- and low-risk groups. The GIMiSig exhibited better predictive performance than current clinical classification methods, with the survival analysis and ROC curves showing that the GIMiSig was an effective tool for predicting survival. Of note, the GIMiSig was associated with frequency of TTN mutation but it predicted prognosis more accurately than TTN mutation status. Furthermore, the GIMiSig exerted an influence on the immune landscape and was also capable of predicting immunotherapy responses. Finally, GSEA revealed that the GIMiSig was closely related to cancer-associated, immune response, and GI-related pathways.

Six miRNAs (miR-99a, miR-100-5p, miR-363-3p, miR-1275, miR-196b-3p, and miR-125b-5p) in the GIMiSig have previously been reported in GC. miR-99a and miR-100-5p were regarded as tumor suppressors in GC due to their function of inhibiting cell proliferation (21, 22). miR-363-3p expression was downregulated in GC and its expression was associated with invasion and lymphatic metastasis (23). miR-1275 prevented metastasis and was associated with survival time in GC (24). Overexpression of miR-196b-3p was discovered in GC samples and affected the epithelial–mesenchymal transition (25). Upregulation of miR-125b-5p was observed in GC and led to invasion and metastasis by two possible pathways: (i) targeting STARD13 and NEU1 mRNAs and (ii) targeting the PPP1CA-Rb axis (26, 27). Specifically, elevated miR-125b-5p resulted in progression of GC and poor response to trastuzumab (28). Though the remaining 2 miRNAs (miR-548v, miR-380-3p) in the GIMiSig have not been reported in GC, they have been identified as potential biomarkers in multiple cancers, including endometrial cancer (29), lung adenocarcinoma (30), breast cancer (31), and neuroblastoma (32).

Interestingly, high somatic mutation and UBQLN4 expression, which indicated high GI, were observed in low-risk patients, who had longer survival time in our study. In these cases, more mutations of tumor driver genes occurred, but more neoantigens would also be presented to trigger immune response (33). Additionally, UBQLN4 can suppress the progression of GC by preventing tumor cell proliferation (34). This may partly explain why high somatic mutation number and UBQLN4 expression was detected in patients with longer survival time. TTN was known as a gene associated with familial hypertrophic cardiomyopathy (35). Yang et al found that TTN mutation was effective for predicting prognosis, TMB, and immunotherapy response in GC (36). Thus, we investigated the relationship between the GIMiSig and TTN mutation status. Our study suggested that the high- and low-risk groups (based on the GIMiSig) exhibited a significant difference in the frequency of TTN mutation, indicating that the GIMiSig was able to capture the TTN mutation status. More importantly, when the TTN gene was used to classify patients into wildtype and mutation groups, there were no significant differences in OS between TTN mutation and wildtype TTN among the high-risk patients or between TTN mutation and wildtype TTN among the low-risk patients. Conversely, the GIMiSig was capable of distinguishing distinct clinical prognoses among patients with TTN mutation and among patients without TTN mutation. This suggests that the GIMiSig is a better predictor of prognosis in GC than TTN mutation status.

In our study, the two groups divided based on the GIMiSig exhibited distinct immune landscapes. After extensively reviewing the literature, we found that miR-99a, miR-100-5p, miR-125b, and miR-363-3p in the GIMiSig were involved in the regulation of immune status, mainly through influencing the function of immune cells. Jaiswal et al showed that miR-99a downregulated TNF-α, resulting in an increased ratio of M2/M1 macrophages (37). Additionally, miR-99a promoted Treg differentiation and decreased cytotoxic T lymphocytes (38). Thus, miR-99a impaired the ability of immune cells to kill tumor cells. Similarly, miR-100-5p expression was also associated with macrophage polarization and Treg differentiation (39, 40). High miR-125b expression could disturb B cell development and suppress effector T cell bioactivities, while enhancing the pro-inflammatory nature of M1 macrophages (41). Of note, due to the crucial roles of miR-125b in the immune system, it demonstrated a promising ability to predict immunotherapy response in several tumors, such as non-small cell lung cancer (NSCLC) (42), prostate cancer (43), and colorectal cancer (44). As for miR-363-3p, it affected several crucial transcription factors that regulated Th17 cell differentiation (45). Combining previous studies and our findings, we speculate that the GIMiSig is suitable for evaluating the immune landscape in GC.

Another highlight of our study is the associations of the GIMiSig with known immunotherapy responsiveness biomarkers (PD-L1, MSI score, aneuploidy score, TMB, and TIDE score). Although ICIs have exhibited encouraging clinical trial results in several tumors (44, 46), their efficacy is not satisfactory in GC. Evidence shows that ICIs are only beneficial for specific GC subtypes, such as those involving deficiency mismatch repair (47). Therefore, predictive biomarkers are urgently needed to distinguish patients that could benefit from immunotherapy. As PD-L1-positive patients are more sensitive to ICIs, PD-L1 is regarded as an index of immunotherapy efficacy (47). High MSI is also associated with better immunotherapy responses in GC (48) while high aneuploidy score correlate with reduced response to immunotherapy (49). TMB reflects the quantity of somatic mutations and is associated with ICI responses in GC, and it can serve as a biomarker for predicting immunotherapy efficacy (50). Patients with high TMB tend to obtain more clinical benefit from ICI treatment (51). The TIDE score is a novel algorithm based on tumor immune escape, which provides clues for selecting patients that are suitable for ICI treatment. In this study, there was no significant difference in PD-L1 between the high- and low-risk patients. However, the low-risk patients had higher MSI scores, aneuploidy scores, and TMB, along with lower TIDE scores. Thus, we inferred that high-risk patients were more likely to benefit from immunotherapy. Moreover, these findings indicate that the GIMiSig could be exploited as a biomarker of immunotherapy responses.

There are several limitations in our study. First, the mechanisms underlying the associations between the GIMiSig and immune landscape, TMB, and ICI responses remain unclear. Second, the data on the cohort used in our study did not include data on GC patients taking immunotherapy, so the ability of the GIMiSig to predict immunotherapy efficacy remains to be elucidated. Third, biology experiments and clinical data are expected to testify the performance of this GIMiSig before clinical application.

In conclusion, we constructed a GIMiSig in GC based on GI-related miRNAs and investigated its prognostic value and potential function. Our study strongly suggests that the GIMiSig is not only of significance for predicting the prognosis and immune landscape of GC, but it may also provide new insights into future immunotherapy targets.

miRNAs: microRNAs

GI: genome instability

GIMiSig: GI-related miRNA signature

TCGA: The Cancer Genome Atlas

GS: genome stability

TTN: Titin

ICIs: immune checkpoint inhibitors

PD-1: programmed cell death-1

PD-L1: programmed cell death ligand-1

TME: tumor microenvironment

GEO: Gene Expression Omnibus

SMV: somatic mutation value

GO: gene ontology

CC: cellular component

MF: molecular function

BP: biological process

ROC: receiver operating characteristic

HR: hazard ratio

ssGSEA: gene set enrichment analysis

DCs : dendritic cells

aDCs: activated dendritic cells

iDCs: immature dendritic cells

pDCs : plasmacytoid dendritic cells

NK: natural killer

Tfh : follicular helper

Th1: type-1 T helper

Th2: type-2 T helper

TILs : tumor infiltrating lymphocytes

Tregs : regulatory T cells

MSI : microsatellite instability

TMB : tumor mutation burden

TIDE : Tumor Immune Dysfunction and Exclusion

KEGG: Kyoto Encyclopedia of Genes and Genomes

OS: overall survival

Ethics approval and consent to participate

Not applicable.

Consent for publication

All authors read and approved the final manuscript.

Availability of data and materials

The data of this study are available in the TCGA (https://portal.gdc.cancer.gov/) and GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This work was financially supported by the National Natural Science Foundation of China (No. 81603670, 81873169, 81774281, and 81303098), Hunan Provincial Natural Science Foundation of China (No. 2017JJ3459, 2020JJ4803, and 2019JJ40509), and China Scholarship Council (No. 202006370051).

Authors’ Contributions

W.J.P conceived and designed the study. Y.Q.L and W.J.P drafted the manuscript and analyzed the data. L.C collected the data and assisted with bioinformatics analysis. W.H and X.C contributed to interpreting data.

Author’s information

Yaqiong Liu¹, Lin Cheng¹, Wei Huang², Xin Cheng¹and Weijun Peng¹

Author details

Department of Integrated Traditional Chinese & Western Medicine, The Second Xiangya Hospital，Central South University, Changsha, Hunan ,410011,China.
Department of Integrated Traditional Chinese and Western Medicine, Xiangya Hospital, Central South University, Changsha, Hunan, 410008, China

Acknowledgments

The authors are grateful to the support of National Natural Science Foundation of China (No. 81603670, 81873169, 81774281 and 81303098 ), Hunan Provincial Natural Science Foundation of China (No. 2017JJ3459, 2020JJ4803, and 2019JJ40509), and China Scholarship Council (No. 202006370051). The authors also thank Lei Liu (Wuhan, Hubei) for helpful discussions.

Kamangar F, Dores G, Anderson WF. Patterns of cancer incidence, mortality, and prevalence across five continents: defining priorities to reduce cancer disparities in different geographic regions of the world. 2006;24(14):2137-50.
Orditura M, Galizia G, Sforza V, Gambardella V, Fabozzi A, Laterza MM, et al. Treatment of gastric cancer. 2014;20(7):1635.
Ribas A, Wolchok JD. Cancer immunotherapy using checkpoint blockade. 2018;359(6382):1350-5.
Dolcetti R, Re VD, Canzonieri V. Immunotherapy for gastric cancer: time for a personalized approach? 2018;19(6):1602.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. 2011;144(5):646-74.
Ottini L, Falchetti M, Lupi R, Rizzolo P, Agnese V, Colucci G, et al. Patterns of genomic instability in gastric cancer: clinical implications and perspectives. 2006;17:vii97-vii102.
Wang T, Wang G, Zhang X, Wu D, Yang L, Wang G, et al. The expression of miRNAs is associated with tumour genome instability and predicts the outcome of ovarian cancer patients treated with platinum agents. 2017;7(1):1-11.
Ng WL, Yan D, Zhang X, Mo Y, Wang Y. Over-expression of miR-100 is responsible for the low-expression of ATM in the human glioma cell line: M059J. 2010;9(11):1170-5.
Lal A, Navarro F, Maher CA, Maliszewski LE, Yan N, O'Day E, et al. miR-24 Inhibits cell proliferation by targeting E2F2, MYC, and other cell-cycle genes via binding to “seedless” 3′ UTR microRNA recognition elements. 2009;35(5):610-25.
Huang J-W, Wang Y, Dhillon KK, Calses P, Villegas E, Mitchell PS, et al. Systematic screen identifies miRNAs that target RAD51 and RAD51D to enhance chemosensitivity. 2013;11(12):1564-73.
Dinami R, Ercolani C, Petti E, Piazza S, Ciani Y, Sestito R, et al. miR-155 drives telomere fragility in human breast cancer by targeting TRF1. 2014;74(15):4145-56.
Chou J, Lin JH, Brenot A, Kim J-w, Provot S, Werb Z. GATA3 suppresses metastasis and modulates the tumour microenvironment by regulating microRNA-29b expression. 2013;15(2):201-13.
Ma L, Young J, Prabhala H, Pan E, Mestdagh P, Muth D, et al. miR-9, a MYC/MYCN-activated microRNA, regulates E-cadherin and cancer metastasis. 2010;12(3):247-56.
Fonseka P, Pathan M, Chitti SV, Kang T, Mathivanan S. FunRich enables enrichment analysis of OMICs datasets. 2020:166747.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. 2013;4(1):1-11.
Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. 2013;39(4):782-95.
Taylor AM, Shih J, Ha G, Gao GF, Zhang X, Berger AC, et al. Genomic and functional approaches to understanding cancer aneuploidy. 2018;33(4):676-89. e3.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. 2018;24(10):1550-8.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. 2012;16(5):284-7.
Jachimowicz RD, Beleggia F, Isensee J, Velpula BB, Goergens J, Bustos MA, et al. UBQLN4 represses homologous recombination and is overexpressed in aggressive tumors. 2019;176(3):505-19. e22.
Wang Z, Zhao Z, Yang Y, Luo M, Zhang M, Wang X, et al. MiR-99b-5p and miR-203a-3p function as tumor suppressors by targeting IGF-1R in gastric cancer. 2018;8(1):1-12.
Chen J-F, Wu P, Xia R, Yang J, Huo X-Y, Gu D-Y, et al. STAT3-induced lncRNA HAGLROS overexpression contributes to the malignant progression of gastric cancer cells via mTOR signal-mediated inhibition of autophagy. 2018;17(1):1-16.
Li C-Y, Liang G-Y, Yao W-Z, Sui J, Shen X, Zhang Y-Q, et al. Identification and functional characterization of microRNAs reveal a potential role in gastric cancer progression. 2017;19(2):162-72.
Mei J-W, Yang Z-Y, Xiang H-G, Bao R, Ye Y-Y, Ren T, et al. MicroRNA-1275 inhibits cell migration and invasion in gastric cancer by regulating vimentin and E-cadherin via JAZF1. 2019;19(1):1-12.
Pant N, Rakshit S, Paul S, Saha I. Genome-wide analysis of multi-view data of miRNA-seq to identify miRNA biomarkers for stomach cancer. 2019;97:103254.
Chang S, He S, Qiu G, Lu J, Wang J, Liu J, et al. MicroRNA-125b promotes invasion and metastasis of gastric cancer by targeting STARD13 and NEU1. 2016;37(9):12141-51.
Wu J-G, Wang J-J, Jiang X, Lan J-P, He X-J, Wang H-J, et al. MiR-125b promotes cell migration and invasion by targeting PPP1CA-Rb signal pathways in gastric cancer, resulting in a poor prognosis. 2015;18(4):729-39.
Sui M, Jiao A, Zhai H, Wang Y, Wang Y, Sun D, et al. Upregulation of miR-125b is associated with poor prognosis and trastuzumab resistance in HER2-positive gastric cancer. 2017;14(1):657-63.
Widodo, Djati MS, Rifa'i M. Role of MicroRNAs in carcinogenesis that potential for biomarker of endometrial cancer. 2016;7:9-13.
Lin K, Xu T, He B-S, Pan Y-Q, Sun H-L, Peng H-X, et al. MicroRNA expression profiles predict progression and clinical outcome in lung adenocarcinoma. 2016;9:5679.
Zhu H, Dai M, Chen X, Chen X, Qin S, Dai SJ. Integrated analysis of the potential roles of miRNA‑mRNA networks in triple negative breast cancer. 2017;16(2):1139-46.
Swarbrick A, Woods SL, Shaw A, Balakrishnan A, Phua Y, Nguyen A, et al. miR-380-5p represses p53 to control cellular survival and is associated with poor outcome in MYCN-amplified neuroblastoma. 2010;16(10):1134-40.
Vareki SM. High and low mutational burden tumors versus immunologically hot and cold tumors and response to immune checkpoint inhibitors. 2018;6(1):1-5.
Huang S, Li Y, Yuan X, Zhao M, Wang J, Li Y, et al. The UbL-UBA Ubiquilin4 protein functions as a tumor suppressor in gastric cancer by p53-dependent and p53-independent regulation of p21. 2019;26(3):516-30.
Gerull B, Gramlich M, Atherton J, McNabb M, Trombitás K, Sasse-Klaassen S, et al. Mutations of TTN, encoding the giant muscle filament titin, cause familial dilated cardiomyopathy. 2002;30(2):201-4.
Yang Y, Zhang J, Chen Y, Xu R, Zhao Q, Guo W, et al. MUC4, MUC16, and TTN genes mutation correlated with prognosis, and predicted tumor mutation burden and immunotherapy efficacy in gastric cancer and pan‐cancer. 2020;10(4).
Jaiswal A, Reddy SS, Maurya M, Maurya P, Barthwal M. MicroRNA-99a mimics inhibit M1 macrophage phenotype and adipose tissue inflammation by targeting TNFα. 2019;16(5):495-507.
Warth SC, Hoefig KP, Hiekel A, Schallenberg S, Jovanovic K, Klein L, et al. Induced miR‐99a expression represses Mtor cooperatively with miR‐150 to promote regulatory T‐cell differentiation. 2015;34(9):1195-213.
Wang W, Liu Y, Guo J, He H, Mi X, Chen C, et al. miR-100 maintains phenotype of tumor-associated macrophages by targeting mTOR to promote tumor metastasis via Stat5a/IL-1ra pathway in mouse breast cancer. 2018;7(12):1-17.
Eniafe J, Jiang S. MicroRNA‐99 family in cancer and immunity. 2020:e1635.
Wang J, Wang Z, Li G. MicroRNA-125 in immunity and cancer. 2019;454:134-45.
Peng X-X, Yu R-Y, Wu X, Wu S-Y, Pi C, Chen Z-H, et al. Correlation of plasma exosomal microRNAs with the efficacy of immunotherapy in EGFR/ALK wild-type advanced non-small cell lung cancer. 2020;8(1).
Triozzi PL, Achberger S, Aldrich W, Elson P, Garcia J, Dreicer, et al. Differential immunologic and microRNA effects of 2 dosing regimens of recombinant human granulocyte/macrophage colony stimulating factor. 2012;35(7):587-94.
Tanaka H, Hazama S, Iida M, Tsunedomi R, Takenouchi H, Nakajima M, et al. miR‐125b‐1 and miR‐378a are predictive biomarkers for the efficacy of vaccine treatment against colorectal cancer. 2017;108(11):2229-38.
Kästle M, Bartel S, Geillinger‐Kästle K, Irmler M, Beckers J, Ryffel B, et al. micro RNA cluster 106a~ 363 is involved in T helper 17 cell differentiation. 2017;152(3):402-13.
Stinchcombe TE. Unmet needs in squamous cell carcinoma of the lung: potential role for immunotherapy. 2014;31(5):960.
Fuchs CS, Doi T, Jang RW, Muro K, Satoh T, Machado M, et al. Safety and efficacy of pembrolizumab monotherapy in patients with previously treated advanced gastric and gastroesophageal junction cancer: phase 2 clinical KEYNOTE-059 trial. 2018;4(5):e180013-e.
Kim ST, Cristescu R, Bass AJ, Kim K-M, Odegaard JI, Kim K, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. 2018;24(9):1449-58.
Davoli T, Uno H, Wooten EC, Elledge SJ. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. 2017;355(6322).
Samstein RM, Lee C-H, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. 2019;51(2):202-6.
Wang F, Wei X, Wang F, Xu N, Shen L, Dai G, et al. Safety, efficacy and tumor mutational burden as a biomarker of overall survival benefit in chemo-refractory gastric cancer treated with toripalimab, a PD-1 antibody in phase Ib/II clinical trial NCT02915432. 2019;30(9):1479-86.

SupplementaryMaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Genome Instability-related MiRNAs Predict Survival, Immune Landscape, and Immunotherapy Responses In Gastric Cancer

Status:

Version 1

Abstract

Figures

Introduction

Methods

Data collection

Identification of GI-related miRNAs

Construction of GI-related miRNA signature (GIMiSig)

Minimally invasive diagnostic value of the GIMiSig

Relationship between the GIMiSig and Titin (TTN) mutation status

Immune landscape analysis and prediction of immunotherapy response

Biological significance of the GIMiSig

Target genes of the GIMiSig

Statistical analysis

Results

Different genomic alterations in the GS and GI groups

Establishment of a GI-related miRNA signature (GIMiSig)

Evaluation of independent prognostic factors

Survival analysis of GC subgroups

Minimally invasive diagnostic value of the GIMiSig

Relationship between the GIMiSig and TTN mutation status

Relationship between the GIMiSig and immune landscape

Prediction of immunotherapy responsiveness

GSEA

Target genes of the GIMiSig

Discussion

Conclusion

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Conflict of Interest

Funding

Authors’ Contributions

Author’s information

Author details

Acknowledgments

References

Supplementary Files

Status:

Version 1