A long non-coding RNA signature to improve prognosis prediction of colorectal cancer

1As the third leading culprit in cancer incidence worldwide, CRC continues to pose significant diagnostic, prognostic and therapeutic tribulations for clinicians. The American Joint Committee on Cancer (AJCC) TNM staging system is currently the only prognostic classification used in clinical practice to select patients for adjuvant chemotherapy. However, the AJCC stage fails to predict recurrence accurately in many patients undergoing curative surgery for localized CRC. This highlights the need for new biomarkers for a more precise prediction of high-risk patients with CRC recurrence and consequently improved personalized cancer care.

2Many studies have exploited microarray technology to investigate gene expression profiles (GEPs) in CRC in recent years, but only a small subset demonstrates clear prognostic significance. Molecular markers such as mutations in Kirsten ras gene (KRAS) and BRAF as well as chromosome instability (CIN) and microsatellite instability (MSI) have been systematically analysed for prognostic potential in CRC. So far, only KRAS mutation analysis has been used in clinical practice as a predictive marker for the effect of EGFR antibodies in metastatic disease. Currently, the roles of dysregulated functional long non-coding RNAs (lncRNAs) in human cancers have received considerable attention. LncRNAs are mRNA-like transcripts ranging in length from 200 nucleotides (nt) to~100 kilobases (kb) that lack significant protein-coding abilities. Increasing evidence suggests that these transcripts are frequently aberrantly expressed in cancers, and some of them have been implicated in diagnosis and prognostication. Searching a lncRNA signature might be of concrete predictive and prognostic value in the management of CRC.

3Presently, a large group of lncRNA-specific probes were fortuitously represented on the commonly used microarray platform (Affymetrix HG-U133 plus 2.0), so we initially mined previously published gene expression microarray data from the Gene Expression Omnibus (GEO), and conducted lncRNA profiling on large cohorts of CRC patients. By using the sample-splitting method and Cox regression analysis, we identified a prognostic, six-lncRNA signature from the GSE39582 test series patients, and validated it in the internal validation series and another two independent GEO cohorts. To further confirm the reliability of the new signature we conducted a series of experiments to investigate three of six lncRNAs' biological behavior in CRC cell lines.

4CRC data sets and corresponding clinical data were downloaded from the publicly available GEO databases. After removal of the samples without survival status, a total of 895 patients were analyzed, as detailed in Supplementary Table S1. These included 553 patients from GSE39582 (436 patients from the test series and 117 from the validation series), 197 patients from GSE14333 and 145 patients from GSE17536.

5The 553 GSE39582 CRC patients were assigned to a test series (N=436) or a validation series (N=117). The test data set was used for the detection of prognostic lncRNAs. By subjecting the lncRNA expression data of the test series to univariable Cox proportional hazards regression analysis using the BRB-Array Tools, we identified a set of six lncRNAs that were significantly correlated with patients' DFS (P<0.0001; Table 1). Of these, a positive coefficient indicated that a higher expression level of two genes (AK024680 and AK026784) was associated with shorter survival. The negative coefficients for the remaining four genes (AK123657, CR622106, BX649059 and BX648207) indicated that their higher levels of expression were associated with longer survival.

6We created a risk-score formula according to the expression of these six lncRNAs for DFS prediction (see Material and methods), as follows: Risk score=(0.07826*expression level of AK024680)+(-0.14300* expression level of AK123657)+(-0.19355* expression level of CR622106)+(-0.00172* expression level of BX649059)+(-0.20855* expression level of BX648207)+(0.24326* expression level of AK026784). We then calculated the six-lncRNA signature risk score for each patient in the test series, and ranked them according to their risk score. As such, patients were divided into a high-risk group (N=218) or a low-risk group (N=218) using the median risk score of the test series as the cutoff point. Patients in the high-risk group had significantly shorter median DFS than those in the low-risk group (log-rank test P<0.0001) (Figure 1A). In details, the survival of patients in the low-risk group was 83.71% at 50 months, 78.09% at 100 months, and 78.09% at 150 months; which compared with 58.76%, 54.35% and 51.49%, respectively, in the high-risk group. The association of the six-lncRNA risk score and DFS was also significant when it was evaluated as a continuous variable in the univariable Cox regression model (Table 2).

