In the present study, two mRNA gene expression datasets (GSE69164 and GSE77509) and one miRNA microarray dataset (GSE76903) were obtained from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo), an important database that contains records from chips, as well as second-generation sequencing and other high-throughput sequencing data. The raw microarray count data from the GEO datasets were normalized by a size factor, and probes were translated to the corresponding gene symbol based on platform annotation information. GSE69164 contains 11 HCC tissue samples and 11 adjacent normal liver tissue samples. GSE77509 and GSE76903 contain 20 HCC tissue samples and 20 adjacent normal liver tissue samples. RNA sequencing raw count data, containing lncRNA sequencing data (from 375 HCC tissues and 49 adjacent normal liver tissues), mRNA sequencing matrix data (from 375 HCC tissues and 49 adjacent normal liver tissues), and miRNA sequencing expression results (from 372 HCC tissues and 50 adjacent normal liver tissues), together with clinical information for 309 HCC patients, were downloaded from TCGA-LIHC (https://portal.gdc.cancer.gov/). Our research conforms fully with TCGA and GEO publication requirements; no ethics committee approval was needed.

The “DESeq2” and “limma” R packages for difference analysis were used in our study. “DESeq2” package applies to RNA sequencing raw count data and only accept integers, and “limma” R package is mainly used for difference analysis of chip data with small sample sizes. In addition, it is generally considered that RNA sequencing raw count data do not conform to a normal distribution, while chip data do; thus, difference analysis using the “limma” package results in large errors for RNA sequencing count data. Therefore, the differentially expressed lncRNAs (DELs) and differentially expressed mRNAs (DEMs) between HCC tissues and adjacent normal liver tissues in TCGA-LIHC were identified using the “DESeq2” package. The screening criteria included two conditions: a false discovery rate (FDR)<0.01 and |log2 fold change (log2 FC)|≥2. The “limma” package was used to identify DEMs and differentially expressed miRNAs (DEMIs) between normal and HCC tissues in the GSE69164, GSE77509, and GSE76903 datasets. The screening criteria for DEMs were |log2 FC|≥2 and adjacent P-value (adj. P)<0.01. For DEMIs, |log2 FC|≥1.5 and adj. P<0.05 were considered to indicate statistical significance.

Overlapping mRNAs were selected as DEMs from the GSE69164, GSE77509, and TCGA-LIHC datasets. miRTarBase (http://mirtarbase.cuhk.edu.cn/php/index.php) (14), a database of experimentally validated miRNA–target interactions, was used to predict miRNA–mRNA interactions based on DEMIs from GSE76903. The target mRNAs were intersected with DEMs. Then, we used the starBase database (http://starbase.sysu.edu.cn/) (15) to identify the interactions between miRNAs and lncRNAs. The predicted target lncRNAs were intersected with DELs from the TCGA-LIHC dataset. Based on the interactions of these DEGs, we constructed a preliminary and unverified lncRNA–miRNA–mRNA regulatory network associated with HCC, which was visualized using the Cytoscape software (version 3.6.1).

PPI networks, which reflect the functional interactions between proteins, are used for the identification of key genes in disease generation and development. In this study, a PPI network was constructed using the STRING (Search Tool for the Retrieval of Interacting Genes; http://string-db.org) (version 11.0) online database. An interaction with a combined score ≥0.4 was considered statistically significant. The PPI network was visualized using Cytoscape version 3.6.1. Nodes of degree ≥35 in the PPI network were selected as hub genes. In addition, a ceRNA sub-network based on hub genes was screened from the whole ceRNA network, and corresponding interactions were visualized with Cytoscape.

To identify relapse-related genes, we first compared expression levels of DEGs between HCC patients with and without recurrence within 3 years. Boxplots were drawn using the “beeswarm” R package. Subsequently, correlations between the expression of DEGs and clinical features (tumor grade and stage) were analyzed based on the clinical information of 309 HCC patients from TCGA. Disease-free survival (DFS) curves were plotted based on the log-rank test using GEPIA (Gene Expression Profiling Interactive Analysis; http://gepia.cancer-pku.cn/) (16), a website for analyzing RNA sequencing data in TPM format from TCGA and the GTEx project.

Correlation analysis was performed to identify a relapse-related lncRNA–miRNA–mRNA network using RNA sequencing data in FPKM format from TCGA-LIHC. The Pearson method was used to determine correlation coefficients. P<0.05 was considered statistically significant. In addition, the expression levels of DEMIs between HCC tissue and normal liver tissue in the network were verified in the TCGA-LIHC dataset. The final relapse-related lncRNA–miRNA–mRNA network was visualized using Cytoscape (17).

Two HCC cell lines, HepG2 and HuH-7, were obtained from the American Type Culture Collection (ATCC). All cells were cultivated in Dulbecco’s Modified Eagle’s Medium (DMEM; Gibco, USA) containing 10% fetal bovine serum (FBS; Gibco, USA) in an incubator containing 5% CO₂ at a constant temperature of 37°C.

HCC cells were inoculated in six-well plates. Seed cells to be 70–90% confluent at transfection according to the instructions of Lipofectamine 3000 (Invitrogen, USA). The transfection plasmids used in this experiment included short interfering RNA (siRNA)-ASF1B, a siRNA negative control (NC), miR-214-3p mimics, NC mimics and SNHG3 overexpression plasmid (SNHG3 pcDNA). Cells were incubated for 48 h after transfection. Transfection efficiency was examined through the subsequent qRT-PCR.

We first used the TRIzol method (Invitrogen, Carlsbad, CA, USA) to extract total RNA from cells. Extracted RNA was reverse transcribed into complementary DNA (cDNA) for qRT-PCR. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and U6 were used as endogenous controls. The qRT-PCR experiments followed the instructions of the CHAMQ SYBR qPCR Master Mix kit (Vazyme, Nanjing, China). The primer sequences used in this study are shown in Table S1 . The mRNA expression levels were analyzed by the 2^-ΔΔCt method (18).

The binding sites of miRNA-214-3p in SNHG3 and ASF1B sequences were first predicted using the starBase platform and then expanded by PCR. Subsequently, the amplified fragment was inserted into a vector to construct wild-type (WT) ASF1B (5′-agUGCCUGUCAAGGCUCCAGUCCUGCUGa-3′) and WT miRNA-214-3p (3′-ugacggacaGACACGGACGACa-5′) plasmids. Then, we altered partial nucleotides by a gene mutation technique and constructed mutant-type (MT) ASF1B (5′-agACGGACAGAAGCCACCACUGGACGACa-3′) and MT miRNA-214-3p (3′-ugacggacaCUCAGCCUGCUGa-5′) plasmids. Cells were co-transfected with the ASF1B-WT or ASF1B-MT plasmid and miR-214-3p mimics, and with the miRNA-214-WT or miRNA-214-MT plasmid and SNHG3 overexpression plasmid. Finally, the luciferase activities in each group were detected.

TIMER is a comprehensive database for analyzing the levels of immune infiltrates in different tumors (https://cistrome.shinyapps.io/timer/) (19). Correlations between relapse-related genes and infiltration levels of different immune cells were estimated using TIMER. Then, the correlations of ASF1B with gene markers of B cells, CD8+ T cells, neutrophils, macrophages, dendritic cells, monocytes, natural killer cells, and regulatory T (Treg) cells in HCC were verified, and tumor purity and patient age were used to correct the P-values. Subsequently, to further verify the relationship between ASF1B and the tumor immune microenvironment, we tested the expressions of partial immune gene markers (COR>0.3 after purity correction) in ASF1B-knockout HCC cells by RT-qPCR.

The HCC cells were divided into 3 groups. The first group was that a blank vector plasmid control. The second group was that SNHG3 was overexpressed in HepG2 cell. The final group was that HepG2 cells were co-transfected with overexpressed plasmid targeted SNHG3 and knockdown plasmid targeted ASF1B. All Cells were incubated for 72 h after transfection. Then, the cells were harvested and incubated with anti-human CD279 (PD-1) Monoclonal Antibody (MIH4) linked with phycoerythrin (PE) fluorochrome (12-9969-42; eBioscience: Thermo Fisher Scientific, Inc.) for 30 min at 4°C in dark. Cells were resuspended in 0.5 ml phosphate buffer saline (PBS) for flow cytometry.

GraphPad Prism (version 8.0) and R statistical software (version 3.6.1) were used for statistical analysis. Experimental data were expressed as mean ± standard deviation. The standard t-test was used to compare the differences between the two groups. Correlation analysis was performed using the Spearman rank correlation test. P<0.05 was considered statistically significant.

A total of 804, 1384, and 1389 mRNAs were screened from the GSE69164, GSE77509, and TCGA-LIHC datasets, respectively, according to the criteria |log2 FC|≥2 and FDR ≤ 0.01 or adj.P<0.01. The 290 overlapping mRNAs among the three datasets were considered to be DEMs. Similarly, 707 upregulated and 16 downregulated DELs were identified from TCGA-LIHC based on the above procedure, as well as 44 upregulated and 38 downregulated DEMIs in GSE76903 based on the criteria of |log2 FC|≥1.5 and adj.P<0.05. Volcano plots of DEMs, DEMIs, and DELs, and a Venn diagram of DEM intersection are shown in Figure 1 .

We employed 290 DEMs, 82 DEMIs, and 723 DELs to build a ceRNA network. The lncRNA–miRNA and miRNA–mRNA interactions were predicted based on DEMIs using the starBase and MiRTarBase online databases. Finally, a preliminary and unverified ceRNA network was constructed from 124 lncRNAs, 35 miRNAs, and 80 mRNAs ( Figure S1 ). There were 314 connections between lncRNAs and miRNAs, and 152 connections between miRNAs and mRNAs in this network.

To determine the interactions between ceRNA-related mRNAs and explore their functions at the protein level, a PPI network was constructed using the STRING database ( Figure 2A ). The PPI network consisted of 59 nodes and 702 edges. No interactions were found for 21 mRNAs in the PPI network, implying that these mRNAs might be unimportant in the development of HCC; thus, we discarded these 21 mRNAs. A total of 22 genes were considered as hub genes with degree ≥35, and their interactions were obtained using Cytoscape for further study ( Figure 2B ).

We first compared expression levels of hub genes between HCC patients with and without recurrence within 3 years. The results showed that the expression of ASF1B (P=0.023), AURKB (P=0.026), CCNB1 (P=0.041), CDKN3 (P=0.042), and DTL (P=0.015) was significantly higher in patients with recurrence compared with those without recurrence ( Figures 3A–E ). The relationships between these five mRNAs and HCC grade and stage were further analyzed. As shown in Figures 3F–O , higher expression levels of these mRNAs were significantly associated with higher tumor grade and more advanced tumor stage. A subsequent DFS analysis using GEPIA showed that HCC patients with high expression of ASF1B, AURKB, CCNB1, CDKN3, and DTL had poor DFS ( Figure 4 ). These findings further indicated that ASF1B, AURKB, CCNB1, CDKN3, and DTL could be prognostic markers of HCC relapse. In addition, an lncRNA–miRNA–mRNA sub-network based on five relapse-related DEMs was constructed using Cytoscape ( Figure 5 ).

According to the lncRNA–miRNA–mRNA sub-network, miR-183-5p, miR-1307-3p, and their corresponding lncRNAs were all upregulated in HCC; however, according to ceRNA theory, there should be a negative correlation between lncRNAs and miRNAs. Therefore, we removed these two miRNAs and their corresponding lncRNAs from the network, and the remaining lncRNAs were used to find relapse-related DELs.

The results showed that nine lncRNAs were highly expressed in HCC patients with relapse compared with those without relapse ( Figure S2 and Figures 6A, D ), and the expression levels of SNHG3 and LINC00205 were closely correlated with HCC grade and stage ( Figures 6B, C, E, F ). Survival analysis further indicated that patients with high expression of SNHG3 and LINC00205 had shorter DFS than those with low expression ( Figures 6G, H ).

As miR-214-3p is the common target of SNHG3 and LINC00205, its expression levels in HCC tissues and normal liver tissues were also verified using TCGA-LIHC data. The results confirmed that it was significantly downregulated in HCC tissues ( Figure 7A ). Correlation analysis was performed to acquire results that were more reliable. The results showed that there was no correlation between LINC00205 and miR-214-3p (R=-0.039, P=0.43; Figure 7B ). MiR-214-3p-ASF1B (R=-0.14, P=0.003) and SNHG3-miR-214-3p (R=-0.11, P=0.02) had negative correlations ( Figures 7C, D ), whereas ASF1B was significantly positively correlated with SNHG3 (R=0.4; P<2.2e-16; Figure 7E ). Finally, a potential lncRNA–miRNA–mRNA regulatory axis associated with HCC recurrence was identified. The interactions among its members are shown in Figure 7F .

First, the binding sites of miRNA-214-3p in the SNHG3 and ASF1B sequences were predicted using the starBase platform. To verify the regulatory effect of miR-214-3p on ASF1B, we constructed WT and MT sequences of ASF1B. ASF1B-WT and ASF1B-MT were co-transfected with miR-214-3p mimics into HepG2 cells. As expected, the dual-luciferase reporter gene assay results indicated that overexpression of miR-214-3p significantly reduced the luciferase activities of cells transfected with vectors containing ASF1B-WT but not ASF1B-MT (P<0.01), and expression of ASF1B was markedly downregulated by overexpression of miR-214-3p (P<0.01; Figure 8A ). These results confirmed that miR-214-3p could directly bind to ASF1B and negatively regulate ASF1B expression. Next, vectors containing miR-214-3p-WT or miR-213-3p-MT were constructed and co-transfected with a SNHG3-overexpression plasmid into HepG2 cells. Luciferase activity was significantly inhibited (P<0.01) in cells containing the miR-214-3p-WT plasmid but not significantly altered in the miR-214-3p-MT group, whereas the expression level of miR-214-3p declined dramatically in HepG2 cells transfected with the SNHG3-overexpression plasmid (P<0.01; Figure 8B ). The results of the dual-luciferase reporter gene assay suggested that lncRNA SNHG3 could directly sponge miR-214-3p. Moreover, SNHG3 could inhibit miR-214-3p expression.

Notably, ASF1B, a key gene associated with HCC recurrence in this study, has not been previously reported in the HCC-related literature. Previous studies have confirmed that tumor infiltration is associated with the recurrence of HCC (8, 20). Therefore, we verified whether expression levels of ASF1B in HCC were correlated with immune infiltration, using the TIMER database. The results showed that expression levels of ASF1B were significantly positively correlated with immune infiltration of B cells (COR=0.489, p<0.001), CD8+ T cells (COR=0.345, p<0.001), CD4+ T cells (COR=0.336, p<0.001), macrophages (COR=0.442, p<0.001), neutrophils (COR=0.357, p<0.001), and dendritic cells (COR=0.464, p<0.001) in HCC ( Figure 9A ).

Subsequently, we further explored the relationships between the expression of ASF1B and gene markers of B cells, CD8+ T cells, neutrophils, dendritic cells, monocytes, Treg cells, macrophages, and natural killer cells in LIHC using TIMER. The results showed that ASF1B was positively correlated with CD19 and CD86 in B cells and with CD8B in CD8+ T cells. ASF1B in LIHC was also positively related to SIGLEC5, KIR2DL4, and CSF3R in neutrophils; ITGAX, STAT4, and STAT1 in dendritic cells; PD1 and CCR8 in Treg cells; and CD68 in monocytes (COR˃0.2). The results remained unchanged after adjustment for tumor purity and patient age ( Table 1 ). As tumor purity has an important influence on immune infiltration levels (21), partial immune genetic markers with COR greater than 0.3 after tumor purity adjustment were selected, and RT-qPCR was used to verify whether their expression was affected by ASF1B. The results showed significantly decreased expression of CD86, CD8, STAT1, STAT4, CD68, and PD1 in HepG2 and Huh-7 cells following ASF1B knockdown ( Figures 9B–G ). These results were consistent with those of the bioinformatics analysis and further confirmed that ASF1B expression in HCC was closely related to immune infiltration. Therefore, we concluded that ASF1B affects recurrence partly via immune infiltration.

In addition, as an important target of immunotherapy, PD-1 was closely related to ASF1B in our study. To further validate the role of SNHG3/miR-214-3p/ASF1B axis on PD-1, HCC cell was transfected with SNHG3 overexpression plasmid and si-ASF1B. The result of flow cytometry showed that PD-1 was significantly upregulated after transfection of SNHG3 overexpression plasmid (p<0.001), and ASF1B knockdown partially reversed the effect of lncRNA SNHG3 on PD-1 expression in HCC cell (p<0.001) ( Figure 9H ). Thus, these data indicated that lncRNA SNHG3 promoted the expression level of PD-1 by regulating ASF1B in HCC.

To test whether there was an interaction between immune infiltration and ASF1B, we performed a prognostic analysis based on the expression levels of ASF1B in different immune cell subgroups using Kaplan–Meier plotter. The results showed that HCC patients with high expression of ASF1B had higher recurrence rates in the decreased B cells [hazard ratio (HR)=1.77, P=0.016], decreased CD8+ T cells (HR=1.84, P=0.012), and decreased neutrophils (HR=1.93, P=0.01) subgroups ( Figures 10A, C, E ). However, there was no significant correlation between high ASF1B expression and recurrence in subgroups enriched in these immune cells ( Figures 10B, D, F ). On the contrary, high expression of ASF1B was closely associated with worse DFS in the enriched CD4+ T cells subgroup (HR=1.95, P=0.0084; Figure 10H ), but there was no such correlation in the decreased CD4+ T cells cohort (HR=1.21, P=0.47; Figure 10G ). These results indicate that high ASF1B expression may influence HCC recurrence partly owing to immune infiltration.