Genome-wide methylation profiling identify hypermethylated HOXL subclass genes as potential markers for esophageal squamous cell carcinoma detection

Data preparation and preprocessing

GSE149608 [27] and GSE52826 [29] datasets were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) as both consisted of tumor and paired normal samples. The two datasets were generated by whole-genome bisulfite sequencing (WGBS) and Illumina HumanMethylation450k, covering over 18 million and 480,000 CpG sites. The methylation value of each CpG site was represented by the percentage of methylated reads in total reads that covered this site for WGBS data and by the signal intensity value of β (range: 0-1) for 450k microarray data. TCGAbiolinks tools [30] was used to download the level 3 methylation data, the corresponding clinical features of ESCC, and paired normal samples from The Cancer Genome Atlas (TCGA) Program (https://portal.gdc.cancer.gov/). CpG sites were removed during the data preprocessing if the methylation values were 0 in all samples of WGBS dataset or the β values were missing in more than 90% of the samples in 450k dataset. Then, KNN algorithm was used to fill the missing values.

The preprocessed datasets were shown in Table 1, containing 29 normal samples and 109 tumor samples. The clinical information of GSE149608 was collected from GEO database simultaneously but was not consistent with the reference paper. Here, we determined their exact clinical characteristics (Supplementary table 1) according to the information provided by the reference, which included two early-stage patients (patient3 and patient10), five intermediate-stage patients (patient2, patient5-8), and two late-stage patients (patient1 and patient4) stages. All cases in GSE52826 dataset were early-stage ESCC (Supplementary table 2).

Human tissue and blood samples

Twenty formalin-fixed paraffin-embedded esophageal squamous carcinoma and adjacent normal sample were collected from the department of pathology of the First Affiliated Hospital of Zhengzhou University in the study. Meanwhile, twenty healthy blood samples were obtained from the blood transfusion department, where the residual blood samples were collected from healthy individuals after blood donation. The healthy blood samples were included because they represented the methylation status of candidate targets on healthy individuals and will also facilitate investigating their potential role as blood diagnostic markers. Twenty plasma samples were also collected from ESCC patients, but only 13 samples obtained adequate cell-free DNA amounts. Clinical features of the samples are displayed in Supplementary table 3. Referring to the NCCN guidelines, in this study, we defined Tis (high-grade dysplasia), T1 (T1a and T1b), and some T2 (without lymph node metastasis and distal metastasis) as early-stage ESCC. All individual identifiers have been removed. The Ethics Committee of the First Affiliated Hospital of Zhengzhou University approved this study (approval number 2020-KY-0152).

Identification of differentially methylated genes

Since the vast majority of CpG sites in the whole genome were covered in WGBS data represented the comprehensive epigenetic information, the GSE149608 dataset was selected to identify differentially methylated CpG sites (DMC) between tumor and paired normal samples using paired student t-test, with P-value <0.05 and fold change ≥2 as the significance threshold. We defined hypermethylated DMGs in normal and tumor samples as NC-DMCs and ESCC-DMCs respectively. The methylation status of adjacent CpG sites are usually highly coupled with each other, and they tend to be located in a small genomic region [31], therefore we further identified the differentially methylated regions (DMR) based on DMCs using a modified sliding window approach described in our previous study [32]. Briefly, to facilitate developing methylation specific PCR (MSP) assay, the maximum length of DMRs was set to 100 bp, and each DMR contained at least three DMCs, and the distance between two adjacent DMCs was less than 50bp. NC-DMRs and ESCC-DMRs were then identified separately, and the genes covered by NC- or ESCC-DMRs were defined as differentially methylated genes (NC-DMG or ESCC-DMG), except for those overlapped with both NC- and ESCC-DMRs. Since ESCC-DMRs are more eligible for MSP, only ESCC-DMGs were selected for subsequent analysis.

Function enrichment analysis

Gene ontology (GO) and KEGG pathway enrichment analysis was performed using the ‘clusterProfiler’ R package for NC-DMGs and ESCC-DMGs. All human genes were used as the background with q-value < 0.05 as the significance threshold. Gene family enrichment analysis was implemented by integrating gene family terms provided by HGNC (HUGO Gene Nomenclature Committee) [33] and the chi-square test. Briefly, the human gene family data were downloaded from HGNC (www.genenames.org), which records the human genes and their corresponding families. For each family, we counted the number of DMGs belonging to this family. The overrepresented gene families were assessed with all human genes and DMGs as background. The chi-square test P-value< 0.05 was selected as the significance threshold to determine enriched families for NC-DMGs and ESCC-DMGs.

DNA extraction and BS-treatment

DNA of tissue and blood samples were extracted using paraffin-embedded tissue DNA Rapid Extraction Kit (TIANGEN^®, Beijing) and TIANamp Genomic DNA Kit (TIANGEN^®, Beijing), respectively. BS-treatment and purification for the extracted DNA were subsequently carried out using the Nucleic Acid Purification Kit (Ammunition^®, Wuhan). The basic principle of BS-treatment is that unmethylated cytosines in denatured DNA can be converted to uracil by bisulfite ions, while methylated cytosines are remained. The methylation status then can be determined by methylation-specific PCR.

Methylation-specific PCR and Sanger sequencing

Methylation-specific PCR experiments and Sanger sequencing were performed to determine the methylation status of candidate targets. Two genes, HOXD1 and GSX1, were selected for MSP because they showed good abilities to discriminate cancer samples from normal samples. Besides, high densities of differentially methylated CpG sites were found within the two genes, which allowed us to design appropriate MSP primers. The designed primers for HOXD1 and GSX1 were listed in Table 2. The amplified regions covered DMR4, DMR5, and DMR6 of HOXD1, and DMR16 and DMR17 of GSX1 (Supplementary table 4). The fully unmethylated and methylated DNA fragments for the two targets were synthesized as negative and positive controls. The MSP amplification system was 20ul, including 7ul of dd-water, 10ul of 2 ×T5 Fast qPCR Mix (SYBRGreenI), 0.5ul of forward and reverse primers (10 uM), and 2ul of BS-converted DNA. PCR reaction was pre-denaturation at 95°C for 1 min, denaturation at 95°C for 10 s, annealing and extension at 60°C for 45s. The melt curve reaction was 95°C ∼ 15s, 60°C ∼ 1 min, and 95°C ∼ 15s. Quantitative real-time PCR was performed on a 7500 device (Thermo Fisher, USA). The MSP products were then used for Sanger sequencing. After BS-treatment, the methylated CpG site would remain C, and the un-methylated C would be T. Therefore, the methylation status of target CpG sites can be determined according to the result of Sanger sequencing.

ROC analysis

For the candidate target genes, we view the sample types (normal or tumor) as the response variable of their methylation values to develop a classifier for ESCC using logistic regression. The logistic probability, herein defined as risk scores for samples, was then estimated. Typically, samples would be recognized as positive if the risk scores were > 0.5 and vice versa as negative. This study used ROC curve to determine the appropriate threshold and AUC value to assess the classifier performance. The optimal threshold was locked when Youden index reached the maximum. Samples were subsequently divided into a positive group if the scores were higher than this threshold or a negative group if the scores were less than this threshold. All samples were classified into four categories, true positive (TP), true negative (TN), false positive (FP), and false negative (FN) according to their sample types and predicted types. The sensitivity and specificity were calculated using the following formulas:

Statistical analysis

Data preprocessing, statistical analysis, and other analysis in this study were implemented in R software (v3.6.1). For continuous variables, paired student t-test was used for comparisons of paired normal and tumor samples, and the wilcoxon rank-sum test for unpaired datasets. The Kruskal-Wallis test was performed for comparisons of multiple groups. For category variables, chi-square test was used to estimate the difference between groups. Logistic regression was implemented using the R package ‘glmnet’ with the parameter ‘family=binomial’. ROC analysis was conducted using the package ‘pROC’ [34] with default parameters.

Landscape of the aberrantly methylated CpGs between normal and ESCC

The flowchart of this study was presented in Figure 1A. The following analyses were performed to survey the landscape of ESCC methylome. Firstly, differentially methylated CpG sites between normal and tumor samples were identified using WGBS data. Due to the highly coordinated CpG sites being often tightly coupled with each other, we further identified differentially methylated regions (DMR) by a modified sliding window method. Differentially methylated genes (DMG) were then identified based on the genomic coordinates of DMRs and genes (see the method). Here we defined DMR or DMG as NC-DMR and NC-DMG if their methylation levels were higher on normal controls than tumor samples and vice versa as ESCC-DMR and ESCC-DMG. GO and KEGG pathway enrichment analyses were conducted to explore the potential functions of these DMGs. We also performed gene family analysis to identify significantly enriched gene families. The HOXL subclass homeobox family was selected to validate on two independent datasets and by Sanger sequencing in our custom samples. We further investigated the potential utility of the family genes as markers for cancer detection.

• Download figure
• Open in new tab

Figure 1. Methylation profiling of ESCC and paired normal samples.

A: An overview of this study. TCGA pan-cancer included the 9 most common cancer types: esophageal carcinoma (ESCA), esophageal squamous cell carcinoma (ESCC), stomach adenocarcinoma (STAD), colorectal cancer (CRC), liver hepatocellular carcinoma (LIHC), pancreatic adenocarcinoma (PAAD), breast invasive carcinoma (BRCA), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC). B: The identified DMCs across the whole genome. The red and gray points represented the ESCC-DMCs and NC-DMCs respectively. Chromosome X and Y were excluded in Figure 1B. WGBS: whole-genome bisulfite sequencing, DMC: differentially methylated CpG, DMR: differentially methylated region, DMG: differentially methylated gene.

