Clinical specimens were acquired from 6 patients who underwent hepatic resection at the First Affiliated Hospital of Weifang Medical University. The cohort comprised 3 patients with HBV-negative hepatic hemangioma (HHE) and 3 patients with HBV-related HCC (positive HBsAg, serum HBV DNA $>$500 IU/mL). All enrolled patients signed informed consent forms, and ethical clearance for this study was obtained from the institutional ethics committee. Strict exclusion criteria were applied prior to patient enrollment, as follows: (1) history of malignancy or concurrent multiple tumors; (2) substantiated evidence of alcohol abuse or a background of nonalcoholic steatohepatitis (NASH); (3) intrahepatic malignancies other than HCC; and (4) previous administration of chemotherapy or local transcatheter arterial chemoembolization (TACE) treatment before surgical intervention. Ultimately, three patients with HBV-related HCC were included in this investigation, with their respective tumor tissues and adjacent liver tissues classified as the HBV-positive tumor group and the HBV-positive liver group (paracancerous group), respectively. Additionally, three patients with HBV-negative hepatic hemangioma were included, and adjacent normal liver tissues were used as the normal group (Supplementary Table 1). Thus, a total of 9 specimens obtained from 6 patients were categorized into 3 groups for comparison (Fig. 1). To investigate the association between HBV and the m7G modification of lncRNAs, comparative analysis was conducted specifically focusing on the paracancerous group and the normal group. Simultaneously, comparison was made between the tumor group and the paracancerous group to assess the plausible correlation between tumorigenesis and anomalous m7G modification of lncRNAs.

Fig. 1.

Flow chart for inclusion and analysis of clinical tissues. HBV, Hepatitis B virus; HHE, hepatic hemangioma; HCC, hepatocellular carcinoma; lncRNAs, long noncoding RNAs; RT-PCR, reverse transcription-PCR; MeRIP-RT-qPCR, methylated RNA immunoprecipitation-reverse transcription-PCR.

