Deciphering the Divergent Gene Expression Landscapes of m6A/m5C/m1A Methylation Regulators in Hepatocellular Carcinoma Through Single-Cell and Bulk RNA Transcriptomic Analysis

Introduction

Hepatocellular carcinoma (HCC) is the most common primary liver cancer worldwide, and is distinguished by its rapid development, high invasiveness, and high mortality rate.¹ Apart from well-explored epigenetic modifications such as DNA methylation, RNA modifications have emerged as significant contributors to the regulation of gene expression and cellular functions. These post-transcriptional modifications dynamically affect various aspects of RNA metabolism, including RNA processing, export, translation, and stability, and play crucial roles in cancer behaviors.^2–8 Research has demonstrated that the development of HCC is governed by intricate genetic and epigenetic mechanisms, in which abnormal RNA modifications play a pivotal role in causing variations in gene expression, thereby influencing the occurrence and progression of HCC.^9–15

Methylation involves a complex and dynamically reversible series of structural alterations induced by various enzymes. The discovery of N6-methyladenosine (m6A) has sparked interest in understanding its biological significance.^16,17 Research has revealed that m6A plays essential roles in various cellular processes, including RNA splicing and alternative splicing, RNA translation and stability, and regulation of gene expression.^18,19 Additionally, the presence of 5-methylcytidine (m5C) is markedly lower than that of m6A, with levels approximately three to ten times lower, accounting for only 0.03–0.1% of cytosines.^20,21 Furthermore, N1-methyladenosine (m1A) is 10 times less abundant than m6A.^21,22 Dysregulation of modification patterns of m5C^23–27 and m1A^28–30 in RNA molecules has been linked to cancer development and progression.

Several studies have demonstrated the importance of the m6A, m1A, and m5C regulatory genes in RNA modifications^31–38 and cancer progression.^30,38–40 To regulate methylation, writers modify DNA through various chemical changes, such as methyl group moieties. These molecular decorations can then attract several proteins known as “readers” that identify these moieties. A collection of “erasers” is used to remove the deposited change. Writers, readers, and erasers collaborate to shape the epigenetic landscape.⁴¹ A closely grouped subset of regulator genes for m6A/m5C/m1A methylation (hereafter referred to as m6A/m5C/m1A genes) was found to be associated with an unfavorable prognosis and immune microenvironment in HCC. These genes have the potential to serve as innovative tools for evaluating the prognosis of HCC patients³¹ and aid in the selection of appropriate immunotherapy strategies.³²

However, the results of numerous studies on m6A/m5C/m1A genes are based on the analysis of bulk sequencing of HCC samples and do not consider the complex and various compositions of cells in the tumor environment. Therefore, this study aimed to explore the expression of regulator genes for m6A/m5C/m1A methylation at the single-cell level and to validate their expression in cancerous and adjacent para-cancerous liver tissues of adult patients with HCC who underwent tumor resection.

Materials and Methods

The Processing of Bulk RNA-Seq Data from TCGA

This study collected a cohort of 371 HCC tissues and 50 adjacent para-cancerous liver tissues from TCGA database. The purpose of this data collection was to assess the expression levels of m6A/m5C/m1A genes in HCC. Based on previously published literature,^{5,30,32,37,38,42–44} 35 m6A, 16 m5C, and 13 m1A regulator genes were selected (Table S1). The FPKM data were first transformed into transcripts per million and normalized. The analysis of differentially expressed genes (DEGs) for the m6A/m5C/m1A genes between HCC and adjacent para-cancerous liver tissues in TCGA cohort was conducted using the “limma” R package (Table S2). The cutoff thresholds employed for identifying DEGs were |log2 fold change|>1 and a false discovery rate of <0.05.

Source of scRNA-Seq and Data Processing

The scRNA-seq data, obtained from the GEO database (GSE149614, GSE156625, and GSE125449) with 10× scRNA-seq data, was processed using the “Seurat” R package. This involved transforming the 10× scRNA-seq data into Seurat objects and filtering low-quality cells. Quality control was performed by calculating the percentage of mitochondrial or ribosomal genes. To identify highly variable genes, the “FindVariableFeatures” program was used to select the top 2000 genes, while an additional 2000 genes were employed for cell subpopulation identification through PCA. The “SingleR” R package was utilized for the annotation of different cell types. To mitigate any batch effects among between GSE samples and integrate them based on annotations, the canonical correlation analysis implemented as part of Seurat software package was applied. Cell clusters were visualized by utilizing a UMAP map for dimensionality reduction of the data. The cell cycle status was analyzed using the “CellCycleScoring” in the Seurat R package.

Differential Expression of scRNA