We identified 438558 DMCs, including 361341 NC-DMCs and 69278 ESCC-DMCs (Figure 1B, supplementary table 5). Chromosome 8 has the most NC-DMCs, while the most ESCC-DMCs were in chromosomes 1 and 2 (Supplementary figure 1A). The distance between ESCC-DMCs was smaller than that of NC-DMCs (Supplementary figure 1B). Further investigations indicated that the methylation levels between adjacent DMCs were strongly correlated, and the correlation coefficients of adjacent ESCC-DMCs were much higher than NC-DMCs (0.73 vs. 0.70, P < 0.05, Supplementary figure 1C).

Identification of DMGs

We obtained 6422 NC-DMRs and 9040 ESCC-DMRs based on the DMCs. Similarly, most NC-DMRs were found in chromosome 8, while the majority of ESCC-DMRs were in chromosome 1 and 2 (Supplementary figure 2A). A smaller distance was also observed for adjacent ESCC-DMRs than that of NC-DMRs (Supplementary figure 2B). The average DMC count of ESCC-DMRs was higher than that of NC-DMRs (Figure 2A), and ESCC-DMRs tended to be located in inner-genic regions, whereas NC-DMRs more often located in intergenic regions (Figure 2B). According to the genomic coordinates of DMRs and genes, we identified 733 NC-DMGs and 906 ESCC-DMGs, of which 71 genes were shared by both (Supplementary table 6). Both NC- and ESCC-DMGs were mostly distributed in chromosome 19, which was different from DMCs and DMRs. The methylation values were able to separate ESCC-DMGs and NC-DMGs clearly, and ESCC-DMGs showed more concentrated than NC-DMGs (Figure 2C, 1st-3rd quantile: [29.32-43.64] vs. [13.86-34.14]). Functional enrichment analysis revealed that neuroactive ligand-receptor interaction was the most significantly enriched pathway for NC-DMGs (Figure 2D). For ESCC-DMGs, calcium signaling pathway was the most significantly enriched pathway (Figure 2E). The term of CD molecules was the major enriched gene family for NC-DMGs, whereas NKL subclass (Figure 2F, Supplementary table 7), HOXL subclass, and Zinc fingers C2H2-type were the 3 top-ranked enriched families for ESCC-DMGs (Figure 2G, Supplementary table 8).

