The graphical abstract presents the flow chart and online resources for this study. The original data of RNA-seq (FPKM format), miRNA-seq, corresponding clinical characteristics, including gender, age, overall survival (OS), and stage, and simple nucleotide variation (SNV) were downloaded from the TCGA-COAD project (December 24, 2020). The tumor mutational burden (TMB) and MSI data were retrieved from cBioportal (Cerami et al., 2012). The clustering method was used to detect and exclude the outlier RNA-seq samples. Eventually, 453 RNA-seq samples were acquired; however, only 446 patients with both matching RNA-seq data and corresponding clinical features were included in this study. All biological and clinical samples in this study are publicly available. The data availability policies of the open-accessed databases were strictly followed. All code utilized in this study can be acquired via email 201800413040@mail.sdu.edu.cn with reasonable grounds.

The R package “limma” was used to extract differential expression genes (DEGs). All tumor RNA-seq samples were sorted according to their SNV numbers, from largest to smallest. The first 25% and last 25% of samples were named the genome-unstable (GU) group and genome-stable (GS) group. Differentially expressed lncRNAs between the GS and GU groups were considered GIRlncRs. The cutoff value was adjusted a p-Value (adj. p) < 0.05, and | log₂ fold-change| > 1. Similarly, genome instability-related mRNA (GIRmR) and miRNA (GIRmiR) were extracted.

The “sparcl” package was implemented to perform unsupervised hierarchical clustering based on GIRlncRs (k = 2). The clusters with higher and lower SNV numbers were named GU-like and GS-like clusters, respectively. We compared the tumor mutational burden (TMB), MSI, expression, and mutation rates of four MMR genes, expression of six immune-related genes, and the mutation rates of POLE (a novel biomarker for immunotherapy outcomes (Wang et al., 2019)) between the two clusters. The MMR and immune-related genes included in this study are listed in Table 1. Single-sample gene set enrichment analysis (ssGSEA) was used to compare tumor microenvironment (TME) infiltration and functions between the two clusters (Supplementary Table 1). Survival analysis using the log-rank test between the GS-like and GU-like clusters was performed to evaluate the prognostic value of clustering. By calculating the Pearson’s correlation coefficients, we ranked the relevance of associations between each mRNA and GIRlncR expression. The first 10 ranked mRNAs were considered coexpressed with GIRlncR. Cytoscape, a biological network modifying software, was used to visualize the lncRNA-mRNA coexpression network (Shannon et al., 2003). To explore the potential functions of GIRlncRs, functional gene enrichment analysis based on the Gene Ontology (GO) (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) databases was performed.

Multivariate Cox regression is the most widely utilized regression model for analyzing medical survival time and survival status data (Bradburn et al., 2003), and it was used to explore the association with event incidence (surviving/deceased in this study). The hazard ratio (HR) represents the probability of event occurrence under the currently observed feature patterns (gene expression pattern in this study). First, we divided 446 patients equally into the training and test cohorts. Chi-square tests were used to detect selection bias for each feature during patient division. Second, univariate Cox regression analysis of each GIRlncR was used for feature selection. Only statistically significant features (p < 0.05) were used for model construction. Third, multivariate Cox regression was performed to construct the genome instability-related lncRNA prognostic signature (GIlncPS) using the training cohort. The constructed model is as follows:

The model was applied to the test cohort. In both training and test cohorts, patients who had a risk score higher than the median value of the risk score were classified as the high-risk group and vice versa. Survival analysis and receiver operating characteristic (ROC) curves and the area under the curve (AUC) were utilized to evaluate the reliability of the GIlncPS. Moreover, given that BRAF is a frequently mutated gene with prognostic value (Sanz-Garcia et al., 2017), we tested whether our signature has better predictive performance than this event. Survival analysis among BRAF mutation-type/high-risk, BRAF mutation-type/low-risk, BRAF wild-type/high-risk, AND BRAF wild-type/low-risk groups was performed using the log-rank test. To demonstrate that our signature has a better predictive performance than other lncRNA signatures, we found two proposed lncRNA signatures for COAD, namely, Li’s signature (Li Z. et al., 2020) and Jin’s signature (Jin et al., 2020). AUC was used to rank the predictive performance of these signatures according to OS at 1, 3, and 5 years for 446 patients. To demonstrate that our signature has prognostic ability independent of general factors, we performed multivariate Cox regression, including clinical features (age, gender stage, T (tumor), M (metastasis), N (lymph node), KRAS mutation type, TP53 mutation type, BRAF mutation type, and POLE mutation type) and the risk score. The nomogram, including our signature, was plotted for clinical reference, which was evaluated using a calibration curve.