7To confirm our findings, we validated our six-lncRNA signature in the internal validation series. By using the same risk score formula, we classified patients into a high-risk (N=44) and a low-risk group (N=73) using the median score of the test series as the cutoff point as for the validation series. In consistence with the findings described above, patients in the high-risk group had significantly shorter median DFS than those in the low-risk group (log-rank test P=0.0212) (Figure 1B). Risk score-based classification of the entire GSE39582 cohort (i.e. combined test and validation series) also yielded similar results (Figure 1C). The DFS for patients with low-risk scores was 82.15% at 50 months, 78.40% at 100 months, and 78.40% at 150 months; which compared with 58.80%, 55.17%, 52.77% in patients with high-risk scores, respectively. In the univariable Cox regression model that the six-lncRNA risk score was evaluated as a continuous variable, similar correlation could be observed (Table 2). We also showed the distribution of lncRNA risk score, the survival status of the CRC patients and the lncRNA expression signature. As shown in Figure 2, in the entire GSE39582 series patients, we found that patients with high-risk scores tended to express high levels of risky lncRNAs (AK024680 and AK026784) in their tumors, whereas patients with low-risk scores tended to express high levels of protective lncRNAs (the remaining four).

8We further validated our six-lncRNA signature in two independent CRC data sets obtained from GEO: GSE14333 and GSE17536. The clinical characteristics of these two cohorts are listed in Supplementary Table 1. These two datasets confirmed the ability of our model in predicting survival. As shown in Fig. 1D, the six-lncRNA model could effectively predict the DFS in patients from GSE14333 (log-rank test P=0.0019). For the GSE17536 set, patients could similarly be segregated into a high-risk group and a low-risk group (DFS: log-rank test P=0.0025 Figure 1E; DSS: log-rank test P=0.0260 Figure 1F). In the univariable Cox regression model, the lncRNA risk score was again significantly associated with DFS as a continuous variable in both GSE14333 and GSE17536 cohorts (Table 2).

9We tested whether the prognostic value of the six-lncRNA signature was independent of AJCC stage. For this, we first performed multivariable Cox regression analysis that included lncRNA risk score, age, and other clinical characteristics such as gender and AJCC stage (when available) as covariables. The results showed that the six-lncRNA risk score remained to be significantly associated with DFS when adjusted by AJCC stage and other variables in every cohort (Table 2). Data stratification analysis was then performed, which stratified the two GEO (GSE39582 and GSE17536) patients (N=698) into a stage I and II stratum (without lymph nodes metastasis) or a stage III and IV stratum (with lymph nodes metastasis). This analysis showed that within each stage stratum, the six-gene risk score could further subdivide the patients into those likely to have longer survival and those likely to have shorter survival (Figure 3A-C), in a manner similar to what we observed in the entire group of patients, results of separate series see Supplementary Figure S1.

10Prognostic value of the signature for the patients with or without postoperative chemotherapy was also assessed. In the combined GSE39582 and GSE14333 series (when chemotherapy information available), high-risk score significantly correlated with an unfavourable DFS in those with and without postoperative chemotherapy (Figure 3D-F), results of separate series see Supplementary Figure S1& Figure S2.

11We performed GSEA to identify associated biological processes and signaling pathways using the lncRNA signature based risk score for classification. Significant gene sets (FDR < 0.01, p <0.005) were visualized as interaction networks with Cytoscape and Enrichment Map (Figure 4a, Supplementary Table S2). The risk score was accompanied with up-regulation of several cancer-related networks, namely Integrin pathway, Extracelluar matrix pathways, Focal adhesion. These related pathways were reported to affect cancer metastasis, and thus the signature might be involved with. Since cancer metastasis is a major influential factor for DFS, we compared the risk score of patients with and without distant metastasis in GSE39582 and GSE17536 series (when the information of AJCC stage available), for which the AJCC stage IV stands for distant metastasis, otherwise no distant metastasis occurs. Surprisingly, patients with distant metastasis tended to get higher risk score than patients without distant metastasis (Figure 4b, P<0.0001), results of separate series see Supplementary Figure S3&Figure S4.

