Twenty protein-coding genes that functioned as m6A regulators[28] were collected, i.e., 11 readers (YTHDF2, YTHDF3, YTHDC1, YTHDC2, YTHDF1, RBMX, HNRNPC, HNRNPA2B1, IGF2BP1, IGF2BP2, and IGF2BP3), seven writers (METTL3, METTL14, RBM15, RBM15B, WTAP, VIRMA, and ZC3H13), and two erasers (ALKBH5 and FTO). We obtained RNA-Seq expression data (including FPKM and read count) and clinical data on 622 CRC (including colon cancer and rectal cancer) patients from the TCGA project (https://xenabrowser.net/datapages/). To validate the prognostic model, we additionally obtained six CRC datasets from the Gene Expression Omnibus (GEO), i.e., GSE17538 (210 patients), GSE39582 (557 patients), GSE33113 (89 patients), GSE31595 (33 patients), GSE29621 (53 patients), and GSE17536 (145 patients), totaling 1,077 CRC patients. They were from the GPL570 platform of U133 plus 2 arrays, which was suitable for probe annotation to obtain lncRNA expression. Gencode.v34 was used for lncRNA annotation.

The differentially expressed lncRNAs were identified by comparing the expression profiles between tumor and normal samples. For the expression data (read count) detected by RNA-seq, we performed the differential expression analysis using R package DESeq2 [29] with FDR≤0.05 and fold change ≥2 or ≤1/2. In order to make the lncRNAs with enough expression and detectable by array, we only kept the differentially expressed lncRNAs with high expression (median FPKM>1) and with probe annotation for the GPL570 platform. Then, m6A-related lncRNAs were determined based on M6A2Target [30] and expression correlation by using four criteria as follows:1) lncRNAs were methylated or demethylated by m6A writers or erasers; 2) or lncRNAs were binding to m6A readers; 3) or the expression level of lncRNAs was influenced by over-expression or knock down of m6A regulators recorded in M6A2Target; 4) and lncRNAs were co-expressed with at least one m6A regulator in the TCGA CRC dataset (p value < 0.05 and Pearson’s coefficient >0.2 or < −0.2).

For m6A-related lncRNAs, we utilized univariate Cox regression analysis to determine candidate factors for PFS. Based on the candidate lncRNAs, we performed LASSO analysis to get succinct and effective prognostic lncRNAs. LASSO analysis was implemented with functions cv.glmnet and glmnet in R package glmnet. The lncRNAs with LASSO regression coefficient not equal to 0 were retained. CRC patients were stratified based on lncRNA expression above or below the median. The survival curves were plotted using the Kaplan-Meier method, and the survival difference of two patient groups was estimated with the log-rank test (p value < 0.05).

The m6A-lncRNA signature model was established with a formula: m6A-LncScore = 0.32* SLCO4A1-AS1 expression +0.41* MELTF-AS1 expression +0.44* SH3PXD2A-AS1 expression +0.39*H19 expression +0.48* PCAT6 expression, where the figures before lncRNAs represent regression coefficients in univariate Cox regression analysis.

CRC patients were stratified into two groups based on whether m6A-LncScore was above or below the median. Receiver operating characteristic (ROC) curve analysis and Area Under Curve of ROC (AUC) was utilized to show prediction power according to m6A-LncScore and other factors. Multi-variate Cox regression analysis was employed to determine the independent prognostic factors for PFS with adjustment for other potential clinicopathologic factors, i.e., age, gender, tumor stage, AJCC-T, AJCC-N, and AJCC-M. A nomogram and calibration plot were adopted to display the predictive ability and power of multiple features using R package rms. The model selection for the nomogram was performed by a backward step-down selection process using a threshold of p value < 0.05. Calibration curves were used to assess the calibration of the nomogram, accompanied by the Hosmer-Lemeshow test.

Our in-house CRC cohort included 55 pairs of fresh specimens from CRC patients (tumor and matched adjacent normal tissue) without radiotherapy or chemotherapy, which were immediately stored in liquid nitrogen after surgery (Supplementary Table S1). All specimens were collected from Zhengzhou Central Hospital affiliated with Zhengzhou University between 2019 and 2020 and this study was approved by the Zhengzhou Central Hospital affiliated with Zhengzhou University. All subjects underwent rigorous screening and provided informed consent.

