Identifying Differentially Expressed MicroRNAs Between Cirrhotic and Non-Cirrhotic Hepatocellular Carcinoma and Exploring Their Functions Using Bioinformatic Analysis

Hepatocellular carcinoma (HCC) is the sixth most common cancer and the second leading cause of cancer mortality worldwide [1]. HCC exists as a high-incidence disease and is usually delayed in diagnosis. As such, HCC inflicts serious disease burden, especially for people in Asia, as the incidence rate is very high because of the high incidence rate of viral hepatitis [2, 3]. The mechanisms underlying the incidence and development of HCC have been found to be very complicated [4, 5]. Although significant effort has been placed towards understanding HCC, the mechanisms have not been well clarified. Consequently, HCC treatments are currently highly limited. Obviously, exploring the exact mechanisms underlying the incidence and development of HCC is highly significant.

To date, many studies, especially clinical studies, have shown that hepatic cirrhosis is an important risk factor in HCC. HCC patients with or without cirrhosis may present with distinctive conditions, as reflected by their tumor sizes, recurrence rates, and patients’ prognoses [6-10]. Prognostic evaluation has been recognized to fall short when only the tumor-related factors are considered because hepatocarcinogenesis is also considerably influenced by the status of the background liver, such as hepatic cirrhosis [11]. Most previous studies have identified this risk factor in the incidence of HCC, but the exact mechanisms underlying this phenomenon remain unknown and require consideration [12]. However, the patients with non-cirrhotic HCC have not been adequately studied [13].

MicroRNAs (miRNAs) are non-coding RNAs 21 to 23 nucleotides in length shown to help regulate the expression of many key genes. MiRNA sequences are relatively conservative in different species and associated with their important roles [14]. On account of the wide distribution and various categories, miRNAs have been recognized to serve numerous biological functions [15, 16]. Even so, the functions of most miRNAs are still unknown [17, 18].

To differentiate the characteristics of cirrhotic and non-cirrhotic HCC in terms of their miRNAs, the present study combined and analyzed the RNA-seq data and clinical data in The Cancer Genome Atlas (TCGA) database. Using bioinformatic analysis, this study also explored the possible roles played by miRNAs in both HCC subtypes. With the data from other public datasets (from Gene Expression Omnibus (GEO) datasets) and the data available at our research center, the expression and clinical significance of the selected four miRNAs were further validated. We successfully identified the differentially expressed miRNAs and their corresponding clinical significance in cirrhotic and non-cirrhotic HCC.