The analysis of DEGs for the m6A/m5C/m1A genes between HCC and adjacent para-cancerous liver tissues in the scRNA was conducted using the “FindMarkers” in the Seurat R package. The cut-off thresholds employed for identifying DEGs were |log2 fold change|>1 and an adjusted p<0.05.

Expression of m6A/m5C/m1A Genes in Collected Human Liver Samples

The study included patients with HCC who underwent tumor resection and were aged 20 years or older. Patients with any other cancer, immunocompromised diseases, or autoimmune diseases were excluded. Informed consent was obtained from all patients prior to sample collection. This study included 20 patients with HCC. Tissue specimens (0.5 × 0.5 cm) were collected from both HCC and adjacent para-cancerous liver tissues. Total RNA was extracted from each sample using an RNeasy Mini Kit (Qiagen, Venlo, Netherlands), and RNA quantity and quality were assessed using an SSP-3000 NanoDrop spectrophotometer (Infinigen Biotech, City of Industry, CA, USA). A real-time PCR was used to measure the expression of dysregulated m6A/m5C/m1A genes identified from the scRNA-seq and selected representative m6A/m5C/m1A genes. cDNA was amplified using Power SYBR Green PCR Master Mix with the specific primers listed in Table S3. The expression of methylation-related genes was quantified using the 2-Ct method with normalized cycle threshold values. All results were presented as the mean ± standard error. One-way analysis of variance was conducted for the overall analysis of the differences between group means, followed by a post-hoc Fisher’s least significant difference test. Statistical significance was set at p<0.05. Kaplan–Meier (K–M) survival analysis of significantly expressed methylation-related genes in the real-time PCR was performed using Timer2.0.⁴⁵

Correlation of Dysregulated m6A/m5C/m1A Genes with Immune Cell Infiltration

Timer2.0⁴⁵ was used to analyze the infiltration relationship between dysregulated m6A/m5C/m1A genes identified by scRNA-seq and specific types of immune cells. A correlation between the dysregulated genes and infiltrated immune cells was considered significant if the absolute value of the correlation coefficient (|Rho|) was >0.1 and p<0.05.

Results

Dysregulated m6A/m5C/m1A Genes Between HCC and Adjacent Para-Cancerous Liver Tissues in the the Cancer Genome Atlas (TCGA) Cohort

From the bulk sequencing of the Cancer Genome Atlas (TCGA) database, 15 upregulated and 6 downregulated m6A genes, 4 upregulated and 2 downregulated m1A genes, and 7 upregulated and 5 downregulated m5C genes were identified in HCC, compared to that in adjacent para-cancerous liver tissues (Figure 1). The clustered heatmap of the dysregulated m6A/m5C/m1A genes revealed their heterogeneous in the adjacent para-cancerous liver tissue compared to that in HCC (Figure 2).

Figure 1 Dysregulated m6A/m5C/m1A genes identified from the bulk sequencing of The Cancer Genome Atlas (TCGA) cohort. *Indicates a significant change (p<0.05) between the hepatocellular carcinoma (HCC) and the adjacent para-cancerous liver tissue.

Figure 2 The clustered heatmap of the dysregulated m6A/m5C/m1A genes in HCC and in adjacent para-cancerous liver tissue.

Single-Cell Transcriptome Profile of HCC and Adjacent Para-Cancerous Liver Tissue

After quality control filtering, 26,708 unique genes were identified in >73,916 cells. Principal Component Analysis (PCA) was conducted using 2000 variable genes to reduce dimensionality and to classify the cells into 20 clusters. After batch effect correction, Uniform Manifold Approximation and Projection (UMAP) analysis (Figure 3) was used to visualize the results based on HCC vs adjacent para-cancerous liver tissue, Gene Series Accession (GSE) categories in the Gene Expression Omnibus (GEO) database, and cell cycles. Unsupervised dimensionality reduction and graph-based clustering analyses were used to identify nine major cell clusters. The following cell types are described: T-cells, CD8+ T-cells, B-cells, endothelial cells, fibroblasts, hepatocytes, monocytes, macrophages, and NK cells. The association of dysregulated m6A/m5C/m1A genes identified from single-cell RNA sequencing (scRNA-seq) with different cell types is shown in Table 1 and the bubble (Figure 4A) and river (Figure 4B) plots in Figure 4. A summary of the expression of scRNA-seq-identified dysregulated m6A/m5C/m1A genes in bulk sequencing from TCGA cohort is provided in Figure 4C. The clustered heatmap of the dysregulated m6A/m5C/m1A genes in HCC and in adjacent para-cancerous liver tissue according to bulk sequencing from TCGA cohort revealed heterogeneous expression of the m6A/m5C/m1A genes within the cancer when compared to a more homogenous expression of these genes in the adjacent para-cancerous liver tissues.

