Transcriptome-Wide 5-Methylcytosine Functional Profiling of Long Non-Coding RNA in Hepatocellular Carcinoma

Introduction

Posttranscriptional RNA modifications have emerged as an additional layer of epigenetic regulation.¹ Over one hundred different chemical modifications in RNA have been described, including pseudouridine, N6-methyladenosine (m6A),² N1-methyladenosine (m1A),³ and 5-methylcytosine (m5C).⁴ Accumulating evidence suggests that RNA modifications play critical roles across a variety of biological functions that can contribute to disease pathogenesis. For instance, m6A modification of Wilms tumor 1-associated protein (WTAP) positively correlates with liver cancer oncogenesis through the HuR-ETS1-p21/p27 axis,⁵ and m1A-related gene dysregulation is associated with tumorigenesis in gastrointestinal cancers,⁶ with ErbB2 and mTOR signaling pathways likely modulated by m1A regulators.

In particular, m5C is one of the most well-known RNA modifications that was first identified in stable and highly abundant ribosomal RNAs (rRNAs) and transfer RNAs (tRNAs).⁷⁻⁹ m5C sites in tRNAs usually occur around the variable region and anticodon loop. As a heavily m5C modified RNA class, tRNA modifications stabilize tRNA secondary structure, confer metabolic stability, and affect aminoacylation and codon recognition.^10,11 m5C sites identified in rRNA play an important role in rRNA processing, structure, and translation.¹² With the development of effective techniques like high-throughput sequencing, it is now possible to identify and quantify m5C modifications in low-abundance RNA species, such as mRNA and long non-coding RNAs (lncRNAs).¹³ A recent study reported that m5C sites were mainly enriched in CG motifs of mRNA, showing conserved, tissue-specific, and dynamic properties in the mammalian transcriptome to regulate mRNA export.¹⁴ To date, the distribution and function of m5C across diverse classes of RNA have been explored.^15,16 However, knowledge surrounding the prevalence and transcriptome-wide distribution of m5C in lncRNA is still very limited.

lncRNAs are a class of RNA molecules that are more than 200 nucleotides long and that derive from non-coding regions of the genome.¹⁷ The functions and mechanisms of most lncRNAs are unclear. In recent years, lncRNAs have obtained widespread attention as key regulators of gene expression in various physiological and pathological processes.¹⁸ Epigenetic regulation is one of the main mechanisms controlling lncRNAs expression and tissue specificity. As an important epigenetic modification, RNA methylation has been identified as an essential marker of oncogenesis.

Hepatocellular carcinoma (HCC), the fourth leading cause of cancer-related deaths globally, is primarily caused by viral hepatitis infections and exposure to aflatoxin and aristolochic acid.^19–21 Although recent decades have yielded breakthroughs in HCC prevention, early detection, diagnosis, and treatment, the global disease burden remains significant due to malignancy and heterogeneity. However, several recent studies have revealed that methylation modifications participate in HCC development and progression, suggesting that methylation could be a vital therapeutic target for HCC.^22–24

In this study, we sought to obtain a deeper understanding of m5C methylation in HCC lncRNA. To this end, we globally mapped m5C using RNA MeRIP-seq in human HCC cells and normal adjacent non-tumor tissues, comparing its prevalence and distribution in both cell/tissue types and cellular compartments. We found the degree of m5C modification in HCC to be much greater than that of adjacent tissue. This difference can be characterized by intra-tissue consistency and inter-tissue differences involving all chromosomes. Our findings may provide new insights into the m5C epigenetic regulation of lncRNA in HCC toward the development of new therapies.

Materials and Methods

Sample Collection and RNA Sequencing Samples Preparation

Specimens of resected HCC and adjacent non-tumor tissues controls were collected from HCC patients. In total, six pairs of HCC and adjacent non-tumor tissue samples were obtained. RNA from each sample was isolated from cells using TRIzol reagent (Invitrogen Corporation, Carlsbad, CA) following manufacturer’s instructions. A NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) was used to quantify RNA concentration. RNA with OD260/OD280 between 1.8 and 2.2 was marked as qualified. LncRNA was fragmented with fragmentation buffer into short fragments (about 100 bp) by GenSeq m5C RNA IP Kit (GenSeq Inc, China).

cDNA Library Construction

Based on the MeRIP-seq technique, we obtained m5C peaks from each sample. RNA immunoprecipitation (RIP) was performed by the GenSeq m5C RNA IP Kit (GenSeq Inc, China). The NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs, Inc.) was applied to construct the cDNA library. The quality and purity of the cDNA library for each sample were validated using an Agilent Bioanalyzer 2100 system. The cDNA libraries were sequenced on a HiSeq 4000 platform (Illumina, San Diego, CA, USA) with 150 bp paired-end reads.

