EPDR1 up-regulation in human colorectal cancer is related to staging and favours cell proliferation and invasiveness

The finding of novel molecular markers for prediction or prognosis of invasiveness in colorectal cancer (CRC) constitutes an appealing challenge. Here we show the up-regulation of EPDR1 in a prospective cohort of 101 CRC patients, in a cDNA array of 43 patients and in in silico analyses. EPDR1 encodes a protein related to ependymins, a family of glycoproteins involved in intercellular contacts. A thorough statistical model allowed us to conclude that the gene is significantly up-regulated in tumour tissues when compared with normal mucosa. These results agree with those obtained by the analysis of three publicly available databases. EPDR1 up-regulation correlates with the TNM staging parameters, especially T and M. Studies with CRC cell lines revealed that the methylation of a CpG island controls EPDR1 expression. siRNA knocking-down and overexpression of the gene following transient plasmid transfection, showed that EPDR1 favours cell proliferation, migration, invasiveness and adhesion to type I collagen fibres, suggesting a role in epithelial to mesenchymal transition. Both statistical and functional analysis correlated EPDR1 overexpression with invasiveness and dissemination of tumour cells, supporting the inclusion of EPDR1 in panels of genes used to improve molecular subtyping of CRC. Eventually, EPDR1 may be an actionable target.