12The analyses of CGH arrays revealed that CRC samples shared a typical DNA copy alteration pattern (Supplementary Figure S5). The signature lncRNAs AK024680, CR622106 and AK026784 showed positive correlations between their somatic copy-number alterations (SCNAs) and expression levels (Figure 5 A-C). While, for AK123657, BX648207 and BX649059, the DNA copy alterations information was unavailable, thus we selected AK123657, BX648207 and BX649059 among the six lncRNAs for further experimental validation and analyzed their expression in 25 pairs of CRC and corresponding nontumor colorectal tissues, using quantitative real-time polymerase chain reaction (qRT-PCR). These data confirmed that AK123657 (P<0.0001), BX648207 (P=0.0188) and BX649059 (P<0.0001) expression were decreased in CRC cancer (Figure 5D-F).

13To evaluate the effects of AK123657, BX648207 and BX649059 on cell biological behaviors, specific small interfering RNAs (siRNAs) were employed to knockdown AK123657, BX648207 and BX649059 expression in human CRC cell lines, HCT116 and SW1116. In case of target-off effect, two small interfering RNAs were used for each lncRNA (Figure 5G&Figure 5H). Cell-counting kit-8 assays indicated that cell proliferation was increased when AK123657, BX648207 and BX649059 was knocked down (Figure 5I&Figure 5J).These results suggested that AK123657, BX648207 and BX649059 played a physiological role in regulating cell proliferation. We further used Transwell assay to monitor the effect of manipulating AK123657, BX648207 and BX649059 expression on cell invasiveness. Knockdown of these lncRNAs significantly increased the number of CRC cells that penetrated the Transwell filter, which demonstrated a substantial gain of cell invasion ability (Figure 5K).

14The conventional view of gene regulation in biology has centered around protein-coding genes until the discovery of thousands of lncRNAs. Numerous reports of dysregulated lncRNA expression across numerous cancer types suggest that abnormal lncRNA expression may be a major contributor to tumorigenesis. The aberrant expressions of specific lncRNAs in cancer could mark the spectrum of disease progression and these lnRNAs may serve as independent biomarkers for diagnosis and prognosis. More recently, lncRNAs have been implicated in CRC pathogenesis. However, the prognostic values of lncRNAs in CRC have not yet been investigated. To explore the prognostic lncRNA genes, we profiled lncRNA by mining the existing microarray gene expression data since a great number of lncRNAs were interrogated on many commonly used commercial arrays. By analyzing the associations between lncRNA expression profiles and clinical outcome of CRC patients in GEO date set, we identified a six-lncRNA signature that was significantly associated with the DFS. Base on our knowledge, this is the first report that relates lncRNA expression patterns with patients' prognosis in CRC.

15By applying the six-lncRNA signature to the GSE39582 test series patients, a clear separation was observed in the survival curves between patients with low- or high-risk signatures. Patients with a low-risk six-lncRNA signature in their tumor specimens tended to have prolonged overall survival, whereas patients with a high-risk signature tended to have shortened survival. The association between the lncRNA signature and survival was significant regardless of the former was evaluated as a continuous variable or category variable (dividing by the median cutoff). The usefulness of this lncRNA signature could be internally validated in the non-overlapping GSE39582 patients (the validation series) and two independent cohorts of GSE14333 and GSE17536 that were profiled through the same platform of GSE39582 patients, indicating good reproducibility of this lncRNA signature in CRC.

16Further analysis revealed that the prognostic value of the six-lncRNA signature was independent of one of the main prognostic factors-AJCC stage. Currently, AJCC stage has been widely accepted as a powerful predictor of treatment response and survival in CRC. In general, stage I and II shows better survival, whereas stage III and IV is correlated with worse survival in a treatment-independent manner. In accordance with other studies, AJCC stage was a significant prognostic factor in our study when assessed in the univariable Cox regression analysis (Table 2). Therefore, it will be important to evaluate whether the prognostic value observed on our six-lncRNA signature is independent of this known strong prognostic factor or not. Here, by performing multivariable Cox regression analysis and stage stratification analysis, we showed the stage-independent prognostic values of the six-lncRNA signature in CRC patients. The exception was GSE14333, in which we did not have the AJCC stage information. Nonetheless, combining all these results together, we could conclude that the prognostic value of the lncRNA signature observed in our study was independent of AJCC stage. Moreover, it was of interest to find that the six-lncRNA signature had a similar survival predictive ability as AJCC stage in the receiver operating characteristic (ROC) analysis. Six-lncRNA signature combine AJCC stage has a stronger power for DFS prediction in the ROC analysis (Supplementary Figure S6). The ability of our six-lncRNA signature in identifying subgroups of CRC patients with identical AJCC stage implies that lncRNA signatures may be used to refining the current prognostic model and facilitating further stratification of patients in the future clinical trials.

