DNA methylation and gene expression profiling highlight retinol metabolism's vital role in Hepatocellular carcinoma development

doi:10.21203/rs.3.rs-3105819/v1

Download PDF

Research Article

DNA methylation and gene expression profiling highlight retinol metabolism's vital role in Hepatocellular carcinoma development

https://doi.org/10.21203/rs.3.rs-3105819/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Purpose

Abnormal DNA methylation patterns play a critical role in the development of hepatocellular carcinoma (HCC). However, the molecular mechanisms associated with these aberrantly methylated genes remain unclear. In this study, we aimed to comprehensively investigate the methylation-driven deregulated gene expressions in HCC using a large cohort of patients with diverse clinical characteristics.

Methods

Whole-genome bisulfite sequencing (WGBS) and RNA sequencing techniques were used to assess the methylation and gene expression profiles of HCC tissue and normal adjacent tissues (NATs). The potential function of candidate genes was then investigated using single cell RNA-seq (sc-RNA seq) data.

Results

We identified 132,773 differentially methylated regions (DMRs) and 4,322 differentially expressed genes (DEGs) between HCC and NATs. Integrated DNA methylation and RNA-seq data obtained 987 methylation-driven candidates, including 970 upregulated and 17 downregulated genes. Four retinol metabolic pathway genes, ADH1A, CYP2A6, CYP2C8, and CYP2C19, were found as hyper-downregulated genes, whose expressions could stratify HCC into three subgroups with distinct survival outcomes, immune cell infiltration, and tumor microenvironments. Validation of these findings in an independent dataset confirmed the high concordance and potential prognostic value of these genes. Sc-RNA seq data revealed the low expression of these genes in immune cells and highlighted their role in promoting malignant cell proliferation and migration.

Conclusion

This study provides an insight into the molecular characteristics of HCC, unraveling the involvement of retinol metabolism-related genes in HCC development and progression. These findings have implications for HCC diagnosis, prognosis prediction, and the development of therapeutic targets.

Hepatocellular carcinoma

DNA methylation

gene expression

retinol metabolic related genes

DNA methylation is a common and critical epigenetic modification of the genome. Abnormal methylation patterns are tightly related to many diseases development, including Hepatocellular carcinoma (HCC) [1, 2]. Early research has demonstrated that advanced HCC exhibits increased methylation of CHFR and SYK, potential tumor suppressor genes, as well as decreased methylation of LINE-1 and satellite 2 repeat elements, which are associated with the progression of chronic hepatitis and cirrhosis to HCC [3]. The widespread application of high-throughput technology has identified several abnormal methylation events in genomic DNA that can predict a higher risk of HCC in individuals [4–6]. In addition, a subgroup named the CpG island methylator phenotype (CIMP) is defined in HCC based on 450k microarray methylation data and exhibited a worse prognosis than the non-CIMP group [7]. These findings suggest that DNA methylation plays a vital role in the initiation and progression of HCC.

However, the molecular mechanisms associated with these aberrantly methylated genes are rarely reported. Unlike SNVs that alter gene structure, methylation modifications occurring in the regulatory regions of target genes are generally believed to silence or downregulate gene expression. Integrating DNA methylation and RNA expression data provides an excellent opportunity to unravel the molecular mechanisms associated with aberrantly methylated genes. A few attempts have obtained valuable findings in this regard, such as methylation-driven deregulated expressions of genes related to retinol metabolism [8] and tyrosine and monocarboxylic acid metabolism [9]. While these findings are not sufficient, which is partly because some studies have focused on discovering diagnostic markers rather than gene function, and other studies are limited to small populations with specific characteristics.

Understanding the molecular pathogenesis of HCC is essential for developing efficient diagnostic and therapeutic approaches because the illness is varied and complex, with various underlying molecular pathways. Therefore, we designed this study to comprehensively investigate the methylation-driven deregulated expressions of genes in HCC using a large cohort of patients with diverse clinical characteristics and, meanwhile, to identify the potential biological pathways affected by these genes. Our findings present an overview of the molecular characteristics of HCC and may have practical implications for HCC diagnosis, prognosis prediction, and treatments.

Samples collection

This study collected 12 fresh HCC and adjacent normal tissues from the First Affiliated Hospital of Zhengzhou University. All samples were frozen and preserved in liquid nitrogen immediately after surgical resection. The clinical characteristics of the 12 patients are shown in supplementary table 3. All participants were obtained signed informed consent forms and were informed of the purpose of use. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University with approval number 2022-KY-0631-002 and had been performed in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki.

Data preparation