Table 1 The Dysregulated m6A/m5C/m1A Genes Identified from the scRNA-Seq in Different Cell Types in the HCC vs the Adjacent Para-Cancerous Liver Tissue

Figure 3 Single-cell RNA-sequencing analysis identified based on the HCC vs adjacent para-cancerous liver tissue, GSE categories, cell cycles, and major clusters of the cells.

Figure 4 (A) The bubble plot conducted to visualize the distribution of m6A/m5C/m1A genes among different cell types. To evaluate the statistical significance of each interaction score, a one-sided Wilcoxon rank-sum test was employed. (B) The river plot shows the association of dysregulated m6A/m5C/m1A genes to different cell types. GSE, Gene Series Accession; UMAP, Uniform Manifold Approximation and Projection; G1, Gap 1; S, Synthesis of DNA; G2 M, Gap 2 and Mitosis. (C) Summary of the expression of scRNA-seq-identified dysregulated m6A/m5C/m1A genes in the bulk sequencing from TCGA cohort. (D) The clustered heatmap of these dysregulated m6A/m5C/m1A genes in HCC and in adjacent para-cancerous liver tissue according to the bulk sequencing from TCGA cohort. *Indicates a significant change (p<0.05) between the hepatocellular carcinoma (HCC) and the adjacent para-cancerous liver tissue.

The Expression of Methylation-Related Genes Determined Using Real-Time PCR

A real-time polymerase chain reaction (real-time PCR) was used to validate the expression of dysregulated m6A/m5C/m1A genes identified by scRNA-seq and revealed significant overexpression of DNMT1, NSUN5, TRMT6, IGF2BP1, and IGFBP3 (Figure 5A). Although YTHDC1 expression was lower in HCC than in the adjacent para-cancerous liver tissue, as shown by scRNA-seq, the difference was not significant. The upregulation of these genes was consistent with that observed by scRNA-seq. Quantification of the selected representative m6A/m5C/m1A genes from the writers, readers, and erasers categories demonstrated that there was no significant difference in the expression levels of these genes in HCC vs adjacent para-cancerous liver tissue, except that IGF2BP2 expression in HCC was considerably elevated compared to its expression in the adjacent para-cancerous liver tissue (Figure 5B). Kaplan–Meier survival analysis (Figure 5C) through the Timer2.0 website demonstrated that out of the six dysregulated genes studied, IGF2BP3 expression exhibited a significant correlation with the survival of the patients with HCC (hazard ratio [HR]= 1.34, p=0.000226).

Figure 5 (A) A real-time PCR was performed to measure the expression of dysregulated m6A/m5C/m1A genes identified from the scRNA-seq. (B) A real-time PCR was performed to measure the expression of dysregulated m6A/m5C/m1A genes identified from selected representative dysregulated m6A/m5C/m1A genes from TCGA. (C) Kaplan–Meier (K–M) survival analysis of six dysregulated methylation-related genes using Timer2.0. *Indicates a significant change (p<0.05) between the hepatocellular carcinoma (HCC) and the adjacent para-cancerous liver tissue.

Correlation of Dysregulated m6A/m5C/m1A Genes with Immune Cell Infiltration

The connection between dysregulated genes involved in the m6A/m5C/m1A processes, as discovered by scRNA-seq, and the infiltration of the corresponding immune cells was examined using Timer2.0.¹¹ The analysis revealed a significant correlation between the five dysregulated genes (ALKBH5, DNMT3A, PPRC2A, WTAP, and YBX1) and the presence of NK cells (Figure 6). The correlation between ALKBH5, DNMT3A, PPRC2A, and NK cell infiltration was negative, whereas the correlation between WTAP and YBX1 was positive. Nevertheless, no statistically significant correlation was found between the expression levels of IGF2BP2 and METTL3 and endothelial cell infiltration. Similarly, the expression levels of RBM15B, YBX1, and ALKBH5 were not significantly correlated with CD8+ T cell infiltration.

Figure 6 Correlation of dysregulated m6A/m5C/m1A genes with immune cell infiltration using Timer2.0. *Indicates a significant change (p<0.05) of the type of infiltrated cells between the hepatocellular carcinoma (HCC) and the adjacent para-cancerous liver tissue.

Discussion