To investigate the relationship between the clusters and risk groups, Sankey plots were constructed. Waterfall plots were generated to explore the diversity of mutation profiles between GS- and GU-like clusters. To extract lncRNAs that had the closest relationship with clustering, LASSO and SVM algorithms were applied. The lncRNAs included in LASSO-screened lncRNAs, SVM-screened lncRNAs, and prognostic signatures were considered hub lncRNAs.

The hub-lncRNA and target mRNA interactions were predicted using an online comprehensive RNA database, ENCORI¹ (Li et al., 2014). The overlapping mRNAs in GIRmRs and hub-lncRNA target mRNAs require further research. Only the statistically significant lncRNA-mRNA pairs (p < 0.05) in COAD remained for further research. We then retrieved the miRNAs mediating the relationships between the lncRNAs and mRNAs. The retrieved miRNAs that were differentially expressed between the GS and GU groups were retained for ceRNA network construction. Figure 4G illustrated the screening process of ceRNA network.

It is hypothesized that the genes included in the ceRNA network play a pivotal role in COAD. A comprehensive single-gene analysis was performed. Differential expression analysis was performed between normal and tumor tissues from COAD patients by employing the Mann-Whitney U test. Thorsson et al. (2018) identified six immune subtypes across 33 TCGA cancer types (C1–C6) and compared the expression of each gene among these six subtypes using the Kruskal-Wallis test. Spearman correlation analysis was performed to detect the correlation between each gene and TME infiltration, the expression of MMR genes, immune-related genes (Table 1), and two previously discovered crucial lncRNAs regulating genome instability (NOARD (Munschauer et al., 2018), GUARDIN (Hu et al., 2018)). The TME was evaluated by the estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) immune data, stromal score (Yoshihara et al., 2013), and cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) and specific immune cell types (Newman et al., 2019). Given that genome instability contributes to tumor heterogeneity, which leads to drug resistance (Morel et al., 2017; Dagogo-Jack and Shaw, 2018; Sansregret et al., 2018), we performed cancer stemness and drug sensitivity analysis for these genes. Malta et al. (2018) identified a novel index to evaluate cancer stemness features based on DNA methylation patterns (DNA stemness score (DNAss)) and mRNA expression patterns (RNA stemness score (RNAss)). We calculated the Spearman correlation between each gene and cancer stemness. Drug sensitivity analysis was also performed using an open-access database, specifically National Cancer Institute (NCI)-60. The Kolmogorov-Smirnov test was applied to verify the normality of the included indexes, including TMB/MSI scores, NCI-60 indexes, TME indexes, and cancer stemness indexes (Supplementary Table 2). We found that the drug sensitivity index was normally distributed. Therefore, the method of Pearson’s correlation to explore the association between gene expression and drug sensitivity is appropriate.

Total RNA was isolated from COAD cell lines named Caco-2, Lovo, HCT116, and HT29 using the RNeasy Mini Kit (QIAGEN, Hilden, Germany), according to the manufacturer’s instructions. Complementary DNA (cDNA) was produced by RNA using the PrimeScript^TM Reverse Transcription Kit (TakaRa, Maebashi, Japan) in an ABI 7500 System (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, United States). The primers specific for target genes were designed and synthesized by Tianyi Huiyuan Biotech (Beijing, China). The following procedures were performed: activation of enzymes at 50°C for 2 min and then 95°C for 2 min, 45 cycles of denaturation at 95°C for 15 s, annealing at 58°C for 20 s, and extension at 72°C for 30 s. The relative expression levels of the target genes were calculated using the 2^−ΔΔCT method. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the internal control. All qPCR reactions were performed in triplicate. The primers used in this study are listed below.