• Download figure
• Open in new tab

Figure 2 Functional enrichment analysis of differentially methylated genes.

A: The number of DMCs in ESCC-DMRs and NC-DMRs. B: The genome distribution of ESCC-DMRs and NC-DMRs. C: The methylation levels of DMGs in normal and ESCC samples. D&E: The enriched KEGG pathways of NC-DMGs (D) and ESCC-DMGs (E). F&G: The top 10 enriched gene families of NC-DMGs (F) and ESCC-DMGs (D). The point size in Figure 2D-F represented the mapped gene numbers, and color indicated the enriched p value.

The aberrant methylation of HOXL subclass homeoboxes

We selected the second-ranked HOXL subclass homeobox family for further analysis, as more than half (n=24, 66.67%) of the genes were identified as ESCC-DMGs. K-means clustering showed that the 36 genes could be clustered into 3 groups (Figure 3A). Group 1 consisted of 11 genes that were hypermethylated on both normal and ESCC samples. Twenty-four genes, all of which were ESCC-DMGs, were in group 2. Furthermore, group 2 genes were divided into three sub-groups, subgroup 1 (n=10), subgroup 2 (n=12), and subgroup 3 (n=2). Only one gene, HOXA7, was in group 3 and showed low methylation levels on normal and ESCC samples. Overall, genes of subgroup 2 showed the lowest methylation levels on normal samples than sub group 1 and sub group 3 (Figure 3B). Correlation analysis revealed significant positive correlations between group 2 genes, except for MNX1 and GBX1 (Figure 3C).