Overall, our findings revealed a discrepancy in the expression of m6A/m5C/m1A methylation regulator genes between bulk sequencing and scRNA-seq. Cellular diversity inside tumors is a key factor explaining the discrepancy in assessments. The clustered heatmap of the scRNA-seq-identified dysregulated m6A/m5C/m1A genes in TCGA cohort revealed heterogeneous expression of these methylation regulators within the cancer, whereas their expression in the adjacent liver tissues was more homogeneous. The dysregulation of methylation genes from bulk sequencing may not reflect those changes inside the hepatocytes but were derived from other cells in the tumor environment, particularly considering that the expression of some regulator genes was correlated with the immune infiltration of certain immune cells. The expression of regulatory genes for m6A/m5C/m1A methylation at the single-cell level may adhere more to the real situation, considering the complex composition of cells and heterogeneity in the tumor environment deterring the evaluation.

This study found five regulator genes (DNMT1, NSUN5, TRMT6, IGF2BP1, and IGFBP3) to be significantly overexpressed in HCC compared to adjacent para-cancerous liver tissue, and one regular gene (IGF2BP) was found to be significantly overexpressed in the real-time PCR, but not in scRNA-seq. DNMT1 is primarily responsible for the preservation of DNA methylation patterns during replication. It detects hemi-methylated DNA, which occurs when one strand of DNA is methylated but the other is not, and methylates the unmethylated strand to ensure correct methylation inheritance. The upregulation of DNMT1 in HCC leads to aberrant patterns of DNA methylation, characterized by reduced methylation levels across the genome and increased methylation levels in specific genomic regions. Global hypomethylation can potentially induce genomic instability and activate transposable elements, whereas localized hypermethylation of gene promoters can silence tumor suppressor genes that play crucial roles in cell cycle regulation, DNA repair, and apoptosis. DNMT1 overexpression frequently results in the silencing of tumor suppressor genes and the activation of oncogenes, such as c-MET, CDH1, and TERT, thereby increasing tumor development and metastasis.

The main role of NSUN5 is in the post-transcriptional modification of RNA, specifically in the methylation of cytosine at the C5 position of the tRNA. Methylation is crucial for maintaining the stability and functionality of tRNA. tRNA methylation, which is mediated by NSUN5, plays a significant role in ensuring precise and efficient translation by facilitating accurate codon–anticodon interactions and ribosome activity. Abundant NSUN5 expression is observed in HCC.⁴⁶ In HCC, NSUN5 is associated with poor prognosis and contributes to disease progression.⁴⁷ However, the exact role of NSUN5 in the development and progression of HCC has not yet been thoroughly investigated.

TRMT6 catalyzes the 2’-O-methylation of the cytosine residue at position 34 within the anticodon loop of specific tRNAs. This modification is essential in ensuring accurate and efficient protein translation.²⁹ Elevated TRMT6 expression has been detected in HCC tissues and is correlated with the staging of Tumor-Node-Metastasis.⁴⁸ Furthermore, TRMT6-mediated m1A methylation is essential for liver tumorigenesis,⁴⁸ promotes HCC progression via the PI3K/AKT signaling pathway,⁴⁹ and functions as a prognostic biomarker of poor clinical outcomes in patients with HCC.^50,51 Furthermore, targeting TRMT6 inhibition shows promise as an effective therapeutic approach for HCC.⁴⁸

Insulin-like growth factor 2 mRNA-binding proteins (IGF2BP1, IGF2BP2, and IGF2BP3) are RNA-binding proteins that are crucially involved in various biological activities encompassing both developmental and pathological contexts. They regulate the stability and translation of crucial regulators involved in cell division and metabolism.^52–54 Specifically, IGF2BP1 and IGF2BP3 are widely recognized oncofetal proteins that are predominantly expressed in cancer cells rather than in normal tissues. In cancer, these proteins play crucial roles in initiating and promoting tumor growth.^55,56 In patients with HCC, these IGF2BPs are dysregulated, suggesting their involvement in tumor development and progression.^57–59 Disrupting the interactions between IGF2BP1, IGF2BP2, IGF2BP3, and IGF1R inhibits HCC growth.⁶⁰

Overall, this study suggests that the expression of methylation regulator genes in traditional bulk sequencing should be inspected more cautiously because of the associated expression of methylation regulator genes in the complex cellular composition of the tumor environment. Among the many regulator genes, five (DNMT1, NSUN5, TRMT6, IGF2BP1, and IGFBP3), which were identified in the scRNA-seq to be significantly overexpressed in HCC compared to the adjacent para-cancerous liver tissue, should be the focus of further studies regarding the function of m6A/m5C/m1A in HCC. Further evaluation of the impact of these five methylation regulators on the total levels of m6A, m1A, and m5C methylation in the cancer cells would be valuable. Such validation would be valuable in the experiments of cultured cells or hepatoma organoid.