Cell culture was conducted using complete medium comprising DMEM (SH30022.01B, HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (S711-001S, Lonsera, Ciudad de la Costa, Carnelones, Uruguay). Normal liver L02 cells and HBV-related hepatoma Hep3B2.1-7 cells were sourced from Duolaimi Biotechnology (Wuhan, Hubei, China) and Pricella Cell Bank (Procell Life Science & Technology Co., Ltd. Wuhan, Hubei, China), respectively. The normal hepatocyte cell line WRL68 and HBV-related hepatoma HepG2.2.15 cells were purchased from Fenghui Biotechnology (Changsha, Hunan, China). All cell lines were validated by STR profiling and tested negative for mycoplasma. Cells were all cultured in a humidified incubator at 37 °C and 5% CO${}_2$. pcDNA3.1 plasmids expressing METTL1 and WDR4 and their small interfering RNAs (siRNAs) were synthesized by Duolaimi Biotechnology (Wuhan, Hubei, China). The pCH-9/3091 plasmid, which contains an HBV 1.1$\times$ genome controlled by the CMV-IE promoter [17], was utilized for HBV transfection. Transfection was carried out using Lipofectamine 2000 (11668019, Thermo Fisher, Waltham, MA, USA) following the manufacturer’s instructions.

After transfecting L02 cells and WRL68 cells with HBV for 72 hours, the cell culture supernatant was collected for analysis. The presence of hepatitis B e antigen (HBeAg) and hepatitis B surface antigen (HBsAg) was determined using diagnostic ELISA kits (202301012 and 202212402, KHB, Shanghai, China) in accordance with the manufacturer’s instructions. To detect HBsAg, 75 µL of the cell culture supernatant sample was added to each well of a reaction plate and incubated at 37 °C for 60 minutes. Subsequently, 50 µL of enzyme conjugate was added to the positive and negative control wells. The plate was gently shaken for 10 seconds and then incubated at 37 °C for 30 minutes. Following a wash step with wash buffer, a chromogenic substrate was added, and the plate was incubated at 37 °C for 30 minutes. The cut-off value (COV) for HBsAg detection was determined as the average optical density (OD) of the negative control plus 0.1. For HBeAg detection, 50 µL of the cell culture supernatant sample was added to each well of the reaction plate and incubated at 37 °C for 30 minutes. After washing, a chromogenic substrate was added, and the plate was incubated at 37 °C for 15 minutes. Finally, optical density (OD) values were measured at 450 nm. The COV for HBeAg was calculated as either 0.05 (if OD 0.05) $\times$ 2.1, or the measured OD (if OD $\ge$0.05) $\times$ 2.1. Samples with OD values greater than the COV were considered positive for HBeAg and HBsAg.

RNA extraction was carried out using RNA Rapid Extraction Kit (220011, Fastagen, Shanghai, China) following established protocols. The concentration of total RNA was determined using a NanoDrop2000 instrument (N50, Implen, LA, USA), and samples with an A260/A280 ratio $\ge$1.8 were deemed suitable for downstream analyses. Reverse transcription (RT) was performed using ReverTra Ace qPCR RT Kit (FSQ-101, Toyobo-Shanghai, Shanghai, China) at 37 °C for 15 minutes. Subsequently, quantitative PCR was conducted using the Prism® 7900HT instrument (ABI, Carlsbad, CA, USA) with SYBR® Premix (QPK-212, Toyobo-Shanghai, Shanghai, China), specific primers (Supplementary Table 2), and ddH${}_2$O. $\beta$-Actin was selected as the internal reference, and quantitative comparisons were obtained using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ calculation method.

Total RNA underwent decapping using an RNA decapping enzyme (M0608S, New England Biolabs, Ipswich, MA, USA) and was subsequently immunoprecipitated using GenSeq® Low Input Whole RNA Library Prep Kit (GS-ET-004, GenSeq, Shanghai, China) following the manufacturer’s instructions. In brief, RNA was randomly fragmented to approximately 200 nt using the RNA fragmentation reagents provided in the kit; the m7G antibody was coupled to Protein A/G beads by rotating at room temperature for 1 hour. The RNA fragments were then incubated with bead-linked antibodies and rotated at 4 °C for 4 hours. After the incubation period, the RNA/antibody complexes underwent washing, and the captured RNA was eluted from the complexes and purified. The libraries were validated using an Agilent 2100 Bioanalyzer (2100, Agilent, Santa Clara, CA, USA) before being sequenced using the NovaSeq platform (Illumina 6000, Illumina-China, Shanghai, China). The m7G-Seq service was provided by CloudSeq Inc. (Shanghai, China).

Post sequencing, MACS software (v1.4.2, Baylor College of Medicine, Houston, TX, USA) [18] was employed to identify methylation sites (peaks) on RNA. Subsequently, diffReps software (v1.55.6, Shen Lab at Mount Sinai, NY, USA) was used to identify and normalize differentially methylated sites based on peak calling with a fold change (FC) $>$2 (or 0.05) and p 0.00001 [19]. Genes of interest were visualized using Integrative Genomics Viewer (IGV) software (v2.4.10, Broad Institute, Berkeley, CA, USA) [20]. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed separately on associated protein-coding genes with differentially m7G methylated peaks derived from lncRNA transcripts. Additionally, sequence motifs associated with m7G were identified within each group [21]. Orthogonal analysis of hyper- or hypo-m7G between HBV+ vs. HBV- and paracancerous vs. tumor samples was conducted to identify potential m7G-related lncRNAs.

Cells were separately transfected with the HBV-1.1 plasmid and overexpressed/knockdown METTL1 and WDR4, followed by total RNA extraction and m7G immunoprecipitation as described above. After RNA immunoprecipitation, differential peaks related to m7G were identified as potential modification regions for each respective gene. Subsequently, we utilized different primers (Supplementary Table 2) in quantitative RT‒PCR to compare differential expression of peak-associated fragments between groups, confirming changes in m7G levels.

Data from three or more independent experiments are presented as the mean $\pm$ standard deviation (SD). Statistical analysis was performed using Student’s t test, and a p value of 0.05 was considered statistically significant. The data were processed and presented using GraphPad Prism software (v9.3, GraphPad Software, San Diego, CA, USA).

The 9 clinical samples obtained from 6 patients were categorized into three groups, as described in Fig. 1. Raw m7G sequencing data were subsequently uploaded to the Gene Expression Omnibus (GEO) database (GSE237003). After summarizing the peaks in each group and conducting intersection analysis, we found that only 97 peaks were shared among all tissues (Fig. 2A). However, corresponding to these peaks, there were 3797 lncRNAs that exhibited the potential for m7G modification across all three groups (Fig. 2B). This indicates the potential presence of multiple peaks available for m7G modifications across different lncRNAs. Furthermore, a consistent overall chromosomal origin of the peaks and lncRNAs was observed (Fig. 2C,D), with lncRNAs primarily being derived from exon sense overlapping regions and intergenic gene regions in all tissues (Fig. 2E). Based on the peak sequences, motif analysis unveiled partially coherent clusters of motifs among the groups. This points toward potential shared enrichment of m7G modification sequences across the different groups (Fig. 2F).

Fig. 2.

The N7-methylguanosine (m7G) modification landscape of lncRNAs within clinical samples. Venn diagrams present (A) detected m7G peaks and (B) matched lncRNAs in different tissues. The distribution of (C) m7G peaks and (D) matched lncRNAs in chromosomal origins. (E) Pie charts depict categorized sources of the detected lncRNAs across various clinical specimens. (F) The top 5 enriched motifs enriched in m7G peaks identified from different clinical tissues.

To investigate variations in m7G methylation among different lncRNA samples, comparative analysis of internal m7G modifications was initially conducted. Compared to the normal group, 856 lncRNAs underwent m7G modification in the paracancerous group, with 76.1% exhibiting hypermethylation (hyper-m7G) and only 23.9% displaying hypomethylation (hypo-m7G) (Fig. 3A). The majority of these m7G-modified lncRNAs originated from exon sense overlapping regions and intergenic regions (Fig. 3B). We observed the same pattern in comparison between the tumor group and the paracancerous group (Fig. 3C,D). This suggests that during progression from HBV infection to HCC, lncRNAs primarily undergo hyper-m7G. Moreover, the chromosomal origin distribution of the differentially m7G-modified lncRNAs was generally consistent with that of the lncRNAs with m7G-related peaks in different tissues (Fig. 3E,F). Interestingly, shorter lncRNAs exhibited a higher frequency of internal m7G modification changes across all groups, regardless of hyper- or hypo-m7G methylation (Fig. 3G,H). This suggests a greater propensity for m7G methylations to occur within shorter lncRNAs.

Fig. 3.

Differential analysis of m7G modifications of lncRNAs across clinical tissues. (A) The distribution and (B) categorization of lncRNAs exhibiting m7G hypermethylation (hyper-m7G) and hypomethylation (hypo-m7G) compared between the paracancerous group and the normal group. (C) The distribution and (D) classification of lncRNAs with alternative m7G modifications compared between the tumor group and paracancerous group. (E) Chromosomal origins of lncRNAs with hyper-m7G and hypo-m7G compared between the paracancerous group and the normal group. (F) Chromosomal origins of lncRNAs with aberrant m7G modification compared between HCC tissue and paracancerous tissue. The length summary of lncRNAs with aberrant m7G modification (G) compared between the paracancerous group and normal group, as well as (H) between HCC tissue and paracancerous tissue; N, normal; P, paracancerous; T, tumor.

We subsequently conducted functional analysis based on the relationship between m7G methylated peaks and the associated functional genes involved via GO and KEGG analyses. When comparing with the normal group, we observed distinct occurrences of hyper-m7G in 689 peaks sourced from different genes in the paracancerous tissue group, whereas 211 peaks displayed hypo-m7G (Fig. 4A). GO and KEGG analyses were then performed using the associated genes of these aberrant peaks. The results revealed that the genes associated with hyper-m7G peaks were enriched in multiple pathways (Fig. 4B–E, Supplementary Fig. 1A–D). Notably, the analysis suggested a potential correlation between the enriched genes and regulation of the Notch pathway (Fig. 4C,E, Supplementary Fig. 1A,D). Multiple pathways were also enriched among the genes associated with hypo-m7G peaks. However, cross-analysis of GO and KEGG did not reveal significant functional pathway convergence (Fig. 4F–I, Supplementary Fig. 1E–H). In subsequent analysis, we compared the HCC group to the paracancerous group. Relative to the paracancerous tissue, the HCC tissue displayed 1086 hyper-m7G peaks and 799 hypo-m7G peaks (Fig. 5A). GO and KEGG analyses revealed that genes associated with both hyper-m7G peaks (Fig. 5B–E, Supplementary Fig. 2A–D) and hypo-m7G peaks were enriched in multiple pathways (Fig. 5F–I, Supplementary Fig. 2E–H), including those associated with cancer regulation (Fig. 5E, Supplementary Fig. 2D,H).

Fig. 4.

Functional analysis of associated functional genes of differentially methylated lncRNAs in paracancerous tissue compared to normal liver tissue. (A) The volcano plot displays the peaks of lncRNAs undergoing m7G methylation changes. For genes associated with hyper-m7G peaks, the top 10 terms of GO analysis according to (B) enriched gene count and (C) differential significance, as well as the top 10 KEGG pathways regarding (D) enriched gene count and (E) differential significance, are summarized. For genes associated with hypo-m7G peaks, the top 10 GO analysis terms based on (F) enriched gene count and (G) lowest p value, as well as the top 10 KEGG pathways regarding (H) enriched gene count and (I) lowest p value, are summarized. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Fig. 5.

Functional analysis of associated genes of differentially methylated lncRNAs in HCC tissue compared to paracancerous tissue. (A) The volcano plot presents the peaks of lncRNAs with alternative m7G modification. For hyper-m7G-associated genes, the top 10 enriched GO analysis terms are summarized based on (B) the number of enriched genes and (C) the extent of differential expression. Similarly, the top 10 enriched KEGG pathways are summarized based on (D) the number of enriched genes and (E) the extent of differential expression. For hypo-m7G-associated genes, the top 10 enriched terms of GO analysis are summarized according to (F) the enriched gene count and (G) the differential significance. Similarly, the top 10 enriched KEGG pathways are summarized regarding (H) the enriched gene count and (I) the lowest p value.

To identify lncRNAs with dynamic m7G modifications from HBV infection to hepatocarcinogenesis, we conducted comprehensive screening based on the extent of m7G alterations across different tissues using a FC threshold of $>$100 or 0.01. It was observed that 470 peaks displayed hyper-m7G and 13 peaks showed hypo-m7G modifications in paracancerous tissue compared to normal liver tissue (Fig. 6A); 352 peaks exhibited FC $>$100 and 338 peaks had FC 0.01 in comparison between the HCC group and the paracancerous group (Fig. 6B). We then matched these peaks with their corresponding lncRNAs and discovered similar proportional changes, and the top 10 FCs of m7G modifications in different group comparisons are presented in Supplementary Tables 3,4. Further intersection analysis revealed that under these screening criteria, 6 lncRNAs (RP11-640N11.2, TEKT4P2, PARP1P1-201, DPP6-210, RP11-989E6.11, DNM1P41) consistently exhibited significant hyper-m7G across intergroup comparisons (Fig. 6C,D). Moreover, we demonstrated the presence of these 6 lncRNAs by visualizing the peaks with significant m7G alterations (Fig. 6E). Therefore, these 6 lncRNAs may serve as potential, dynamically evolving molecular markers during the process from HBV infection to tumorigenesis.

Fig. 6.

Selection of lncRNAs with dynamically changing m7G modification. Heatmaps summarizing m7G alterations of peaks (FC $>$100 or FC 0.01) in the comparison (A) between the paracancerous tissue and the normal tissue and (B) between the HCC tissue and the paracancerous tissue. (C) Venn diagram illustrating the intersection of lncRNAs with significant m7G changes (FC $>$100 or FC 0.01) across different intergroup comparisons. (D) Heatmap showing the occurrence of m7G alterations across the 6 screened lncRNAs across different intergroup comparisons. (E) Visualization of the peak regions with m7G alterations in comparison between the paracancerous tissue and normal tissue (upper 6 panels) and between the HCC tissue and paracancerous tissue (lower 6 panels). MeRIP and Input were overlaid in the same track. N, normal; P, paracancerous; T, tumor; FC, fold change.

To further substantiate the identified 6 m7G-related lncRNAs with dynamic alterations, we transfected the HBV 1.1 plasmid into normal liver L02 cells and WRL68 cells (Supplementary Fig. 3A,B) and subsequently evaluated expression levels of the 6 selected lncRNAs. The results showed significant upregulation exclusively in TEKT4P2 and DNM1P41, but the remaining 4 lncRNAs exhibited no significant changes in expression (Fig. 7A,B). Moreover, HBV transfection led to an increase in m7G levels in both TEKT4P2 and DNM1P41 (Fig. 7C,D). We then investigated alteration of the m7G methyltransferases METTL1 and WDR4 after HBV transfection. The findings revealed significant upregulation of METTL1 and WDR4 in response to HBV in both L02 and WRL68 cells (Fig. 7E,F). Therefore, we hypothesize that expression levels of TEKT4P2 and DNM1P41 are regulated by HBV-mediated m7G modifications. To test our hypothesis, we overexpressed METTL1 and WDR4 in L02 and WRL68 cells (Supplementary Fig. 3C,D) and observed significant upregulation of TEKT4P2 and DNM1P41 expression levels (Fig. 8A,B). In addition, we overexpressed METTL1 and WDR4 in HBV-related hepatoma cells (Supplementary Fig. 3E,F), Hep3B2.1-7 and HepG2.2.15 cells and found upregulation of TEKT4P2 and DNM1P41 (Fig. 8C,D). To ascertain the association between this upregulation and m7G modification, we overexpressed METTL1 and WDR4 in the aforementioned cells and assessed changes in m7G modifications of TEKT4P2 and DNM1P41. The results indicated that the m7G levels of TEKT4P2 and DNM1P41 were indeed influenced by regulation of METTL1 and WDR4, thereby affecting their own expression levels (Fig. 8E–H). Finally, to confirm the significant regulatory roles of HBV and m7G methyltransferases on TEKT4P2 and DNM1P41, we initially examined differences in m7G modifications of TEKT4P2 and DNM1P41 among different cell lines. The results showed that the m7G levels of TEKT4P2 and DNM1P41 were significantly higher in hepatoma cells than in normal hepatocytes (Fig. 9A,B). We thus transfected HBV 1.1 into L02 and WRL68 cells while simultaneously knocking down m7G methyltransferases (METTL1 or WDR4) (Supplementary Fig. 3G,H). The results demonstrated that both HBV and m7G methyltransferases have remarkable impacts on TEKT4P2 and DNM1P41 m7G modification (Fig. 9C–F).

Fig. 7.

Preliminary validation of identified m7G-related lncRNAs. qRT-PCR was performed to assess expression of the 6 screened potential m7G-related lncRNAs in (A) L02 cells and (B) WRL68 cells following transfection with the HBV 1.1 plasmid. Using (C) L02 cells and (D) WRL68 cells, MeRIP-PCR was carried out to examine changes in m7G levels of TEKT4P2 and DNM1P41 after HBV transfection. Levels of the m7G methyltransferases METTL1 and WDR4 were assessed by qRT-PCR in (E) L02 cells and (F) WRL68 cells following HBV transfection. The mean and SD from three independent experiments are depicted. **p 0.01, ***p 0.001.

Fig. 8.

The expression pattern of m7G modification-dependent dynamic lncRNAs. Expression of TEKT4P2 and DNM1P41 was assessed by qRT‒PCR in normal hepatocytes (A) L02 cells and (B) WRL68 cells and hepatoma (C) Hep3B2.1-7 cells and (D) HepG2.2.15 cells after transfection of METTL1 and WDR4. MeRIP-PCR was carried out to detect m7G alteration of TEKT4P2 and DNM1P41 in (E) L02 cells, (F) WRL68 cells, (G) Hep3B2.1-7 cells, and (H) HepG2.2.15 cells following overexpression of METTL1 and WDR4. ***p 0.001.

Fig. 9.

The impact of HBV and methyltransferase on m7G modification of dynamic lncRNAs. m7G levels of (A) TEKT4P2 and (B) DNM1P41 in hepatocytes and hepatoma cells were determined by MeRIP-PCR. MeRIP-PCR was performed to detect changes in m7G modification of TEKT4P2 and DNM1P41 after HBV transfection and METTL1 knockdown in (C) WRL68 cells and (D) L02 cells. Alteration of m7G levels of TEKT4P2 and DNM1P41 was assessed by MeRIP-PCR following transfection of HBV and knockdown of WDR4 in (E) WRL68 cells and (F) L02 cells. *p 0.05, **p 0.01, ***p 0.001.