• Download figure
• Open in new tab

Figure 3 Methylation patterns of HOXL subclass homeoboxes.

A: The methylation levels of HOXL subclass homeoboxes genes in normal and ESCC samples. B: The methylation levels of subgroup 2 genes between normal and ESCC samples. C: The methylation correlation of HOXL subclass homeobox genes. The correlation coefficients were estimated using Pearson’s method. The β values of each gene in Figure 3A-C were calculated using TCGA ESCC dataset.

Validation of sub-group 2 genes of HOXL subclass homeobox

Using GSE52826 dataset, we verified the methylation characteristics of subgroup 2 genes (n=10) in ESCC and paired normal samples (Figure 4A). After the exclusion of outlier samples (GSM1276746 and matched GSM1276750) (Supplementary figure 3), six genes were validated hypermethylated in tumor samples (Figure 4B). Since HOXC12 showed high methylation levels in normal samples, we further investigated methylation patterns of the rest five genes on TCGA ESCC dataset (Figure 4C). The methylation of HOXC11, HOXC10, and HOXD1 in normal samples were lower than GSX1 and CDX2 (Figure 4D), suggesting their low methylation background. ROC analysis indicated that methylation values of HOXC10 and HOXD1 showed the best performance in discriminating ESCC from normal samples, with both AUC reached 0.85 (Figure 4E). Further investigations revealed that four DMRs located in HOXC10, and HOXC10-R3 presented the largest delta methylation value (Figure 4F). HOXD1 contained seven DMRs, with HOXD1-R4 showed the largest delta methylation value (Figure 4G). When stratified TCGA ESCCs by different pathological stages, both genes showed hypermethylated across all stages, including stages IA and IB (supplementary figure 4).

• Download figure
• Open in new tab

Figure 4 Methylation patterns of sub-group 2 genes in validation sets.

A: The methylation levels of sub-group 2 genes in normal and ESCC samples in GSE52826 dataset. B: The differentially methylated six genes between ESCC and paired normal samples. Paired student t-test was used to estimate significant p values. C: The methylation levels of five genes in normal and ESCC samples in TCGA-ESCC. D: The differentially methylated four genes between ESCC and normal samples. The Wilcoxon rank sum test was used to estimate significant p values. E: ROC curves of the four genes in classifying TCGA-ESCC and normal samples. F-G: The identified ESCC-DMRs of HOXC10 and HOXD1 represented by R1-R4 and R1-R7 respectively.

The performance of HOX-score for ESCC classification

We used logistic regression to develop an ESCC classification model based on HOXC10/HOXD1 methylation using the TCGA ESCC cohort. The risk score, defined as HOXL-score, was then estimated for each sample. The normal samples showed the lowest risk scores than ESCC and ESCA (Figure 5A, median: 0.28, 0.88, and 0.99). According to the mechanism of logistic regression, the HOXL score represented the probability that a sample was classified as cancer. Therefore, we attempted to classify ESCC and normal samples using HOXL scores. ROC curve analysis indicated that the AUC reached 0.96 (95% CI: 0.91-0.99) for ESCC (Figure 5B). The optimal threshold determined by the Youden index was 0.72, with a sensitivity of 94.8% and specificity of 87.5%. Besides, for ESCA, the AUC was 0.83 (95% CI: 0.72-0.93), with sensitivity and specificity of 83.1% and 87.5% at the optimal threshold (Figure 5B). No significant variations were observed for HOXL-score in detecting both ESCC and ESCA stratified by gender, age, and stage at the threshold of 0.72 (Table 3).

• Download figure
• Open in new tab

Figure 5 The performance of HOXL score for ESCC classification.

A: The HOXL scores of normal, ESCC and ESCA samples. For two group comparisons, the Wilcoxon rank sum test was used. For three group comparison, the Kruskal test was used. The three lines from top to bottom in each category indicated the median and 90th percentile of scores. B: ROC curves of HOXL-score for ESCC and ESCA classification. The points indicated the best cutoff of HOXL-scores, and the percentages were the best specificity and sensitivity estimated by Youden’s index.

