The ‘limma’ package was obtained from Bioconductor (bioconductor.org/biocLite) and was employed to analyze the data via R software (version 4.0.0). RNA sequencing (-seq) data was imported into the R environment and converted into a data structure supported by ‘DESeq2’ package (bioconductor.org/biocLite). Original data was preprocessed by standardizing the data, removing low-expression genes and batch effects. R language ‘limma’, ‘dplyr’ and ‘tidyr’ packages (bioconductor.org/biocLite) were used for differential expression gene analysis. R package ‘ggplot2’ was applied to visualize the results. Heatmap results of the miRNAs and differential gene expression were produced with TBtools (version 1.09876) (34).

The miRNA expression database was established by searching GEO Data Sets (ncbi.nlm.nih.gov). The inclusion criteria were samples must be pathologically confirmed as HCC patients and miRNA expression profiles are detected in tumor and adjacent normal tissues. Through screening, two datasets containing samples from tissues of patients with HCC were included: GSE67882 and GSE69580. A total of 13 HCC and nine non-tumor tissues were included in the present study. These patients with HCC had diagnosis confirmed through pathological examination. The ‘limma’ R package was applied to compare expression levels of miRNAs between HBV-related HCC and the corresponding non-tumor tissues from the two series of matrix files. After data normalization, expression levels of the miRNAs in the two datasets were analyzed.

The prognostic values of DE-miRNAs in HCC were analyzed using the Kaplan-Meier Plotter database (kmplot.com/analysis/index.php?p=service&cancer=liver_mirna). Through comparing the relationship between miRNA expression levels and the survival rate of HCC patients, the DE-miRNAs with an impact on the prognosis of patients with HBV-related HCC were selected for further analysis.

miRNet (mirnet.ca/miRNet/home.xhtml), an easy-to-use analysis tool based on an online network, was used to predict the miRNA target genes (35). miRNet integrates data from 11 high-quality miRNA-target interaction databases and displays the results through a network visualization system. Gephi software (version 0.9.2, gephi.org/users/), a network visualization software based on the JAVA working environment, displayed the regulatory association between miRNAs and their target genes.