Specific somatic mutations drive cancer development (Koch, 2017). Therefore, the mutation analysis of the target mRNAs was applied, including somatic mutations and copy number variation (CNV) analysis. The R package “RCircos” was used to explore the chromosome location and CNV number for each gene. We then explored the alteration rate, type, and site in the COAD datasets and examined the association between the alteration status and gene expression. We calculated the Spearman’s correlation TMB and MSI scores and assessed the expression of each gene to explore whether the target mRNAs were directly associated with genome instability. Furthermore, the association between target mRNA alteration status and MSI status was explored. All mutation analyses were performed using the online database cBioportal.

SPSS Statistics (version 26.0; IBM, Armonk, NY, United States) and Microsoft Excel (Microsoft Corporation, Redmond, WA, United States) were used to present the retrieved data and perform the Chi-square test. R language (version 4.0.3) was used to perform machine learning algorithms and correlation analysis (R Core Team, 2013). The Benjamini and Hochberg method was used to adjust the p-Value for DEG extraction.

In total, 123 GIRlncRs were extracted between the GS and GU groups (Figure 1A). We constructed a mRNA and GIRlncR coexpression network (Figure 1B) to show the top 10 close interactions between each GIRlncR and corresponding mRNAs, indicating that the GIRlncRs regulate genome instability and complex tumor cell biological processes. Functional analysis further illustrated the relationship between GIRlncRs and tumor immunity (Figures 1C,D). Almost all GO- and KEGG-enriched pathways involved immune responses, such as antigen processing and T-cell activation. The unsupervised hierarchical clustering based on GIRlncRs was employed to divide COAD patients into two clusters with various genome instability statuses and immune therapy response rates. In total, 453 samples were clustered into GS-like (302 samples) and GU-like groups (151 samples; Figure 1E). Surprisingly, the TMB, MSI score, and expression of all immune checkpoint-related genes were higher in the GU-like cluster (p < 0.05), whereas expression of the three MMR genes (MLH1, PMS2, MSH2) was lower (p < 0.05; Figures 1F–I). The mutation rates of all four MMR genes were higher in the GU-like group (p < 0.05). The ssGSEA showed that almost all immune cells and immune functions were increased in the GU-like group (Figures 1J,K). Furthermore, the mutation rate of the novel immunotherapy response rate biomarker POLE was also significantly higher in the GU-like group (p < 0.05; Figure 1L). These findings strongly suggest that our clustering based on GIRlncRs could predict TME infiltration and ICI response rates in COAD patients. However, as shown in Figure 1M, the OS of the two clusters did not differ, indicating that the clustering method could not fulfill the purpose of prognosis.

Since our clustering could not predict OS in COAD patients, we constructed a novel GIlncPS. The basic characteristics of the training and test cohorts are shown in Supplementary Table 3. The Chi-square test showed that the division of training and test cohorts was proper without selection bias (all p > 0.05). The regression results are shown in Supplementary Table 4. An eight-lncRNA GIlncPS was acquired as follows: risk score = 0.239 × LINC01807 + 0.235 × AC009237.14 + 0.310 × LOXL1-AS1 + 0.028 × AC005392.2 + 0.751 × AP003555.1 + 0.384 × BOLA3-AS1 – 0.330 × PTPRD-AS1 + 0.258 × AC004009.1. The high- and low-risk groups were divided according to median risk score. Survival analysis indicated that OS was significantly longer in the low-risk group in the training, test, and entire TCGA cohorts (p < 0.05; Figures 2A–C). Notably, all the included lncRNAs, except PRPRD-AS1, were risk factors for poor prognosis in COAD patients. The heatmap shows that all of these lncRNAs, except PRPRD-AS1, were overexpressed in high-risk groups because only PRPRD-AS1 was a protective factor (Figures 2D–F). The AUCs of the ROC in the three cohorts were acceptable for the lncRNA prognostic signature. Only the AUC at 1 year in the test cohort did not reach 0.7 (Figures 2G–I). Compared with the prognostic models of Li and Jin, our GIlncPS had the best performance (AUC = 0.735, 0.745, and 0.741, respectively; Figures 2J–L). Independent prognostic analysis indicated that old age, advanced stage, and occurrence of metastasis were significantly related to poor prognosis, and a low-risk score was significantly related to good prognosis (p < 0.05; Figures 3A,B). Figure 3C shows the nomogram for clinical use. The calibration curve demonstrated the reliability of the nomogram (Figure 3D). We also compared TMB/MSI and TME infiltration and functions between the high- and low-risk groups. However, the results differed from those of hierarchical clustering. As shown in Supplementary Figures 1A–H, all ICI outcome predictive biomarkers, except MSI, were not significantly different between the two groups, indicating that our prognostic signature behaved poorly in predicting genome instability, TME infiltration, and ICI response rates for COAD patients. Furthermore, Supplementary Figure 1I shows that the mutation rate of BRAF was higher in the high-risk group (19% vs. 11%). The prognostic ability of the GIlncPS was better than that of BRAF (p < 0.01; Supplementary Figure 1J). The OS in the high-risk groups was always lower than that in the low-risk groups, regardless of the BRAF mutation type. To conclude, we constructed a GIlncPS with strong prognostic ability, but it could not predict the TME infiltration and ICI response rates for COAD patients.