The diagnostic performance of 13 potential methylation markers, CDKN2A, CDKN2B, TFF1, MGMT, MLH1, DAPK1, SCGB3A1, TFPI2, DACH1, SOX17, CHFR, CDH1, APC that have been reported in ESCC were also evaluated using the same approach in TCGA ESCC dataset (supplementary figure 5). We observed the highest AUC value (0.96 [95%CI: 0.91-0.98]) for DAPK1 methylation, which is comparable to the HOXL-score. The AUC values of CDKN2A and TFF1 were 0.85, which equaled HOXD1 but lower than HOXC10. These results suggested that the two genes and their combination showed great potential for ESCC detection.

The methylation patterns of HOXL subclass homeoboxes in pan-cancers

We observed extensive hypermethylation of HOXL subclass homeobox genes on ESCC and then investigated their methylation patterns on nine other cancer types, including STAD (n=393), COAD/READ (n=379), LIHC (n=374), PAAD (n=183), PRAD (n=495), BRCA (n=778), LUAD (n=456) and LUSC (n=364). Forty (83.33%) and eight genes (16.67%) were found hypermethylated on cancer and normal samples, respectively. HOXD10 showed the highest delta β between normal and tumor samples, with widespread hypermethylated on most cancer types, except for LIHC (Figure 6A). The similarity was observed for HOXD12, HOXD1, and HOXC5 (Supplementary table 9). GSX1 had the smallest standard error (sd) of β values on both cancer and normal samples, indicating its good uniformity across the ten cancer types (Figure 6B). Therefore, we further evaluated the potential utility of GSX1 methylation in detecting pan-cancers. Pan-cancer samples were randomly divided into training and testing sets by 1:1. AUC was 0.89 (95% CI: 0.87-0.91) for training set with a sensitivity and specificity of 84.1% and 85.4% at the optimal threshold of 0.88. For the testing set, the AUC was 0.88 (95% CI: 0.86-0.90), and the sensitivity and specificity were 84.2% and 82.7%, respectively (Figure 6C). In addition, HOXC10 showed the highest methylation levels on ESCC than other cancer types and normal samples (Figure 6D), suggesting that it could be used as a marker specifical for ESCC.

• Download figure
• Open in new tab

Figure 6 The potential of HOXL subclass homeoboxes for pan-cancers detection.

A: The methylation levels of HOXL subclass homeoboxes genes in normal samples and 9 cancer types. The gene methylation level of each category was represented by the average β values of all samples in each cancer. B: The β values of GSX1 in normal samples and multiple cancer types. C: ROC curves showed the classification performance of GSX1 methylation in training and test sets. The red points indicated the optimal cutoff of probabilities estimated by logistic regression. D: The β values of HOXC10 in normal samples and multiple cancer types. For panel B and D, the red lines represented the median and 90th percentile β values.

Validation of the methylation status of HOXD1 and GSX1 by sequencing

Among all collected samples, amplification products of HOXD1 were obtained from 4 normal samples, 19 ESCC tissues and 19 blood samples. For GSX1, products were obtained from 14 normal samples, 19 cancerous tissues and 19 blood samples. The results of Sanger sequencing indicated that methylation events of the target regions on both genes occurred more frequently on ESCC samples (Table 4). For all CpG sites, the methylation frequency on ESCC samples was significantly higher than that on normal samples (Supplementary table 10). Sanger sequencing was also conducted for 13 cfDNA samples collected from ESCC patients’ plasma. In the 13 cfDNA samples, we observed that 90.77% (118/130) of HOXD1 CpG sites were methylated, and 92.31% (108/117) of GSX1 CpG sites were methylated (Figure 7A&B, Table 4). The methylation ratios of cfDNA samples were a little lower than tissue samples (99.47% for HOXD1 and 99.42% for GSX1), which might be attributed to the low cfDNA amount in plasma resulting in failed detection. Besides, we observed lower methylation levels of HOXD1 than GSX1 on normal samples in TCGA dataset, which was also revealed by the sequencing results that fewer methylation events occurred on normal samples for HOXD1 than GSX1 (Figure 7A&B).

• Download figure
• Open in new tab

Figure 7 The methylation status of HOXD1 (A) and GSX1 (B) between normal, ESCC-tissue and -plasma samples.

Normal and ESCC samples were indicated by green and red side bars. BN, EN, EA, and cfDNA represented the blood samples, adjacent normal samples, ESCC-tissue and plasma samples, respectively. The methylated and unmethylated status of each CpG site were presented by black and white color in the heatmap. P and N indicated the positive and negative controls.