ALKBH5 suppresses malignancy of hepatocellular carcinoma via m6A-guided epigenetic inhibition of LYPD1

Two HCC cohorts were included in this study, which was approved by Institutional Ethics Committee in First Affiliated Hospital of Zhejiang University. Cohort one contained 80 HCC patients who had undergone curative surgery from 2015 to 2018 in our hospital. Specimens of tumor and adjacent tissues were collected from these patients. RNA (70 pairs) or proteins (10 pairs) were isolated from frozen tissues for quantitative real-time PCR (qPCR) or western blotting assay to assess the expression of ALKBH5 in HCC. Cohort two consisted of 90 HCC patients which were the source of commercial tissue microarrays (TMA) supplied by Shanghai Outdo Biotech (LivH180Su07, Shanghai, China) together with integrated follow-up and clinical information data (illustrated in Table S). This TMA cohort was employed to evaluate the role of ALKBH5 in HCC prognosis and construct the correlation of ALKBH5 and LYPD1 expression. Written informed consents were acquired from each patient relying on guidelines of the Declaration of Helsinki. 1

The human HCC cell lines Huh7, MHCC97H, HCCLM3, HepG2, Hep3B, PLC/PRF/5, SMCC7721 and BEL7402 were obtained from the Shanghai Institutes of Biological Sciences (Shanghai, China). STR finger printing authentications of all employed HCC cells are available upon request. And here we present STR certificates for three of them (Huh7, MHCC97H and HCCLM3), which are mainly investigated in our study (Additional file). These cells were incubated at 37 °C in a 5% CO2 incubator (ThermoFisher, USA) with the humidified environment. And they were cultured with Minimum Essential Media (MEM, BI, Israel), which were routinely supplemented with 10% fetal bovine serum (FBS, BI), penicillin (100 units/ml) and streptomycin (100 μg/ml). 10

