This article is mentioned in:

Introduction

Liver cancer is one of the most common malignant tumor types worldwide, with its mortality rate ranking second among malignant tumors, and its incidence rate increasing (1). Every year, liver cancer cases in China account for >50% of worldwide cases (2) and the median survival time of patients with late liver cancer is <5 months (3). Although the pathogenesis of liver cancer has been studied in detail, identifying therapeutic targets for liver cancer remains urgent.

RNA N6-methyladenosine (m6A) methylation modification is a common modification in mRNA and noncoding RNA (ncRNA). It affects RNA cutting, translation, stability and epigenetic effects in certain ncRNAs. Studies have indicated that m6A gene modification can regulate the occurrence and development of liver cancer. For instance, it has been shown that methyltransferase 14, N6-adenosine-methyltransferase subunit (METTL14) upregulates ubiquitin specific peptidase 48 to alleviate liver cancer by regulating Sirtuin 6 stability (4). In addition, it was found that, compared to the normal group, the expression level of vir like m6A methyltransferase associated in the liver cancer group was also higher, and that KIAA1429 was modified by reinforcing ID2 m6A-modified ID2 mRNA expression, promoting the invasion and metastasis of tumor cells (5). It has also been reported that m6A methylation can modify long (l)ncRNA and affect its stability and metabolism (6). To date, only a small number of studies have explored how m6A methylation-modified lncRNA can promote the occurrence and development of liver cancer. Therefore, further study on m6A methylation-modified lncRNA can provide new ideas for the study of lncRNA.

In the present study, based on the The Cancer Genome Atlas (TCGA) dataset, bioinformatics and statistical analyses were conducted on the relevant data of cancer patients to determine the prognostic significance of m6A-related lncRNAs. Using bioinformatics analysis, a risk model consisting of nine genes was established. RP1-228H13.5/hsa microRNA (miR)-205/zinc finger protein interacting with K protein 1 (ZIK1) was selected from the constructed competitive endogenous RNA (ceRNA) network for subsequent research. Finally, the specific effects and molecular mechanisms of RP1-228H13.5 on liver cancer cells were studied using a cell assay.

Materials and methods

Data collection