The raw data and clinical information were downloaded from the TCGA Research Network (https://cancergenome.nih.gov/). Only the samples with both miRNA-seq data and clinical information were included in our study. Patients were categorized into two groups (cirrhotic and non-cirrhotic HCC groups) based on the Ishak score [19]. With normalized data, the analysis to acquire differentially expressed miRNAs was conducted using the limma package in R. The fold changes (FCs) and P values of the expression of each miRNA were calculated. Considering the relatively small FCs in differentially miRNAs, this study identified the target miRNAs mainly by restricting the P value to lower than or equal to 0.05.

In our study, we collected 27 cases of HCC samples (including 14 non-cirrhotic HCC and 13 cirrhotic HCC) from the Second Affiliated Hospital of Chongqing Medical University. The follow-up data were collected by phone every two months until the patients died or could not be contacted. The median follow-up time was 36 months. The general characteristics (age, gender, HBV infection, and tumor size) of these patients between the cirrhotic group and non-cirrhotic group were comparable (data not shown). The entire study design followed appropriate medical ethics principles.

Chi-square tests, Unpaired t-tests, or Mann-Whitney tests were used to assess the relationship between differential miRNA expression and the patients’ clinical features. Kaplan-Meier survival analysis (Log-rank method) and univariate/multivariate Cox proportional hazard regression analysis were performed to explore the relationship between each parameter and the patients’ overall survival.

The target genes of the miRNAs were predicted using TargetScan [20], miRDB [21], DIANA [22], miRanda [23], and starBase [24] online analysis tools. The target genes of each miRNA were predicted using at least four of these tools. To further enhance the reliability of the bioinformatic analysis, Venn diagrams were used to determine the overlapping target genes. Given these overlapping genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were then performed using DAVID Bioinformatics Resources [25]. After setting the cutoff criteria of P≤0.05 for the results of gene enrichment, the genes most associated with the behavioral characteristics of the tumor were selected and employed to predict the function of these miRNAs. The interactions of the target genes were analyzed by String online tools [26]. Then, the interaction network between target genes was built using Cytoscape software (version 3.5.1) [27].

Total RNA was isolated from HCC samples using TRIzol reagent. iScript cDNA synthesis kit, SsoFast^TM EvaGreen Supermix, and the CFX96^tm Real-Time PCR system (Bio-Rad, Hercules, CA) were used for qRT-PCR reactions. MiRNA levels were normalized to U6 gene expression, and relative gene expression levels were analyzed using the ∆∆t method. The U6 primer was purchased from Guangzhou RiboBio Co., Ltd. The primers of other miRNAs were designed by our teamer and synthesized by Sangon Biotech. The primer sequences of the miRNAs were summarized in Table 1.

Statistical analysis was performed using IBM SPSS Statistics software (version 22.0) and R software. Considering the relatively small sample size in some comparisons, P values of less than or equal to 0.10 were considered statistically significant. The column diagram was graphed with Graphpad Prism software (version 7.0a).

In the TCGA database, a total of 377 HCC patients were found. Among these patients, 159 were excluded because of the lack of clinical information from the Ishak score. The retained patients (n=218) were then assigned to two groups on the basis of their different Ishak scores. The patients with scores lower than or equal to 4 were included in the non-cirrhotic HCC group, and the patients with scores over 4 were included in the cirrhotic HCC group. In the TCGA database, 372 of 377 patients had information on miRNA detection. After choosing patients with both miRNA data and enough clinical information, we finally included 215 cases (135 non-cirrhotic HCC patients and 80 cirrhotic HCC patients) for statistical analysis. The detailed clinical characteristics, including the age at initial diagnosis, gender, presence or absence of cirrhosis, family cancer history, AFP level, pathological stage, T stage, N stage, M stage, and existence of vessels in the tumor are shown in Table 2.

From the TCGA databases, we noted that 1046 miRNAs were analyzed for each patient. After restricting the cut-off criteria (P≤0.05), we found 41 differentially expressed miRNAs between the cirrhotic and non-cirrhotic HCC. These miRNAs included 32 upregulated and 9 downregulated miRNAs (non-cirrhotic compared with cirrhotic HCC, Fig. 1 and Table 3).

Herein, on the basis of the miRNA list sorted by P value, the top ten miRNAs were primarily selected for next-step analysis, and their differential expression is presented in Fig. 2A. After screening for the miRNAs with stable high expression in most of the HCC samples, we initially obtained five miRNAs (hsa-mir-3676, hsa-mir-1296, hsa-mir-23c, hsa-mir-149, and hsa-mir-95). However, the miRNA hsa-mir-3676 was excluded because the deep sequencing read mapping to the predicted hairpin did not support this miRNA processing. Therefore, 4 miRNAs were selected as the potential miRNA signatures for downstream analysis.

To identify the clinical significance of miRNAs, we evaluated the association between miRNA expression level and the clinical information of the patients with non-cirrhotic (Table 4) and cirrhotic HCC (Table 5). In the non-cirrhotic HCC, the expression of hsa-mir-149 was associated with the gender and T stage. Moreover, hsa-mir-23c was found to be linked to various clinical traits, such as age at initial diagnosis, gender, family cancer history, and the N and M tumor stages. Additionally, the miRNA hsa-mir-1296 was associated with age at initial diagnosis and family cancer history. For cirrhotic HCC, hsa-mir-95 was associated with AFP level. Hsa-mir-23c was also found to be strongly associated with the characteristics of cirrhotic HCC, such as pathologic stage and the T and N tumor stages. All of these relationships were statistically significant (P≤0.05).

The survival data were extracted from the TCGA databases. During this process, we found that 1 patient in the cirrhotic group lacked follow-up information. Thus, the survival analysis was conducted with only 214 patients (including 135 noncirrhotic and 79 cirrhotic patients). The survival rates of the patients grouped by different parameters, such as miRNA expression levels, age at initial diagnosis, gender, and AFP level, were compared. We noted that increased hsa-mir-1296 (Hazard ratio (HR), 2.12) and hsa-mir-149 (HR, 3.74) expression, as well as T (HR, 1.85), N (HR, 2.35), and M (HR, 2.57) tumor stages, were risk factors, whereas increased hsa-mir-23c (HR, 0.46) expression was a protective factor in the survival rate in the subgroup of non-cirrhotic HCC patients (Table 6, P≤0.05). We also found increased age at initial diagnosis (HR, 3.09) and existence of vessels in the tumor (HR, 3.62) were significant risk factors of the survival rate in the subgroup of patients with cirrhotic HCC (Table 7, P< 0.05). The survival curves describing the relationship of the different miRNA expression levels to the survival rate in both HCC types are shown in Figures 3 and 4. Additionally, high expression of hsa-mir-149 (HR, 3.06) was found to be an independent risk factor of prognosis in non-cirrhotic HCC, and the age at initial diagnosis (HR, 2.75) and existence of vessels in the tumor (HR, 3.16) were noted to be independent risk factors of prognosis in cirrhotic HCC (P< 0.05).

The target genes of the four miRNAs were predicted using at least four of the online databases we introduced above. A total of 45 overlapping genes of miR-1296, 313 overlapping genes of miR-23c, 235 overlapping genes of miR-149, and 28 overlapping genes of miR-95 were identified (Fig. 2B). Then, enrichment analysis of all these overlapping genes was carried out to elucidate the function of the potential target genes. Given our analysis, we found that the GO-biological process terms were mainly enriched in angiogenesis, regulation of cell adhesion, negative regulation of cell proliferation, negative regulation of translation, negative regulation of epithelial cell proliferation, and positive regulation of protein secretion (Fig. 5A). In addition, the KEGG pathways were significantly enriched in the pathways of cancer, focal adhesion, Ras signaling, miRNAs in cancer, PI3K-Akt signaling, and hepatitis B (Fig. 5B). Interaction analysis showed that CDC23, FASLG, and BCL2L11 genes presented with the most links with surrounding genes, and they were thought to be potential hub genes of the target genes (Fig. 5C).

Validation of the expression of the miRNAs. In one dataset (GSE10694) that included 38 non-cirrhotic HCC and 40 cirrhotic HCC patients, we analyzed the different miRNA microarray expression profiles of both types of HCC using 121 miRNA probes [28]. We found there were eight differentially expressed miRNAs (6.61%), which supported our conclusion of the existence of miRNA expression differences in cirrhotic and non-cirrhotic HCC (Fig. 6A). However, the analysis of this dataset did not provide the expression level of the selected four miRNAs, possibly because of the shortage of a gene microarray. The hsa-miR-149 expression in HCC tissue with cirrhosis was shown to be lower than that in HCC tissue without cirrhosis (Fold Change = 0.66, P = 0.0109), which was presented in a report using datasets (GSE31384) including 166 HCC patients [29] and was in agreement with our findings. To further present the result of miRNA expression in cirrhotic and non-cirrhotic HCC, another dataset (GSE69580) was used for analyzing the selected miRNAs (has-mir-1296, hsa-mir-149, and hsa-mir-95) [30]. Using this dataset, we found that only the expression of hsa-mir-95 was consistent with the results we found with the TCGA database (Fig. 6B). The inconsistent results might come from the relatively small sample size (including two non-cirrhotic and three cirrhotic HCC patients).

Considering this issue, we decide to conduct a qRT-PCR for evaluation of the miRNA expression in non-cirrhotic HCC (n=14) and cirrhotic HCC (n=13) patients. Similarly, we used the Ishak score to group patients. However, to improve the reliability of liver cirrhosis scoring, we set the diagnostic threshold at 6 points. As we found, the expression of mir-1296, mir-149, and mir-95 in non-cirrhotic HCC was higher than that in cirrhotic HCC (Fig. 6C, P≤0.10). Although the expression of mir-23c was lower in non-cirrhotic HCC than that in cirrhotic HCC, no significant difference was found (Fig. 6C, P> 0.10). Anyway, all of these results validated the differentially expressed miRNAs we found in the analysis of the TCGA database.

Validation of the clinical significance of the miRNAs. Here, we attempted to use another dataset (GSE21362) to analyze the expression of the selected four miRNAs because this dataset recruited patients with HCC with mild cirrhosis [31]. It was found that hsa-mir-1296, hsa-mir-149, and hsa-mir-95 expressions were not different between cirrhotic HCC tissue and corresponding normal tissue (Fig. 6D), which could be an indication of the lack of importance of these miRNAs in cirrhotic HCC. Furthermore, in another reported dataset (GSE6857) including 244 HCC patients [32], hsa-miR-149 expression in HCC tissue with metastasis was found to be higher than that in HCC tissue without metastasis (Fold Change = 1.279, P = 0.002), which supported the negative role found in our study. In the GSE12717 dataset (including 6 low grade cases and 2 high grade cases), the expression of hsa-mir-95 in the low grade group was higher than that in the high grade group [33] (Fig. 6E), which is not consistent with the results of our study and indicates the positive effect of hsa-mir-95 in HCC. Even so, this result is limited by the small sample size. Furthermore, the analysis of another dataset (GSE36915) that included 26 patients with T1 grade, 26 patients with T3a grade, and 16 patients with T3b grade [34] showed that the high expression of hsa-mir-1296 and hsa-mir-149 was positively associated with tumor invasion of the portal vein or hepatic vein (P< 0.05). The expression of hsa-mir-95 also trended to be positively associated with tumor invasion of the portal vein or hepatic vein, although it was not statistically significantly different between groups (Fig. 6F). All of these results suggested the clinical significance of the selected four miRNAs in HCC.

Validation of the prognostic significance of the miRNAs. To evaluate the prognostic significance of miRNA, we detect the miRNA expression differences between the patients who were dead or alive in the follow-up data. The analysis of the GSE6857 dataset showed that high expression of mir-149 was positively correlated with lower survival time only when the tumor had not metastasized (Fig. 6G, P< 0.05) [32]. Furthermore, mir-95 expression was not different between patients who were dead or alive at follow-up, which indicates the poor predictive ability of mir-95 for patients survival.

For the clinical samples collected in our study, we found that patients with cirrhotic HCC had lower postoperative recurrence-free survival (P=0.02) and overall survival (P=0.01) when compared with patients with non-cirrhotic HCC (Fig. 6H). Furthermore, we attempted to evaluate the correlations among miRNA expression and HCC patient survival using subgroup analysis. In the non-cirrhotic HCC subgroup, one recurrence and no deaths were found in the postoperative follow-up. Therefore, we did not evaluate the correlations among miRNA expression and non-cirrhotic HCC patient survival. In the cirrhotic HCC subgroup, we found no correlation between miRNA expression and postoperative recurrence-free survival or overall survival (Fig. 6I, P> 0.10). All of these prognostic indicators strongly supported our findings.

Recently, accumulating evidence has shown that hepatic cirrhosis is an independent risk factor of HCC development. Patients with or without cirrhosis may reasonably present unique clinical characteristics. Some patients with HCC lack a background of hepatic cirrhosis. In this case, investigating the differences between these HCC types would be meaningful. However, previous studies were not able to determine the origins of such differences. The present study is the first study to analyze the differences in miRNA expression and explore the relationships of such differences to corresponding clinical traits through miRNA-seq results from the TCGA database.

After analysis of the TCGA database, we found 41 differentially expressed miRNAs between both types of HCC. Both HCC types may involve distinct regulatory modes for miRNA. Furthermore, after selecting four representative miRNAs for the next-step study, we noted that the distinct expression patterns of these miRNAs had different associations with the extent of the patients’ clinical traits. Two miRNAs were found to be closely associated with the clinical traits, such as T, M, and N tumor stages, in patients with non-cirrhotic HCC. However, only hsa-mir-23c was associated with the T and N tumor stages of the patients with cirrhotic HCC. Besides, non-cirrhotic HCC owned more postoperative survival-associated miRNAs than cirrhotic HCC. Hence, we concluded that miRNA expression was associated more closely with non-cirrhotic than with cirrhotic HCC clinical traits. In short, the non-cirrhotic HCC traits appeared more dependent on the miRNA regulation network. Survival analysis also showed that both HCC types involve separate risk factors. Increased T, N, and M tumor stage were found to be risk factors for the prognosis of non-cirrhotic HCC, whereas the existence of vessels in the tumor is a risk factor for the prognosis of cirrhotic HCC. Most importantly, our study attained new information that can guide further studies. First, although we used the same variables for survival analysis of the two HCC types, no shared risk factors between the two types were observed. This result means that the differences in survival between these two tumor types may originate from different risk factors. Second, hsa-mir-149 expression was shown to be an independent risk factor of prognosis in non-cirrhotic HCC but not in cirrhotic HCC, and the age at initial diagnosis was noted as an independent risk factor in cirrhotic HCC but not in non-cirrhotic HCC.

Additionally, in order to validate the significance of miRNAs and their clinical correlation, we widely searched external datasets (from GEO datasets), including cirrhotic and non-cirrhotic HCC patients. A total of 7 extra datasets were analyzed. The results validated the selected four miRNAs and their clinical significance in cirrhotic and non-cirrhotic HCC. Furthermore, these results were also validated by our experiments.

After reviewing previous articles, we found that hsa-mir-149 has been well-studied in many different cancers. The gene serves as an anti-oncogene in colorectal cancer [35], HCC [36], laryngeal squamous cell carcinoma [37], and renal cell carcinoma [38]. In contrast, the gene functioned as an oncogene gene in T-cell acute lymphoblastic leukemia [39] and prostate cancer [40]. Hence, we may infer that hsa-mir-149 plays completely different roles across different tumor types. Intriguingly, hsa-mir-149 was considered an anti-oncogene for HCC in previous studies [36]; this finding contradicts our results. After reviewing the miRNA databases, we found that hsa-mir-149 contains two segments, named hsa-miR-149-5p (also named hsa-miR-149) and hsa-miR-149-3p (also named hsa-miR-149*). In fact, it was inferred that the previous study might regard hsa-miR-149-3p as an object of study. Meanwhile, our study treated hsa-miR-149-5p as a study object. In other words, this opposite function was possibly followed by investigation of different genes. This information should be paid adequate attention through further investigation. Although the important role of hsa-miR-149 has been well-studied in many tumor types, previous works did not describe its role in the two different HCC types explored in our research. In this study, we found that hsa-mir-149 (also called hsa-miR-149-5p) was an independent risk factor for prognosis in non-cirrhotic HCC but not in cirrhotic HCC. Therefore, additional studies exploring the association of hsa-mir-149 and HCC development should be conducted mainly in the non-cirrhotic group. Furthermore, the notion that age at initial diagnosis is an independent risk factor in cirrhotic HCC was supported by another study [41]. Even so, this risk factor may not substantially help clinical decision making because the patients commonly diagnosed at an old age may have a shorter survival time than younger patients.

A currently well-known miRNA mechanism involved in human cancer is that miRNA operates by combining with the 3′ untranslated regions of target complementary paired mRNA and then promotes their degradation or inhibits their expression [42-45]. Therefore, a definite relationship exists between miRNAs and their target mRNAs. Hence, in our study, four representative miRNAs were used to predict target genes. To obtain a concise prediction, we used at least four databases for each miRNA, and the obtained target genes were collected and intersected. Afterwards, the regulated target gene was found. Functional analysis of the GO and KEGG pathways was then conducted to search for the function of the four miRNAs. We found that these miRNAs were significantly enriched in various interesting signal pathways, such as angiogenesis, regulation of cell adhesion, negative regulation of cell proliferation, negative regulation of translation, pathways in cancer, focal adhesion, and Ras signaling. Differentially expressed miRNAs between both HCC types may target the mRNAs by changing these pathways and then result in different clinical and pathological traits between HCC types. After an interaction analysis of their target genes, we found that these differentially expressed miRNAs serve various functions by targeting certain mRNAs. Additionally, the interaction analysis showed that the CDC23, FASLG, and BCL2L11 genes are potential hub genes among these target genes. Thus, additional studies are necessary to investigate this regulatory relationship.

In fact, miRNAs may play many essential roles in the incidence and development of HCC, and these molecules may be used for diagnosing, treating, and prognosticating or therapeutically evaluating HCC patients. Recently, serum miRNAs have been recognized as significant indicators of a diagnosis of HCC [46]. One study showed that miR-203 overexpression may inhibit the augmented proliferation and lung metastasis of residual HCC [47]. Another study proposed that serum miR-10b-3p levels may be valuable for HCC diagnosis and survival prediction in sorafenib-treated patients [48]. MiR-135a-5p was also suggested to act as a mechanistic regulator of hepatic protein tyrosine phosphatase receptor delta (PTPRD) expression in patients with HCV and HCC [49]. In addition, HCC surveillance may help detect early incidence or recurrence after liver resection and consequently improve patients’ survival outcomes. However, HCC surveillance is not properly conducted and requires improvement [50, 51]. Selecting a target patient group instead of all patients for surveillance may be a suitable strategy, but few studies have been designed to obtain such a criterion [52]. In our study, we successfully determined that HCC patients with different expressions of miRNAs would present with different prognosis. We also showed this correlation would be different in cirrhotic and non-cirrhotic HCC types. Therefore, HCC patients with differentially expressed miRNAs and with or without cirrhosis would be part of the criterion for accepting extra surveillance to obtain timely therapy and improved prognosis.

The major limitation in this study is that many unreported miRNAs were not evaluated experimentally despite being differentially expressed between the two HCC types. Even so, our findings remain helpful for guiding subsequent experimental studies.

In conclusion, our study successfully identified that the clinical characteristics of cirrhotic and non-cirrhotic HCC patients differ in relation to the expression of various miRNAs. When compared with cirrhotic HCC, non-cirrhotic HCC appears more dependent on the regulation of the miRNA network. Differentially expressed miRNAs may play different roles between the two HCC types through different regulatory pathways. Even so, further experimental studies are still needed to validate these claims.

The authors confirm that there is no conflict of interests.