Expression profiles and prognostic significance of RNA N6-methyladenosine-related genes in patients with hepatocellular carcinoma: evidence from independent datasets

Introduction

HCC represents a major public health problem, with an incidence currently rising in almost all countries.¹ It is generally admitted that liver tumorigenesis is a consequence of the multistep process that involves complicated interplays between genetic, epigenetics, and transcriptomic alterations.^2,3 However, till now, the detailed mechanics regulating progression of HCC are still not fully elucidated. The recent discovery of reversible N6-methyladenosine (m6A) of RNA provided a novel layer of post-transcriptional gene regulation.⁴ Modifications of m6A are dynamic and reversible in mammalian cells, which have been recognized as the most abundant internal modification in eukaryotic mRNA. m6A-dependent mRNA regulation is essential in diverse key biological processes, including embryonic development, stem cell differentiation, neural development, and stress responses.^5–8 Consistent with its significant roles, the effects of m6A mRNA methylation on cancer progression are beginning to emerge.^9,10

The m6A modifications are mediated through a complex interplay among “writer”, “eraser”, and “reader” proteins.^4,11 Writers and erasers regulate the abundance, prevalence, and distribution of m6A, whereas readers modulate m6A modification-related functions, such as mRNA splicing, RNA export, RNA stability, and translation. A series of proteins function as writers, erasers or readers have been uncovered, such as methyltransferase-like 3 (METTL3), methyltransferase-like 14 (METTL14), members of YT521-B homology (YTH) domain family, alkB homologue protein family, and so on.^12–15 Although previous studies have reported some m6A-associated proteins were involved in the progression of several cancers,^16,17 the expression pattern as well as the prognosis value have not been fully elucidated in HCC. In this study, we investigated the expression pattern and prognostic value of writer”, “eraser”, and “reader” proteins in patients with HCC through analyzing a HCC cohort of our center, the TCGA, and a Gene Expression Omnibus (GEO) dataset.

Methods

Patients and clinical samples

A total of 162 frozen tumor specimens and the paired para-tumor non-tumor tissues were collected from patients with hepatocellular carcinoma (HCC) who underwent surgery at Guangdong General Hospital (Guangzhou, China) between February January 2013 and May 2014. All fresh tissues were immediately frozen in liquid nitrogen and stored at −80°C until required. Patients were staged according to the seventh edition of the International Union against Cancer TNM classification system. The 162 patients were classified to the Guangdong General Hospital (GDGH) cohort in this study. All protocols were approved by the Research Ethics Committee of Guangdong General Hospital (No. GDREC2015151H(R1)), and written informed consent was obtained from all participants. All study procedures were carried out in accordance with the ethical standards of the Helsinki Declaration.

Data resources and reprocessing

The datasets GSE14520 and GSE63898 were downloaded from expression database GEO (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/).^18,19 GSE14520 included a total of 488 samples, 241 samples were paired non-tumor samples, while the other 247 samples were HCC samples. Platform Information was [HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array for 43 samples, and [HT_HG-U133A] Affymetrix HT Human Genome U133A Array for the other 445 samples. GSE63898 included 228 HCC and 168 cirrhotic samples, and platform was [HG-U219] Affymetrix Human Genome U219 Array for all samples. Limma package in Bioconductor of R language was used to preprocess the normalization of transcriptome data from different platforms. If a number of probe sets correspond to the same gene, then the mean value of probe set was considered as the expression level of this gene. Heatmap was generated by using the pheatmap package of R language. The preprocessed level 3 transcriptome expression data and corresponding clinical information of HCC patients were collected from The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/) (as of September 2018). The gene expression profiling data of 374 HCC tumor samples and 50 normal samples were downloaded from the TCGA database, and a total of 370 HCC patients with follow-up information were enrolled in the survival analysis.

RNA isolation and qRT-PCR