17Our signature predicted an unfavourable prognosis of those with or without postoperative chemotherapy in combined cohorts, indicating that some signature lncRNAs might involve in chemotherapy response and contribute to therapeutic outcome, thus needed to be studied further.

18Examination of associated molecular pathways revealed that the six-lncRNA signature was more likely to involve with Integrin pathway, extracelluar matrix pathways and focal adhesion. Extracellular matrix via its receptors, the integrins, and kinases downstream of β1 integri-focal adhesion kinase (FAK) have emerged as a major pathway leading to several cellular responses including migration, differentiation, and proliferation. Furthermore, patients with distant metastasis were validated in our study to achieve higher risk score than patients without distant metastasis. Studied have revealed that distant metastasis played a significant role in DFS, cases without distant metastasis were more likely to get prolonged DFS and vice versa. Therefore, the significant signaling pathways might support that lncRNA signature has DFS prediction power and suggest possible avenues for future targeted therapies.

19However, gene signatures attempted to serve as prognostic and/or predictive markers for CRC hardly achieve widespread use despite accumulating signatures have been explored. This may partly result from ignoring of biologic rationale during gene selection, which may derive from computer-based algorithms that incorporate a cluster of unrelated genes. SCNAs is an important form of somatic genetic alteration in cancer, and within that a genomic region is either amplified or deleted. Some of the genes within amplified (or deleted) regions exhibit increased (or decreased) expression levels resulting in altered activity in cancer cells. As lncRNAs do not have proteins encoding abilities, their functions are closely related with their transcript abundance. The signature lncRNAs AK024680, CR622106 and AK026784 showed positive correlations between their SCNAs and expression levels, which we reasoned could possibly lead to altered lncRNA activity in CRC. The other three lncRNAs from the gene signature were functional studied by further experiments. We focused on the lncRNAs' biological function and analyzed their relation with clinicopathological features. AK123657, BX649059 and BX648207 were significantly down-regulated in CRC tissues compare to normal colorectal tissues, suggesting a protective role in CRC biogenesis. Cell proliferation and invasion ability were enhanced after decreasing expression of these lncRNAs, which was in accordance with the gene signature that patients suffered distant metastasis achieve higher risk score than patients without distant metastasis.

20The limitations should be acknowledged for this study. In particular, treatment response was not included because this information was not available for a substantial proportion of cases. Thus, the significance and robustness of the signature as a prognostic classification requires further confirmation, ideally with large prospective patient cohorts included in adjuvant trials.

21The present study demonstrated the associations between the expressions of these lncRNAs and DFS of CRC patients. The findings may identify high-risk patients for more intensive adjuvant therapy in addition to the standard regimen; low-risk patients, on the other hand, may not need to receive intensive, and potentially toxic, therapies.

22Biological behavior of three lncRNAs from the gene signature was studied because we supposed gene signature consisted of functional lncRNAs tended to have better prediction. However, the roles of the other lncRNAs in CRC pathogenesis are presently unclear, and our findings suggest that they deserve further studied. Additional functional investigations of these lncRNAs on cancer cell lines and xenograft models may increase our outstanding of their roles in determining CRC prognosis. Last but not least, all the lncRNAs were derived through the re-annotation algorithm in this study and should be validated in the further studies.

23In conclusion, this study presents a powerful lncRNA signature by probing and integrating currently available microarray data. This innovative lncRNA signature showed independence of one of the main prognostic factors-AJCC stage and the signature might help personalize prediction of CRC prognosis. The GSEA analysis suggested that the signature might involve with several cancer pathways, most likely with cancer metastasis, which provided support for DFS predictive ability of the signature. These identified signature lncRNAs played vital roles in the CRC cell lines, which might serve as alternative biomarkers and therapeutic targets for CRC.

24The raw CEL files were downloaded from GEO database and background were adjusted using Robust Multichip Average. GATExplorer was used to process microarrays on a local computer for gene expressions of lncRNAs. This GATExplorer provides a series of R packages, designed to be used with BioConductor tools, that allow to apply in a simple way the probe mapping data included in GATExplorer. A type of files called ncRNA Mapper were also obtained from GATExplorer, which include the probes that do not map to any coding region but that were mapped to a database for non-coding RNA of human and mouse derived from RNAdb. A customized R scripts was used to perform a microarray expression calculation according to the re-mapping data(file ncrnamapperhgu133plus2cdf_3.0). Each LncRNA should include at least a minimum of 3 probes mapping in the corresponding ncRNAs entity.