The public data of PRJNA762641 project comes from the study by Huang P et al. [9], consisting of whole genome bisulfite sequencing (WGBS) and RNA-seq data. We downloaded the raw data of 33 HCC samples from the Sequence Read Archive (SRA) website (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA762641). The information for HCC patients is listed in Supplementary table 1. The GSE70090 dataset is from Li X et al.'s study [10], which includes WGBS data for 28 normal and lung/liver cancer tissues (Supplementary table 2). Only the raw read files of eight HCC and matched adjacent normal tissues (NATs) in the GSE70090 dataset were used in this study.

The GSE151530 [11] dataset is single-cell RNA-seq data containing 46 samples from 37 patients with liver disease (25 HCC, 12 CC). We removed the samples of CC patients and obtained 50023 cells. Then, we excluded the samples with cells less than 200 and obtained 49094 cells, which came from 23 samples of 17 HCC patients.

Bisulfite sequencing and methylation calling

Genomic DNA from fresh tissue samples was extracted (Qiagen, Germany) and treated by bisulfite using the EZ DNA Methylation Gold Kit (Zymo Research), followed by purification to obtain the BS-DNA. The BS-DNA was then randomly interrupted, and added adapters to construct sequencing libraries. Subsequently, whole genome sequencing was performed on the libraries using the high-throughput sequencing platform DNBSEQ. Next, raw reads were filtered using the SOAPnuke [12] tool to remove low-quality reads, adapter contamination, and other artifacts, with the parameter of '-n 0.001 -l 20 -q 0.4 --adaMR 0.25 --ada_trim --polyX 50'. At last, clean reads were obtained, and their base quality was set to Phred + 33.

The first step in analyzing BS-seq data is to align the short sequence reads to the reference genome and then count the number of C-to-T conversions in the aligned reads at each CpG site. This study used the Bismark [13] tool to carry out methylation calling. The reference genome sequence (hg38) was first indexed using the Bismark genome preparation tool, bismark_genome_preparation. Then, clean reads were mapped to the reference genome using bowtie2. The Bismark's methylation extractor tool was used to extract methylation values of each CpG. Methylation values were calculated as the percentage of methylated cytosines at each CpG site across all reads covering that site.

Differential methylation analysis

Three datasets, PRJNA762641, GSE70090, and our custom 24 samples, were integrated into one dataset for differential methylation analysis. One normal sample (SRR2074679) in GSE70090 was excluded because the complete fastq file was unavailable. Therefore, we obtained WGBS data from 97 samples, including 48 NATs and 49 HCCs. Differentially methylated CpGs (DMC) were identified by comparing the methylation levels between NATs and HCCs using a rank-sum test followed by multiple testing corrections. Significant DMCs were defined if the false discovery rate (FDR) was less than 0.05. Abnormally methylated regions were determined using a circular binary segmentation algorithm (CBS) [14], which has been reported in previous studies for differentially methylated regions (DMR) identification [15, 16]. In brief, the average methylation value of each CpG site in HCC and NAT samples was calculated, respectively. Then, delta (△_value) value of each CpG is calculated using the formula: T _average – N _average, where T _average and N _average represent the average methylation value of the CpG in HCC and NAT samples. DMRs were determined using the CBS algorithm with a maximum size of less than 1 Mb and at least covering 3 CpGs.

DMRs were classified into hypermethylated DMRs (hyper-DMR) and hypomethylated DMRs (hypo-DMR) according to the direction of methylation events. Hyper-DMRs were defined as those with △_value ≥ 10, while hypo-DMRs were defined as those with △_value ≤ -15. These two thresholds were selected for the following reason. We observed that the △_values of DMRs exhibited multimodal distribution patterns, with three peaks around − 15, 0, and 10, respectively. The threshold of 10 represents the majority of DMRs that are hypermethylated on HCC (accordingly, hypomethylated in NATs). And the threshold of -15 represents most of DMRs that are hypermethylated in NATs (accordingly, hypomethylated in HCC). DMR locations were annotated as Upstream (within 2000bp of the gene TSS), Upstream-Body, Inner-genic, Body, Body-Downstream, Downstream (within 200bp downstream of the gene), and Intergenic according to the genomic annotation information (GRCh38 gtf).

Differential expression analysis

RNA-seq data of 33 HCCs and matched NATs from the PRJNA762641 project were used to identify differentially expressed genes (DEG) between normal and tumor samples. Transcripts per million (TPM) values for all genes in both HCC and NATs were calculated using RSEM (RNA-Seq by Expectation-Maximization). Rank-sum test was performed to determine significant DEGs (FDR < 0.05). Low-expressed genes with an average TPM < 10 across all samples were excluded. Methylation-driven DEGs (Methy-DEGs) were those hyper-DMGs with low expression values in HCC or hypo-DMGs with high expression values in HCC. The identified DEGs were further analyzed using gene ontology (GO) and pathway enrichment analyses to determine the biological processes and molecular pathways disturbed in HCC. It should be noted that hyper-DMGs with high expressions and hypo-DMGs with low expressions in HCC were not in the scope of this study.