Quantitative RT-PCR (qRT-PCR) was employed to detect the RNA expression of the five lncRNAs (SLCO4A1-AS1, MELTF-AS1, SH3PXD2A-AS1, H19, and PCAT6). In brief, total RNAs of 55 pairs of tissue specimens were extracted using the Trizol method. After testing for concentration, purity, and integrity, an equal number of RNAs was used to synthesize cDNA. Finally, the SYBR Green Quantitative Kit (DBI, Germany) and 7500 Fast Quantitative PCR System (AB, United States) were used for detection. The housekeeping gene GAPDH was used as an internal reference, and the relative gene expression was expressed as 2^−ΔΔCt. Primer sequences are shown in Supplementary Table S2.

Compared with 51 normal adjacent samples, 3452 differentially expressed lncRNAs were identified in 622 CRC tumor samples, which comprised 2212 up-regulated and 1240 down-regulated lncRNAs (Supplementary Figure S1). Only 157 lncRNAs with high expression (median FPKM>1) remained. In order to enable lncRNAs to be verified by other datasets, we only focused on 43 recurrent lncRNAs (Figure 1A and Supplementary Table S3), whose expression could be also detected by the GPL570 platform. Interestingly, 18 of the 43 lncRNAs could interact with m6A regulators in NPInter V4 [31] (Supplementary Table S4); 32 of 43 lncRNAs and 41 of 43 lncRNAs could act as miRNA sponges and indirectly regulate m6A regulators via miRNAs in starBase V3 [32] and DIANA-LncBase V3 [33], respectively. For example, lncRNA PCAT6 bound to RBM15 and IGF2BP3 in several cancer cell lines, which were determined by CLIP and eCLIP technology.

Considering that these lncRNAs may act as the targeting genes of m6A modification, we further identified 24 m6A-related lncRNAs. They could receive m6A modification and bind to m6A readers, or their expression could be influenced by over-expression or knock down of m6A regulators in the M6A2Target database (Supplementary Tables S5–S7). Meanwhile, they were significantly co-expressed with at least one m6A regulator in the TCGA dataset (Supplementary Table S8). Notably, PCAT6 was significantly co-expressed with 12 m6A regulators (one positive and 11 negative relationships), suggesting its role in m6A RNA methylation (Supplementary Figure S2).

For PFS, univariate Cox regression analysis identified five prognostic lncRNAs (SLCO4A1-AS1, MELTF-AS1, SH3PXD2A-AS1, H19, and PCAT6) (Table 1). LASSO analysis suggested they could form the simplest and most effective combination for predicting PFS (Supplementary Figure S3, Supplementary Table S9). Compared with normal samples, the five lncRNAs in CRC tumors had obviously higher expression (Figure 1B). We subsequently detected their expression status in our in-house CRC cohort by doing a qRT-PCR assay. Compared with matched adjacent normal tissues, their RNA expression was obviously up-regulated in most of the 55 tumors, and the difference was highly statistically significant (Figure 1C, Supplementary Table S10). The patients with high expression had a significantly shorter PFS time than other patients (Figure 2A).

The m6A-Lnc signature was established with a formula: m6A-LncScore = 0.32* SLCO4A1-AS1 expression +0.41* MELTF-AS1 expression +0.44* SH3PXD2A-AS1 expression +0.39*H19 expression +0.48* PCAT6 expression. We calculated m6A-LncScore of 622 patients and divided the patients into two groups based on whether they scored above or below the median. Patients at high risk had a significantly shorter PFS time than those at low risk (Figure 2B). The high-risk patient group had higher lncRNA expression than the other group (Figure 2C). We also observed that m6A-LncScore was significantly co-expressed with 11 m6A regulators (two positive and nine negative relationships, Supplementary Figure S4).

Furthermore, we obtained the gene expression data of 1,077 CRC patients, including 210 patients from GSE17538, 557 patients from GSE39582, 89 patients from GSE33113, 33 patients from GSE31595, 53 patients GSE29621, and 145 patients from GSE17536, and validated the prognostic model in these six independent datasets (Figures 2D,E). Using the lower or upper quartile as the threshold, we also observed the statistical significance in most datasets (Supplementary Figures S5, S6), suggesting the robustness of the m6A-Lnc signature using different thresholds to classify patients as high or low risk. In addition, we found that the m6A-Lnc signature was also suitable for predicting overall survival (OS) in the TCGA CRC and COAD datasets (Supplementary Figure S7).