25GSEA was performed by the JAVA program (http://www.broadinstitute.org/gsea) using MSigDB C2 CP: Canonical pathways gene set collection(1320 gene sets available).The GSEA, visualized in Cytoscape (version 2.8.2), and the Enrichment Map software, was used to determine if the members of a given gene set were generally associated with risk score, and was therefore performed on all mRNA genes on the HG-U133 Plus 2.0 ranked by enrichment score from most negative to most positive. 1000 random sample permutations were carried out, and the significance threshold set at FDR<0.01. If a gene set had a positive enrichment score, the majority of its members had higher expression accompanied with higher risk score, and the set was termed “enriched”.

26The association between the lncRNA expression and patient's disease free survival (DFS) and disease specific survival (DSS) was assessed by univariable Cox regression analysis along with a permutation test using Biometric Research Branch-Array (BRB-Array) Tools. Genes were considered statistically significant if their permutation P values were less than or equal to 0.01. To construct a predictive model, the selected genes were fitted in a multivariable Cox regression model in the test series as described. A risk score formula was then established by including each of these selected genes, weighted by their estimated regression coefficients in the multivariable Cox regression analysis. With this risk score formula, patients in each set were classified into high-risk or low-risk group by using the median risk score of the test series as the cutoff point. Survival differences between the low-risk and high-risk groups in each set were assessed by the Kaplan–Meier estimate, and compared using the log-rank test. To test whether the risk score was independent of AJCC stage, multivariable Cox regression analysis and data stratification analysis were performed. We used ROC curves to compare the sensitivity and specificity of the survival prediction based on the lncRNA risk score and AJCC stage. In the log-rank test, Cox regression analysis and ROC analysis, the significance was defined as P values being less than 0.05.

27A total of 402 of the 553 primary CRC samples from the entire GSE39582 cohort could be analyzed for array-based comparative genomic hybridization (CGH). To exam whether lncRNAs within amplified (or deleted) regions have increased (or decreased) expression levels, Mann-Whitney U (MWU) test was used to determine the significance of the comparisons.

28Cell viability was assessed by the Cell Counting Kit 8 (CCK-8; Dojindo) as described previously. Briefly, control and treated HCT116 and SW1116 cell lines were seeded into 96-well plates at an initial density of 5000cells/well. At each time points, 10 μl of CCK-8 solution was added to each well and incubated for 2 h. The absorbance was measured by scanning with a microplate reader at 450 nm. Data were expressed as the as follows: relative viability= A450 (treated) − A450 (blank) or (A450 (control) − A450 (blank).

29AK123657 siRNA1 CCUCCAGACUGUGAGUAAUTT AUUACUCACAGUCUGGAGGTT

30siRNA2 GGAGGUGCAUGACUAACAATT UUGUUAGUCAUGCACCUCCTT

31BX648207 siRNA1 GGCCUGAAUUUGGUUACAUTT AUGUAACCAAAUUCAGGCCTT

32siRNA2 GCCACUUGCAAGUGGAAUATT UAUUCCACUUGCAAGUGGCTT

33BX649059 siRNA1 GGCUUAAAGUAGGUAUUUATT UAAAUACCUACUUUAAGCCTT

34siRNA2 CCUGAAGAGUACAGAUAAATT UUUAUCUGUACUCUUCAGGTT

35Tumor cell invasion assays were performed using Boyden chambers with filter inserts (pore size, 8-μm, Millipore) coated with Matrigel (40μm; BD Biosciences) in 24-well dishes (Corning) as described previously. Briefly, 3×105 cells after transfected with lncRNA siRNAs or control siRNA were seeded in the upper chamber, while the RPMI-1640 medium (Invitrogen) supplemented with 30% fetal bovine serum was placed in the lower chamber. The plates were incubated for 48h. Then the cells were fixed in 4% formaldehyde and stained with 0.05% crystal violet for 20min at room temperature. Cells on the upper side of the filters were removed by cotton-tipped swabs, and the filters were washed with PBS. The cells on the lower side of the filters were defined as invasive cells.