Statistical analysis

Both data processing and analysis are implemented in R software (version 4.1.1). We used the Kaplan-Meier method to estimate the survival curves of multiple groups and performed log-rank tests to assess the differences. ROC curve analysis was implemented using the R package 'pROC', and the area under ROC curve (AUC) value was calculated to evaluate the model's performance. Continuous variables were compared using rank-sum tests or Kruskal tests, depending on the number of groups being compared and the data distribution. The category variables were compared using chi-square or Fisher's exact tests depending on the expected cell counts. A p-value or FDR < 0.05 was considered statistically significant for all analyses.

Differential methylation analysis based on WGBS data

After excluding one abnormal sample, QC results showed that the average sequencing depth of WGBS data for the 24 samples was 21.42 [minimum 15.85, maximum 26.2]. The average base quality of WGBS data was 31.68, and the average mapped quality was 34.40. These results reflect the high quality of WGBS data from 24 fresh tissue samples. We combined the PRJNA762641 and GSE70090 datasets with the 24 samples of this study as an integrated dataset. The MDS plot showed that the tumor and normal samples of the three datasets clustered clearly in two separate groups, indicating good consistency among the datasets (Fig. 1A). Overall, CpGs in tumor samples showed lower methylation levels than normal samples, as seen in the density plot (Fig. 1B).

We use the circular binary segmentation (CBS) algorithm to identify DMRs between HCC and NATs, as described in the methods. We obtained 132773 DMRs, including 12317 hypermethylated and 120456 hypomethylated regions in HCC (Fig. 1C). We found that the hyper-DMRs displayed a significantly smaller length than hypo-DMRs and non-DMRs (Fig. 1D), suggesting that CpG in hyper-DMR tends to accumulate in shorter regions (e.g., CpG islands). The top 100 DMRs with the most significant differences are shown in Supplementary table 5. Some of these DMRs, like the hypermethylation of the MGMT gene promoter region (Fig. 1E) and DAPK1 (Fig. 1F), have been linked to HCC. In addition, we investigated the methylation status of 27 genes reported hypermethylated in HCC by previous studies [8, 17] and found that ten were covered by these hyper-DMRs (Supplementary table 6).

Integration of DNA methylation and gene expression data

In addition to WGBS data, the PRJNA762641 project also included transcriptomic data from the same patient, allowing us to investigate the potential impact of aberrant DNA methylation on gene expression. We identified 4,322 DEGs, including 4662 upregulated and 264 downregulated genes in HCC (Fig. 2A). Interestingly, we found a high proportion of upregulated DEGs (48.58% vs. 2.75%) (Fig. 2B), suggesting extensive transcriptional activation events during HCC development, consistent with the frequent hypomethylation events observed in tumor tissues. The functional enrichment results showed that the upregulated DEGs in HCC were mainly enriched in biological processes related to RNA metabolism (including RNA splicing, ribonucleoprotein complex biogenesis, etc.) (Fig. 2C). In contrast, down-regulated DEGs were significantly enriched in entirely different biological processes, mainly in GO terms like small molecule and organic acid metabolism (Fig. 2D).

By integrating the results of DMRs with gene expression, we attempted to identify the genes whose expression was tightly correlated with methylation levels. DEGs were first linked with the nearest DMRs (within 2kb) and then grouped into four categories: hyper-upregulated, hyper-downregulated, hypo-upregulated, and hypo-downregulated. We focused on hyper-downregulated and hypo-upregulated genes based on the hypothesis that methylation is negatively regulated with expression, as mentioned in the methods. Consequently, we obtained 987 methylation-driven candidates, including 970 upregulated and 17 downregulated genes (Fig. 2E). The methylation-driven DEGs accounted for 20.1% of the total DEGs, suggesting a significant impact of DNA methylation on gene expression.

Further investigations revealed that hypo-upregulated genes are enriched in the biological processes of RNA biosynthesis and transport (Fig. 2F). In contrast, hyper-downregulated genes are mainly associated with the metabolism of retinol, cellular hormones, and olefinic acid (Fig. 2G). Of note, four of the 17 hyper-downregulated genes are retinol metabolic pathway genes, namely ADH1A, CYP2A6, CYP2C8, and CYP2C19. Correlation analysis further confirmed that methylation levels of CpGs within the DMRs of these four genes significantly correlate with gene expression negatively, while no significant correlation was observed for CpGs outside of the DMRs (Supplementary Fig. 1). In addition, we evaluated the abilities of four genes’ DMRs to discriminate HCC from NATs. AUC values for the four DMRs were 0.93 [0.88 − 0.97], 0.72 [0.59 − 0.84], 0.91 [0.86 − 0.97] and 0. 90 [0.81 − 0.95] (Supplementary Fig. 2), indicating their potential for HCC diagnosis.