Total RNA isolated from HCC cells or subcutaneous tumors was mixed in three times volume of incubation buffer and denatured at 65 °C for 5 min. Samples (400 ng, 200 ng or 100 ng) dissolved in SSC buffer (Sigma-Aldrich, Germany) were deposited on an Amersham Hybond-N+ membrane (GE Healthcare, USA) which was settled on the Bio-Dot Apparatus (Bio-Rad, USA). Then the membrane was crosslinked by UV light for 5 min, followed by the staining with 0.02% Methylene blue (Sangon Biotech, China). Scanning of blue dots were performed to show the input RNA content. And the membrane was hatched with m⁶A antibody (1:5000, Synaptic System, #202003) overnight at 4 °C. Dot blots were visualized by the imaging system after incubation with secondary antibody.

MeRIP-sequencing and following data analyses were mainly supported by Genesky Biotechnologies Inc. (Shanghai, China). HCCLM3 cells with stable overexpression of ALKBH5 and control cells transfected with an empty vector were collected (two replications; labeled as ALKBH5 and Vector, respectively). More than 250 μg of total RNA was extracted from each group, and mRNA was further purified with the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB#E7490) using oligo (dT) beads. Concentration and integrity of RNA were evaluated using NanoDrop Spectrophotometer and Agilent 2100 Bioanalyzer (Agilent, USA). The mRNA was then chemically fragmented into ~ 150 nts nucleotides with fragmentation buffer. After 10% of fragmented mRNA was saved as input, m⁶A-modified mRNA was immunoprecipitated with anti-m⁶A antibodies (Synaptic System, #202003) and eluted. RNA sequencing libraries for input mRNA (RNA-seq) and m⁶A-enriched mRNA (MeRIP-seq) were simultaneously constructed with the VAHTS Total RNA-seq (H/M/R) Library Prep Kit for Illumina (Box2&3, Vazyme#NR603), followed by the sequencing on Novaseq sequencer (Illumina, USA) with PE150 strategy. The MeRIP-seq data was analyzed based on the published standardized pipeline [24]. In brief, the raw data was aligned to human genome GRCh37/hg19 by the HISAT2 software (v2.0.5). Then m⁶A peaks were determined by the ExomePeak software (v2.6.0) and annotated according to the Ensembl database. Integrative Genomics Viewer (IGV) software was applied to present the visualization of the m⁶A peaks distribution. On the other hand, RNA-seq reads of input samples were normalized with Cufflinks (v2.2.1) [25], and Cuffdiff was employed to determine differentially expressed genes [26]. And major R codes during analyses were provided in Additional file 11.

RIP assay was conducted with Magna RIP Kit (Millipore, Germany) according to manufacturer’s illustrations. Briefly, magnetic beads were mixed with 5μg anti-ALKBH5 (Sigma-Aldrich, Germany) or IGF2BP1 (Abclonal, China) and anti-rabbit IgG (Millipore, Germany) before the addition of cell lysates (approximately 2*10⁷ cells for each sample). After the treatment of proteinase K, interested RNAs were eluted from immunoprecipitated complex and purified for further analysis using qPCR. Relative enrichment was normalized to the input: %Input =1/10 × 2^{Ct [IP] – Ct [input]}.

MeRIP assay was performed with the Magna MeRIP™ m⁶A Kit (Millipore, Germany) to determine the m⁶A modification on individual transcripts. In brief, 150 μg total RNA was isolated from pretreated cells and randomly fragmented into a size of 100 or less nucleotides. RNA samples were then immunoprecipitated with magnetic beads pre-coated by 10 μg anti-m⁶A antibody (Millipore, Germany) or anti-mouse IgG (Millipore). And N6-methyladenosine 5′-monophosphate sodium salt (6.7 mM) were applied to elute the m⁶A-modified RNA fragments. Based on MeRIP-seq results, we focused on the sites of LYPD1 transcript where differential m⁶A peak was identified between ALKBH5-overexpressing cells and empty control cells (Fig. 5a). Specific primers were designed for MeRIP-qPCR analysis according to the information from MeRIP-seq and a motif-dependent m⁶A site predictor SRAMP (http://www.cuilab.cn/sramp↗) (Forward: AGCAGAATTGGCTGGTTTCG; reverse: AGCCCCAGTCTAAGTCCCA). Relative enrichment of m⁶A was normalized to the input: %Input =1/10 × 2^{Ct [IP] – Ct [input]}.

Statistical analysis was performed with the GraphPad Prism 8.0 (GraphPad, Inc., USA) and SPSS 22.0 (SPSS, Inc., USA) software. Experiments were independently repeated for at least three times. Representative data was exhibited as the means ± SD. Quantitative data was compared using two-tail Student t test, while qualitative data was evaluated by Chi-Square test. The overall and recurrence-free survival were analyzed with Kaplan–Meier method and log-rank test. And univariate and multivariate Cox regression models were employed to investigate independent prognostic factors. In addition, correlational analysis of gene expression was conducted with linear regression. P-values for every result were labeled on figures, and P < 0.05 was reckoned as statistically significant (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

More detailed methodology could be obtained in Supplementary Materials and Methods.

To investigate the expression profile of ALKBH5 in HCC, we analyzed the mRNA and protein levels of ALKBH5 in HCC and matched adjacent tissues, and found that ALKBH5 was significantly down-regulated in HCC (Fig. 1a-c). And subsequent immunohistochemistry (IHC) staining with TMA from an independent HCC cohort confirmed these results (Fig. 1d, e). Besides, HCC patients with lower ALKBH5 expression obtained shorter overall survival (OS) and recurrence-free survival (RFS) (Fig. 1f). This finding was further validated by results from the Cancer Genome Atlas (TCGA) database (Additional file 5: Figure S1a). The expression of ALKBH5 also seemed to perform well in survival prediction of early-stage HCC patients (Additional file 5: Figure S1b). Moreover, loss of ALKBH5 was identified as an independent prognostic factor for HCC patients (HR = 2.24, P = 0.007) (Fig. 1g). It implies that dysregulation of ALKBH5 may be involved in the progression of HCC.

To evaluate the functional roles of ALKBH5 in HCC, we firstly examined the expression of ALKBH5 in HCC cell lines (Additional file 5: Figure S1c, d). Huh7 or MHCC97H and HCCLM3 cells were chosen to establish ALKBH5-silencing and ALKBH5-overexpressing models, respectively. And the transfection efficiency was validated by qPCR and western blotting (Additional file 5: Figure S1e-h). As were indicated by CCK-8 and colony formation assays, knockdown of ALKBH5 enhanced the proliferation capability of HCC cells, while up-regulation of ALKBH5 showed the opposite effect (Fig. 2a). Similarly, EdU assay also demonstrated that ALKBH5 could suppress cell growth in vitro (Fig. 2b-e). Moreover, the re-expression of ALKBH5 was introduced into stable ALKBH5-knockdown cells (Additional file 6: Figure S2a, b). As expected, the restoration of ALKBH5 could rescue the enhanced proliferation capabilities induced by ALKBH5 loss in both Huh7 and MHCC97H cells (Additional file 6: Figure S2c, d), which was supported by EdU results as well (Additional file 6: Figure S2e-g).

To further address the anti-oncogenic role of ALKBH5 in HCC, we conducted in vivo experiments with subcutaneous tumor models. When ALKBH5 was silenced (Additional file 5: Figure S1i, j), volumes and weights of xenografted tumors increased compared with control group (Fig. 2f). On the contrary, ALKBH5 overexpression retarded tumor growth with considerably diminished tumor volumes and weights (Fig. 2g). And attenuation of PCNA, a marker of proliferation, was also detected along with the elevated expression of ALKBH5 (Additional file 6: Figure S2h). These results suggest that ALKBH5 exerts an inhibitory effect on HCC tumor growth in vitro and in vivo.

We performed transwell assays and noticed that inhibition of ALKBH5 promoted both migration and invasion abilities of HCC cells, while overexpression of ALKBH5 impaired these phenotypes (Fig. 3a). Then the wound healing assay also indicated that ALKBH5 tended to attenuate migration of HCC cells (Fig. 3b). Interestingly, we always observed the altered cell morphology under microscope when ALKBH5 was silenced (data not shown). To verify whether this phenomenon was due to the reshaping of the cytoskeleton, phalloidin staining was performed subsequently. As expected, knockdown of ALKBH5 led to a looser and more divergent pattern of cytoskeleton through the rearrangement of microtubules and microfilaments (Fig. 3c, d), which denoted a more active migrating form.

To clarify the effects of ALKBH5 on HCC metastasis in vivo, ALKBH5-overexpressing and negative control HCCLM3-luc cells were implanted into BALB/c mice via tail vein injection, followed by the bioluminescence imaging. It seemed that activation of ALKBH5 damaged the metastatic potential of HCC cells with lower luciferase activity and less pulmonary metastasis (Fig. 3e, f), which was confirmed by HE staining results (Fig. 3g). In contrast, silencing of ALKBH5 promoted the metastasis of HCC (Additional file 6: Figure S2i). Therefore, ALKBH5 suppresses the migration/invasion abilities of HCC cells in vitro and their metastatic capabilities in vivo.

We firstly applied dot blot assays to examine the role of ALKBH5 in modulating m⁶A modification. Loss of ALKBH5 led to a convincingly increased m⁶A level in both Huh7 and MHCC97H cells, while ALKBH5 overexpression generated the opposite result (Fig. 4a). And the same conclusion could be obtained from the analysis of subcutaneous tumors (Additional file 6: Figure S2j).

To find out the precise mechanisms underpinning the observed ALKBH5-dependent phenotypes, an integrated approach combining MeRIP-seq and RNA-seq was employed using stable ALKBH5-overexpressing and vector-transfected HCCLM3 cells. MeRIP-seq revealed 1538 differential m⁶A peaks with reduced abundance (1344 corresponding transcripts) when ALKBH5 was up-regulated. Meanwhile, RNA-seq uncovered 481 down-regulated transcripts upon ALKBH5 overexpression.

We attached more importance to oncogenes whose methylation patterns and expression levels were regulated by ALKBH5. Therefore, merely those transcripts owning both hypo-m⁶A-peaks and decreased expression upon ALKBH5 overexpression were selected for following investigations (Fig. 4b). To further narrow down the scope of candidates, we focused on the top 10 genes from the overlap, namely COCH, LYPD1, ADAMTS14, ABCA4, TP53I11, COLCA2, TMED7, CYP4F3, IL17RB and VCAN, listed in ascending order of expression fold change. They were subject to preliminary validation in ALKBH5-silencing or -overexpressing cells by qPCR (Fig. 4c). Intriguingly, only LYPD1 was consistently found to be inversely regulated by ALKBH5 in all three HCC cells (Fig. 4d-g; Additional file 7: Figure S3a-e), which was further confirmed by western blotting results (Fig. 4h). Taken together, LYPD1 may be the direct downstream target of ALKBH5.

Our MeRIP-seq analysis suggested that m⁶A peak of LYPD1 in 3’UTR shrank remarkably with the overexpression of ALKBH5 (Fig. 5a). To substantiate this result, we first conducted RIP assays using the anti-ALKBH5 antibody in Huh7 and HCCLM3 cells. We observed that ALKBH5 could enrich LYPD1 mRNA (Fig. 5b), implying that LYPD1 may be regulated in RNA level upon interaction with ALKBH5. Then MeRIP-qPCR assays with specific primers aiming at potential m⁶A sites revealed that knockdown of ALKBH5 could promote m⁶A modification of LYPD1 in 3’UTR, while activation of ALKBH5 led to a decreased m⁶A level in this site (Fig. 5c). To further demonstrate the essential role of m⁶A in the regulation of LYPD1, we designed a luciferase reporter inserting a wild-type (WT) LYPD1–3’UTR sequence or mutant (Mut) counterpart whose putative m⁶A sites were mutated (Fig. 5d). As expected, the luciferase activity of cells transfected with LYPD1-WT plasmid tended to increase when ALKBH5 was silenced, while that of mutant group seemed to be unaffected. And analogous results could be verified in ALKBH5-overexpressing cells (Fig. 5e). In addition, we found that ALKBH5 deficiency induced a slower degradation rate of LYPD1 mRNA, whereas ALKBH5 overexpression abolished the stability of LYPD1 (Fig. 5f).

Now that “readers” were crucially responsible for the direct effect on m⁶A-modified transcripts, we investigated potential effectors participating in the process illustrated above. As YTHDFs and IGF2BPs were extensively involved in the modulation of RNA stability [11], YTHDF1–2 and IGF2BP1–3 were knockdown successively in two HCC cells to examine the alterations of LYPD1 expression. We noticed that intervene of YTHDF1/2 and IGF2BP2/3 hardly impacted LYPD1 expression (Additional file 7: Figure S3f-i). Nevertheless, LYPD1 was significantly inhibited when IGF2BP1 was impaired (Fig. 5g), which was consistent with the knowledge that IGF2BP1 intended to promote the transcription of its targets [10]. And the interaction between IGF2BP1 protein and LYPD1 mRNA was confirmed by RIP assays (Fig. 5h). Moreover, knockdown of IGF2BP1 counteracted the accumulation of LYPD1 caused by ALKBH5 loss (Fig. 5i). In summary, LYPD1 is governed by ALKBH5-mediated m⁶A modification and recognized by IGF2BP1 which enhances its stability.

To illustrate the role of LYPD1 in HCC, we established LYPD1-knockdown Huh7 and MHCC97H cell lines (Fig. 6a, b; Additional file 8: Figure S4a). CCK-8 and colony formation assays indicated that silencing of LYPD1 suppressed cell growth and viability (Fig. 6a, b), which was consistent with results of EdU (Fig. 6c, d). Moreover, loss of LYPD1 led to the inhibition of migration and invasion abilities of HCC cells (Fig. 6e, f). To evaluate the role of LYPD1 in vivo, lentiviruses carrying shRNA targeting at LYPD1 were transfected into Huh7 and MHCC97H cells with verified efficiency (Additional file 8: Figure S4b). Subsequently, subcutaneous implantation experiments were conducted in nude mice. As expected, knockdown of LYPD1 markedly impaired the growth of xenografted tumors (Fig. 6g, h).

Furthermore, bioinformatics analysis was carried out to explore the clinical association of LYPD1. The analysis of HCC cohort from TCGA and three other cohorts from Gene Expression Omnibus (GEO) datasets demonstrated that LYPD1 was up-regulated in tumorous tissues compared with normal tissues (Fig. 6i; Additional file 8: Figure S4c, d). Moreover, the up-regulation of LYPD1 was frequently detected in HCC patients who suffered nodal metastasis or belonged to higher tumor grades/stages (Additional file 8: Figure S4e-g). And pan-cancer analysis manifested that expression of LYPD1 was widely elevated across numerous cancers (Additional file 8: Figure S4h). Besides, Kaplan-Meier analysis implied higher LYPD1 expression correlated with poorer OS and disease-free survival (DFS) in HCC (Fig. 6j). Taken together, LYPD1 is activated during HCC development and promotes the oncogenesis of HCC.

To confirm that the observed phenotypes were mediated by the dysregulation of ALKBH5-LYPD1 axis, we conducted several functional rescue assays. As CCK-8 and colony assays showed, knockdown of ALKBH5 led to the enhanced proliferation capacity in two HCC cells, which could be reverted by LYPD1 silencing (Fig. 7a-d). Knockdown of LYPD1 also significantly abolished the increased mobility ability induced by ALKBH5 loss (Fig. 7e-g). Besides, wound healing assays suggested that the inhibition of ALKBH5 failed to promote cell migration in LYPD1-silenced Huh7 and MHCC97H cells (Fig. 7h). To sum up, dysfunction of LYPD1 may account for the ALKBH5-mediated proliferation or mobility signatures of HCC cells.

To further explore the correlation between expression of ALKBH5 and LYPD1 in HCC tissues, IHC staining of these two proteins were performed on TMA from the second cohort. As expected, approximately 62.2% of specimens with lower expression of ALKBH5 presented stronger LYPD1 staining, while nearly 66.7% of those with higher ALKBH5 expression exhibited weaker LYPD1 dyeing (Fig. 8a, b). Furthermore, the analysis of two independent GEO datasets revealed that ALKBH5 negatively interrelated with LYPD1 in RNA level (Fig. 8c). In conclusion, the expression of ALKBH5 and LYPD1 are inversely correlated in HCC samples.

All data created and analyzed during this current work are involved in this published article (and its supplementary information files) or available on published databases (TCGA or GEO). And GEO accession number of our sequencing is GSE149510↗.