Cytoscape software (version 3.7.1, https://github.com/cytoscape/cytoscape) was used to evaluate the miRNA hub genes. For all DE-miRNA target genes, Molecular Complex Detection plug-in (version 3.3.10, http://apps.cytoscape.org/) was used for scoring according to the interaction between the genes. The higher the score, the more key the role the hub genes play in the entire gene network.

Database for Annotation, Visualization and Integrated Discovery (DAVID; david.ncifcrf.gov/home.jsp) was applied for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DE-miRNA target genes. For GSE101728 data, containing the expression profiling of lncRNA and mRNA in HCC (36), the enrichment of KEGG pathway according to the top 25% of the highest and lowest GRB2 expression patients' genes expression profiling was assessed.

HBV-associated experimental data were screened through the GEO database. In the GSE100400 dataset, the researchers combined circularized chromosome conformation capture (4C) with RNA-seq and chromatin immunoprecipitation (ChIP)-seq to illuminate the nuclear landscape associated with HBV episomes and HBx (37). To the best of our knowledge, the aforementioned study is the first to use 4C to identify regional virus-host genome interactions. By using R programming language tool, expression matrix of RNA-seq data in GSE100400 was obtained and the difference in gene expression was calculated.

UALCAN database (ualcan.path.uab.edu/) is an easy-to-use interactive web portal that for in-depth analysis of TCGA gene expression data. UALCAN includes TCGA RNA-seq and clinical data from 31 types of cancer (38). This was used to select genes with significant differences in expression and prognosis between tumors and corresponding adjacent tissues.

The Tumor Immune Estimation Resource (TIMER) database (cistrome.shinyapps.io/timer/) is designed to investigate the molecular characteristics of tumor-immune interactions and analyze the correlation between gene expression and the abundance of immune infiltrates (39). This was used to screen genes related to immune infiltration.

Human hepatoblastoma cell lines (HepG2 and HepG2.2.15) were purchased from Guangzhou Cellcook Biotech Co., Ltd. and were authenticated via STR profiling. The hepatoblastoma cells were grown in Minimal Essential Medium supplemented with non-essential amino acids and 10% fetal bovine serum (all WISENT, Inc.). Cells were placed in a cell incubator (Thermo Fisher Scientific, Inc.) at 37˚C with 5% carbon dioxide. miRNA-93 mimics and inhibitors and negative controls (Guangzhou RiboBio Co., Ltd.) were prepared for cell transfection. The confluence of the cells was maintained at ~60% before transfection. Cells were starved in serum-free Minimal Essential Medium (WISENT, Inc.) at 37˚C with 5% carbon dioxide for 2 h and transfected with 1 µg miRNA-93 mimics and inhibitors and their respective blank vectors using Lipofectamine 2000 reagent (Invitrogen; Thermo Fisher Scientific, Inc.) for 36 h at 37˚C. Subsequent experiment was conducted within 72 h. The sequences of the miRNA-93 mimics/inhibitors are displayed in Table SI. The sequence of mature miRNA-93 was 5'-CAAAGUGCUGUUCGUGCAGGUAG-3'. The sequence of mimics negative control was 5'-UUGUACUACACAAAAGUACUG-3' and inhibitor negative control was 5'-CAGUACUUUUGUGUAGUACAA-3'.

Total RNA of the cells was extracted using TRIzol (Invitrogen; Thermo Fisher Scientific, Inc.). RT was performed using a PrimeScript^™ RT reagent Kit (Takara Bio, Inc., No RR047A) according to the manufacturer's instructions. A SYBR Premix Ex Taq II kit (Takara Bio, Inc.) was used for RT-qPCR with thermocycling conditions as follows: initial denaturation 95˚C, 30 sec; 40 cycles of denaturation 95˚C, 30 sec, annealing at 60˚C 30 sec and extension 72˚C, 30 sec. StepOne Plus system (Applied Biosystems; Thermo Fisher Scientific, Inc.) was employed. Target gene expression was normalized using the 2^-ΔΔCq method (40). The sequences of mature miRNAs are presented in Table SII; RT-qPCR primers (Shanghai GenePharma Co., Ltd.) are displayed in Table SIII.

All experimental data with at least 3 times are presented as the mean ± standard deviation using GraphPad Prism (version 5.0; Dotmatics) and the statistical analysis was performed using SPSS 20.0 (IBM Corp.). Differences between the two groups was analyzed using paired t tests. For >2 groups in TCGA database, One-way ANOVA was used followed by Fisher's least significant difference post hoc test. The correlation analysis between two genes was evaluated using the Spearman rank correlation coefficient. All plots were drawn with GraphPad Prism software. For the results of miRNA microarray, P<0.05 and |log₂fold-change (FC)|> 1 were considered to define the DE-miRNAs. The association between expression of DE-miRNAs and the prognostic value was analyzed using log-rank test. P<0.05 was considered to indicate a statistically significant difference.

The original data of miRNA arrays GSE67882 and GSE69580 were obtained and normalized (Fig. S1). A total of 13 HCC and nine non-tumor tissues were included for further analysis. The ‘limma’ R package was used to generate volcano plots of miRNA expression (Fig. S2). The top 10 up- and downregulated genes are displayed as heatmaps (Fig. 1A and B). Intersection analysis showed that there were 63 up- and 22 downregulated genes (Fig. 1C and D). The candidate DE-miRNAs were defined according to two indicators |log₂FC|>1 and P<0.05. Following screening, nine up- and one downregulated DE-miRNA were obtained (Table I).

To evaluate candidate DE-miRNAs, their prognostic value was analyzed. Kaplan-Meier Plotter database showed that among nine upregulated miRNAs, only four miRNAs, namely miRNA-93, miRNA-106b, miRNA-210 and miRNA-374b, had a significant effect on prognosis of patients with HCC (Fig. 2A-D). The prognostic analysis of miRNA-125b showed that the lower the expression of miRNA-125b, the worse the prognosis (Fig. 2E). Finally, five miRNAs whose expression levels and prognosis were significantly different between HCC and normal tissue were identified as candidate DE-miRNAs for further study.

The target genes of DE-miRNAs were predicted in miRNet database. miRNA-93 had 1,220, miRNA-106b had 1,091, miRNA-210 had 124, miRNA-374b had 234 and miRNA-125b had 432 potential target genes (Fig. S3A and B; Table II). Protein-protein interaction network analysis was performed on predicted target genes of DE-miRNA using the STRING database and Cytoscape software (version 3.7.1; Fig. 3A and B). The top 20 resulting hub genes are displayed, as well as their scores (Table SIV).

The present study evaluated the terms involved in these target genes through GO analysis. DAVID (version 6.8) database was used to perform GO enrichment analysis of 2,669 predicted target genes of up- and 432 predicted target genes of downregulated DE-miRNAs. The results were divided into three categories: Biological process, cellular component and molecular function. The top ten terms are listed according to P-value (Fig. 4A and B). KEGG pathway enrichment analysis of target genes was performed using the DAVID database to evaluate the signaling pathways. The top 20 terms are displayed according to P-value (Fig. 4C and D). Notably, three signal pathways (‘pathway in cancer’, ‘hepatitis B’ and ‘microRNAs in cancer’) were co-enriched, which were associated with HBV-related HCC. To identify potential target genes regulated by HBV, these signaling pathways were selected for intersection analysis, resulting in 11 and 8 sets of intersecting genes for up- and downregulated genes, respectively (Fig. 4E and F). For these predicted target genes, experimental data from GEO database were further investigated. In GSE100400 dataset, RNA-seq, ChIP-seq and 4C methods were combined to identify target genes regulated by HBV (37). Since this dataset provided a reliable and well-analyzed differential gene expression profile, it was cross-analyzed with the present findings to predict the potential miRNA-mRNA regulatory relationship network in HBV-related HCC. An intersection analysis of 11 up- and 8 downregulated target genes with the results of GSE100400 was conducted. A total of seven HBV-associated intersecting genes was obtained: CCND1, MYC, E2F1, E2F2, E2F3, GRB2 and CCNE1. The sequencing expression profiles of these genes in the GSE100400 dataset is shown in Fig. 4G. The expression levels of the DE-miRNA target genes were analyzed via UALCAN cancer database (Fig. S4). The present study summarized the potential miRNA-mRNA regulatory association network model in HBV-associated HCC (Fig. 5).

To validate the results of the bioinformatics analysis, the expression levels of DE-miRNAs and target genes in HepG2 and HepG2.2.1.5 (a cell line that stably expresses HBV) (41) cells were detected. Compared with those in HepG2 cells, miRNA-93, miRNA-106b, miRNA-210 and miRNA-374b showed significantly higher expression levels in HepG2.2.1.5 cells, while miRNA-125b showed significantly lower expression (Fig. 6A and B). Similarly, compared with HepG2 cells, expression levels of E1F1, E2F3, GRB2, CCND1 and CCNE1 were higher than those in HepG2.2.15 cells, while expression levels of E2F2 and MYC were not significantly different (Fig. 6C).

Base pairing analysis demonstrated that miRNA-93 had binding sites not only in the 3' but also the 5' UTR of GRB2. A total of nine potential binding sites was predicted between miRNA-93 and 5' UTR of GRB2 (Fig. 6D). As aforementioned, in HepG2.2.15 cells with HBV overexpression, the expression levels of miRNA-93 and GRB2 were upregulated. According to the predicted binding site, miRNA-93 has a regulatory effect of activation by targeting the GRB2 5' UTR. Therefore, more attention was paid to the active regulation between miRNA-93 and GRB2.

Transfection was performed to detect the effect of miRNA-93 on GRB2 expression levels in vitro. Mimics of miRNA-93 were transfected into HepG2 cells and changes in GRB2 were detected. With the increase in miRNA-93, the expression of GRB2 also increased (Fig. 6E). In addition, transfection of miRNA-93 inhibitors in HepG2.2.15 cells confirmed that the expression of GRB2 decreased with the knockdown of miRNA-93 (Fig. 6E). Thus, GRB2 was a target gene of miRNA-93 and was positively regulated by miRNA-93.

TCGA database demonstrated that expression levels of GRB2 in HCC were significantly higher than those in non-tumor tissues (Fig. 7A). Expression of GRB2 was associated with lymph node metastasis in patients with HCC (Fig. 7B). The higher the GRB2 expression, the lower the overall survival rate of patients with HCC (Fig. 7C).

The present study analyzed the role of GRB2 in occurrence and development of HCC. After GSE101728 data normalization, the present study assessed enrichment of the KEGG pathway between the top 25% of the highest and lowest GRB2 expression (Fig. 7D-I). Primary immunodeficiency signaling pathway had the smallest P-value in the low GRB2 expression group.

The association between GRB2 expression and immune infiltration levels in HCC were assessed using the TIMER database. Spearman's correlation was assessed based on microarray expression levels. High GRB2 expression had positive correlations with the infiltration levels of B (R=0.324) and CD4⁺ T cells (R=0.415), macrophages (R=0.483), neutrophils (R=0.429) and dendritic cells (R=0.386) in HCC (Fig. 7J). GRB2 expression had a weakly positive correlation with the infiltration of CD8⁺ T cells (R=0.222) and no significant correlation with tumor purity (R=0.026). These findings suggest that GRB2 may serve a specific role in immune infiltration in HCC, especially for CD4⁺ T cells, macrophages and neutrophils.