The Sankey plot showed the relationship among the GS-and GU-like clusters, risk groups, and survival status (Figure 4A). The results showed no significant difference between the clusters and risk groups, and the death proportion in the high-risk group was higher. The top three genes with the highest mutation frequencies in the GS-like cluster were APC (86%), TNN (40%), and TP53 (65%) (Figure 4B), whereas those in the GU-like cluster were APC (53%), TTN (69%), and TP53 (31%) (Figure 4C). In general, the mutation frequency in the GU-like cluster was higher than that in the GS-like group, and the frequencies of multi-hit and frame-shift mutations were much higher. LASSO and SVM algorithms were utilized to extract hub-lncRNAs in GS- and GU-clusters for further investigation. We selected the minimum value of λ, a penalty parameter of 0.1, and the L1 norm of the default (Figures 4D,E). Here, 43 lncRNAs were screened using LASSO regression. The average weight of each lncRNA acquired from the 10-fold test SVM algorithm is provided in Supplementary Table 5. The cut-off value was set at | 0.28|. The Venn plot of the selection process of hub-lncRNA among SVM, LASSO, and the prognostic signature is presented in Figure 4F. In total, 55 hub lncRNAs were identified. The screening process of the lncRNAs, miRNAs, and mRNAs included in the ceRNA network is presented in Figure 4G.

The ENCORI database was utilized to predict the interactions between hub lncRNAs and their target mRNAs. The constructed ceRNA network is shown in Figure 5A and Table 2. Seven lncRNA-miRNA-mRNA pairs were found in the ceRNA network. To comprehensively understand these three lncRNAs (BOLA3-AS1, AC124067.4, and AL161772.1) and four target RNAs (NINL, SALL4, TNFSF11, PHYHIPL), we performed differential expression, cancer stemness, immune subtype, TME infiltration, and drug sensitivity analyses.

First, we compared the expression of seven genes between normal and tumor tissues using the Mann-Whitney U test in COAD. As shown in Figures 5B–H, AC124067.4, AL161772.1, SALL4, and TNFSF11 were highly expressed in COAD, whereas BOLA3-AS1 and PHYHIPL were expressed at low levels (p < 0.05). These results indicate that AC124067.4, AL161772.1, SALL4, and TNFSF11 might promote COAD development and that BOLA3-AS1 and PHYHIPL might inhibit its development (Figure 5I). As expected, the four lncRNA-mRNA expression patterns met the ceRNA criterion, especially AC124067.4-PHYHIPL (r = 0.309, p = 6.88e-12) and BOLA3-AS1-SALL4 (r = 0.377, p = 2.56e-17). It has been reported that genome instability is closely related to tumor immunity. ANOVA was used to detect whether the expression of the seven genes was significantly related to the immune subtypes in COAD (Figure 5J). The results showed that six genes (BOLA3-AS1, AC124067.4, AL161772.1, SALL4, PHYHIPL, and TNFSF11) were significantly associated with the COAD immune subtype (p < 0.05), indicating that these seven genes could modulate tumor immune functions at the genetic level. These seven genes were also closely related to genome instability and immune-related genes (Supplementary Figure 2). All seven genes were positively related to MLH1, PMS2, and NORAD (p < 0.01). In addition, most genes (4/7, 57.1%) were negatively correlated with PD-L1 (CD274; p < 0.05). These findings suggest that the high expression of these seven genes contributes to genomic integrity and predicts poor ICI outcomes. Moreover, these seven genes are closely related to COAD TME infiltration and function (Figure 5K). Most genes were positively related to stromal cell infiltration (5/7, 71.4%) and immune cell infiltration (4/7, 57.1%) and negatively related to tumor purity (5/7, 71.4%). According to the Spearman’s correlation between cancer stemness features and these seven genes, we found that most genes were negatively related to cancer stemness (5/7 (71.4%) for RNAss and 4/7 (57.1%) for DNAss; Figure 5K). Finally, drug sensitivity analysis was performed to explore the relationship between the seven genes and drug resistance. The results showed that the expression of these seven genes was positively associated with the sensitivity of drugs regulating DNA replication and synthesis. For example, TNFSF11 expression was positively associated with sensitivity to oxaliplatin, a common platinum-based antineoplastic drug for colorectal cancer, which is believed to function by blocking the duplication of DNA (r = 0.429, p < 0.001). Meanwhile, some protein kinase inhibitors, such as perifosine and cobimetinib, were negatively related to these genes. In conclusion, we proposed a novel ceRNA network that modulates genome instability and found that these seven genes are involved in TME infiltration, cancer stemness, and drug resistance.