HCC can be stratified by the genes of retinol metabolism

To validate these findings, we investigated the methylation and expression characteristics of the four genes in the TCGA-HCC dataset, which consists of 365 HCCs and 50 matched NATs. We obtained 12 probes associated with the four genes (Fig. 3A) which all showed high concordance with the results of WGBS, except for cg19645066 (CYP2A6). For example, two probes on ADH1A showed hypomethylated in HCC, which was in line with the results of WGBS (Supplementary Fig. 3A). Similar results were observed for the probes of the other three genes (Supplementary Fig. 3B-D). However, significantly different methylation levels were only found for six probes between HCC and NATs (FDR < 0.05, Supplementary table 7), indicating the shortcomings of 450k microarray technology in determining the methylation patterns of these genes.

We then investigated the expression profiles of these four genes in the TCGA HCC dataset. Three of them had expression values, showing significant downregulation on tumor samples (Fig. 3B), consistent with our findings based on 33 tumor samples and matched NATs. Interestingly, we found that the expression levels of the three genes were not uniform across the HCCs, with some HCCs exhibiting higher expression levels and others significantly lower expression levels (Fig. 3B), suggesting a high heterogeneity of the three genes expression in HCC. Unsupervised clustering of the HCC samples using K-means method indicated that they could be clustered into three groups (Fig. 3C). Subgroup 1 (C1, n = 214) displayed the highest expression of the three genes, while subgroup 2 (C2, n = 118) and subgroup 3 (C3, n = 33) showed lower expression levels (Supplementary Fig. 4A). UMAP analysis also showed similar results: HCC could be stratified into distinct subtypes by the three genes (Fig. 3D).

Additionally, survival analysis demonstrated that patients in C1 had better overall survival rates than those in C2 and C3 (p < 0.05, Fig. 3E). To exclude bias in a single dataset, we verified these results in an independent liver cancer cohort from the ICGC-JP project. Using the same method, we obtained a remarkable agreement between the two cohorts for the K-means results (Fig. 3F) and the UMAP analysis (Fig. 3G). Moreover, the prognosis for C1, C2, and C3 subgroups also showed significant differences in the ICGC’s LIRC cohort, and the trend was consistent with the TCGA HCC cohort (Fig. 3H). These findings suggest that the genes related to retinol metabolism can be used to stratify hepatocellular carcinoma into different subtypes, with potential prognostic value for patient outcomes assessment.

Clinical features and genomic characteristics of the three subgroups

To better understand the differences between the three subgroups, we calculated the average expression values of the three genes as the retinol scores of HCC patients. We observed significant differences in the distribution of several clinical features between the three subgroups in TCGA-LIHC cohort (Supplementary Fig. 4B-H). Overall, retinol scores were lower in patients with advanced, high-grade HCC. For instance, patients in T1 stage showed higher scores than those in T2 and T3 (Supplementary Fig. 4B). However, the degree of lymph node metastasis showed an opposite trend with patients’ retinol scores (Supplementary Fig. 4C). While no significant difference in retinol scores was found between M0 and M1 stage patients, which may be related to the small number of cases in M1 stage (n = 3, Supplementary Fig. 4D). Interestingly, retinol scores showed a decreasing trend over the AJCC stage (from I to II and III) and grade (from G1 to G2, G3, and G4) (Supplementary Fig. 4E-F). Patients with HBV or HCV infections did not display significantly different retinol scores (Supplementary Fig. 4G). It should be noted that younger patients (age < 50) showed significantly lower retinol scores than older patients (age > = 50). This may be attributed to the higher malignancy of younger people in the TCGA HCC cohort because we found that patients under 50 were more likely to be AJCC stage III/IV, T3/T4, and G3/G4 (data not supplied).

We then investigate the mutation profiles of the three subgroups (Fig. 4A). The most commonly mutated genes in HCC patients, such as TP53, CTNNB1, TTN, MUC16, and PCL, were shared across the three groups (Supplementary Fig. 5A-C). Tumor mutation load also did not show significant differences between the three groups (Supplementary Fig. 5D). Nevertheless, mutational enrichment analysis showed that MDN1, NLGN4X, and COL11A1 were significantly enriched on C1, while CTNNB1, PREX2, and AXIN1 were significantly enriched on C2, and OBSCN, BAP1, and TSC2 were significantly enriched on C3 (Supplementary Fig. 5E). In addition, we found that the three groups shared the mutation signature of SBS22 (exposure to aristolochic acid) (Supplementary Fig. 6). These results suggest that three subgroups shared similar mutation patterns.