Sequencing Data Analysis

Libraries were sequenced and visualized on a HiSeq 4000 (Illumina). ‘Q30’ identifies the mapping quality, and reads higher than 30 were considered as reliable mapping and retained. Library mapping was drawn through STAR software (http://gingeraslab.cshl.edu/STAR).²⁵ Sequencing datasets were aligned to the reference genome (//http.ncbi.nlm.nih.gov/genomes/Helicoverpa-armigera/) using HISAT software (2.0.4). Methylation peaks and their differences between HCC and adjacent non-tumor tissue were identified by MACS software and DiffReps, respectively.^26,27

GO and KEGG Pathway Databases Analysis

GO (http://www.geneontology.org) and KEGG (https://david.ncifcrf.gov) pathway enrichment analyses were performed using DAVID to reveal the functions of these differentially methylated genes. GO analysis was performed based on differentially methylated genes to analyze the biological processes, cellular components, and molecular functions of these genes. P < 0.05 was considered as the cutoff for significant GO terms. KEGG pathway enrichment for the differentially methylated and expressed genes was also performed to investigate their biological pathway and functions, as previously described.²⁸ P-values <0.05 were considered significant.

Ethical Approval

The study has been approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University, and the patient’s written informed consent was obtained before the study began. And the study was conducted in accordance with the Declaration of Helsinki.

Results

Transcriptome-Wide m5C Methylation of lncRNA in Human HCC

To gain insight into transcriptome-wide m5C methylation, we performed RNA MeRIP-seq of lncRNA and identified 21,960 m5C peaks in HCC, and 15,325 m5C peaks in adjacent non-tumor tissues (Figure 1A and B). Up to 19,935 annotated genes (HCC tissues) and 14,127 annotated genes (adjacent tissues) were mapped, respectively. Among them, 12,928 m5C peaks were detected within both HCC and its adjacent non-tumor tissues. m5C sites with fold changes (FCs) >2 are considered to be specific to a certain tissues. Comparing HCC with adjacent non-tumor tissue revealed that 39% of identified sites were specific to HCC, whereas 14% were specific to adjacent non-tumor tissues. This showed that the overall occurrence of m5C methylation is higher in HCC than adjacent non-tumor tissue.

Figure 1 Transcriptome-wide m5C methylation and characters of lncRNA in human HCC. (A) Venn diagram of m5C methylation sites identified in lncRNA from human HCC and adjacent non-tumor tissues. (B) Venn diagram of m5C genes in HCC and adjacent tissues. (C) The sequence motif of m5C sites in human HCC and adjacent non-tumor tissues. (D) Percentage of lncRNA harboring different numbers of m5C peaks in the two tissue types, with the majority harboring only one m5C peak.

UARUCCCA is a Conserved Sequence Motif for m5C-Containing Regions

To determine the presence of an m5C motif, we scanned the sequences of m5C-methylated peaks. UARUCCCA (R = purine) was determined to be the most reliable motif in HCC (Figure 1C). Interestingly, the most reliable motif for m5C peaks in the adjacent non-tumor tissues was also UARUCCCA, indicating that the same sequence motif may be necessary for m5C methylation in HCC lncRNA. We speculate that m5C in HCC and adjacent non-tumor tissues may be caused by the same methylase.

Abundance of m5C Peaks in lncRNA

We next determined the abundance of peaks in lncRNA according to tissue type. Notably, 91.4% of the lncRNA contained one m5C peak, which was slightly lower than the 92.9% of single peak lncRNA in adjacent non-tumor tissues (Figure 1D). Percentages for two peaks (7.7% HCC vs 6.4% adjacent non-tumor tissues), three peaks (0.79% HCC vs 0.58% adjacent non-tumor), and more than three peaks (0.17% HCC vs 0.16% adjacent non-tumor tissues) were also determined.

Heatmap of m5C Levels in Human HCC and Adjacent Non-Tumor Tissues

Quantitatively, the methylation of m5C sites varied between HCC and adjacent non-tumor tissue based on unsupervised hierarchical cluster analysis. A heatmap of methylation levels showed the apparent expression difference between the two tissue types (Figure 2A). Specifically, cluster analysis revealed that the degree of m5C methylation had intra-tissue consistency and inter-tissue differences between HCC and adjacent non-tumor tissue.

Figure 2 m5C levels and distribution of chromosome in human HCC and adjacent non-tumor tissue. (A) A cluster tree at the top indicates the relatedness of the two tissues (x-axis) or methylation levels (y-axis). The heatmap uses a color scale to indicate the relative methylation level at each locus. Each row of colored lines (N = 117) represents the methylation level for each CpG locus: red for hypermethylated and yellow for hypomethylated. Each column (n = 12) corresponds to each tissue. (B) Circos plot showing the distribution of m5C methylation sites on each chromosome. Red represents HCC, blue represents adjacent non-tumor tissue.

Chromosome Visualization of m5C in lncRNA

In HCC and adjacent non-tumor tissues, we investigated the distribution of m5C methylation sites throughout the whole genome (Figure 2B). We found that m5C sites were dispersed across all chromosomes, which may indicate the extent to which m5C functions in cells. m5C methylation levels and distribution were different on each chromosome between the two tissue types.

The Differences in lncRNA Expression According to m5C Methylation

To explore differences in lncRNA expression according to m5C methylation, we researched lncRNA expression in HCC and adjacent non-tumor tissues using the RNA sequencing data. We found that more lncRNA were up-regulated by methylation and tended to have corresponding higher expression levels in HCC tissues (Figure 3A). On the contrary, more lncRNA were down-regulated by methylation in adjacent non-tumor tissues than HCC. Interestingly, we found that some lncRNA with minimal expression in adjacent non-tumor tissues were highly expressed in HCC and up-regulated by methylation (Table 1). These genes may be related to the pathogenesis of HCC and are worth further investigation. Among lncRNAs with m5C modifications, more were down-regulated than up-regulated in the HCC tissues (Figure 3B). Together, the data indicated that m5C tends to correlate with gene expression in a large fraction of transcripts in HCC tissues.

Table 1 Top 10 Highly Expressed lncRNAs in HCC

Figure 3 The differences in lncRNA expression according to m5C methylation. (A) Differentially expressed lncRNAs in HCC and adjacent non-tumor tissues. Genes up-regulated by methylation are red, and genes down-regulated by methylation are blue. (B) Cumulative distribution of lncRNA expression changes between HCC and adjacent non-tumor tissue for m5C up-regulated genes (red) and m5C down-regulated genes (green), whereas blue represents others.

Gene Ontology (GO) Enrichment Analysis

To examine the function of m5C methylation in HCC and adjacent non-tumor tissues, GO enrichment analysis was performed based on biological processes (BP), cellular components (CC), and molecular functions (MF). For BP, genes with up-methylated m5C sites were significantly enriched in cell development, system development, multicellular organismal development, and neurogenesis (Figure 4A), whereas genes with down-methylated m5C sites were highly enriched in cellular processes such as nervous system development, neurogenesis, and generation of neurons and synaptic transmission (Figure 4D). For MF, genes with up-methylated m5C sites were notably enriched in enzyme binding, sequence-specific DNA binding transcription factor activity, nucleic acid binding transcription factor activity, and sequence-specific DNA binding (Figure 4B). Furthermore, genes with up-methylated m5C sites were enriched in gated channel activity, substrate-specific channel activity, ion channel activity, ligand−gated channel activity, and others (Figure 4E). For CC, genes with up-methylated m5C sites were mainly enriched in cell projection, neuron projection, neuron part, and others, as was the function of genes with down-regulated methylation (Figure 4C and F).

Figure 4 Gene Ontology (GO) term enrichment analysis of m5C genes in HCC lncRNAs. (A-C) The top 10 GO terms for (A) biological processes, (B) molecular functions, and (C) cellular components were significantly enriched for up-methylated m5C genes in HCC. (D–E) The top 10 gene GO terms of (D) biological processes, (E) molecular functions, and (F) cellular components were significantly enriched for the down-methylated m5C genes in HCC.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

Notably, KEGG pathway analysis indicated that hypermethylated genes were significantly related to endocytosis, thyroid hormone signaling pathway, pathways in cancer, and dilated cardiomyopathy (Figure 5A), whereas down-methylated genes were associated with glutamatergic synapse, glycosaminoglycan biosynthesis-heparan sulfate, transcriptional misregulation in cancer, and retrograde endocannabinoid signaling (Figure 5B).

Figure 5 KEGG pathway analysis of m5C genes in HCC lncRNAs. (A) Bar plot showing the top 10 enrichment scores of significantly enriched pathways for up-methylated m5C genes in HCC. (B) Bar plot showing the top 10 enrichment scores of significantly enriched pathways for down-methylated m5C genes in HCC.

Discussion

HCC has become the second leading cause of cancer-related deaths, which accounts for approximately 90% of all cases of primary liver cancer.²⁹ HCC pathogenesis is a complex process, with epigenetic alterations, heritable gene expression changes, and chromatin organization shown to drive tumor initiation and progression.³⁰ Recently, epigenetic alterations, in particular methylation, have come into focus for cancer pathogenesis and prognosis. Epigenetic alterations to DNA, RNA methylation, and histone deacetylation are the most frequent epigenetic events that play numerous roles in tumorigenesis.³¹

Characteristically, DNA methylation could regulate tumor metastasis through up-regulation of DNA methyltransferases and a series of molecular mechanisms in HCC.³² Besides, CPG island hypermethylation in promoter regions was found in HCC tissues, which could silence the expression of tumor suppressor genes and promote tumor development.³³ In addition to DNA methylation, RNA modifications were revealed to moderation of RNA structure, function, and stability. Previous study reported that RNA m6A methyltransferase-like 3 was up-regulated in HCC samples by m6A-YTHDF2-dependent posttranscriptional silencing of suppressor of cytokine signaling 2.³⁴ Another research found that ALYREF, a m5C-related regulatory gene, was associated with cell cycle regulation and mitosis in HCC progression.⁴

m5C is ubiquitous in nature and one of the most important RNA modifications identified in tRNA, rRNA, and mRNA.³⁵ m5C has been implicated in protein translational regulation, RNA processing, and stress responses, and been closely correlated with cell proliferation and differentiation.^36,37 The levels of m5C regulation are significantly different according to tissue type.³⁸ For instance, m5C levels in metastatic neoplasms were much lower than that of benign neoplasms or normal tissue. It is worth mentioning that in humans, aberrant regulation of m5C causing CpG events is underrepresented.³⁹ Over one-third of germline point mutations are closely correlated with human genetic diseases, and many somatic mutations leading to cancer involve the loss of CpG.

It has been reported that genes have a high occurrence of methylation in gastric cancer, breast cancer, and Hodgkin’s and non-Hodgkin lymphomas, ranging from 60% to 90%.⁴⁰ Xue et al demonstrated that TRDMT1, an important member of the m5C family, plays a crucial role in inhibiting the proliferation and migration of HEK293 cells.⁴¹ Li et al showed that m5C enhanced the extent of overall DNA modifications, which lays a foundation for somatic mutation.⁴² However, m5C expression levels and its function in HCC remains poorly defined. We estimated the degree of cellular methylation by measuring the ratio of m5C in HCC tissues to adjacent non-tumor tissues.

For this reason, we focused our research on m5C methylation in HCC. Our results show that the number and distribution of m5C on HCC and adjacent non-tumor tissues are markedly different. In short, the m5C sites of lncRNA in HCC and the genes mapped by m5C sites far exceed those in adjacent tissues. These differences in methylation are spread across all chromosomes. Cluster analysis results show that the degree of methylation can clearly distinguish HCC from adjacent tissues, which further confirms the potential relationship between m5C and HCC. Interestingly, we found that some genes minimally expressed in adjacent non-tumor tissues are not only highly expressed in HCC, but also have an increased degree of methylation. The role of methylation of these genes in the pathogenesis of HCC deserves further discussion.

As a potentially critical layer of bioregulation, lncRNAs were reported to be aberrantly expressed through histone modification or DNA methylation, or transcription factors, and posttranscriptional modification in HCC.⁴³ Dys-regulated lncRNAs might change the expression and activities of mRNAs, miRNA and proteins, resulting in epigenetic alteration in HCC. For instance, the up-regulated lncRNA DANCR in HCC could bind to 3 UTR of CTNNB1 mRNA, blocking the inhibitory effects of mir-214 miR320a and miR199a, thus increasing the stemness-like characteristics of HCC.⁴⁴ Besides, lncRNA PVT1 could promote cell proliferation, cell cycle progression, and increased stem-like characteristics through stabilizing the nucleolar protein NOP2 in HCC.⁴⁵ Transcriptional expression levels of lncRNAs in HCC have been widely studied. However, the mechanism of m5C methylation in lncRNAs promoting HCC progress is unclear, and deeper exploration would be helpful for understanding the pathogenesis of HCC.

Conclusion

These findings have significant implications for understanding the mechanisms of m5C regulation in HCC. Further exploration of epigenetic alterations in the pathology of HCC is still needed.