To ascertain these potential prognostic biomarkers, three lncRNAs (AC124067.4, AL161772.1, BOLA3-AS1) with significant expression alterations and one mRNA SALL4 were further performed by qRT-PCR validation. The results showed that AC124067.4, AL161772.1, and SALL4 gene levels were significantly increased in COAD cell lines (Figures 6A–D) whereas BOLA3-AS1 (Figure 6) was significantly reduced, which was consistent with the above analysis. Taken together, the expressions of AC124067.4, AL161772.1, and SALL4 were significantly upregulated, and their expression level was also associated with the OS rate in patients with COAD. These findings indicate the possibility of using the panel as a prognostic biomarker for COAD.

First, we explored the rates of alteration and types with respect to the four target mRNAs. As shown in Figures 7A–D, the rates of NINL, PHYHIPL, SALL4, and TNFSF11 alteration were 5%, 0.8%, 6%, and 2%, respectively. For NINL and SALL4, somatic mutations and CNVs accounted for almost half of these for each. Regarding PHYHIPL alterations, somatic mutations were predominant, whereas CNVs formed the majority of TNFSF11 alterations. Missense mutations accounted for most somatic mutations among the four gene alterations, and amplification accounted for most CNVs among the NINL, SALL4, and TNFSF11 alterations. The alteration sites of each target gene were also explored. Notably, all four gene alteration sites were dispersed across the entire gene, rather than clustering on specific domains, indicating that the four-gene alterations mainly affect mRNA expression but not the coded protein activity (Figures 7E–H). Therefore, we investigated whether alterations in the four genes were associated with their expression. As expected, three genes (NINL, SALL4, and TNFSF11) were expressed at low levels for mutated types, and the CNVs of the three genes were positively related to mRNA expression. However, the rate of PHYHIPL alterations was only 0.8%, and PHYHIPL expression was not related to its alteration type.

Microsatellite instability has become a recognized biomarker for immunotherapy and has been widely used in clinical practice. It is speculated that the seven genes proposed in this study could serve as clinically novel TMB and MSI biomarkers. As shown in Figures 8A–D, we explored the association between each gene and TMB/MSI in the pan-cancer database (the abbreviations are provided in Supplementary Table 6). Surprisingly, all the expression of all seven genes was negatively related to TMB and MSI scores in COAD (p < 0.05), indicating that these seven genes could maintain genomic integrity in COAD, which was particularly apparent for AC124067.4 (correlation coefficient R reached approximately –0.4 for TMB and –0.5 for MSI) and SALL4 (R reached approximately –0.3 for TMB and –0.35 for MSI; Figures 8A–D). Moreover, we compared MSI status between patients with mutated and wild-type disease using cBioPortal. The MSI-H proportions were greater in patients with the mutated type than in those with the wild type according to all four genes (Figure 8E). Especially, the MSI-H proportion reached almost 60% in patients with mutated SALL4, whereas only 10% of wild-type patients had MSI-H. We next explored the association between the seven CNVs and the MSI status. As shown in Figure 8F, the CNVs of NINL, SALL4, and TNFSF11 were fewer in patients with MSI-H COAD. In conclusion, the seven genes proposed in this study were closely related to genome instability in COAD and could serve as novel immunotherapy outcome biomarkers.