Next, we explored the tumor microenvironment (TME) profiles of the three subgroups. Cell proliferation and wound healing scores gradually increased from C1 to C3, suggesting that the subgroups have distinct tumor cell aggressiveness and blood vessel formation (Fig. 4B-C). We found that IFN-γ and TGF-β response scores exhibited opposite trends among the three groups (Fig. 4D-E), with C3 having the lowest IFN-γ response scores but the highest TGF-β response scores, indicating a lower immune response in C3 than C1. In addition, we investigated the 24 immune cell scores between the three groups (Supplementary Fig. 7). The top 4 significantly altered immune cells are macrophage M0, memory B cells, resting mast cells, and Th17 cells. We observed the highest scores of macrophage M0, memory B cells, and resting mast cells in C3 (Fig. 4F–H) and the highest score of Th17 cells in C1 (Fig. 4I). These results highlight a suppressed immune response in C3 compared to C1 and C2.

Integration of single-cell RNA-seq data revealed potential functions of retinol metabolism-related genes

In order to elucidate the biological processes in which these genes may be involved, we integrated single-cell RNA-seq data from HCC patients and investigated the expression patterns of the three retinol metabolism-related genes across different cell types. The scRNA-seq dataset, GSE151530, was used in this study, which contained 49432 cells from 23 tumor tissues after quality control. These cells were classified into six categories, namely T cells, B cells, cancer-associated fibroblasts (CAFs), tumor-associated macrophages (TAMs), tumor-associated endothelial cells (TECs), and epithelial cells, according to the annotation of Gao et al. (Fig. 5A). We calculated the retinol scores (average expression) for all cells using the expression values of the three genes (CYP2C19 was excluded due to lack of expression). Malignant cells have the highest scores and exhibit huge variation, while retinol scores are pretty low for other cells (almost equal 0) (Fig. 5B). Therefore, we focused on malignant cells and re-clustered them. Malignant cells clustered into 15 subpopulations (Fig. 5C), with the highest scores in subpopulation 3 and the lowest scores in subpopulation 5, 14, 9, and 12 (Fig. 5D). Cell types were annotated using Human Primary Cell Atlas Data, and most of the cells in subpopulation 5, 14, 9, and 12 were annotated as hepatocyte cells (Supplementary Fig. 8). These findings suggest that a subpopulation of cells, differentiated from malignant cells, exhibits significantly varied retinol scores. We then identified marker genes for the 15 subpopulations and conducted a GO-enrichment analysis. The results showed that marker genes of subpopulation 14, 9, and 12 were predominantly involved in immune system (Fig. 5E). Moreover, we found marker genes of subpopulation 5 significantly enriched in wound healing, endothelium/vasculature development, and cell migration (Fig. 5E). In addition, we compared the DEGs between TCGA HCC-C1 (highest retinol scores) and HCC-C3 (lowest retinol scores) with marker genes of subpopulation 5 (Supplementary Fig. 9A). There were 91 genes shared by both, which significantly enriched in similar biological processes such as cell-substrate adhesion, regulation of actin cytoskeleton, and wound healing (Supplementary Fig. 9B).

This study comprehensively assessed the methylation and gene expression profiles of human HCC tissue and non-cancer liver tissue through WGBS and RNA-seq techniques. Integrating the methylation patterns of regulated regions with gene expression allowed us to identify a panel of retinol metabolism-related genes whose expression showed negative correlations with their methylation levels. Based on the expression profiles of these genes, HCC can be divided into three subgroups that showed significant variations in survival, immune cell infiltration, and tumor microenvironment. In addition, results of single-cell data identified a subpopulation of cells with low expression of retinol metabolism-related genes, characterized by highly active proliferation and migration behavior, further confirming the functional implications of these genes in malignant cells.

Our findings reveal extensive hypo-methylated events in HCC compared to normal tissue, evidenced by the bimodal distribution of CpG methylation values in HCC, and approximately 90% of DMRs were determined to be hypo-DMR. In general, genome-wide hypomethylation, as well as site-specific hypermethylation, is one of the characteristics of cancer cells [2], which is consistent with the findings of this study. However, we observed a pretty low proportion of upregulated DEGs in HCC, implying more gene expression turn-on and activation of oncogenic pathways than inactivation of suppressors during the development of HCC. Hypomethylation events contribute to the activation of gene expression [18] and chromosome instability [19, 20] that leads to cancer development. Recent studies have also highlighted the importance of cancer-associated CpG island hypomethylation in driving the overexpression of oncogenic factors.

