ENPP2 Methylation in Health and Cancer

Autotaxin (ATX) encoded by Ectonucleotide Pyrophosphatase/Phosphodiesterase 2 (ENPP2) is a key enzyme in Lysophosphatidic Acid (LPA) synthesis implicated in cancer. Although its aberrant expression has been reported, ENPP2 methylation profiles in health and malignancy are not described. We examined in silico the methylation of ENPP2 analyzing publicly available methylome datasets, to identify Differentially Methylated CpGs (DMCs) which were then correlated with expression at gene and isoform levels. Significance indication was set to be FDR corrected p-value < 0.05. Healthy tissues presented methylation in all gene body CGs and lower levels in Promoter Associated (PA) regions, whereas in the majority of the tumors examined (HCC, melanoma, CRC, LC and PC) the methylation pattern was reversed. DMCs identified in the promoter were located in sites recognized by multiple transcription factors, suggesting involvement in gene expression. Alterations in methylation were correlated to an aggressive phenotype in cancer cell lines. In prostate and lung adenocarcinomas, increased methylation of PA CGs was correlated to decreased ENPP2 mRNA expression and to poor prognosis parameters. Collectively, our results corroborate that methylation is an active level of ATX expression regulation in cancer. Our study provides an extended description of the methylation status of ENPP2 in health and cancer and points out specific DMCs of value as prognostic biomarkers.


Introduction
ATX encoded by ENPP2 is a secreted lysophospholipase D (lysoPLD) and belongs to the ENPP (1-7) protein family [1]. ATX is responsible for the synthesis of the majority of extracellular LPA in blood. LPA acts locally upon increased ATX levels through at least six G protein-coupled receptors [2]. Increased ATX activity and levels have been correlated with several inflammatory [3] and fibroproliferative conditions [4], as well as with cancer [5]. In particular, increased expression of ATX in blood and the subsequent increase of LPA have been correlated with cancer invasiveness [6]. In addition, it has been shown that ATX expression is upregulated in cancerous [7,8] and fibrotic tissues [9].
ENPP2 contains 26 introns and 27 exons and is located in the human chromosomal region 8q24 [10], a region with frequent genetic alterations in many cancers [11]. ENPP2 is characterized by alternative splicing of mRNA. The best-known splice variants of ENPP2 are isoforms alpha, beta and gamma; between them, differences in the stability and expression pattern have been documented among several tissues [12].
A study in Biliary Atresia (BA) showed hypomethylation at four CpGs of ENPP2 promoter in the blood and liver of BA patients in relation to normal tissue and was correlated to increased ATX expression [21]. ENPP2 promoter hypermethylation and gene underexpression was found in primary invasive breast carcinomas [13]. Similarly, in breast cancer cell lines, a promoter-associated CpG (cg02156680) of ENPP2 was found highly methylated [22].
In the present study, we studied in silico the methylation of ENPP2 in health and several malignancies and correlated it with gene and isoform expression, aggressiveness and prognosis. Cancer types included in our study were chosen based on their high world incidence, mortality and prevalence [23], as well as access to readily available suitable highthroughput datasets. We examined publicly available methylation datasets from readings by the Illumina methylation bead-chip arrays found in Gene Expression Omnibus (GEO), to identify Differentially Methylated CpGs (DMCs) of ENPP2 between health and disease. Lung, prostate and liver cancer presented a greater number of Promoter Associated (PA) DMCs for ENPP2 and were further pursued using large datasets retrieved from The Cancer Genome Atlas (TCGA), which allowed DMC correlation to clinicopathological parameters and gene expression. A workflow of our study is presented in Figure 1. aberrant gene-specific methylation has been correlated with many pathologies, such cancer [15][16][17][18][19][20]. However, data on the methylation profile of ENPP2 in health and path ogy are fragmented. A study in Biliary Atresia (BA) showed hypomethylation at f CpGs of ENPP2 promoter in the blood and liver of BA patients in relation to normal tis and was correlated to increased ATX expression [21]. ENPP2 promoter hypermethylat and gene under-expression was found in primary invasive breast carcinomas [13]. Si larly, in breast cancer cell lines, a promoter-associated CpG (cg02156680) of ENPP2 w found highly methylated [22].
In the present study, we studied in silico the methylation of ENPP2 in health a several malignancies and correlated it with gene and isoform expression, aggressiven and prognosis. Cancer types included in our study were chosen based on their high wo incidence, mortality and prevalence [23], as well as access to readily available suita high-throughput datasets. We examined publicly available methylation datasets fr readings by the Illumina methylation bead-chip arrays found in Gene Expression Omni (GEO), to identify Differentially Methylated CpGs (DMCs) of ENPP2 between health a disease. Lung, prostate and liver cancer presented a greater number of Promoter Asso ated (PA) DMCs for ENPP2 and were further pursued using large datasets retrieved fr The Cancer Genome Atlas (TCGA), which allowed DMC correlation to clinicopatholog parameters and gene expression. A workflow of our study is presented in Figure 1.