Total RNA was extracted from tissue samples using TRIzol (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions. Total RNA was converted to complementary DNA by reverse-transcription using oligodT primers and SuperScript II reverse transcriptase (Invitrogen). For qRT-PCR, three replicates per sample were amplified and analyzed. Reactions were carried out in a 20-μL volume using SYBR Green Reaction Mix (Qiagen Science, GER) with 0.5 mM primer. To detect the mRNA expression levels in the tissues, GAPDH was used as the normalization control. Primers aree listed in Table S1. The relative gene expression levels were determined using the comparative threshold cycle (2^−ΔΔCT) method.

Immunohistochemistry

Tissue samples were fixed with 10% formaldehyde in PBS, embedded in paraffin, and cut into serial 4-μm sections. The tissue sections were deparaffiized in xylene and rehydrated with graded ethanol and then treated with 0.3% hydrogen peroxide for 10 mins, followed by antigen retrieval in 0.01 M citrate buffer (p H: 6.0, 10×). After 20 min cooling, sections were incubated with the primariy METTL3 antibody (rabbit monoclonal; no. ab195352, Abcam Inc., USA) or primarily YTHDF1 antibody (rabbit monoclonal; no. ab230330, Abcam Inc., USA) at 4°C overnight. The sections were then incubated with the secondary antibody and then were visualized.

Gene set enrichment analysis (GSEA)

GSEA was performed using normalized data by GSEAv3.0 tool (http://software.broadinstitute.org/gsea/index.jsp).^20,21 To explore the differences in potential biological functions in the low- and high-expression sets of prognostic genes, GSEA was used using the Molecular Signatures Database (MSigDB) of KEGG gene sets (c2.cp.kegg.v6.2.symbols).

Statistics

All quantitative data are presented as the mean ± standard deviation from at least three independent experiments. Unless otherwise noted, continuous variables were analyzed using the Student’s t-test, and the Mann–Whitney U test was used for independent samples when the population could not be assumed to be normally distributed. Kaplan–Meier curves of overall survival (OS) were compared with the log-rank test. Associations between the variables and survival were also evaluated by using univariate and multivariate analyses with the Cox proportional hazard model. All tests were two-sided, and p-values <0.05 were considered statistically significant. All statistical analyses were conducted using SPSS 17.0 (SPSS Inc., Chicago, IL, USA).

Results

Identification of differentially expressed m6A-related genes in hepatocellular carcinoma

To acquire the transcriptional profiles of m6A-related genes in HCC, HCC tissues (n=162) and paired adjacent non-tumor tissues (n=162) were analyzed using qRT-PCR. We found that two m6A writers (METTL3, WTAP) and six readers (EIF3, YTHDC1, YTHDF1, YTHDF2, YTHDF3, and HNRNPA2B1) were significantly upregulated in cancer tissues compared with paired non-tumor tissues (Figure 1A). To verify our result, the expression of m6A-related genes were also examined in the TCGA dataset and two GEO datasets (GSE14520, GSE63898). The heatmaps of genes in TCGA, GSE63898, and GSE14520 are shown in Figure 1B–D, respectively. As summarized in Figure 1E, a total of five genes, including METTL3, YTHDF1, YTHDF2, YTHDF3, and EIF3, showed consistent trends of upregulation in HCC tissues across all four datasets.

Figure 1 Expression profiles of m6A-related genes in four independent datasets. (A) Expression levels of twelve m6A-related genes in 162 paired HCC tissues and corresponding adjacent non-tumor tissues (GDGH cohort) was examined via quantitative reverse transcription-polymerase chain reaction (qRT-PCR). GAPDH was used as internal control. Relative gene expression was determined using the comparative delta-delta CT method, and data are presented as △△Ct. (B, C, D) Expression heatmap plotting of m6A-related genes in TCGA dataset (B), GSE63898 dataset (C), and GSE14520 dataset (D). (E) The change profiles of m6A-related genes in the four datasets are summarized and compared. The red block indicates the corresponding gene was significantly upregulated in HCC tissues compared with the non-tumor control tissues; the green red block indicates the corresponding gene was significantly downregulated in HCC tissues compared with the non-tumor control tissues; the black red block indicates the expression of corresponding gene was not significantly changed in HCC tissues compared with the non-tumor control tissues; the grey block indicates the expression level of corresponding gene was not available.

Prognosis significance of m6A related genes in hepatocellular carcinoma

To further investigate the clinical significance of m6A-related genes, we assessed survival rates in three independent datasets with prognosis information (TCGA, GSE14520, and GDGH cohort). The clinical characteristics of patients in each cohort are summarized in Table 1. In each dataset, patients were divided into two groups by the median value of the candidate gene, and OS curves for each gene were plotted by the Kaplan–Meier method. Analysis of these datasets revealed that among all genes examined (Figures S1, S2, S3), only two genes (METTL3 and YTHDF1) showed significant associations with survival throughout all three datasets. As shown in Figure 2, upregulation of METTL3 and YTHDF1 was correlated with shorter OS rates in the HCC patients of all datasets. The protein expression of METTL3 and YTHDF1 was further validated in HCC tissues of 10 patients by immunohistochemical staining. All patients showed positive staining of both METTL3 and YTHDF1. The representative images are shown in Figure 3. We found that METTL3 was mainly expressed in nuclear (Figure 3A,B) and YTHDF1 was mainly expressed in cytoplasm (Figure 3C,D).

Table 1 Clinical characteristics of patients in the TCGA, GDGH, and GSE14520 cohorts

Figure 2 Correlation between METTL3 or YTHDF1 expressions and survival rates in patients with HCC. (A, B) Kaplan–Meier OS curve based on METTL3 expression (A) and YTHDF1 expression (B) in TCGA dataset. (C, D) Kaplan–Meier OS curve based on METTL3 expression (C) and YTHDF1 expression (D) in GDGH dataset. (E, F) Kaplan–Meier OS curve based on METTL3 expression (E) and YTHDF1 expression (F) in GSE14520 dataset. Patients were grouped into low- vs high-expression levels based on the median values of METTL3 or YTHDF1 expressions in each dataset. Log-rank test was used to compare differences in survival times.

Figure 3 Immunohistochemical staining of METTL3 (A, B) and YTHDF1 (C, D) in HCC tissues Scale bars: 200 μm.

Cox regression analysis for the prognosis significance of METTL3 and YTHDF1

Since survival curves showed that both METTL3 and YTHDF1 were correlated with survival, we performed multivariate Cox proportional-hazards analyses to further investigate the prognostic value of METTL3 and YTHDF1. According to the results of multivariate Cox regression in the GDGH cohort, METTL3 expression (HR =1.78, 95%CI =1.09–2.91, P=0.021) and YTHDF1 expression (HR =1.53, 95%CI =1.08–2.16, P=0.017) were both independent prognostic factors for patients with HCC (Table 2).

Table 2 Univariate and multivariate analyses of factors associated with survival

The two-gene signature predicts the survival of patients with HCC

Due to the importance of m6A modifications in biological processes, as well as the significant prognostic value of METTL3 and YTHDF1 in HCC, we developed a risk stratification method based on these two genes. Patients were divided into three groups based on median expression levels of METTL3 and YTHDF1 (Figure 4). Kaplan–Meier curves showed that patients in low-risk groups had significantly increased survival than those in other risk groups (log-rank test, P=0.0002 in TCGA, P=0.0249 in GSE14520, P＜0.0001 in GDGH cohort) (Figure 4A–C). Overall, after we integrated all patients in one survival analysis, we observed that patients in different groups had distinctively different outcomes (P＜0.0001) (Figure 4D).

Figure 4 Kaplan–Meier OS curves of HCC patients stratified into low-, medium-, and high-risk groups in GDGH dataset (A), TCGA dataset (B), GSE14520 dataset (C), and the integration of three databases (D). Low-risk group was defined as patients with both low METTL3 expression (<50% patients) and low YTHDF1 expression (<50% patients); medium-risk group was defined as patients with either high METTL3 expression (≥50% patients) or high YTHDF1 expression (≥50% patients). High-risk group was defined as patients with both high METTL3 expression (≥50% patients) and high YTHDF1 (≥50% patients) expression. Log-rank test was used to compare differences in survival times.

GSEA and functional annotation of METTL3 and YTHDF1

To explore the potential function of METTL3 and YTHDF1 in HCC, we applied GSEA on the TCGA and GSE14520 datasets using the KEGG gene sets (c2.cp.kegg.v6.2.symbols). The results of GSEA are listed in Supplementary files S1-S4. By using the P-value ＜ 0.05 and the FDR q-value ＜ 0.05 as filter condition, the most enriched KEGG pathways were screened. As shown in Figure 5A, seven enriched KEGG pathways according to METTL3 expression appeared in GSEA result based on TCGA dataset and that based on GSE14520. The seven pathways are summarized in Table 3, which were spliceosome, cell cycle, DNA replication, basal transcription factors, homologous recombination, base excision repair, and RNA degradation. Meanwhile, eight enriched KEGG pathways associated with high YTHDF1 were also screened (Table 4). The eight pathways were associated with ribosome, cell cycle, spliceosome, DNA replication, base excision repair, RNA degradation, nucleotide excision repair, and mismatch repair. The GSEA plots of the above pathways are attached in Figure S4.

Table 3 Gene set enrichment analysis for METTL3 with KEGG pathway gene sets gene symbols using the TCGA and GSE14520 datasets

Table 4 Gene set enrichment analysis for YTHDF1 with KEGG pathway gene sets gene symbols using the TCGA and GSE14520 datasets

Figure 5 GSEA results of METTL3 and YTHDF1 in TCGA and GSE14520 patients. (A) Flow chart showing the screening process to identify KEGG pathways associated with METTL3 expression, YTHDF1 expression, or both METTL3 and YTHDF1 expression. (B) GSEA plots displayed the significant association between five KEGG pathways, including spliceosome, cell cycle, DNA replication, base excision repair, and RNA degradation, and both METTL3/YTHDF1 expression.

Because we found that patients with high expression of both METTL3 and YTHDF1 had the worst prognosis, we further screened the enriched KEGG pathways associated with both METTL3 expression and YTHDF1 expression (Figure 5A), they were spliceosome, cell cycle, DNA replication, base excision repair, and RNA degradation. The GSEA plots are shown in Figure 5B. This result indicated that both METTL3 and YTHDF1 play important roles in the above five biological processes.

Discussion

RNA can be modified by more than one hundred distinct post-transcriptional modifications, and m6A is now recognized the most abundant internal modification in eukaryocyte since it was discovered in the 1970s.^22,23 However, due to the lack of effective techniques, functional characterizations of m6A was just beginning to emerge since the antibody-based high-throughput sequencing technology appeared.²⁴ Recent work has uncovered that the m6A modifications were controlled by the complex interplay of m6A deposition, removal, and recognition factors (writers”, “erasers”, and “readers”). Although the role and mechanism of these factors in regulating m6A were conservative, their pathological implications can be different depending on the disease. For example, a major m6A “writer”, METTL3, is significantly up-regulated in HCC and serve as an oncogene,¹⁷ but a recent study revealed that knockdown of METTL3 in endometrial cancer cell lines promoted cell proliferation, anchorage-independent growth, colony formation, migration, and invasion,¹ which was totally opposite to its role in HCC.

To explore the role of m6A-related factors in HCC, we evaluated their expression pattern and prognosis significance in four HCC cohorts. By comparing results from TCGA data, two GEO datasets, and our HCC cohort, we found that the expression trends of five genes (METTL3, YTHDF1-3, EIF3) were consistent throughout all datasets. We also identified that overexpression of METTL3 or YTHDF1, but not other genes, was associated with poor prognosis in patients with HCC. Moreover, patients with co-overexpression of METTL3 and YTHDF1 showed the poorest prognosis. These findings suggested important roles of METTL3 and YTHDF1 as oncogenes in the progression of HCC.

METTL3 is a key component of m6A methyltransferase complex, and it determines the levels and distribution of target-specific m6A modifications through a complex interplay with other components.²⁵ A recent study reported that knockout of METTL3 remarkably suppressed HCC tumorigenesis and progression in vitro and in vivo via a YTHDF2-dependent mechanism, suggesting a close collaboration between m6A “writers” and “readers” in HCC.² The YT521-B homology (YTH) domain family of proteins (YTHDF1, YTHDF2, YTHDF3, YTHDC1, and YTHDC2) are direct readers of m6A and have a conserved m6A-binding pocket.^26–28 Although YTHDF1, YTHDF2, YTHDF3, YTHDC1, and YTHDC2 all belong to YTH domain family, their functions are different. For example, YTHDC1-2 locate in the nuclear compartment and YTHDF1-3 locate in the cytoplasmic compartment.^29,30 The nuclear reader YTHDC1 regulates alternative splicing by binding m6A directly and recruiting the splicing factors,²⁹ while cytoplasmic reader YTHDF1 mediates translation initiation of m6A-containing transcripts,²⁶ and YTHDF2 promotes mRNA decay by facilitating the recruitment of the CCR4–NOT complex to target mRNA.³¹ Intriguingly, our present study showed that all of the three cytoplasmic m6A readers (YTHDF1, YTHDF2, and YTHDF3) were overexpressed in HCC patients throughout the four datasets. Therefore, overexpression of YTHDF1, YTHDF2, and YTHDF3 might be a common feature of the expression profile of m6A-related genes in HCC. Moreover, among the overexpressed m6A readers, YTHDF1 showed a significant and consistent prognostic value among all four datasets, suggesting a potential important role of YTHDF1 in HCC progression. Further investigations will be required to address the role of YTHDF1 in HCC.

All together, we reported that both METTL3 and YTHDF1 were independent prognostic factors in patients with HCC. Moreover, patients with different co-expression levels of METTL3 and YTHDF1 showed diverse survival outcomes. Future studies could focus on the molecular mechanisms of the potential interplay between METTL3 and YTHDF1 in HCC progression.