We obtained 17 potential methylation-driven DEGs, four of which are retinol metabolism-related, whose regulatory regions are hypermethylated with significantly downregulated expressions in HCC. Correlation analysis also indicated that the methylation levels of those CpGs located in regulatory regions were negatively correlated with gene expression. These findings suggest that methylation events may play a pivotal role in inhibiting these genes expression. In the TCGA HCC dataset, the methylation profiles of all probes associated with the four genes showed high agreement with the results of WGBS. However, we did not observe significant negative correlations between the methylation values and gene expression. Further investigations showed that these probes were not located in DMRs, providing a possible explanation for this result. Besides, very limited probes were available for these genes, suggesting the weaknesses of methylation array technique in some critical genomic regions.

This study also found that HCC could be stratified into three subgroups based on the four genes' expression. HCC patients with low expression values of the four genes had a significantly worse prognosis and showed increased proliferation and wound healing activity, suggesting that abnormal expression of these genes may affect HCC prognosis by altering cell proliferation and expansion. Moreover, we also found that HCC-C3, which had the worst prognosis, exhibited immunosuppressive characteristics compared to HCC-C1, which had the best prognosis, indicating that their prognostic impact also involves remodeling the immune system. Interestingly, single-cell RNA seq data also revealed extremely low expression levels of these genes in immune cells and high expression in malignant cells. Overall, these findings suggest potential molecular mechanisms for the association of retinol metabolism-related genes with the development and progression of HCC.

Several studies have demonstrated that abnormal expression of retinol metabolism-related genes can be used as diagnostic and prognostic markers for HCC [21–23]. In addition, the blood retinol level has been used for HCC diagnosis and prognosis [24]. Researchers found that retinol metabolism can distinguish HCC from non-tumor tissues and is an independent prognostic predictor for HCC, with low expression associated with poor survival time. The current study also revealed their downregulation in HCC, and the low expression indicates a worse prognosis. Moreover, we observed hypermethylation events in the regulatory regions of these genes, and their methylation levels were significantly negatively correlated with expression, suggesting that methylation events may play a critical role in driving altered gene expression. Meanwhile, we found that these DMRs are able to discriminate HCC from NATs, highlighting their diagnostic potential for HCC.

ScRNA-seq results reveal that these four genes are only expressed in malignant cells with high variation. In contrast, they are rarely expressed in other cell types (T cells, B cells, CAFs, TAMs, TECs and epithelial cells). Further investigations showed the lowest expressions in a subpopulation of malignant cells characterized by activated biological processes of wound healing, endothelium/vasculature development, and cell migration, which was then validated in the TCGA HCC dataset. A recent study reported that the overexpression of CYP2C8 inhibits the proliferation, clonality, migration, invasion, and cell cycle of HCC cells through the PI3K/Akt/p27Kip1 axis [25]. Meanwhile, in vitro experiments also revealed that a synthetic retinoid, ACR, acts as an agonist for RARs and RXRs, inducing apoptosis and inhibiting cell proliferation in human hepatocellular carcinoma cell lines [26]. Interestingly, our findings also suggest that these genes may act as good prognostic factors by inhibiting biological processes such as cell proliferation and migration.

The current study has some drawbacks. First, we only focused on the four hypermethylated genes related to retinol metabolism. Other genes identified in this study may also play important roles in HCC development and progression, and more research is still needed to unveil their potential functions. Second, the findings of this study show that the four gene methylation can distinguish HCC from NATs well, indicating their diagnostic value. However, the current evidence is insufficient to support their good performance in clinical practice, which needs to be evaluated in a larger population. Third, gene expression is regulated in a complex way. In addition to abnormal epigenetic modification, other factors, such as miRNA and copy number variation, also affect gene expression. This study found that the methylation values of four gene DMRs showed significant negative correlations with their expression. However, no direct evidence suggested that the downregulated expression was due to hyper-methylation events. More molecular and cellular experiments need to be conducted to investigate their relationships in the future.

We observed extensive demethylation events in HCC as well as activation of gene expression. Four genes related to retinol metabolism are able to classify HCC into three subgroups, showing different prognoses and immune microenvironments. Integration of single-cell data allowed us to identify a class of cells that exhibit low retinol scores and are characterized by active wound healing and cell proliferation. Overall, our study provides new insight into HCC development and progression from a multidimensional perspective. Moreover, the current study obtains a panel of aberrant methylation and expression candidates that might serve as potential methylation markers, therapeutic targets, and prognostic markers. However, further research is needed to uncover the functions of other identified genes and evaluate the diagnostic and prognostic performance of the identified biomarkers in larger populations.