Methylation and Statistical Analysis
Methylation analysis was carried out using normalized beta values ranging between 0 (no methylation) and 1 (full methylation) representing methylation levels of each CpG site (Level 3 data). The Kolmogorov-Smirnov test was applied to check for normality in distribution. Statistical analysis was performed using IBM SPSS 19.0 statistical software (IBM Corp. 2010. IBM SPSS Statistics for Windows, Version 19.0. Armonk, NY, USA). One-Way ANOVA tests followed by Bonferroni post-hoc or Kruskal-Wallis tests were used for comparisons of continuous variables among three or more different groups. In the case of binary variables, t-test or Mann-Whitney tests were also applied. Pearson or Spearman correlation was applied to compare two continuous variables. Differentially Methylated CpGs (DMCs) for ENPP2 were identified based on the False Discovery Rate (FDR-adjusted p-value < 5.00 × 10 −2 ).

In Silico Determination of Transcription Factor (TF) Binding
In order to examine if the DMCs identified were correlated to ENPP2 expression, we further analyzed promoter regions to locate TF binding sites. Hence, PROMO (http://alggen. lsi.upc.es/, accessed on 18 March 2021) [37] tool was used in order to define possible TFs binding in identified DMCs inside ENPP2 promoter. Only human factors and human sites were considered for a TFs search.

Expression and Methylation Correlation Analysis Using TCGA Datasets
Normalized (gene and isoform level) RNA-seq (Illumina HiSeq), level 3 methylation legacy data (Infinium Human Methylation 450 K bead-chip) and corresponding available clinical data were retrieved from prostate adenocarcinoma, lung adenocarcinoma and liver hepatocellular carcinoma TCGA projects representing PC, LC and HCC cases, using the TCGAbiolinks R package [38]. In total, 212 LC, 235 PC and 241 HCC cancer cases along with adjacent healthy lung (15 cases), prostate (35 cases) and liver (42 cases) tissues were obtained. More specifically, cases were chosen to include both mRNA expression (gene and isoforms) and methylation data, from which all matched control and tumor cases were retrieved along with 200 additional tumor samples per cancer type. In the rare case of a case ID being represented by more than one methylation or expression file, the weighted average of the respective values was used for downstream analysis (all weights sum up to the unit). Spearman correlation was calculated per cancer type using the cor.test function. The cutoff level of significance was set to be FDR corrected p-value < 0.05. Last, p-values of linear models fitted between methylation and expression levels (lm R function) were used to test and establish the importance of small correlation coefficients. Analyses were performed using R version 4.0.4.

Expression, Methylation and Survival Analysis Using the UALCAN Database
In order to further verify our results, we used the UALCAN database (http://ualcan. path.uab.edu/, accessed on 10 September 2021) [39] that enables researchers to analyze cancer archived omics data. We performed expression, methylation and survival analysis of ENPP2 in the three types of cancer used in our analysis (PC, LC, HCC) along with the corresponding controls. According to UALCAN, different beta value cut-offs have been considered to indicate hyper-methylation [beta value: 0.7 − 0.5] or hypo-methylation [beta-value: 0.3 − 0.25]. For mRNA expression, methylation and survival, we used TCGA gene analysis, and the screening conditions were as follows: gene "ENPP2", TCGA dataset "Prostate adenocarcinoma", "Lung adenocarcinoma", "Lung Squamous cell carcinoma", "Liver hepatocellular carcinoma". We then used "expression", "methylation" and "survival" as links for analysis of each cancer. Protein expression data were available only for lung adenocarcinoma, and Clinical Proteomic Tumor Analysis Consortium (CPTAC) datasets were used for analysis. For protein analysis, Z-values represent standard deviations from the median across samples for the given cancer type. Log2 Spectral count ratio values from CPTAC were first normalized within each sample profile and then normalized across samples.

Analysis of ENPP2 Methylation from GEO Datasets
In silico methylation analysis of ENPP2 was performed using methylome data retrieved using the GEO. The results are described below.

ENPP2 Methylation in Normal Tissues
In order to examine the methylation profile of ENPP2 across different human healthy tissues, we analyzed methylome data from 17 healthy tissues included in the GSE50192 study. We observed a consistent methylation pattern across all studied tissues ( Figure 2), with methylation being increased in all 7 CGs in the gene body region and decreased in 5 CGs in the Transcription Start Site (TSS) and 1 in the 1st exon.

ENPP2 Methylation in Normal Tissues
In order to examine the methylation profile of ENPP2 across different human healthy tissues, we analyzed methylome data from 17 healthy tissues included in the GSE50192 study. We observed a consistent methylation pattern across all studied tissues ( Figure 2), with methylation being increased in all 7 CGs in the gene body region and decreased in 5 CGs in the Transcription Start Site (TSS) and 1 in the 1st exon.

ENPP2 Methylation in Tissues from Cancer Patients
In order to unravel aberrant ENPP2 methylation in cancer, we compared methylomes of malignant vs. respective benign tissues from 7 different cancer types, using 10 GEO datasets (GSE113017, GSE113019, GSE120878, GSE27130, GSE63704, GSE76938, GSE98534, GSE46306, GSE134772, GSE97686) ( Table 1). In total, 13 DMCs were identified in 5 cancers, i.e., HCC (12 DMCs), PC (10 DMCs), LC (9 DMCs), melanoma (7 DMCs), CRC (1 DMC) ( Table 2), most of which were common between them (Table 3), whereas no DMCs were identified in Precancerous Interepithelial Cervical Neoplasia (CIN) and Cancer (CC) and in Gastric Cancer (GC). With two exceptions, all gene body DMCs showed decreased methylation in cancer in relation to their controls. Most importantly, all CGs located in the promoter-associated region and the 1st exon, regions known to hold an important role in transcriptional regulation [38,40] were DMCs across different cancer types, all presenting increased methylation. These results demonstrate aberrant methylation of ENPP2 in the majority of cancer types studied, following a specific pattern pointing to down-regulation of expression.

ENPP2 Methylation Was Correlated to Aggressiveness in Cancer Cell Lines
In order to study any relation of ENPP2 methylation to cancer aggressiveness, we compared cell lines from hepatocellular and prostate cancer presenting a more (SKHEP1 and PC3 respectively) or less (HEPG2 and LNCAP respectively) invasive phenotype, using the GSE71627 study dataset [33]. In total, 12 DMCs were identified ( Table 4), 6 of which were common in both cancer types. All 8 DMCs identified in HCC cell lines were also found in the liver tumor methylomes, whereas the common DMCs between prostate tissues and cell lines were 6/8. With two exceptions, all DMCs across the whole gene (1st exon, TSS, body) showed higher methylation in the more aggressive hepatocellular and prostate cell lines. These observations suggest an involvement of ENPP2 methylation in cancer aggressiveness. Interestingly, treatment of colon cancer cell lines with the DNA methylation inhibitor 5-aza-2 -deoxycytidine (GSE51815 study) caused a decrease of methylation in all 14 DMCs located throughout ENPP2 (Supplementary Table S1), implying that methylation could present a potential therapeutic target to reverse the aggressive phenotype.

In Silico Analysis of TF Binding on the ENPP2 Promoter
Regulation of gene expression via DNA methylation occurs mainly by disturbing TF and RNA polymerase binding to sites known to be necessary for initiation of transcription [41]. To support that the identified DMCs on ENPP2 may actually play a role in regulating ATX expression, we examined if they are located within TF binding promoter regions that could initiate transcription. Analysis using the PROMO tool predicted 39 putative TFs that could bind to the ENPP2 promoter (Figure 3), 4 of them (TFIID, GR, GR-beta, C/EBPbeta) on or in proximity to cg04452959, 7 TFs (TFII-I, GR-alpha, GATA-1, E2F-1, Pax-5, p53, Sp1) in cg02709432, 7 TFs (C/EBPbeta, C/EBPalpha, Pax-5, p53, ENKTF-1, YY1, GR-beta) in cg02156680 and 3 TFs (PEA3, GATA-1, XBP-1) in cg06998282. Interestingly, those 4 CGs located 200 nucleotides upstream of and up to the TSS (TSS200) (first 2), or 200 to 1500 nucleotides upstream of the TSS (TSS1500) (last 2) were identified as DMCs in most of the malignancies examined and between more and less aggressive cell lines. Collectively, these findings show that DMCs identified in the ENPP2 promoter in cancer are found in sites significant for TF binding, and therefore, altered methylation is likely to affect transcription and expression of ATX. Interestingly, treatment of colon cancer cell lines with the DNA methylation inhibitor 5-aza-2′-deoxycytidine (GSE51815 study) caused a decrease of methylation in all 14 DMCs located throughout ENPP2 (Supplementary Table S1), implying that methylation could present a potential therapeutic target to reverse the aggressive phenotype.

In Silico Analysis of TF Binding on the ENPP2 Promoter
Regulation of gene expression via DNA methylation occurs mainly by disturbing TF and RNA polymerase binding to sites known to be necessary for initiation of transcription [41]. To support that the identified DMCs on ENPP2 may actually play a role in regulating ATX expression, we examined if they are located within TF binding promoter regions that could initiate transcription. Analysis using the PROMO tool predicted 39 putative TFs that could bind to the ENPP2 promoter (Figure 3), 4 of them (TFIID, GR, GR-beta, C/EBPbeta) on or in proximity to cg04452959, 7 TFs (TFII-I, GR-alpha, GATA-1, E2F-1, Pax-5, p53, Sp1) in cg02709432, 7 TFs (C/EBPbeta, C/EBPalpha, Pax-5, p53, ENKTF-1, YY1, GR-beta) in cg02156680 and 3 TFs (PEA3, GATA-1, XBP-1) in cg06998282. Interestingly, those 4 CGs located 200 nucleotides upstream of and up to the TSS (TSS200) (first 2), or 200 to 1500 nucleotides upstream of the TSS (TSS1500) (last 2) were identified as DMCs in most of the malignancies examined and between more and less aggressive cell lines. Collectively, these findings show that DMCs identified in the ENPP2 promoter in cancer are found in sites significant for TF binding, and therefore, altered methylation is likely to affect transcription and expression of ATX.

ENPP2 Methylation and Expression Analysis from TCGA Datasets
An important objective of our study was to address if and how aberrant methylation of ENPP2 is related to alterations in gene expression. Based on our findings, three cancer entities presenting the highest number of DMCs were selected for further study, in order to confirm altered ENPP2 methylation in larger cohorts and correlate them with expression at gene and isoform levels. For this purpose, several available datasets including PC, LC and HCC readings of mRNA expression (gene and isoforms) and methylation along with the available clinical and demographic data were downloaded from TCGA.

ENPP2 Methylation and Expression Analysis from TCGA Datasets
An important objective of our study was to address if and how aberrant methylation of ENPP2 is related to alterations in gene expression. Based on our findings, three cancer entities presenting the highest number of DMCs were selected for further study, in order to confirm altered ENPP2 methylation in larger cohorts and correlate them with expression at gene and isoform levels. For this purpose, several available datasets including PC, LC and HCC readings of mRNA expression (gene and isoforms) and methylation along with the available clinical and demographic data were downloaded from TCGA.

ENPP2 Methylation and Expression Analysis in Prostate Cancer
Comparisons were performed between 235 prostate adenocarcinoma tumors and 35 healthy prostate tissues (Table 5). In general, results confirmed those from the GEO datasets. In total, 12 DMCs were identified between cancer and control tissues (5 in TSS, 1 in 1st exon and 6 in the gene body), 10 of which were common to those found in the GSE76938 dataset. All DMCs in TSS and the 1st exon presented increased methylation in PC in relation to controls, whereas decreased methylation was noticed in 3 out of 6 DMCs in the gene body area. Following this, we examined correlation of DMCs to clinicopathological patient characteristics, to reveal associations with prognosis. Methylation analysis in relation to available patient data (age, race, nodal status, relapse, tumor size and treatment response) showed a significant correlation with tumor size, as increased methylation of 3 CGs, namely, cg02534163 (1st exon), cg02709432 (TSS200) and cg23725583 (gene body), was found in larger tumors in relation to smaller tumors (p = 0.032). Furthermore, non-response to pharmacotherapy was correlated with increased methylation of cg01243251 in the gene body region (p = 0.023). No other correlations were found in relation to age, race, nodal status and the event of relapse. mRNA expression analysis in the same samples showed decreased levels in PC in relation to normal tissues (LogFC: −0.379, FDR: 3.70 × 10 −2 ), indicating that the increased methylation of ENPP2 in PA regions is correlated with the decreased expression of ENPP2 in PC. Spearman correlation of mRNA ENPP2 expression (at gene and isoform level) per CG site revealed statistically significant correlations shown in Figure 4A and Table 6. Between gene body CGs, a tendency towards positive correlation of mRNA expression to cg01243251 and cg20162626 methylation was observed and a negative to cg07236691. TSS CG sites cg02156680, cg02709432, cg06998282, cg14409958 and 1st exon cg02534163 methylation showed a negative correlation with expression, showing that the increased methylation at these regions is associated with decreased expression. Interestingly, although Spearman's coefficient is relatively small for TSS and 1st exon CGs, successful fit of a linear model further supports the existence of an expression-methylation relationship ( Table 6, coefficient p-value column). No significant correlations emerged between ENPP2 expression and methylation in control prostate tissue.     In order to unfold the impact of CG methylation on ENPP2 isoform expression [12,42], we downloaded mRNA expression data from ENPP2 isoforms, i.e., isoform alpha (uc003yos.1), isoform beta (uc003yor.1 and uc003yot.1) and isoform gamma (uc010mdd.1). Uc003yot.1, uc003yos.1 and uc003yor.1 isoform expression showed statistically significant correlation with the methylation of several CGs. In specific, although they were all characterized by small effect sizes, uc003yor.1 expression is linearly related to the methylation levels of TSS and 1st exon CGs cg02156680, cg06998282, cg14409958 and cg02534163, respectively, further strengthening the observed correlation (Supplementary Figure S1 and Table S2). Finally, no relation emerged between expression of any of the ENPP2 isoforms and methylation of CGs in healthy prostate samples, consistent with the observations at the gene level.

ENPP2 Methylation and Expression Analysis in Lung Cancer
Analysis was performed between 212 LC adenocarcinoma tumors and 15 healthy lung tissues, and results presenting statistically significant correlations are shown in Table 7. Findings confirmed those from the GEO datasets. Eight DMCs were identified between cancer and control tissues (3 in the TSS, 1 in the 1st exon, 4 in the gene body) and 6 of them were common to those found in the GSE76938 dataset. DMCs of ENPP2 showed upregulation of methylation in TSS (cg04452959, cg06998282, cg14409958) and the 1st exon (cg02534163) and downregulation in the gene body (cg07236691, cg09444531, cg20048037, cg20162626). Methylation was also correlated to the available clinicopathological characteristics of LC and normal lung tissue samples (gender, age, nodal status, distance metastasis, relapse, tumor size and treatment response and stage). In LC samples, increased methylation of cg14409958 (TSS) was significantly correlated with advanced cancer stage (p = 0.035). Differential mRNA expression analysis in the same samples showed decreased levels in LC in relation to normal tissues (LogFC: 1.285, FDR: < 1.00 × 10 −2 ) similarly to PC, indicating that in cancer the increased methylation of PA CGs is correlated to decreased autotaxin expression. The impact of ENPP2 methylation on its expression was examined in LC and healthy lung tissue samples. Spearman correlation of mRNA expression per CG resulted in a single statistically significant correlation ( Figure 4B and Table 6). A reverse correlation was noticed between methylation of cg14409958 (TSS) and mRNA expression, suggesting again the DNA methylation role in repressing expression. Fit of a linear model once again reinforced the observed correlation. On the other hand, control samples did not show any statistically significant correlation after p-value correction, and only the methylation of body site cg7236691 showed a significant correlation coefficient along with a linear relationship to ENPP2 expression levels. Last, no significant correlations were witnessed between methylation and expression levels of all isoforms examined, yet large rho values and significant linear model fit propose the existence of such a relationship.

ENPP2 Methylation and Expression Analysis in Hepatocellular Carcinoma
Analysis was performed between 241 HCC tumors and 42 control liver tissues. Statistically significant correlations are presented in Table 8. In total, 13 DMCs were identified between cancer and control (5 in the TSS, 1 in the 1st exon, 7 in the gene body) and 12 were common to those found in GSE113017 and GSE113019. Again, downregulation of methylation was noticed in all gene body CGs (cg00320790, cg23725583 and cg01243251) and upregulation of methylation in all TSS and 1st exon related CGs (cg02156680, cg02534163, cg02709432, cg04452959, cg06998282 and cg14409958) in HCC. Methylation of ENPP2 was also correlated to available clinical and demographic characteristics of the HCC cohort. Interestingly, in the tumor samples, increased methylation of the majority of the ENPP2 CGs (cg00320790, cg01243251, cg02156680, cg02709432, cg07236691, cg09444531, cg14409958, cg20048037, cg20162626, cg23725583) (all p < 0.05) was noticed in women in relation to men. In addition, a negative correlation was found between age and methylation of cg00320790, cg01243251, cg07236691, cg09444531, cg20048037, cg20162626 and cg23725583 (all p < 0.001), i.e., younger people presented increased methylation in relation to older. Finally, increased methylation of cg04452959 was correlated to tumors with macro invasion in relation to those with no or micro invasion (p = 0.044). No correlation was noticed between methylation and BMI, hepatic inflammation, Ishak fibrosis, relapse, family history, grade, stage or tumor size. Analysis in normal samples showed a gender correlation only for one CG (cg20048037) which presented increased methylation in females (p = 0.033) in relation to males. Negative correlation was also noticed between cg01243251 methylation and age (p = 0.037). Finally, no relationship was found between BMI and ENPP2 methylation in normal samples. mRNA expression analysis in the same samples showed increased levels in HCC in relation to normal tissues (LogFC: 0.710, FDR: 1.00 × 10 −2 ), i.e., the opposite of LC and PC observations, suggesting a methylation-independent and a cancer type-specific regulation of ENPP2 in HCC. Spearman correlation of mRNA expression (at gene and isoform levels) per CG site revealed the most numerous statistically significant correlations, compared to PC and LC samples, shown in Figure 4C and Table 6. Between gene body CGs, a positive correlation of mRNA expression to cg00320790, cg01243251, cg07236691, cg09444531, cg20048037 and cg20162626 methylation was observed. Apart from the significant correlations established, methylation of the aforementioned CGs was characterized by a linear relationship to ENPP2 expression, further supporting dependence of the latter on the former. Control samples also showed positive correlation between ENPP2 mRNA expression and methylation of 2 gene body CGs (cg20048037, g20162626). Finally, cg23725583 of body and cg02709432, cg04452959, cg06998282 and cg02156680 of TSS regions showed negative correlation of methylation in relation to ENPP2 mRNA expression. Isoform analysis for the control tissues showed similar correlation patterns (Supplementary Figure S2 and Table S2).

ENPP2 Methylation, Expression and Survival Analysis via UALCAN
In order to further verify our findings, we conducted expression, methylation and survival analysis of ENPP2 in PC (all adenocarcinoma cases), LC (adenocarcinoma and squamous cell carcinoma cases) and HCC using the UALCAN database. Analysis confirmed the above results as ENPP2 mRNA was under-expressed in PC (p = 9.31 × 10 −3 , Figure 5A) and LC (adenocarcinoma, p = 1.68 × 10 −3 and squamous cell carcinoma, p = 4.52 × 10 −3 , Figure 6A,C) and upregulated in HCC (p = 2.38 × 10 −10 , Figure 7A). Protein expression analysis was available only for LC adenocarcinoma cases, showing downregulation in primary tumor tissues in relation to normal tissues ( Figure 6E, p = 1.78 × 10 −4 ). Next, methylation analysis revealed upregulation in all cancer types in relation to normal tissues (PC, p = 1.62 × 10 −12 , HCC, p = 1.11 × 10 −16 and LC, p = 1.62 × 10 −12 for both types) as depicted in Figures 5-7, in accordance with our previous observations. Methylation and expression results via UALCAN strengthen our findings, showing that the ENPP2 gene is methylated in LC, HCC and PC and this is related to under-expression in LC and PC, suggesting a causative relationship in these two cancer types and a cancer-specific regulatory mechanism in HCC. Finally, survival analysis did not reveal any statistical significance for any of the studied cancers, as depicted in Supplementary Figure S4A-D.
lationship to ENPP2 expression, further supporting dependence of the latter on the former Control samples also showed positive correlation between ENPP2 mRNA expression and methylation of 2 gene body CGs (cg20048037, g20162626). Finally, cg23725583 of body and cg02709432, cg04452959, cg06998282 and cg02156680 of TSS regions showed negative corre lation of methylation in relation to ENPP2 mRNA expression. Isoform analysis for the control tissues showed similar correlation patterns (Supplementary Figure S2 and Table S2).

ENPP2 Methylation, Expression and Survival Analysis via UALCAN
In order to further verify our findings, we conducted expression, methylation and survival analysis of ENPP2 in PC (all adenocarcinoma cases), LC (adenocarcinoma and squamous cell carcinoma cases) and HCC using the UALCAN database. Analysis confirmed the above results as ENPP2 mRNA was under-expressed in PC (p = 9.31 × 10 −3 Figure 5A) and LC (adenocarcinoma, p = 1.68 × 10 −3 and squamous cell carcinoma, p = 4.52 × 10 −3 , Figure 6A,C) and upregulated in HCC (p = 2.38 × 10 −10 , Figure 7A). Protein expres sion analysis was available only for LC adenocarcinoma cases, showing downregulation in primary tumor tissues in relation to normal tissues ( Figure 6E, p = 1.78 × 10 −4 ). Next methylation analysis revealed upregulation in all cancer types in relation to normal tissues (PC, p = 1.62 × 10 −12 , HCC, p = 1.11 × 10 −16 and LC, p = 1.62 × 10 −12 for both types) as depicted in Figures 5-7, in accordance with our previous observations. Methylation and expression results via UALCAN strengthen our findings, showing that the ENPP2 gene is methylated in LC, HCC and PC and this is related to under-expression in LC and PC, suggesting a causative relationship in these two cancer types and a cancer-specific regulatory mecha nism in HCC. Finally, survival analysis did not reveal any statistical significance for any of the studied cancers, as depicted in Supplementary Figure S4A-D.

Discussion
ATX encoded by ENPP2 is a secreted glycoprotein that forms LPA [42]. The ATX-LPA axis is related to many physiological processes, including embryonic development and wound healing. Dysregulation of ATX expression is connected with various pathological conditions such as cancer, inflammatory diseases and fibrosis [3][4][5]. The exact mechanism by which ENPP2 expression is regulated is still not fully understood, whereas recently, it has been proved that ENPP2 is prone to epigenetic alterations [13]. Still, very little information is available about its DNA methylation pattern and the consequent impact in gene expression in health and human pathology.
In the present study we adopted a bioinformatic in silico approach using publicly available datasets from healthy tissues and different cancer tissues and cell lines to analyze methylation patterns of ENPP2. Our analysis showed a consistent methylation pattern throughout the gene's regions across human tissues, i.e., increased methylation in the gene body and decreased methylation in TSS and the 1st exon. Given the fact that ENPP2 is expressed in almost all tissues and biological fluids [12,43,44], we can postulate that the decreased methylation in the TSS and 1st exon is associated with the active transcription of the gene in most human tissues.
Analysis of cancer datasets revealed aberrant ENPP2 methylation, showing a malignant-specific profile throughout different cancer types. In general, methylation was increased in the TSS and 1st exon, regions known to hold an important role in gene expression, and decreased in the gene body region. A large number of DMCs were identified between malignant and respective benign tissues. Most importantly, all six DMCs of ENPP2 located at TSS in the promoter or at the 1st exon showed increased methylation across different cancer types, including HCC, melanoma, CRC, LC and PC. These results corroborate and expand recent observations showing a hypermethylated ENPP2 promoter in primary tumors of LC and squamous cell carcinoma patients [45] and in breast cancer [13,22,46].
Based on these interesting observations, we next performed in silico analysis of ENPP2 methylation in datasets retrieved from the TCGA, focusing on those cancer types presenting the greatest number of DMCs, i.e., LC, PC and HCC. TCGA datasets are generally larger compared to those of other research efforts, allowing comparisons of stronger statistical relevance, and most significantly, they contain several clinical and demographic parameters of each patient. In addition, the datasets selected included also mRNA expression data and were therefore suitable for addressing an important objective of this study,

Discussion
ATX encoded by ENPP2 is a secreted glycoprotein that forms LPA [42]. The ATX-LPA axis is related to many physiological processes, including embryonic development and wound healing. Dysregulation of ATX expression is connected with various pathological conditions such as cancer, inflammatory diseases and fibrosis [3][4][5]. The exact mechanism by which ENPP2 expression is regulated is still not fully understood, whereas recently, it has been proved that ENPP2 is prone to epigenetic alterations [13]. Still, very little information is available about its DNA methylation pattern and the consequent impact in gene expression in health and human pathology.
In the present study we adopted a bioinformatic in silico approach using publicly available datasets from healthy tissues and different cancer tissues and cell lines to analyze methylation patterns of ENPP2. Our analysis showed a consistent methylation pattern throughout the gene's regions across human tissues, i.e., increased methylation in the gene body and decreased methylation in TSS and the 1st exon. Given the fact that ENPP2 is expressed in almost all tissues and biological fluids [12,43,44], we can postulate that the decreased methylation in the TSS and 1st exon is associated with the active transcription of the gene in most human tissues.
Analysis of cancer datasets revealed aberrant ENPP2 methylation, showing a malignantspecific profile throughout different cancer types. In general, methylation was increased in the TSS and 1st exon, regions known to hold an important role in gene expression, and decreased in the gene body region. A large number of DMCs were identified between malignant and respective benign tissues. Most importantly, all six DMCs of ENPP2 located at TSS in the promoter or at the 1st exon showed increased methylation across different cancer types, including HCC, melanoma, CRC, LC and PC. These results corroborate and expand recent observations showing a hypermethylated ENPP2 promoter in primary tumors of LC and squamous cell carcinoma patients [45] and in breast cancer [13,22,46].
Based on these interesting observations, we next performed in silico analysis of ENPP2 methylation in datasets retrieved from the TCGA, focusing on those cancer types presenting the greatest number of DMCs, i.e., LC, PC and HCC. TCGA datasets are generally larger compared to those of other research efforts, allowing comparisons of stronger statistical relevance, and most significantly, they contain several clinical and demographic parameters of each patient. In addition, the datasets selected included also mRNA expression data and were therefore suitable for addressing an important objective of this study, i.e., if aberrant methylation is correlated to gene expression. Methylation, clinical and expression data were recovered for the three cancer types. Differential methylation analysis of ENPP2 revealed that all emerged DMCs identified in transcription-related (TSS and 1st exon) regions were hypermethylated in all three cancers compared to healthy controls, confirming the analysis of the GEO datasets. In addition, the majority of DMCs located at the gene body were hypomethylated in the studied cancers in relation to controls. mRNA levels were decreased in PC and LC in relation to normal tissues. Collectively, our results indicate that the increased methylation of PA and 1st exon CGs is correlated with decreased expression in lung and prostate cancer. This is in line with previous studies in LC and BC showing that ENPP2 is hypermethylated in tumor tissues in relation to normal, causing down regulation in gene expression [13,45]. In PC, ATX protein was not or was weakly expressed in non-neoplastic epithelial cells and in high-grade intra-epithelial neoplasia, while in cancer cells ATX was only expressed in half of the tumors and was correlated with adverse tumor parameters [47]. A relevant study in LC showed that ATX protein expression and activity was increased in LC tissues and sera [48]. As far as HCC is concerned, our analysis showed upregulation of expression in HCC in relation to normal liver, showing a TSS and 1st exon methylation-independent and a cancer type-specific role of ENPP2 expression regulation. In a previous study, ATX overexpression in HCC tissues was correlated with inflammation and liver cirrhosis. In addition, liver cancer cell lines presented stronger ATX expression in relation to normal hepatocytes [49]. It should be noted that many authors have demonstrated that the relationship between mRNA expression and protein differs in many cancers. It has been reported in lung cancer and glioblastoma that, for many genes, mRNA expression is lower but protein levels are higher compared with the control [50][51][52][53].
In agreement with the above findings, analysis using the UALCAN database showed that ENPP2 is hypermethylated and under-expressed in LC and PC, suggesting that DNA methylation regulates expression in LC and PC. However, no regulatory relation was observed between methylation and expression in HCC, as both were upregulated, pointing again to a cancer-specific methylation-independent ENPP2 regulation. Different mechanisms between cancer types are common. Here, our presented results from the cancer types studied indicate a cancer type-specific profile of ENPP2 methylation rather than a similar pan-cancer dysregulation. Without availability of suitable methylome datasets or targeted methylation studies of ENPP2 in each different cancer type, we cannot extrapolate conclusions between cancers.
The same correlation pattern was noticed for ENPP2 isoforms in all cancer types studied. Interestingly, there was a significant negative correlation between mRNA expression (gene and isoform alpha and beta) and promoter methylation in four CGs (cg02156680, cg02709432, cg04452959 and cg06998282) in PRAD. In LC samples, the methylation of cg06998282 and cg02709432 was negatively correlated with the expression of ENPP2 and also with isoform beta and gamma. Finally, in the case of HCC, only the methylation of cg06998282 was negatively correlated with the expression of ENPP2 and isoform beta. The above findings indicate that the promoter methylation of specific CGs is negatively correlated with ENPP2 and isoform expression differs between cancers, with cg02709432 being a common site in PC and LC but not in the case of HCC. This CG is located at a site that can bind E2F-1 TF, which has been shown to be inhibited by CG methylation [54], and Sp1 TF, which has been found to regulate ENPP2 transcription [55]. Thus, we hypothesized that as the level of methylation increases, methylation of cg02709432 hinders the binding of the TFs to the promoter, thus leading to reduction in ENPP2 gene and isoforms expression.
The expression pattern of isoforms differs between tissues as high expression levels of isoform beta were found in peripheral tissues and plasma, while isoform gamma was mostly found in the brain, and isoform alpha is considered to be the most under-expressed in brain and peripheral tissue in comparison to the other two [56]. According to a relevant study, isoform alpha has a deletion of exon 12, isoform beta a deletion of exons 12 and 21 and isoform gamma a deletion of exon 21 [12], leading to different splice variants. None of the identified DMCs were located at these regions, explaining similar patterns of ENPP2 mRNA and isoform expression. DNA methylation within promoters is known to modulate the binding of TFs to regulatory elements, thus resulting in transcriptional repression [57]. In our study, we predicted 39 TFs which can regulate transcription through binding to ENPP2 promoter's DMCs. Therefore, any aberrant methylation events in these DMCs during pathological transformation may block TF binding and related transcription. This is further supported by reports involving the identified TFs in ENPP2 and ATX expression. Indeed, among the predicted TFs, NF kappaB, AP-2 and E2F have been previously shown to be sensitive to CG methylation with consequent inhibition of their DNA binding activities [54]. Another TF predicted to bind DMCs of ENPP2, NFAT1, has been shown to mediate ATX overexpression in MDA-MB-435 cells [58]. It has also been shown that blocking the expression of NFAT1 results in downregulation of ATX expression, leading to inhibition of melanoma and metastasis [35]. High C-Jun levels seem to enhance ENPP2 expression [59]. Interestingly, SP was found to regulate ENPP2 transcription in neuroblastoma cells by activating a CRE/AP-1-like element at position −142 to −149 and a GAbox at position −227 to −235 near the TSS of ENPP2 [55]. This is in accordance with our finding that Sp1 can bind near the cg02709432 located at TSS200.
In order to assess any correlation of ENPP2 methylation to tumor prognosis, clinical characteristics analysis was performed and showed that increased methylation of some CGs was correlated with poor tumor parameters. Indeed, in PC it was associated with larger tumors and non-response to pharmacotherapy, in LC it was connected to the advanced cancer stage and in HCC it was associated with macro-invasion. Hence, ENPP2 methylation in the identified CGs could be pursued further and be evaluated in clinical cancer samples as biomarkers of cancer progression and poor outcome. In addition, these results corroborate previous data showing that low mRNA expression was associated with worse prognosis in LC [45].
The involvement of ENPP2 methylation in tumor progression and prognosis was also addressed by analyzing methylomes from cell lines presenting a more or less aggressively invasive phenotype, revealing several DMCs. Higher methylation was observed in the more aggressive in relation to less aggressive HCC and PC cell lines, indicating a connection of ENPP2 methylation with worse prognostic behavior, in accordance with our findings in the clinical samples.
Finally, analysis of colon cell lines treated with DNA methyltransferase inhibitors showed that 5-AZA caused a decrease of methylation in all CGs in relation to untreated controls in the three studied cell lines, showing a clear demethylation effect in the ENPP2 gene. Given the contribution of ENPP2 in a variety of pathologies, further studies could assess a methylation-based reprogramming of ENPP2 via a variety of methylation inhibitors. Similarly, previous studies have demonstrated that targeting the ATX-LPA-LPP axis is an attractive strategy for introducing new therapeutic choices [60,61].
In conclusion, healthy tissues presented increased methylation of ENPP2 in the gene body and decreased in the promoter and 1st exon connected to the active transcription of the gene in most human tissues. A different pattern was described in HCC, melanoma, CRC, LC and PC, showing a malignant-specific profile of ENPP2 methylation. Further analysis of independent TCGA datasets confirmed these results as increased methylation of promoter and 1st exon CGs and decreased ENPP2 mRNA expression in PC and LC in relation to healthy tissues were found. Furthermore, increased methylation of ENPP2 was connected to poor prognostic parameters in the same cancers, which was also supported by analysis of cell line datasets. We also found a negative correlation between mRNA expression at gene and isoform levels and methylation of PA CGs that present TF binding sites. In specific, we postulate that the methylation of promoter CGs may hinder the binding of TFs, and thus, the expression of ENPP2 and isoforms may be reduced.
Our findings contribute to the understanding of methylation events and regulatory mechanism of ENPP2 in cancer and provide a full description of DMCs to be further validated in functional and clinical studies.