To estimate whether m6A-LncScore could act as an independent factor in CRC to predict PFS, we performed univariate and multivariate Cox regression analysis. The results showed that m6A-LncScore (p < 0.0001; HR = 1.43, CI = 1.26–1.62) and several clinicopathological factors (tumor stage, AJCC‐T, AJCC‐N, AJCC‐M and CEA level) were significantly relevant to patient PFS in the TCGA dataset (Table 2). Considering that the information on tumor stage was relatively complete in the TCGA dataset and the other three datasets, while other factors in many patients were missing or even unrecorded, we only adopted tumor stage into the multivariate Cox regression analysis. After adjustment by tumor stage, m6A-LncScore was still significant (p = 0.0021; HR = 1.25, CI = 1.08–1.44) (Table 3), indicating its independent prognostic potential. The independent prognostic value of m6A-LncScore was validated in three other datasets (GSE17538, GSE39582, and GSE33133) (Supplementary Tables S11, S12).

The result of nomogram analysis also showed the good predictive ability of m6A-LncScore, as well as clinicopathological factors (Figure 3A). Since tumor stage contained the complete information on AJCC‐T, AJCC‐N, AJCC‐M, the integrated model combining m6A-LncScore with independent prognostic factors (tumor stage and risk score) was further established. We found that the AUC of m6A-LncScore was 0.75, 0.73, 0.76 (for 1‐, 3‐, and 5‐year PFS, respectively), the AUC of tumor stage was 0.7, 0.75, 0.72, while the AUC of the model integrating m6A-LncScore with tumor stage was 0.79, 0.81, 0.82. The results suggested that the integrated model for predicting PFS was superior to m6A-LncScore or tumor stage (Figure 3B, Supplementary Tables S13, S14).

Furthermore, the calibration plot showed good consistency between observation and predictive values for 1‐, 3‐, and 5‐year PFS (Figure 3C). The ROC analysis in GSE17538, GSE39582, and GSE33113, also confirmed that m6A-LncScore had high accuracy in predicting patient PFS (Supplementary Figures S8–S10). In the TCGA dataset, even considering CEA level, we found that the integrated model for predicting PFS was superior to m6A-LncScore or CEA level (Supplementary Figure S11, Supplementary Table S15).

We found some clinicopathological factors had a significant association with m6A-LncScore, especially tumor stage, AJCC‐T, AJCC‐N, and AJCC‐M (Figure 4A). When the patients were stratified by these factors, m6A-LncScore was still statistically significant for patients when comparing the high- and low-risk groups (Figure 4B). The results demonstrated that m6A-LncScore was completely independent of four factors (age, gender, lymph node count, and cancer type) and partially independent of another four factors (tumor stage II, AJCC T3, N0, and M0). Taken together, m6A-LncScore was an independent prognostic biomarker for PFS.

Previous studies revealed three lncRNA-related signatures relevant to PFS in CRC patients [10–12]. Thus, we compared the prognostic potential of m6A-Lnc signature (called m6A-LncSig) to these three lncRNA-related signatures (Zhao-LncSig, Huang-LncSig, and Gu-LncSig). In TCGA patients, the three signatures had a similar tendency for high-risk patients to have a shorter PFS period. They all showed a good ability of predicting PFS (p = 0.032, HR = 1.35; p = 2e-04, HR = 1.82; p = 0.0045, HR = 1.57, log-rank test) (Figures 5A–C). Dependent ROC analysis was conducted to compare the prognostic power of m6A-LncSig and the three signatures in the TCGA dataset. In general, the AUC at 1‐, 3‐, and 5‐ PFS for the m6A-LncSig was 0.75, 0.73, and 0.76, which was significantly higher than that of Zhao-LncSig (AUC = 0.55, 0.53, 0.53), Huang-LncSig (AUC = 0.6, 0.62, 0.63), and Gu-LncSig (AUC = 0.63, 0.62, 0.63) (Figure 5D, Supplementary Tables S16). Even integrated with tumor stage, the predictive power (1‐year, 3‐year, and 5‐year PFS) of m6A-LncSig was comparable with or significantly higher than the other three signatures (Figure 5E, Supplementary Tables S16). These results demonstrated that the prognostic power of the m6A-lncRNA signature was superior to three known lncRNA-related signatures.

Functional annotation was further conducted using gene-set enrichment analysis (GSEA) [34] for cancer hallmarks from MsigDB [35], which was implemented by R package “clusterProfiler” [36]. The differentially expressed genes between the high- and low-risk groups based on m6A-LncScore were enriched in immune-related cancer pathways and hallmarks (Figures 6A,B), such as interferon alpha/gamma response and inflammatory responses. The result indicated that the CRC patients with high m6A-LncScore had poor PFS time, which may be related to the immunosuppression of the tumor microenvironment.