Acknowledgments

We are grateful to Dr. Xu at Zhongnan Hospital of Wuhan University for his help. All authors listed have read the complete manuscript and have approved submission of the paper.

Consent to publish

All participants were obtained signed informed consent forms and were informed of the purpose of use. Their personal identifications were removed in this study.

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Funding

This study was supported by science and technology project of Henan Province (232102310028).

Authors' Contributions

Study design: SC Wang, YT Zhao and KK Wan

Study supervision: QK Yang

Sample collection and preparation: J Wang, SY Wang, DH Zhou and W Zhang

Resource support: LL Zhang and LL Dong

Data analysis and interpretation (e.g., statistical analysis, computational analysis): YT Zhao, KK Wan, and J Wang

Data curation and supervision: ZY Du, YL Chang

Manuscript writing: YT Zhao and KK Wan

Manuscript review, and revision: SC Wang and QK Yang

Ethics statements

This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University with approval number 2022-KY-0631-002.

Data availability

The raw fastq data of 24 samples (12 HCC and 12 NATs) were submitted to the Sequence Read Archive database (PRJNA984754). The methylation data and associated R code used in this study are available on GitHub (https://github.com/amsinfor/HCC-WGBS-analysis).

Schübeler D. Function and information content of dna methylation. Nature. 2015; 517(7534): 321–6.doi: 10.1038/nature14192.
Nishiyama A, Nakanishi M. Navigating the dna methylation landscape of cancer. Trends Genet. 2021; 37(11): 1012-27.doi: 10.1016/j.tig.2021.05.002.
Hattori N, Ushijima T. Epigenetic impact of infection on carcinogenesis: mechanisms and applications. Genome Med. 2016; 8(1): 10.doi: 10.1186/s13073-016-0267-2.
Gonçalves E, Gonçalves-Reis M, Pereira-Leal JB, Cardoso J. Dna methylation fingerprint of hepatocellular carcinoma from tissue and liquid biopsies. Sci Rep-Uk. 2022; 12(1): 11512.doi: 10.1038/s41598-022-15058-0.
Bai Y, Tong W, Xie F, Zhu L, Wu H, Shi Ret al. Dna methylation biomarkers for diagnosis of primary liver cancer and distinguishing hepatocellular carcinoma from intrahepatic cholangiocarcinoma. Aging (Albany NY). 2021; 13(13): 17592–606.doi: 10.18632/aging.203249.
Li R, Shui L, Jia J, Wu C. Construction and validation of novel diagnostic and prognostic dna methylation signatures for hepatocellular carcinoma. Front Genet. 2020; 11: 906.doi: 10.3389/fgene.2020.00906.
Cheng J, Wei D, Ji Y, Chen L, Yang L, Li Get al. Integrative analysis of dna methylation and gene expression reveals hepatocellular carcinoma-specific diagnostic biomarkers. Genome Med. 2018; 10(1): 42.doi: 10.1186/s13073-018-0548-z.
Udali S, Guarini P, Ruzzenente A, Ferrarini A, Guglielmi A, Lotto Vet al. Dna methylation and gene expression profiles show novel regulatory pathways in hepatocellular carcinoma. Clin Epigenetics. 2015; 7(1): 43.doi: 10.1186/s13148-015-0077-1.
Huang P, Xu M, Han H, Zhao X, Li MD, Yang Z. Integrative analysis of epigenome and transcriptome data reveals aberrantly methylated promoters and enhancers in hepatocellular carcinoma. Front Oncol. 2021; 11: 769390.doi: 10.3389/fonc.2021.769390.
Li X, Liu Y, Salz T, Hansen KD, Feinberg A. Whole-genome analysis of the methylome and hydroxymethylome in normal and malignant lung and liver. Genome Res. 2016; 26(12): 1730–41.doi: 10.1101/gr.211854.116.
Ma L, Wang L, Khatib SA, Chang CW, Heinrich S, Dominguez DAet al. Single-cell atlas of tumor cell evolution in response to therapy in hepatocellular carcinoma and intrahepatic cholangiocarcinoma. J Hepatol. 2021; 75(6): 1397 – 408.doi: 10.1016/j.jhep.2021.06.028.
Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li Set al. Soapnuke: a mapreduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience. 2018; 7(1): x120.doi: 10.1093/gigascience/gix120.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-seq applications. Bioinformatics. 2011; 27(11): 1571-2.doi: 10.1093/bioinformatics/btr167.
Olshen AB, Venkatraman ES, Lucito R, Wigler M. Circular binary segmentation for the analysis of array-based dna copy number data. Biostatistics. 2004; 5(4): 557 – 72.doi: 10.1093/biostatistics/kxh008.
Gong B, Purdom E. Methcp: differentially methylated region detection with change point models. J Comput Biol. 2020; 27(4): 458 – 71.doi: 10.1089/cmb.2019.0326.
Jühling F, Kretzmer H, Bernhart SH, Otto C, Stadler PF, Hoffmann S. Metilene: fast and sensitive calling of differentially methylated regions from bisulfite sequencing data. Genome Res. 2016; 26(2): 256 – 62.doi: 10.1101/gr.196394.115.
Zhang C, Li J, Huang T, Duan S, Dai D, Jiang Det al. Meta-analysis of dna methylation biomarkers in hepatocellular carcinoma. Oncotarget. 2016; 7(49): 81255–67.doi: 10.18632/oncotarget.13221.
Good CR, Panjarian S, Kelly AD, Madzo J, Patel B, Jelinek Jet al. Tet1-mediated hypomethylation activates oncogenic signaling in triple-negative breast cancer. Cancer Res. 2018; 78(15): 4126–37.doi: 10.1158/0008-5472.CAN-17-2082.
Kawano H, Saeki H, Kitao H, Tsuda Y, Otsu H, Ando Ket al. Chromosomal instability associated with global dna hypomethylation is associated with the initiation and progression of esophageal squamous cell carcinoma. Ann Surg Oncol. 2014; 21 Suppl 4: S696-702.doi: 10.1245/s10434-014-3818-z.
Sheaffer KL, Elliott EN, Kaestner KH. Dna hypomethylation contributes to genomic instability and intestinal cancer initiation. Cancer Prev Res. 2016; 9(7): 534 – 46.doi: 10.1158/1940-6207.CAPR-15-0349.
Zhao N, Zhang Y, Cheng R, Zhang D, Li F, Guo Yet al. Spatial maps of hepatocellular carcinoma transcriptomes highlight an unexplored landscape of heterogeneity and a novel gene signature for survival. Cancer Cell Int. 2022; 22(1): 57.doi: 10.1186/s12935-021-02430-9.
Awan FM, Naz A, Obaid A, Ali A, Ahmad J, Anjum Set al. Identification of circulating biomarker candidates for hepatocellular carcinoma (hcc): an integrated prioritization approach. Plos One. 2015; 10(9): e138913.doi: 10.1371/journal.pone.0138913.
Wang X, Yu T, Liao X, Yang C, Han C, Zhu Get al. The prognostic value of cyp2c subfamily genes in hepatocellular carcinoma. Cancer Med-Us. 2018; 7(4): 966 – 80.doi: 10.1002/cam4.1299.
Han J, Han ML, Xing H, Li ZL, Yuan DY, Wu Het al. Tissue and serum metabolomic phenotyping for diagnosis and prognosis of hepatocellular carcinoma. Int J Cancer. 2020; 146(6): 1741-53.doi: 10.1002/ijc.32599.
Zhou X, Li TM, Luo JZ, Lan CL, Wei ZL, Fu THet al. Cyp2c8 suppress proliferation, migration, invasion and sorafenib resistance of hepatocellular carcinoma via pi3k/akt/p27(kip1) axis. J Hepatocell Carcino. 2021; 8: 1323–38.doi: 10.2147/JHC.S335425.
Shirakami Y, Sakai H, Shimizu M. Retinoid roles in blocking hepatocellular carcinoma. Hepatobil Surg Nutr. 2015; 4(4): 222-8.doi: 10.3978/j.issn.2304-3881.2015.05.01.

No competing interests reported.

Supplementarymaterials.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

DNA methylation and gene expression profiling highlight retinol metabolism's vital role in Hepatocellular carcinoma development

Status:

Version 1

Abstract

Purpose

Methods

Results

Conclusion

Figures

Introduction

Materials and methods

Samples collection

Data preparation

Bisulfite sequencing and methylation calling

Differential methylation analysis

Differential expression analysis

Statistical analysis

Results

Differential methylation analysis based on WGBS data

Integration of DNA methylation and gene expression data

HCC can be stratified by the genes of retinol metabolism

Clinical features and genomic characteristics of the three subgroups

Integration of single-cell RNA-seq data revealed potential functions of retinol metabolism-related genes

Discussion

Conclusions

Declarations

Acknowledgments

Consent to publish

Competing Interests

Funding

Authors' Contributions

Ethics statements

Data availability

References

Additional Declarations

Supplementary Files

Status:

Version 1