The search term ‘hepatocellular carcinoma’ was used in the TCGA database (https://portal.gdc.cancer.gov) to obtain relevant RNA expression profile data, and the corresponding clinicopathological data were downloaded from the cBioPortal website (https://www.cbioportal.org/), including patient ID, age, survival time, survival status, tumor grade, stage and TNM stage. The deadline for downloading the data was August 9, 2021. In the end, data pertaining to 374 liver cancer samples and 50 control samples were obtained. Samples with incomplete clinical data were excluded from the above samples. Finally, 341 patients with liver cancer samples and 40 normal controls were included in the present study. By reviewing m6A-related literature (7–9), a total of 24 m6A-related genes were extracted from RNA expression profiles, including writers (METTL3, METTL14, METTL16, WTAP, RBM15, RBM15B, ZCCHC4, CBLL1 and ZC3H13), erasers (FTO and ALKBH5) and readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1, RBMX, LRPPRC, FMR1, IGF2BP3, IGF2BP1 and IGF2BP2). Differentially expressed genes (DEGs) between the normal and cancer groups were identified using the DESeq2 package (DOI: 10.18129/B9.bioc.DESeq2), in which P<0.001 indicated extremely significant, P<0.01 very significant and P<0.05 significant. The ‘Pheatmap’ and ‘ggplot 2’ software packages in R language (version 4.3.0; https://r-charts.com/correlation/pheatmap/; R Foundation for Statistical Computing) were used to draw the expression heat and violin maps of DEGs. The lncRNA annotation file was downloaded from the GENCODE website (https://www.gencodegenes.org/) to identify and annotate the lncRNAs in the TCGA dataset.

Correlation analysis between m6A genes and m6A gene-related lncRNA

Pearson's correlation coefficient was used to analyze the correlation between m6A genes and lncRNAs. LncRNAs with an absolute correlation coefficient R>0.5 and P<0.001 were considered to be m6A-related lncRNAs.

Bioinformatics analysis

In order to determine the prognostic value of m6A-related lncRNAs, univariate Cox regression analysis was performed on m6A-related lncRNAs. First, the original data of m6A-related lncRNAs were standardized in the form of log2(x+1). Univariate Cox analysis was performed on m6A-related lncRNAs using the ‘survival’ package in R language to screen prognosis-related lncRNAs. Next, the R software package ‘glment’ was used for Cox regression analysis of Least Absolute Shrinkage and Selection Operator (LASSO), a Cox proportional risk model was built and the risk score was calculated. The calculation formula of the risk score was as follows:

where Coefi refers to the coefficients and Xi is the FPKM value of each m6A-related lncRNA. The risk score of all patients was calculated, the median value of the risk score was selected as the critical value and patients were divided into high- and low-risk groups. The ‘survival’ package was further used to compare the overall survival (OS) rate between high- and low-risk groups through Kaplan-Meier analysis, and the survival curves were drawn. The ‘Survival Receiver Operating Characteristic (ROC)’ package was used to evaluate the predictive value of the LASSO regression Cox proportional hazards construction model through the area under the ROC curve (AUC). In order to study the difference between risk grouping and clinical characteristics, a χ2 test was performed, the difference in clinical characteristics between the high- and low-risk groups was evaluated and a heat map was drawn using the ‘pheatmap’ software package. Univariate and multivariate Cox analysis were performed using the ‘survival’ package to select independent predictors of OS, and visualized through the ‘forest’ software package.

ceRNA network and functional enrichment analysis of genes

The targeted miRNAs of m6A-related lncRNAs were predicted using the Perl programming language through the miRcode database (http://www.mircode.org/index.php miRcode) and the targeted mRNAs of these miRNAs were predicted using the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), miRDB (http://www.mirdb.org/) and TargetScan (http://www.targetscan) databases. The ceRNA network was constructed using ‘Cytoscape’ software (http://www.cytoscape.org/) and the results were visualized. The patients were divided into high- and low-risk groups using the median risk score. The differential gene expression of the high-risk and low-risk groups was analyzed. DEGs were screened using P<0.05, and |log2 fold change (FC)|>2 was used to perform the functional and pathway enrichment analysis of DEGs in the high- and low-risk groups, and 201 target mRNAs in the ceRNA network. Log2FC refers to the difference in gene expression between two samples. This is a commonly used concept in bioinformatics analysis. Specifically, log2FC is an indicator that measures the magnitude of changes in gene expression levels. By calculating the expression ratio between two samples (such as the expression level of sample 1/sample 2) and then taking logarithms (based on 2), the value of log2FC can be obtained. Compared to directly calculating the difference in gene expression levels, using log2FC can more accurately reflect the degree of difference between gene expression levels. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis is a bioinformatics method used to explain and understand the function of genes or proteins and the metabolic pathways involved. Gene ontology (GO) enrichment analysis is a statistical method widely used in bioinformatics research, mainly used to test the degree of enrichment of certain functions or features in a gene set. Using KEGG and GO enrichment analysis to analyze differentially expressed genes.

Cell culture and cell transfection

The THLE-2 human normal liver cell line and the SNU-398, Huh-7 and Hep3B2.1–7 human liver cancer cell lines were purchased from the American Type Culture Collection and cultured in DMEM (Gibco; Thermo Fisher Scientific, Inc.) containing 10% FBS (Gibco; Thermo Fisher Scientific, Inc.). The cells were cultured at 37°C in a saturated humidified atmosphere with 5% CO2. The day before transfection, the cells were seeded into 96-well plates with an inoculation density of 1×106 cells/well. When the cell density reached 50–80%, short hairpin RNA (sh)-RP1-228H13.5 (sequence, 5′-CAGAGAAAGGTGACTGAAGGA-3′), overexpression vector (OV)-RP1-228H13.5, mimics-negative control (NC) (sequence, sense: 5′-UUGUACUACACAAAAGUACUG-3′; antisense: 5′-GUACUUUUGUGUAGUACAAUU-3′) [corresponding sh-NC (5′-CAGUACUUUUGUGUAGUACAA-3′) and OV-NC], miR-205-5p mimics (sequence, sense: 5′-UCCUUCAUUCCACCGGAGUCUG-3′; antisense: 5′-GACUCCGGUGGAAUGAAGGAUU-3′), inhibitor-NC (sequence, 5′-CAGUACUUUUGUGUAGUACAA-3′), miR-205-5p inhibitor (sequence, 5′-CAGACUCCGGUGGAAUGAAGGA-3′), Sh-METTL14 (sequence, 5′-ATGGATAGCCGCTTGCAGGAG-3′) and OV-METTL14 [corresponding Sh-NC (5′-CAGUACUUUUGUGUAGUACAA-3′) and OV-NC] were transfected into the cells using Lipofectamine™ 2000 (Invitrogen; Thermo Fisher Scientific, Inc.), according to the manufacturer's instructions. To achieve overexpression of RP1-228H13.5 and METTL14, pcDNA3.1-RP1-228H13.5 and pcDNA3.1-METTL14 vectors were designed and constructed. The final transfection concentration of the vector in a 6-well plate was 2 µg/ml. The above-mentioned oligonucleotide primers were synthesized by Shenggong Corporation. pcDNA3.1™ plasmid was purchased from Invitrogen (Thermo Fisher Scientific, Inc.).

Reverse transcription-quantitative PCR (RT-qPCR)

The experiment was performed according to the manufacturer's instructions. Total RNA was extracted by the TRIzol® one-step method and 2 µg of total RNA was reverse-transcribed into cDNA using the PrimeScript RT Master Mix RT Kit (Takara Biotechnology Co., Ltd.). A final reaction volume of 20 µl was established in accordance with the SYBR Green PCR Master Mix (Thermo Fisher Scientific, Inc.). The PCR thermocycling parameters were set to 95°C for 5 min for initial denaturation, followed by 35 cycles of denaturation at 94°C for 30 sec, annealing at 60°C for 30 sec and extension at 72°C for 30 sec. Primers were designed by Primer Premier Software 5.0 (Premier Biosoft International) and synthesized by Invitrogen (Thermo Fisher Scientific, Inc.). Real-time qPCR was performed on an ABI 7500 system (Thermo Fisher Scientific, Inc.). GAPDH was used as the internal reference in the experiment. The 2−ΔΔCq method (10) was employed to assess relative expression levels. The primer sequences used are listed in Table I.

Table I. Primer sequences.

Table I.

Primer sequences.

Gene	Sense	Antisense
RP1-228H13.5	ACTCAGGCAGGGAAATACGG	CTGCCCATCAGTCATACAGGA
miR-205-5p	TCCTTCATTCCACCGGAGTCTG	GTGCGTGTCGTGGAGTCG
ZIK1	GCCCCGACTCAGGTTACTGT	CCCACTCGTCCTGTGAGAAG
METTL14	GTAGCACAGACGGGGACTTC	GCCAGCCTGGTCGAATTGTA
GAPDH	TCATGACCACAGTCCATGCC	TTCTAGACGGCAGGTCAGGT
U6	CTCGCTTCGGCAGCACA	AACGCTTCACGAATTTGCGT
18S	GTGGAGCGATTTGTCTGGTT	AACGCCACTTGTCCCTCTAA

[i] METTL14, methyltransferase 14, N6-adenosine-methyltransferase subunit; ZIK1, zinc finger protein interacting with K protein 1.

Cell Counting Kit-8 (CCK-8) assay

Cell proliferative activity was examined using a CCK-8 kit (Dojindo Molecular Technologies, Inc.). One day before the experiment, the cells in the logarithmic growth phase were inoculated into a 96-well plate at a density of 1,000 cells/well and then cultured in an incubator at 37°C. The next day, the cells were transfected with Sh-RP1-228H13.5 and OV-RP1-228H13.5 (NC). Following transfection, the cells were cultured for 24, 48 or 72 h, CCK-8 reagent (10 µl/well; Thermo Fisher Scientific, Inc.) was added to the culture wells, followed by incubation for 2 h at 37°C. The absorbance at 450 nm was read by a microplate reader (Bio-Rad Laboratories, Inc.). A cell proliferation curve was then drawn.

Wound-healing assay

The transfected cells were inoculated in 6-well plates at a density of 1×106/well. Serum-free culture medium was used for the wound-healing assay. Using a 200-µl pipette tip, three scratches perpendicular to the marking line were drawn in the cell culture plate. The floating cells were then washed away with PBS. Three scratch locations were selected and photographed at 0 and 24 h to record cell migration, which was quantified as follows: Migration rate=(0 h wound distance-24 h wound distance)/0 h wound distance ×100%.

Transwell assay

Matrigel® (BD Biosciences) was added into the upper chambers (8 µm pore size) of Transwell inserts (Thermo Fisher Scientific, Inc.) and incubated in a cell incubator for ~2 h for the cell invasion assay. The cells were transfected for 24 h, digested with trypsin (Thermo Fisher Scientific, Inc.) and resuspended with 400 µl serum-free cell culture medium. Cell suspension (5×103 cells in 100 µl) was inoculated into the upper Transwell chambers and 600 µl DMEM with 10% serum was added to the lower chambers. After 24 h of culture, the medium in the upper chamber was removed and the cells in the upper chamber that had not transgressed to the lower side of the membrane were gently wiped off using defatted cotton. Next, the chamber was placed in 4% paraformaldehyde solution, fixed at room temperature for 20 min and then stained with 0.1% crystal violet solution at room temperature for 15 min. A total of five visual fields were randomly observed and photographed by Olympus CX43 Fluorescence Microscope and the number of invaded cells was calculated. In the cell migration assay, no Matrigel® was added to the upper chamber of the Transwell, while the other steps were the same as those in the cell invasion assay.

Flow cytometry

After 48 h of cell transfection, the cells were collected and rinsed with PBS twice. They were then mixed with 500 µl of pre-cooled 1X binding buffer and 5 µl Annexin-V-FITC (cat. no. 88-8102-72; Invitrogen; Thermo Fisher Scientific, Inc.), incubated at room temperature for 15 min, and 2.5 µl propidium iodide was added for 5 min at room temperature. The apoptotic rate was detected using flow cytometry (FACSCalibur; BD Biosciences) and anlysed by BD CellQuest Pro software (version 5.1). The experiment was repeated three times.

Fluorescence in situ hybridization (FISH)

The coverslips which coated with poly-L-lysin were placed on a 24-well plate and an appropriate number of cells were cultured (~6×10−4/well). Next, 4% paraformaldehyde was used to incubated cells for 10 min at 4°C. Pre-cooled (4°C) permeabilization solution (PBS containing 0.5% Triton-X100) was added to each well, incubated at 4°C for 5 min and then discarded. Pre-hybridization solution (200 µl) was added to each well and sealed at 37°C for 30 min. Under dark conditions, 2.5 µl 20 µM RP1-228H13.5 probe mix storage solution (probe sequence: 5′-TGAATTATCTAAAGTAAAGCTATAAATCCAGAATATATACCAGTATCTCTAGAATTTTTCATATTACTACATTTACTTAA-3′) or internal reference FISH Probe Mix storage solution was added to 100 µl of hybridization solution. Probes were synthesized by Shenggong Corporation and the in situ hybridization kit was purchased from Beyotime Institute of Biotechnology. The pre-hybridization solution was discarded from each well, 100 µl probe hybridization solution containing probes was added and hybridization was performed overnight at 37°C while avoiding light. Subsequently, while avoiding light, each well was washed three times at 42°C with hybrid washing solution I for 5 min each time to reduce the background signal. Subsequently, avoiding light, cells were washed once with hybrid detergent II and once with hybrid washing solution III, both at 42°C. Next, avoiding light, cells were washed with 1X PBS at room temperature for 5 min. Staining with DAPI solution was then performed for 10 min at room temperature. Under dark conditions, the coverslips were carefully removed from the well and fixed on a slide with a sealing agent (Beyotime Institute of Biotechnology) for fluorescence detection at room temperature.

Dual-luciferase assay

Starbase 2.0 (http://starbase.sysu.edu.cn/) was used to predict the binding sites of RP1-228H13.5 and miR-205-5p, as well as those of miR-205-5p and ZIK1. The construction of reporter vector pmirGLO-negative control (NC), pmirGLO-RP1-228H13.5-wild-type (wt), pmirGLO-RP1-228H13.5-mutant (mut), pmirGLO-ZIK1-wt and pmirGLO-ZIK1-3′-untranslated region mut (complementary mutations were performed on the sequence combining ZIK1 with miR-205-5p; ATGAAGGA mutates to UACUUCCU), or synthesis of miR-NC mimics or miR-205-5p mimics were performed by Sangon Biotech Co., Ltd. The pmirGLO vector was from Promega Corp. Upon cell confluence reaching 70–80%, each reporter plasmid with miR-NC mimics or miR-205-5p mimics were co-transfected into 293T cells (cat. no. ACS-4500; American Type Culture Collection) using Lipofectamine® 2000 (Invitrogen; Thermo Fisher Scientific, Inc.). At 24 h after transfection, the relative firefly luciferase activity was detected using a Dual-Glo Luciferase Reporter assay system (Promega Corp.).

RNA pull-down assay

RNA pull-down assays were performed using a Pierce Magnetic RNA-Protein Pull-Down kit (cat. no. 20164; Thermo Fisher Scientific, Inc.), according to the manufacturer's instructions. In brief, 1×107 cells were treated with 500 µl RIPA buffer (Thermo Fisher Scientific, Inc.), centrifuged at 4°C for 15 min at 10,000 × g, and the lysates in the supernatant were collected. Subsequently, biotinylated miR-205-5p probe (5′-UCCUUCAUUCCACCGGAGUCUG-biotin-3′) and the control probe (5′-CUCUGAGGCCACCUUACUUCCU-biotin-3′), synthesized by Sangon Biotech Co., Ltd., along with streptavidin magnetic beads, were incubated at 4°C for 2 h with 500 µl cell lysate. Subsequently, the mixture was washed with wash buffer twice and the eluate was then analyzed using RT-qPCR as indicated above.

For RNA chromatin immunoprecipitation, the supernatant of the lysates was collected and IgG or METTL14 antibody was added from the Magna RIP Kit (cat. no. 17-700; Sigma-Aldrich; Merck KGaA). The immunoprecipitation of METTL14-associated long intergenic RNA was verified by RT-qPCR.

Western blot analysis

Cells were inoculated in a 6-well plate, transfected with plasmid and placed in a cell culture incubator with a constant temperature and humidity for 48 h, and cells were then collected. Cells were lysed with radioimunoprecipitation assay buffer (Beyotime Institute of Biotechnology) and a BCA kit (Beyotime Institute of Biotechnology) was used to determine the protein concentration, followed by denaturation, electrophoresis separation, membrane transfer, washing and sealing (11). Next, a ZIK1 primary antibody (dilution, 1:1,000; cat. no. bs-19266R; Shanghai Bioplus Biotech Co., Ltd.) was added for overnight incubation at 4°C, followed by the addition of a secondary antibody (dilution, 1:1,000; cat. no. bs-0295G-HRP; BIOSS Antibodies Co., Ltd.). The blot was finally developed using enhanced chemiluminescence reagents (Amersham; Cytiva) and GAPDH (dilution, 1:10,000; cat. no. ab181602; Abcam) was used as the internal reference.

Statistical analysis

R software (version 4.3.0) was used for statistical analysis and continuous variables were expressed as the mean ± standard deviation or median ± interquartile range. Fisher's exact test or the χ2 test were used to test differences between groups. Pearson's correlation coefficient was calculated to determine the correlation between parameters. Kaplan-Meier curves were drawn and the log-rank test was used to analyze the survival rate. Multivariate Cox regression analysis was used to determine the independent prognostic factors of OS. GraphPad Prism 5.0 software (GraphPad Software; Dotmatics) and SPSS 18.0 statistical software (SPSS, Inc.) were used for statistical data analysis. One-way ANOVA was used to compare the mean differences among three or more groups and the Student-Newman-Keul's test was used as a post-hoc test, while unpaired t-tests were used to calculate the mean differences between two groups. P<0.05 was considered to indicate a statistically significant difference.

Results

The 24 m6A regulatory genes and related lncRNAs identified in patients with liver cancer

The expression profile data of 24 m6A-related genes were extracted from the TCGA dataset, of which 5 (RBMX, RBM15B, IGF2BP1, METTL3 and YTHDF1) were significantly upregulated in tumor tissues and 11 (ZCCHC4, METTL14, ZC3H13, ALKBH5, RBM15, YTHDC1, YTHDF2, YTHDC2, YTHDF3, METTL16 and FTO) were significantly downregulated in tumor tissues (Fig. 1A and B). The lncRNA annotation file was downloaded from the GENCODE website and 14,816 lncRNAs were identified from the TCGA dataset. Through Pearson's correlation analysis, using |Pearson R|>0.5 and P<0.001 as the threshold, the co-expressed lncRNA and mRNA pairs were screened and the co-expression relationship pairs containing 22 m6A genes and 287 lncRNAs were finally obtained.

Figure 1. Prognostic analysis of m6A methylation regulatory factor-related lncRNAs. (A) Heatmap of the correlations between m6A-related genes and the 24 prognostic m6A-related lncRNAs. (B) Violin diagram of mRNA expression distribution of 24 RNA m6A methylation regulators in liver cancer tumor tissues and adjacent normal tissues. Blue represents normal tissues and red represents tumor tissues. (C) Heat map of the association between m6A-related genes and m6A-related lncRNAs. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. lncRNA, long non-coding RNA; m6A, N6-methyladenosine.

Univariate Cox regression analysis of m6A gene-related lncRNA expression and prognosis

Univariate Cox regression was used to screen out prognosis-related lncRNAs from 287 m6A-related lncRNAs (P<0.05). It was found that 89 m6A-related lncRNAs were significantly associated with the prognosis of patients with liver cancer. The association between 24 m6A-related genes and prognosis-related lncRNAs in the TCGA dataset is presented in Fig. 1C (the figure only shows 20 lncRNAs). The results of the univariate Cox analysis of m6A-related lncRNAs are provided in Table II.

Table II. N6-methyladenosine modification-associated genes related prognostic long non-coding RNAs.

Table II.

N6-methyladenosine modification-associated genes related prognostic long non-coding RNAs.

Gene	Hazard ratio	95% CI	P-value
RP4-681N20.5	3.175	1.310–7.697	0.011
WAC-AS1	1.093	1.045–1.144	<0.001
BX322557.10	1.642	1.125–2.396	0.010
RP11-73M18.10	1.789	1.173–2.728	0.007
STX18-AS1	6.91	1.203–39.699	0.030
RP11-527J8.1	2.013	1.350–3.001	0.001
GS1-279B7.2	2.463	1.090–5.563	0.030
CTD-2537I9.12	1.116	1.012–1.231	0.028
AC092614.2	2.626	1.039–6.639	0.041
RP11-332H14.2	2.562	1.354–4.848	0.004
RP11-815I9.4	12.405	2.714–56.702	0.001
CTD-2510F5.4	1.149	1.090–1.212	<0.001
CTD-3099C6.9	1.71	1.015–2.883	0.044
RP11-884K10.7	1.385	1.083–1.770	0.009
BMS1P20	1.793	1.231–2.611	0.002
LENG8-AS1	1.265	1.025–1.562	0.029
ZEB1-AS1	1.442	1.143–1.820	0.002
RP11-435O5.2	1.541	1.172–2.026	0.002
AC073254.1	4.077	1.189–13.978	0.025
CTD-2012J19.3	16.19	4.298–60.981	<0.001
AP001258.4	1.428	1.005–2.030	0.047
AP001469.7	3.064	1.080–8.691	0.035
RP11-480A16.1	2.511	1.022–6.168	0.045
PKD1P6	3.513	1.017–12.143	0.047
CTD-2528L19.6	1.411	1.091–1.826	0.009
AC025335.1	1.712	1.257–2.331	0.001
MCM3AP-AS1	8.044	2.296–28.178	0.001
RP11-261C10.5	5.215	1.262–21.550	0.023
RP1-257A7.4	3.588	1.331–9.673	0.012
SNHG1	1.043	1.015–1.071	0.002
RP11-355O1.11	2.418	1.154–5.068	0.019
RP11-968A15.2	1.387	1.104–1.743	0.005
RP11-15N24.4	5.639	1.172–27.13	0.031
CAPN10-AS1	1.656	1.014–2.706	0.044
LINC00665	1.123	1.042–1.211	0.003
NUTM2B-AS1	2.326	1.178–4.596	0.015
RP11-119F7.5	1.347	1.006–1.803	0.045
AC074117.10	1.411	1.147–1.736	0.001
RP11-1246C19.1	1.438	1.045–1.977	0.026
NCK1-AS1	1.467	1.180–1.823	0.001
RP1-92O14.6	2.485	1.524–4.055	<0.001
RP11-35G9.3	2.145	1.517–3.033	<0.001
RP11-147L13.13	1.177	1.054–1.314	0.004
SNHG21	1.735	1.260–2.388	0.001
RP11-506M12.1	1.864	1.149–3.024	0.012
RP1-39G22.7	1.123	1.022–1.235	0.016
ERVK3-1	1.185	1.021–1.374	0.025
CTD-2574D22.4	2.028	1.299–3.166	0.002
RP11-152N13.5	1.637	1.283–2.088	<0.001
RP4-769N13.6	1.836	1.139–2.959	0.013
BACE1-AS	1.207	1.091–1.335	<0.001
AC010761.8	1.466	1.170–1.839	0.001
ZNF674-AS1	1.486	1.116–1.977	0.007
ZNF529-AS1	1.301	1.038–1.631	0.022
TRAF3IP2-AS1	4.961	1.349–18.238	0.016
DDX11-AS1	4.625	2.197–9.735	<0.001
SNHG4	1.472	1.271–1.705	<0.001
RP11-398K22.12	10.874	1.844–64.126	0.008
RP11-121C2.2	2.033	1.373–3.008	<0.001
RP1-228H13.5	2.27	1.768–2.914	<0.001
AC099850.1	2.164	1.325–3.534	0.002
RP11-15E18.1	2.927	1.252–6.840	0.013
RP11-443B20.1	1.848	1.490–2.291	<0.001
SNHG12	1.142	1.044–1.249	0.004
DYNLL1-AS1	1.257	1.134–1.393	<0.001
RP11-199F11.2	1.476	1.191–1.830	<0.001
RP11-378J18.8	3.329	1.458–7.603	0.004
LINC00205	1.234	1.096–1.389	0.001
RP11-216B9.6	5.889	2.084–16.636	0.001
PTOV1-AS1	1.24	1.040–1.477	0.016
NIFK-AS1	1.274	1.059–1.532	0.010
RP5-1068E13.7	1.792	1.014–3.166	0.045
SBF2-AS1	1.426	1.113–1.828	0.005
CTD-2116N20.1	2.317	1.507–3.563	<0.001
RP5-967N21.11	1.086	1.014–1.163	0.019
RP11-498C9.15	1.755	1.344–2.291	<0.001
RP11-148K1.12	1.54	1.110–2.138	0.01
RP5-821D11.7	1.472	1.160–1.868	0.001
RP11-501C14.5	1.414	1.037–1.927	0.029
RP11-14N7.2	1.533	1.173–2.004	0.002
RP11-923I11.6	1.134	1.064–1.208	<0.001
AP001469.9	2.058	1.355–3.125	0.001
RP5-1112D6.8	2.302	1.129–4.691	0.022
RP11-701H24.8	2.185	1.141–4.183	0.018
CTD-2008L17.1	1.938	1.170–3.212	0.010
RP11-817I4.1	3.997	2.091–7.642	<0.001
HMGN3-AS1	1.988	1.284–3.078	0.002
LINC00094	1.294	1.059–1.582	0.012
AL133243.1	1.844	1.178–2.887	0.007

LASSO regression analysis of m6A regulatory gene-related lncRNA expression and prognosis in patients with liver cancer

LASSO regression analysis was performed on 89 genes screened by single-factor Cox regression analysis; a risk prediction model was established and the risk value was calculated. The risk model includes nine lncRNAs (CTD-2510F5.4, CTD-2012J19.3, SNHG4, RP1-228H13.5, RP11-443B20.1, DYNLL1-AS1 RP11-498C9.15, RP11-923I11.6 and RP11-817I4.1), which are nine m6A-related lncRNAs, and the coefficient for each lncRNA is shown in Fig. 2A and B. For each patient in the TCGA dataset, the risk score was calculated according to the coefficient of each lncRNA and the risk score is shown in Fig. 2C. According to the median risk score, patients with liver cancer were divided into low- and high-risk subgroups and the R language was used for survival analysis. The patient prognosis was made according to their risk scores. The results of the Kaplan-Meier survival curve analysis show risk scores related to patient prognosis. The clinical prognosis of patients with liver cancer with high risk scores was poor (the OS rate was low and the OS was short), while the clinical prognosis of patients with liver cancer with a low risk score was better (higher OS rate, longer OS time), and differences were significant (Fig. 2D). The ROC curve showed that the risk model had a good predictive value for prognosis (1-year OS AUC=0.738; 2-year OS AUC=0.71; Fig. 2E).

Figure 2. Least absolute1 shrinkage and selection operator (LASSO) regression was performed, calculating the minimum criteria (A-C) coefficient risk score corresponding to each lncRNA was calculated using LASSO regression analysis. (D) Kaplan-Meier curve indicating that the overall survival rate of the high-risk subgroup was lower than that of the low-risk subgroup. (E) Receiver operating characteristic curve analysis was used to predict the one- and two-year survival rate. LASSO, Least Absolute Shrinkage and Selection Operator; lncRNA, long non-coding RNA.

Prognostic analysis of nine m6A-related lncRNAs

A total of nine lncRNAs were included in the risk assessment model and single-variable COX regression analysis was used to assess the prognosis of the nine lncRNAs. The forest plot indicated that the nine lncRNAs (CTD-2510F5.4, CTD-2012J19.3, SNHG4, RP1-228H13.5, RP11-443B20.1, DYNLL1-AS1, RP11-498C9.15, RP11-923I11.6 and RP11-817I4.1) are risk factors for patients with liver cancer [hazard ratio (HR)>1; P<0.001; Fig. 3A]. The heat map showed that the expression of CTD-2510F5.4, CTD-2012J19.3, SNHG4, RP1-228H13.5, RP11-443B20.1, DYNLL1-AS1, RP11-498C9.15, RP11-923I11.6 and RP11-817I4.1 increased as the risk score increased. The expression of the nine lncRNAs was also associated with clinicopathological features of liver cancer, such as stage, gender, age and fustat (Fig. 3B). The Kaplan-Meier survival curves showed that in the TCGA dataset, low expression of CTD-2510F5.4, CTD-2012J19.3, SNHG4, RP1-228H13.5, RP11-443B20.1, DYNLL1-AS1, RP11-498C9.15, RP11-923I11.6 and RP11-817I4.1 was associated with favorable prognosis (Fig. 3C-K).

Figure 3. Prognostic analysis of nine lncRNAs associated with the m6A gene. (A) Forest map results of the prognostic ability of nine m6A gene-related lncRNAs. (B) Heatmap of the association between the expression level of nine m6A gene-related lncRNAs and clinicopathological features. (C-H) Kaplan-Meier curves showing survival according to (C) RP11-923I11.6 m6A-related lncRNA expression levels; (D) RP11-817I4.1 m6A-related lncRNA expression levels; (E) RP11-498C9.15 m6A-related lncRNA expression levels; (F) CTD-2510F5.4 m6A-related lncRNA expression levels; (G) DYNLL1-AS1 m6A-related lncRNA expression levels; (H) CTD-2012J19.3 m6A-related lncRNA expression levels; (I) RP11-443B20.1 m6A-related lncRNA expression levels; (J) RP1-228H13.5 m6A-related lncRNA expression levels; and (K) SNHG4 m6A-related lncRNA expression levels. lncRNAs, long non-coding RNA; m6A, N6-methyladenosine.

Association between clinicopathological features and risk score

To determine whether clinicopathological features are associated with risk scores, association analysis was performed. The results showed that the risk score was correlated with TNM stage, gender and age (P<0.05; Fig. 4A-F). In order to better evaluate the prognostic ability of the risk model, hierarchical analysis was conducted to determine whether the risk model can predict the prognosis of different subgroups. The results showed that compared with low-risk patients, high-risk patients with the following clinicopathological characteristics had a poor prognosis: T1, T3, age >65 or age ≤65 years, stage 3 and male gender (Fig. 4G-M).

Figure 4. Association between clinicopathological features and risk score. (A-F) Patients with different clinicopathological features have different levels of risk scores: (A) TNM stage, (B) stage, (C) gender, (D) M stage, (E) N stage and (F) age. (G-M) Prognostic analysis of risk model in multiple subgroups of patients with liver cancer: (G) T1, (H) T3, (I) Stage III, (J) age >65 years, (K) age ≤65 years, (L) males and (M) females.

Enrichment analysis of DEGs in the high- and low-risk groups

The aim of the present study was to determine the potential biological processes and pathways of molecular heterogeneity between the low- and high-risk subgroups. Using P<0.05 and log2F old Change >1 as the screening criteria, 2,352 DEGs between the high- and low-risk groups were identified. KEGG analysis of these DEGs showed that they were mainly enriched in ‘neuroactive ligand-receptor interaction’, ‘calcium signaling pathway’, ‘bile secretion’, ‘metabolism of xenobiotics by cytochrome P450’ and ‘retinol metabolism’ (Fig. 5A). GO enrichment analysis showed that they were mainly enriched in biological process terms of ‘small molecule process’, ‘embryonic organ morphogenesis’ and ‘hormone metabolic process’, cellular component terms of ‘collagen-containing extracellular matrix’, ‘apical part of cell’ and ‘synaptic membrane’ and molecular function terms of ‘channel activity’, ‘passive transmembrane transporter activity’ and ‘ion channel activity’ (Fig. 5B); These results may provide some clues for us to understand the cell biology-related effects of different genes in the high- and low-risk groups.

Figure 5. Enrichment analysis of differentially expressed genes in low- and high-risk subgroups, and univariate and multivariate Cox regression analysis of risk score. (A) Gene Ontology functional analysis of 2,352 differentially expressed genes. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis of 2,352 differentially expressed genes. Both (C) Univariate and (D) multivariate regression analysis showed that the risk score of patients was a risk factor.

The risk score model is an independent prognostic factor for patients with liver cancer

To assess whether the risk score is an independent prognostic factor in patients with liver cancer, the data of liver cancer patients from the TCGA dataset were analyzed. The results of the univariate Cox regression analysis indicated that the risk score model was significantly associated with the OS hazard ratio (HR: 1.318; P<0.001), multivariate Cox regression analysis further showed that the risk scoring model was an independent predictor of OS (HR: 1.280, P<0.001; Fig. 5C and D). These results further suggested that the risk scoring model can be used as an independent prognostic factor, which may be helpful for the evaluation of clinical prognosis.

Construction of ceRNA network and gene function enrichment analysis

In order to further clarify the m6A-related lncRNAs, sponge miRNAs and regulated mRNA in hepatocellular carcinoma, a ceRNA network was constructed based on m6A-related lncRNAs. The miRNAs targeted by nine lncRNAs were predicted from the miRcode database; 3/9 lncRNAs were extracted (DYNLL1-AS1, RP1-228H13.5 and RP11-923I11.6), and 91 pairs of interactions between 3 lncRNAs and 74 miRNAs were identified. After intersecting the predicted miRNA with the differential miRNAs screened by sequencing, 6 intersecting miRNAs were obtained (hsa-miR-205, hsa-miR-216b, hsa-miR-217, hsa-miR-137, hsa-miR-182 and hsa-miR-183). Eight pairs of interactions between 3 lncRNAs and 6 miRNAs were identified. The target genes predicted by miRNA and the differential mRNAs screened by sequencing were crossed and 201 genes were obtained. Finally, a ceRNA network containing 3 lncRNAs, 6 miRNAs and 201 mRNAs was obtained (Fig. 6A). Due to the large number of genes regulated by RP1-228H13.5, which may be more important in liver cancer, the RP1-228H13.5/hsa-miR-205/ZIK1 signaling axis was selected for further research. GO enrichment analysis indicated that the 201 target genes were mainly concentrated in ‘response to forskolin’, ‘membrane transporter activity’, ‘membrane organization’ and ‘cell junction’; regarding KEGG pathways, target genes were mainly concentrated in the ‘cGMP/PKG signaling pathway’, ‘cellular sensitivity’, ‘colorectal cancer’ and ‘apoptosis’ (Fig. 6B and C).

Figure 6. Construction of ceRNA network map and enrichment analysis of target genes. (A) LncRNA-miRNA-mRNA ceRNA network map. (B) GO enrichment analysis of target genes. (C) KEGG enrichment analysis of target genes. ceRNA, competing endogenous RNA; lncRNAs, long non-coding RNA; miRNA/miR, microRNA; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Expression of m6A methylation-related RP1-228H13.5 in liver cancer cells and its impact on cancer cell proliferation, migration, invasion and apoptosis. RT-qPCR was used to detect the expression levels of RP1-228H13.5, miR-205-5p and ZIK1 in normal and liver cancer cells. The results showed that RP1-228H13.5 and ZIK1 were significantly upregulated in liver cancer cell lines compared with THLE-2 cells (P<0.05), while miR-205-5p was significantly downregulated in liver cancer cell lines (P<0.05; Fig. 7A-C). After overexpressing and knocking down RP1-228H13.5 in cells, the expression levels of RP1-228H13.5, miR-205-5p and ZIK1 were detected. The results showed that after overexpressing RP1-228H13.5, RP1-228H13.5 and ZIK1 were significantly upregulated (P<0.05), while miR-205-5p was significantly downregulated (P<0.01; Fig. 7D-F). The effects of RP1-228H13.5 on the proliferation, migration, invasion and apoptosis of liver cancer cells were further investigated through cell experiments. The results of the CCK-8 assay showed that RP1-228H13.5 overexpression promoted the proliferation of Hep3B2.1–7 cells, while knockdown of RP1-228H13.5 inhibited the proliferation of Hep3B2.1–7 cells (Fig. 7G). The results of the wound-healing assay showed that RP1-228H13.5 overexpression promoted the migration of Hep3B2.1–7 cells, while RP1-228H13.5 knockdown inhibited it (Fig. 7H and I). Transwell assay results showed that RP1-228H13.5 overexpression promoted the migration and invasion of Hep3B2.1–7 cells, while RP1-228H13.5 knockdown inhibited them (Fig. 7J and K). The results of the cell apoptosis assay showed that RP1-228H13.5 overexpression inhibited Hep3B2.1–7 cell apoptosis, while RP1-228H13.5 knockdown promoted Hep3B2.1–7 cell apoptosis (Fig. 7L and M). The above results indicated that RP1-228H13.5 promotes the proliferation, migration and invasion of Hep3B2.1–7 cells and inhibits their apoptosis.

Figure 7. RP1-228H13.5 is upregulated in liver cancer and promotes the proliferation, migration and invasion of Hep3B2.1–7 cells, while inhibiting cell apoptosis. (A) RT-qPCR assay was used to detect the expression levels of RP1-228H13.5 in THLE-2, SNU-398, Huh-7 and Hep3B2.1–7 cells. (B) RT-qPCR assay was used to detect the expression levels of ZIK1 in THLE-2, SNU-398, Huh-7 and Hep3B2.1–7 cells. (C) RT-qPCR assay was used to detect the expression levels of miR-205-5p in THLE-2, SNU-398, Huh-7 and Hep3B2.1–7 cells. (D-F) Following the overexpression and knockdown of RP1-228H13.5, RT-qPCR was performed to detect the expression levels of (D) RP1-228H13.5, (E) ZIK1 and (F) miR-205-5p. (G) Following the overexpression and knockdown of RP1-228H13.5, CCK-8 assay was used to detect the cell proliferation ability. (H and I) Following the overexpression and knockdown of RP1-228H13.5, a wound-healing assay was used to detect cell migration. (H) Representative images (scale bar, 200 µm) and (I) quantitative results. (J and K) Following overexpression and knockdown of RP1-228H13.5, a Transwell assay was used to detect cell invasion ability. (J) Representative images (scale bar, 200 µm) and (K) quantitative results. (L and M) Following the overexpression and knockdown of RP1-228H13.5, flow cytometry was used to detect cell apoptosis. (L) Representative graphs and (M) quantitative results. Values are expressed as the mean ± standard error of the mean. An unpaired Student's t-test was used for two-group comparison and one-way ANOVA for multi-group comparisons. *P<0.05, **P<0.01. RT-qPCR, reverse transcription-quantitative PCR; CCK-8, cell-counting kit; NC, negative control; OV, overexpression; Sh, short hairpin RNA; miR, microRNA; ZIK1, zinc finger protein interacting with K protein 1.

RP1-228H13.5 can target miR-205-5p to regulate ZIK1 expression

FISH is an accurate method for determining the localization of RNA in cells. In this experiment, liver cancer cell lines were collected and cultured for FISH to determine the cell localization of RP1-228H13.5. The results showed that RP1-228H13.5 was localized in the nucleus of liver cancer cell lines (Fig. 8A). The binding sites of RP1-228H13.5 and miR-205-5p were predicted using the Starbase 2.0 database, as presented in Fig. 8B. The results of the dual luciferase assay showed that following miR-205-5p overexpression, the luciferase activity of RP1-228H13.5-wt was significantly reduced (P<0.05), while the relative value of the luciferase activity of RP1-228H13.5-mut showed no significant change (P>0.05; Fig. 8C), indicating that RP1-228H13.5 is able to regulate miR-205-5p. The database predicted the binding site between miR-205-5p and ZIK1, as presented in Fig. 8D. The results of the dual luciferase assay showed that following miR-205-5p overexpression, the activity of ZIK1-wt luciferase significantly decreased (P<0.05), while the relative value of ZIK1-mut luciferase activity did not change significantly (P>0.05; Fig. 8E), indicating that miR-205-5p can regulate ZIK1. An RNA pull-down assay was also conducted to confirm the binding of miR-205-5p to RP1-228H13.5 and ZIK1. The results showed that compared with the control probe group, the fold enrichment of RP1-228H13.5 and ZIK1 in the miR-205-5p probe group was significantly increased, indicating that miR-205-5p can bind to RP1-228H13.5 and ZIK1 (Fig. 8F and G). The RT-qPCR assay results showed that overexpression of miR-205-5p significantly increased the expression level of miR-205-5p; after knocking down miR-205-5p, the expression level of miR-205-5p significantly decreased, indicating successful transfection of miR-205-5p in Hep3B2.1–7 cells (Fig. 8H). To further clarify the regulatory relationship between miR-205-5p and RP1-228H13.5 and ZIK1, OV-RP1-228H13.5 alone or with miR-205-5p mimics were transfected into cells to detect the expression level of ZIK1. The results showed that, compared with the NC control group, the expression of ZIK1 was significantly promoted in the OV-RP1-228H13.5 group, while compared with the OV-RP1-228H13.5 group, co-transfection of OV-RP1-228H13.5 and miR-205-5p mimics significantly reduced the expression level of ZIK1 (Fig. 8I-K). Following the transfection of cells with Sh-RP1-228H13.5 alone or with miR-205-5p inhibitor, the expression level of ZIK1 was detected. The results showed that, compared with the NC control group, the expression of ZIK1 was significantly inhibited in the Sh-RP1-228H13.5 group, while, compared with the Sh-RP1-228H13.5 group, co-transfection of Sh-RP1-228H13.5 and miR-205-5p inhibitor significantly increased the expression level of ZIK1 (Fig. 8L-N). The above results suggested that miR-205-5p can reverse the effect of Sh-RP1-228H13.5 on the expression of ZIK1. To clarify the effect of miR-205-5p on the expression of ZIK1, ZIK1 expression was detected following the overexpression and knockdown of miR-205-5p. The results showed that, as compared with the NC group, the expression of ZIK1 was significantly inhibited in the miR-205-5p mimics group, while it was significantly promoted in the miR-205-5p inhibitor group (Fig. 8O-Q). The above study suggested that Sh-RP1-228H13.5 can target miR-205-5p to regulate the expression level of ZIK1.

Figure 8. RP1-228H13.5-targeted miR-205-5p regulates the expression of ZIK1. (A) Fluorescence in situ hybridization detected the localization of RP1-228H13.5 in cells (scale bar, 200 µm). (B) Prediction of the binding sites of RP1-228H13.5 and miR-205-5p through the Starbase 2.0 database. (C) The combination of RP1-228H13.5 and miR-205-5p was verified in a dual luciferase experiment. (D) Prediction of the binding sites of miR-205-5p and ZIK1 through the Starbase 2.0 database. (E) The combination of miR-205-5p and ZIK1 was verified by a dual luciferase experiment. (F) RNA pull-down assay verified the binding of miR-205-5p and RP1-228H13.5. (G) RNA pull-down assay verified the binding of miR-205-5p and ZIK1. (H) RT-qPCR assay detected the expression of miR-205-5p after overexpression and knockdown of miR-205-5p in Hep3B2.1–7 cells. (I) Following the transfection of cells with OV-RP1-228H13.5 alone or with miR-205-5p mimics, RT-qPCR was used to detect the expression of ZIK1. (J and K) Following the transfection of cells with OV-RP1-228H13.5 alone or with miR-205-5p mimics, western blotting was used to detect the expression of ZIK1. (J) Representative western blots and (K) quantified results. (L) Following the transfection of Sh-RP1-228H13.5 alone or with miR-205-5p inhibitor into cells, RT-qPCR were used to detect the expression of ZIK1. (M and N) Following the transfection of Sh-RP1-228H13.5 alone or with miR-205-5p inhibitor into cells, western blotting was used to detect the expression of ZIK1. (M) Representative western blots and (N) quantified results. (O) Following the overexpression and knockdown of miR-205-5p, the expression level of ZIK1 was detected by RT-qPCR. (P and Q) Following the overexpression and knockdown of miR-205-5p, the expression level of ZIK1 was detected by western blotting. (P) Representative western blots and (Q) quantified results. Values are expressed as the mean ± standard error of the mean. An unpaired Student's t-test was used for two-group comparison. *P<0.05, **P<0.01. ns, no significance; RT-qPCR, reverse transcription-quantitative PCR; NC, negative control; OV, overexpression; Sh, short hairpin RNA; miR, microRNA; ZIK1, zinc finger protein interacting with K protein 1; mut, mutant; wt, wild-type.

m6A methyltransferase METTL14 mediates the regulation of miR-205-5p expression by RP1-228H13.5

RT-qPCR detection of METTL14 expression levels in normal and liver cancer cells showed that METTL14 was significantly downregulated in liver cancer cell lines compared with THLE-2 cells (P<0.05; Fig. 9A). The RT-qPCR assay results showed that overexpression of METTL14 significantly increased the expression of METTL14; after knocking down METTL14, the expression level of METTL14 significantly decreased, indicating successful transfection of METTL14 and its inhibitor in Hep3B2.1–7 cells (Fig. 9B). To investigate the regulatory effect of m6A methyltransferase METTL14 on RP1-228H13.5, METTL14 was overexpressed or knocked down. The results showed that METTL14 overexpression inhibited the expression of RP1-228H13.5, while knocking down METTL14 promoted the expression of RP1-228H13.5 (Fig. 9C). A further RNA pull-down assay confirmed that METTL14 can bind to RP1-228H13.5 (Fig. 9D). The cells were co-transfected with shMETTL14 or shMETTL14 + sh-RP1-228H13.5. RT-qPCR assay was performed to detect the expression of miR-205-5p. The results showed that knockdown of METTL14 impaired the promoting effect of shRP1-228H13.5 on miR-205-5p in cells (Fig. 9E). After co-transfecting cells with OV-METTL14 or OV-METTL14 + OV-RP1-228H13.5, the expression of miR-205-5p was also detected by RT-qPCR assay. The results showed that METTL14 upregulation impaired the inhibitory effect of OV-RP1-228H13.5 on miR-205-5p in cells (Fig. 9F). The above results suggested that m6A methyltransferase METTL14 mediates the regulation of miR-205-5p expression by RP1-228H13.5.

Figure 9. N6-methyladenosine methyltransferase METTL14 mediates the regulation of miR-205-5p expression by RP1-228H13.5. (A) RT-qPCR was used to detect the expression level of METTL14 in normal and liver cancer cells. (B and C) METTL14 overexpression and knockdown were performed, and the expression level of (B) METTL14 and (C) RP1-228H13.5 was detected by RT-qPCR in Hep3B2.1–7 cells. (D) RNA pull-down assay verified the binding of METTL14 with RP1-228H13.5. (E) After co-transfecting cells with sh-METTL14 or sh-METTL14 + sh-RP1-228H13.5, the expression of miR-205-5p was detected using RT-qPCR. (F) After co-transfecting cells with OV-METTL14 or OV-METTL14 + OV-RP1-228H13.5, the expression of miR-205-5p was detected by RT-qPCR. Values are expressed as the mean ± standard error of the mean. An unpaired Student's t-test was used for two-group comparison and one-way ANOVA for multi-group comparisons. *P<0.05, **P<0.01. RT-qPCR, reverse transcription-quantitative PCR; NC, negative control; OV, overexpression; Sh, short hairpin RNA; miR, microRNA; Ab, antibody; METTL14, methyltransferase 14, N6-adenosine-methyltransferase subunit.

Discussion

m6A methylation modification of RNA is a common modification after transcription in eukaryotes (12). It is reported in the literature that ~20% of human mRNA can be modified by conventional methylation. m6A is often located near the stop codon of mRNA. It has been determined that >7,000 different mRNA molecules have m6A methylation modification (13), which means that m6A modification may widely affect gene expression. m6A methylation-modified RNA can change the structure of RNA by weakening base pairing, increasing protein binding recognition sites and recruiting protein complexes to participate in cell biology processes, including mRNA splicing, RNA output, stability and immune tolerance (14,15). It was found that m6A methylation can also modify lncRNA and affect its stability and metabolism (16). m6A regulators can maintain the malignancy of a variety of tumors by modifying specific lncRNAs. METTL3 promotes liver cancer progression through YTHDF2-dependent post-transcriptional silencing of SOCS2 (17); METTL16 promotes the progression of hepatocellular carcinoma by downregulating RAB11B-AS1 in an m6A-dependent manner (18); the ALKBH5/MAP3K8 axis can promote the progression of hepatocellular carcinoma by regulating PD-L1 and macrophage infiltration (19); HIF-1α induces YTHDF1 to promote the translation of ATG2A and ATG14 to drive hypoxia-induced autophagy and malignancy of hepatocellular carcinoma (20); YTHDF2 regulates OCT4 expression through m6A RNA methylation, promoting liver cancer cell phenotype and cancer metastasis (21); YTHDF3 mediates m6A modification and PFKL as a functional loop in glycolysis in hepatocellular carcinoma (22); lncRNA PCAT6 regulates the miR-326/HNRNPA2B1 axis in liver cancer and promotes proliferation and invasion of liver cancer cells (23); IGF2BP1 regulates the stability of MGAT5 mRNA and promotes the phenotype of liver cancer stem cells (24); IGF2BP2 promotes the growth of liver cancer cells through a m6A-FEN1-dependent mechanism (25); ALKBH5-mediated m6A demethylation of the lncRNA RMRP has a carcinogenic role in lung adenocarcinoma (26). ALKBH5 can upregulate the expression of lncRNA NEAT1 through demethylation and inhibit the malignant behavior of colon cancer (27). LCAT3 is a novel m6A-regulated lncRNA that activates c-MYC by binding to FUBP1 and has a carcinogenic role in lung cancer (28). It has been shown that m6A-modified lncRNAs can affect the occurrence and development of tumors, and lncRNAs may act as ceRNAs to target m6A regulators, thus affecting the invasive progression of tumors. In summary, it is indicated that m6A is modified to target lncRNAs, and accordingly, increased attention should be paid to the interactions and functions of lncRNAs and m6A modifications to determine potential tumor prognostic markers or treatment targets.

In the present study, 374 patients with liver cancer from the TCGA dataset were included to explore the prognosis of m6A-related lncRNAs. The prognostic value of 89 m6A-related lncRNAs was confirmed and nine of them were used to establish a risk model to predict the OS of patients with liver cancer, including CTD-2510F5.4, CTD-2012J19.3, SNHG4, RP1-228H13.5, RP11-443B20.1, DYNLL1-AS1 RP11-498C9.15, RP11-923i11.6 and RP11-817I4.1. The established ceRNA network identified the RP1-228H13.5/hsa-miR-205/ZIK1 signaling axis for subsequent research, of which RP1-228H13.5 has not been reported in previous studies. In the present study, it was found that RP1-228H13.5 is significantly upregulated in liver cancer and can promote the proliferation, migration and invasion of liver cancer cells, while inhibiting cell apoptosis. The above results suggest that RP1-228H13.5 is a tumor-promoting gene in liver cancer. The Starbase database confirmed that RP1-228H13.5 can be combined with miR-205-5p. Research has found that miR-205-5p has an important role in the invasion and migration of liver cancer and that is has a promoting role in various malignant tumors, such as gastric cancer (29), colon cancer (30) and non-small cell lung cancer (31), participating in tumor invasion and migration. Previous studies have found that miR-205-5p is a key regulatory factor for VEGFA in liver cancer-related angiogenesis (32) and YB1 regulates the miR-205/200b-ZEB1 axis by inhibiting the maturation of miRNAs in liver cancer (33). DNAJA1 stabilizes EF1A1 and promotes miR-205-5p-mediated cancer cell proliferation and metastasis (34). miR-205-5p regulates the chemotherapeutic resistance of hepatocellular carcinoma cells by targeting the PTEN/JNK/ANXA3 pathway (35). The database confirmed that ZIK1 can be combined with hsa-miR-205. ZIK1 was initially considered a transcription inhibitor that binds to nuclear RNA protein particle K protein (36,37) and exhibits high methylation in gastrointestinal metaplasia (38). Abnormal expression of ZIK1 is associated with the occurrence of renal cell carcinoma (39), diffuse large B-cell lymphoma (39), esophageal cancer (40), colorectal cancer and gastric cancer (41). Previous studies have also found that miR-197-3p may partially regulate hepatocellular carcinoma survival by downregulating ZIK1. Therefore, the miR-197-3p/ZIK1 axis may serve as a new therapeutic target for patients with liver cancer (42). The present study found that m6A methylation-related RP1-228H13.5 can target and bind with hsa-miR-205 to regulate the expression of ZIK1.

Our research team is interested in the METTL14 gene and has previously found that METTL14 may regulate the expression of RP1-228H13.5. Therefore, METTL14 was selected for subsequent analysis. Other m6A-related genes will be our future research endeavor. Research has also found that METTL14 can regulate the expression of RP1-228H13.5. METTL14 is an important m6A methylase that has a key role in the m6A methylation process of RNA, and has also been found to have a role in the occurrence or development of various tumors (43). It was reported that METTL14 inhibits the potential metastatic ability of hepatocellular carcinoma by regulating m6A methylation-dependent miRNA processes (44). Lin et al (45) found a decrease in METTL14 expression in lung cancer tissue, but to date, its function in lung cancer has not been well studied. The present study found that METTL14 can mediate the regulation of miR-205-5p expression by RP1-228H13.5.

In the present study, bioinformatics and statistical analyses were conducted on relevant data of patients with cancer from the TCGA database to determine the prognostic significance of m6A-related lncRNAs. An RP1-228H13.5/hsa-miR-205/ZIK1 signaling axis was identified. Finally, through cell experiments, it was determined that m6A methyltransferase METTL14 may mediate RP1-228H13.5 to promote the proliferation, migration, invasion and inhibition of liver cancer cell apoptosis by targeting hsa-miR-205/ZIK1, providing ideas for the study of targets and mechanisms of liver cancer occurrence and development. However, the present study only analyzed METTL14 and the data analyzed were from 2021, which is a limitation of this study. Further research will be conducted on other m6A-related genes in the future.

Acknowledgements

Not applicable.

Funding

Funding: No funding was received.

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

CL, KQ and FM conceived the study, planned and carried out experiments and were involved in data analysis and manuscript writing-reviewing and editing. JX, CL, KQ and YW conceived the study, planned and carried out experiments and were involved in data analysis and original draft preparation. JYZ and SL carried out experiments and performed data analysis. KQ and FM have checked and confirmed the authenticity of the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References