EPDR1 is up-regulated in human cRc patients. EPDR1 expression was first analysed in a TissueScan cDNA array (OriGene) of human CRC patients (see Material and Methods). The level of whole EPDR1 and of its isoform 2 was determined by RT-qPCR in all the samples. Expression of total EPDR1 is detected in normal mucosa at a very low level, but is clearly expressed in almost all the tumour samples. Only in 5 out of the 43 tumour samples the fold change of EPDR1 expression relative to the normal mucosa is less than 5. As shown in Fig. 1a, the up-regulation of the whole gene is found in all the CRC stages. The analysis of the gene expression using the usual ANOVA approach and multiple comparisons, revealed a significant difference (p = 0.041) between normal and tumour tissues. The behaviour of isoform 2 was somewhat different. First of all, as previously reported 27 , its expression level is low, and it is undetectable in normal tissues. Moreover, independently of the stage, the plots showed that in some CRC patients isoform 2 is not detected or expressed at a negligible level, while there are other patients showing a higher expression (Fig. 1b). Stage II patients show a large dispersion of expression values, both in whole EPDR1 and in its isoform 2.
Then, a prospective cohort of 101 CRC patients who underwent surgery at the University Hospital (Valencia) in the last 3 years (Supplementary Table S1) was studied. The use of paraffin-embedded samples did not allow us determining the level of the minor isoform 2 and only the expression of the whole gene is reported here.
We first determined the means of EPDR1 mRNA expression, measured by RT-qPCR, in 97 paired tumour and normal samples from the 101 patients' cohort. The remaining 4 samples were not apt for processing. A paired t-test showed that the average expression in tumours is significantly higher (p = 0.005) than in normal tissues (Fig. 2a). This may be also observed in the density plot built with the experimental results obtained from the 97 paired samples (Fig. 2b). Most of the expression values of EPDR1 in normal adjacent tissues group as a single skewed peak, with some outliers at high values. However, the expression of the gene in tumour samples results in a wider distribution. The arrows mark the positions of the means, clearly larger than the medians given in the box plot of Fig. 2a. This is related to the skewness to the right observed in the density plot. The plots for single paired samples obviously show that, in most cases, the expression of the gene is higher in tumour tissues than in adjacent mucosa in both non-metastatic and metastatic patients (Fig. 2c,d). Moreover, the up-regulation of the gene is observed irrespective of the AJCC stage of the patients. In the case of stage I patients the difference is close to significant, while the significance is clear for the remaining stages ( Table 1). As a further control, we determined the expression of the gene in the only four benign colon polyps available from patients. The resulting value (0.013 Figure 2. Expression of total EPDR1 in a cohort of 101 CRC patients from our University Hospital. (a) Box plots with whiskers with maximum 1.5 IQR of total EPDR1 expression measured by qPCR in tumour and nontumour adjacent tissues. Outliers with expression values higher than 200 (4 from tumours and 3 from non tumour samples) are not included in the graph, but they are shown in panel d. Statistical analysis was carried out with Student's t-test. (b) Density plot of the expression values for tumour (red) and non-tumour adjacent tissues (blue). The arrows point to the position of the means. The inset shows a magnification of the plot to appreciate the differences between the two tissue types for expression values lower than 150. (c) Up-regulation of EPDR1 in localized tumour samples when compared with paired normal adjacent tissue. (d) As c, but with samples from metastatic patients. Statistical analysis was carried out with tailed Student's t-test.
± 0.006 relative to β-actin gene) was not included in Table 1 because of the small sample size, but it falls in the range of non-tumour samples given in that Table. Next, our statistical analysis focused on the study of the differential expression of EPDR1 with respect to important predictors, using the data derived from our local patients' cohort. It is usual to include marginal descriptive studies of the expression of a gene in tumour and adjacent tissues with respect to these predictors and these are included in the Supplementary Material. Nevertheless, a thorough statistical analysis is described here. Both the expression differences between tumour and adjacent tissues and the expression ratio tumour/adjacent tissue were calculated for all the 97 paired samples and their dependence with respect to the patients' characteristics and other predictors (T, N, M and AJCC stages) was determined by using multiple regression. A stepwise variable selection has been performed as described under Materials and Methods.
When the analysis included the TNM and the AJCC stages, a significant coefficient in the ratio tumour/adjacent tissue expression was observed with T3 and T4 stages (p = 0.0021 and p = 0.0009, respectively), with M1 stage (p = 0.0004) and with AJCC stage II (p = 0.0005). No significant coefficients were observed between the ratios and N stages or with AJCC stage III. These results are summarized in Table 2. The p-value of the test for null coefficients corresponding to the variables T2, T3 and T4 is 0.00038. When the categorical predictors corresponding to N and M stages were removed from the model, a significant correlation with AJCC stages III and IV (p = 0.0013 and p = 4.8·10 −7 ) was also observed (see Supplementary Material). All these results manifest that a strong relation exist between EPDR1 expression in CRC tissues and the invasiveness of tumour cells.
When the analysis was done with the differences of expression between tumour and adjacent mucosa as response rather than with the ratios tumour/adjacent tissue, roughly similar results were obtained, but significant coefficients were also observed with the N stages ( Table 2). The reasons for using this double system of comparing the expression in tumour and adjacent tissues will be commented under the Discussion.
Next we asked if there exist a dependence of the patients' survival with respect to the level of gene expression. To answer this question, we have used a Cox regression analysis (see Supplementary Material, Statistical Analysis). The dependence of the survival time on EPDR1 expression in tumours is close to significance (p = 0.085), but a very significant dependence was observed with respect to the expression of the gene in the histologically normal www.nature.com/scientificreports www.nature.com/scientificreports/ tumour borders (p = 0.0057). When a Cox regression including the AJCC stage was carried out, the survival time was also found significantly dependent on the EPDR1 expression in tissues adjacent to the tumour (p = 0.0344) but not on the expression in the tumour itself (p = 0.1161). Similar results were obtained if the Cox regression included the T stage. In this case significant dependence on the expression in non-tumour adjacent tissue was found (p = 0.0127) but no significant dependence on the expression in tumour was found (p = 0.2188). In other words, the expression of EPDR1 in the environment of the tumour tissues seems to be a better predictor of survival time than the expression in the tumour itself.
The results obtained with the local cohort of CRC patients were compared with those obtained in silico by using TCGA data. This database includes expression values for EPDR1 in a cohort of 622 CRC patients. When a Welch two sample t-test was used to compare the expression means in both groups of samples, tumour and normal tissue, it was found that the expression in tumour is 2.08 fold higher than in normal tissues (p = 9.67·10 −9 ). The comparison of the paired samples included in the TCGA cohort (n = 50) also allowed us to conclude that EPDR1 expression is significantly higher in tumour than in paired adjacent mucosa (Supplementary Table S2). Finally, Kim et al. 28 generated a series of RNA-seq data from 18 CRC metastatic patients in which they included samples from normal mucosa, primary tumours and metastasis. The database is publicly available (https://www. ncbi.nlm.nih.gov/bioproject/218851) and we used it to recover the EPDR1 expression values. The gene expression in primary tumours is 2.15-fold increased over that in normal mucosa and, interestingly, the fold change increase is 2.518 when the expression in metastases is compared with that of normal mucosa (Supplementary Table S3). All the results acquired in silico both from TCGA and from the Kim et al. 28 data agree with those obtained from our local patients' cohort. A complete description of the model used for statistical analysis of the data sets is given in Supplementary Material.
In summary, the analysis of expression data in a large number of patients, both in vivo and in silico, allowed us to conclude that EPDR1 is significantly up-regulated in tumour tissues when compared to adjacent mucosa. Moreover, the statistical analysis indicated that expression in tumour correlated with invasiveness and dissemination of tumour cells. Therefore, we next addressed the questions of the possible causes leading to a differential expression of the gene and of the mechanisms that might relate EPDR1 up-regulation to CRC progression. To answer some of these questions we used in vitro models with human CRC cell lines, as described below.

EPDR1 expression is regulated by DnA methylation in cRc cell lines.
We have previously found that the expression of the whole EPDR1 gene is not detected in RKO, almost negligible in SW48 (0.00012 ± 0.00007 relative to ACTB), low in HCT116 (0.0054 ± 0.0004), and high in DLD1 (0.0252 ± 0.0056) and D-Mut1 (0.0454 ± 0.0055) 27 . We hypothesize that epigenetic changes may cause these differences in EPDR1 expression among the cell lines used. To check that hypothesis, the methylation of a segment of 1241 bp, which includes a CpG island encompassing EPDR1 exon 3, the first transcribed exon except for isoform 3, and the first translatable exon in the major isoform, was determined. A schematic map of the analysed region, including the location of amplicons used, is depicted in Fig. 3a. Interestingly, the methylation level of the analysed CpGs is compatible with the idea that methylation of the CpG island is the cause of EPDR1 repression in SW48 and RKO cells and of the low expression in HCT116 (Fig. 3b). The methylation of the distal CpGs upstream of the island does not vary from line to line and this serves as an internal control of the results. In a comparison of the methylomes of normal colon mucosa, tubular adenomas, and colorectal cancers. Luo et al. described the methylation status of 5 of the CpGs we have analysed 29 and they are identified at the bottom of Fig. 3b. Interestingly, these authors showed that the CpG at location 37,920,084 is hypermethylated irrespective of the nature of the sample, and we found that this CpG is hypermethylated regardless of whether the gene is expressed or not in the cell lines tested. Therefore, the results presented here agree with the constitutive hypermethylation found in that CpG in human colonic tissues.
To further check that the CpG island methylation is related to EPDR1 transcription, the expression of the gene was determined in RKO and SW48 cells grown in the presence of 5-azacitidine. A significant, dose-dependent increase in transcription of the gene was observed in both cell lines (Fig. 3c), showing that demethylation of DNA actually causes the activation of EPDR1 expression. Anyway, the fact that the island is virtually non-methylated in both DLD1 and D-Mut1cells, which show high, but distinct, levels of EPDR1 transcription, suggests that some additional mechanism must exist to account for that difference.

EPDR1 increases cell proliferation in cRc cell lines.
To explore the function of EPDR1 in human CRC progression, its effects on cell growth were first analysed. To do this, the gene was independently silenced with two EPDR1 siRNAs, further referred to as si1 and si2, as mentioned under Material and Methods. Both siRNAs efficiently knocked down the expression of EPDR1 both at the level of mRNA and protein ( Supplementary  Fig. S2). The consequences of this knockdown were evaluated by MTT, colony formation assays and flow cytometry ( Fig. 4 and supplementary Fig S3). The assays show that EPDR1 significantly favours proliferation of DLD1 and HCT116 cells. Silencing the gene with siRNAs affects the progression of the cell cycle in both DLD1 and HCT116 cells (Fig. 4c).
When the results obtained with si1 ( Fig. 4) and si2 (supplemental Fig. S3) are compared some differences were observed, which may suggest that the cytotoxicity of both siRNAs is different in each cell line. For instance, the sub-G0 population at 48 h is higher in si1-treated DLD1 than in HCT116 cells, while an opposite result was obtained with si2. A reliable characterization of the sub-G0 population cannot be achieved by a simple cell cycle analysis. To overcome this difficulty, an analysis in the presence of 7-aminoactinomycin D and Annexin V was carried out to discriminate between apoptotic and necrotized cells. Supplementary Fig. S4 shows the results of an experiment with DLD1 in which untreated cells and cells transfected with scrambled siRNA or with si1 or si2 were analysed after 72 h of transfection. It is clear that the differences in cytotoxicity found between si1 and si2-silenced cells are mainly attributable to a greater necrosis induced by si1 in this cell line. The use of untreated DLD1 cells in the experiment of supplemental Fig. S4 shows that neither si2 nor control siRNA affect the viability (2020) 10:3723 | https://doi.org/10.1038/s41598-020-60476-7 www.nature.com/scientificreports www.nature.com/scientificreports/ of the cells in a significant way. Moreover, both si2 and scrambled siRNA induce apoptosis and necrosis to an extent similar to that of untreated cells. To further explore the differences in cytotoxicity between si1 and si2, we transfected RKO cells, in which EPDR1 is not expressed 27 , with si1 and si2. No consequences in proliferation or cell migration were observed ( Supplementary Fig. S5) and the sub-G0 population was similar between cells transfected with si1 and si2 (12.39 ± 1.66% and 13.29 ± 1.96% respectively). These results allowed us to discard the possibility that the phenotypic effects of si1 and si2 were due to off-target effects. EPDR1 promotes cell migration, invasiveness and adhesion to collagen type I fibres. The effects of silencing EPDR1 on cell migration were next studied in DLD1 and HCT116 cell lines. The results of the transwell analysis (Fig. 5a) showed that EPDR1 knockdown significantly reduces cell migration, and the wound healing assay (Fig. 5b) confirmed these results. When the transwell analysis was carried out through a Matrigel layer, the results indicated that EPDR1 is also involved in migration through extracellular matrix (Fig. 5c). It can be, therefore, concluded that EPDR1 promotes cell migration and invasiveness.
As migration of cells is facilitated by binding of cells to collagen type I fibres 30,31 , the effects of silencing EPDR1 on cell adhesion to collagen-coated plates was next examined. The results show that knocking-down the gene significantly reduces the adhesion of DLD1 and HCT116 cells to collagen type I fibres (Fig. 5d). When si2 was used to knock-down EPDR1 similar results in cell migration, invasiveness and adhesion to collagen assays were obtained (Supplementary Figure S6).
Once studied the effects of knocking-down EPDR1, it can be wondered whether the opposite effects are observed by overexpressing the gene. To answer this question, we used the RKO and SW48 cells, which, as mentioned above, are characterised by a null or negligible expression of EPDR1. Supplementary Figure S7 shows that  www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ transfecting the cells with pEPDR1 resulted in a dramatic increase of gene expression. In both, RKO and SW48, overexpression of EPDR1 resulted in a significant increase in colony formation, migration, invasiveness and adhesion to type I collagen of the cells (Fig. 6). These results confirm the proposed role for EPDR1 in increasing the dissemination of cancer cells. 32 and proposals for CRC have been done 8 . It has been recently emphasized that transcriptional signatures may allow a refinement of the CRC subtyping 9 and, in this sense, finding novel genes dysregulated in the disease, either in tumour itself or in its environment, may aid to construct gene panels useful for diagnosis, prognosis or targeted therapy selection.

Molecular subtyping based on gene expression is gaining interest in cancer stratification
In the present paper, we show that EPDR1 gene is up-regulated in CRC. The first report in this sense was published long time ago but it was limited to three patients 14 . When all the experimental data described in the present paper had been obtained, Chu et al. reported that EPDR1 is down-regulated in CRC tumour cells when compared with adjacent mucosa. They analysed a local cohort of 75 patients, but the mRNA expression of EPDR1 mRNA was only determined in a sub-group of 23 paired samples 33 . The cohort analysed in the present work is large enough and the thorough statistical analysis carried out allows us to draw more significant conclusions. Moreover, the results presented here include as well an additional analysis of 43 patient samples from a cDNA array, in which we have also shown the overexpression of isoform 2, and all of our data agree with the huge number of values recovered in silico from publicly available databases. Nevertheless, the above comments do not imply that every single tumour express higher levels of EPDR1 than normal mucosa cells. Figure 2a shows that the expression level of the gene in some individual tumours is as low as in normal adjacent tissues or in benign colon polyps. A similar comment can be done regarding the EPDR1 expression in cell lines. In a previous paper we reported the analysis of EPDR1 in five CRC cell lines and some of their isogenic derivatives, to show that the expression of the gene is negligible in some of them 27 , as it occurs with the colorectal normal NCM460 cells 14 .
Apart from the bioinformatic analysis of TCGA data reported above, the expression of EPDR1 in other tumours was recovered from that database by using the facility Gene Expression Profile Interactive Analysis (http://gepia.cancer-pku.cn/). In some other cancers, such as hepatocellular carcinoma, kidney carcinoma, cholangiocarcinoma, pancreatic adenocarcinoma and stomach adenocarcinoma, the gene expression profile is also elevated when compared with paired normal tissues. Nevertheless, uterine carcinosarcoma, bladder urothelial carcinoma and lung squamous cell carcinoma, among others, behave in an opposite way.
The results of gene expression in tumour and adjacent normal tissues are usually given as the fold change, i.e., the ratio of expression in tumour and adjacent tissue. Though useful for comparative purposes, some flaws may exist in this procedure, especially when the expression in normal tissues is very low or even below the detection limit. In these cases, the absolute error of the determination may result in large errors in the ratio, abnormally high values in it, or even a senseless infinite value. To prevent these errors, in the present paper, apart from the customary fold change values, we have also statistically analysed the results subtracting the expression in normal tissues from that in tumour samples. Both procedures, fully presented in the Supplementary Material, roughly gave similar results, except in some instances, mentioned under Results. Of note, high-level expression of EPDR1 is associated with T and M parameters and, less significantly, with N. This is clearly in accordance with the results obtained after knocking-down or overexpressing the gene in CRC cell lines. As mentioned above, the in silico analysis of the data obtained by Kim et al. 28 revealed that EPDR1 expression in metastatic tumours is higher than that in the primary CRC tissues. These results are in accordance with the high level expression of the gene in metastatic tumours described here and with the proposed role for EPDR1 in EMT.
In strong disagreement with the data of Chu et al. 33 our results on the effects of silencing EPDR1 clearly show that the gene favours not only cell proliferation, but also invasiveness ( Fig. 4a and Supplementary Fig. 3a). These facts may be the cause of the involvement of EPDR1 with tumour development in vivo. Moreover, our results strongly suggest that the ectodomain of EPDR1 possesses antiadhesive properties as the related ependymins do 12 . The dissemination of cancer cells may be facilitated by the effects of EPDR1 on adhesion to type I collagen fibres, which have been described as the "highways" for tumour cell migration 31 . The interaction of fish ependymins with collagen is calcium-dependent 13 ; should this dependence also take place in EPDR1, an interesting consequence might result. STIM1 has been found to be overexpressed in CRC and promotes cell motility and EMT 34,35 . The gene encodes an endoplasmic reticulum sensor, which regulates the concentration of Ca 2+ , in cancer cells 36 . Therefore, it is tempting to speculate that EPDR1 might be a downstream target of STIM1, which would ultimately cause, in a Ca 2+ -dependent manner, the cell migration and EMT by affecting cell adhesion and/or interaction with type I collagen fibres. Experiments to check this hypothesis are currently being done in our laboratory.
We have also shown that EPDR1 is regulated by DNA methylation at an island encompassing the first translatable exon in the major isoform, although some other regulatory mechanism must exist. Methylation-dependent regulation of EPDR1, using a region shorter than ours, but overlapping with it, has also been reported in the recent paper of Chu et al. 33 . Luo et al., in the previously mentioned methylome analysis, found that in most stage IV patients the methylation of CpGs located at positions 37,920,649 37,920,715 and 37,921,129 is very low, although there is a subgroup of patients in which these CpGs are hypermethylated 29 . Present results show that the highest expression values of EPDR1 are found in some stage IV patients (Fig. 2d) and, therefore, our findings are in agreement with an important, though not unique, role of DNA methylation in the regulation of EPDR1 expression.
It is worth noting that patients' survival depends more on the expression of EPDR1 in non-tumour adjacent mucosa than in tumour itself. The adjacent non-tumour tissues are defined as normal from a histological point of view, but they may be abnormal when their molecular characteristics are taken into account. For instance, it has been shown that histologically normal tissues surrounding breast cancer tissues exhibit some gene expression properties reflecting those of the cancer subtypes, probably causing in some instances the recurrence after conservative therapies. Consequently, the importance of studying biomarkers in adjacent tissues to predict the risk of recurrence has been emphasized 37 . Recently, Aran et al. 38 have conducted an exhaustive transcriptomic study in which compared the tumour, non-tumour adjacent tissues and healthy tissues obtained by post-mortem collection. In some instances, specific gene activation occurs in adjacent tissues. In our case, as EPDR1 overexpression is compatible with a diminution of cell adhesion and in higher cell association with type I collagen fibres, increased expression in adjacent tissues may facilitate dissemination and the risk of metastasis.
At any rate, present results may be useful to consider the inclusion of EPDR1 in panels of genes used to improve molecular subtyping of CRC and perhaps to consider the gene as an actionable therapeutic target.  Fig. S2).

Methods
For transient overexpression of EPDR1, a plasmid containing the full-length cDNA of the canonical isoform 1, under the control of immediate-early cytomegalovirus promoter (Sino Biological, HG13665-ACG), was transfected into RKO and SW48 cells. An empty plasmid was used as control. For transfection, 3·10 5 cells were cultured in 6 wells plate, and grown to 80% confluence. Then, a mixture of 4 μg of plasmid and 10 μl of Lipofectamine 2000 Reagent (Invitrogen, #11668-027), solved in Opti-MEM Medium (Thermo Fisher Scientific, #31985047) were added independently to each three wells. After 6 hours, this medium was removed and standard medium was added. The effects of EPDR1 overexpresion were analysed after 72 h ( Supplementary Fig. S7).
All the experiments described under this heading were repeated twice.
cell proliferation assays. Cell proliferation was measured by the 3-(4,5-dimethylthiazol-2-yl)-2,5-diph enyltetrazolium bromide (MTT) assay. After three independent transfections, either with scrambled siRNA or EPDR1 siRNAs, cells were trypsinized and re-plated in quadruplicate (3, Analysis of DnA methylation. Quantitative DNA methylation analysis was carried out using a MassARRAY platform (Agena Bioscience) as previously described 39 . The primers used, their localisation and the size of defined amplicons are given in the Supplementary Table S4.
Rt-qpcR analysis of the expression of EPDR1. The transcription level of the whole EPDR1 gene was determined by RT-qPCR with the primers given in the Supplementary Table S5, which define an amplicon at the 5′ end of exon 6, common to all the isoforms. The level of isoform 2 was determined with primers from exon 4, unique to this isoform 27 (Supplementary Table S5). The ACTB gene was used as loading control.
Western blot analyses. Western blot analyses were carried out essentially as previously described 40 . Briefly, cultured cells were extracted with RIPA buffer (50 mM Tris-HCl pH 7.5, 150 mM NaCl, 1% Triton X-100, 0.1% SDS, 0.5% deoxycholic acid sodium salt (w/v)) supplemented with 2 μl/ml protease inhibitor cocktail (Sigma). Samples were sonicated and centrifuged at 13,000 × g for 30 min at 4 °C. Total protein was determined by a BCA protein assay kit (Thermo Fisher Scientific, 23225). Electrophoresis was carried out in 12% polyacrylamide gels and the transfer was done overnight. Immunoblots were visualized using either the ECL Western Blotting detection kit reagent (GE Healthcare) or SuperSignal West Femto (Thermo, #34095) and the ImageQuant LAAS 400 system (Healthcare Bio-Sciences). The antibodies used were: anti-EPDR1 (Abcam, ab-197932); anti-β-actin (Abcam, ab-8227). For a semi-quantitative determination of protein in the western blots, four grey values relative to the loading control were measured with ImageJ and averaged.
Statistical analyses. The differential expression between tumour tissue and control was quantified using the difference of the observed expressions or their fold-change between tumour and adjacent tissue. The expression of the gene for both types of tissues was compared using a paired t-test. Multiple regression was used to study the relationship between the difference of the expressions or the fold-change and the predictors (StT, StN, StM and StStage). A stepwise variable selection was performed using the function stepAIC 41 . Survival functions, estimated using the Kaplan-Meier estimators, were plotted using survminer 42 . The relationship between the survival time and the covariables (expressions in tumour and normal tissue and the covariables StT, StN, StM and StStage) was studied using a Cox regression performed using the R package survival 43 . The survival functions were compared using the log-rank test when a single categorical covariable had been used.
The in silico analysis of The Cancer Genome Atlas (TCGA) data included the EPDR1 differential expression analysis using edgeR bioconductor package and a survival analysis using the R package coin 44 . The expressions measured in the OriGene arrays for different stages were compared using a unilateral Kolmogorov-Smirnov test. A detailed description can be found in the Supplementary Material.
All the statistical analyses were performed using the software R 45 . The plots have been made using ggplot2 46 . Other R packages used in the analysis are rcompanion 47 and xlsx 48 . The code with detailed comments can be found in Supplementary Material.

Data availability
Materials, protocols and data reported in this paper are available from the corresponding author on reasonable request.