Up-Regulation of MiR-330-5p is Associated with the Advanced Clinical Stage of Hepatocellular Carcinoma

Purpose: A opposite expression levels and roles of miR-330-5p in hepatocellular carcinoma (HCC) have been reported in previous studies, so the clinicopathological implications and the prospective molecular mechanism of miR-330-5p in HCC require elaboration. Materials and methods: To examine the mRNA expression prole of hsa-miR-330-5p in HCC, a in-house RT-qPCR was performed, and the expression data was extracted from sequencing data and gene microarrays from the datasets of The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), ONCOMINE, ArrayExpress, and Sequence Read Archive (SRA). To have an overview of the clinical value of miR-330-5p, all possible data were integrated to calculate the standard mean difference (SMD) and area under the curve (AUC) from summary receiver operating characteristic curve (sROC). The target genes of miR-330-5p were also predicted and the relative signaling pathways were investigated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Results: According to the integrative analysis, the expression value of miR-330-5p in HCC was higher than in non-tumor tissue (SMD=0.47, 95CI%=0.31-0.63), and high expression of miR-330-5p was related with advanced HCC stage (SMD=0.27, 95CI%=0.08-0.46), poor differentiation (SMD=0.39, 95CI%=0.21-0.58) and poor prognosis (HR=1.91, Log-rank P=0.014). GLUD1, GOT1 and GLS2 were highly likely to be targeted by miR-330-5p, and the progression and prognosis of HCC might be inuenced via the miR-330-5p–GLUD1/GOT1/GLS2 axis. Conclusion: The result showed a high-expressed level of miR-330-5p was in the HCC tissue compared with non-tumor tissue, and the overexpression of miR-330-5p indicated

2.1.2 MiR-330-5p expression in HCC based on miRNA-microarray and miRNA-sequencing data Since the in-house RT-qPCR detection was carried out based on a small sample size with 52 paired HCC cases, we further mined data from miRNA-microarray and miRNA-sequencing to evaluate the miR-330-5p expression in HCC with various detecting methods. Datasets including The Gene Expression Omnibus (GEO), ArrayExpress, Sequence Read Archive (SRA), and peer reviewed publications were used to screen relevant miRNA-microarray and miRNA-sequencing data involving miR-330-5p expression in HCC (33)(34)(35) .
The keywords "hepatocellular carcinoma" and "HCC" were searched in the above datasets. The inclusion criteria were set as: 1) human samples; 2) miR-330-5p expression data being available; 3) non-cancerous liver controls being su cient. Finally, 22 miRNA-microarrays from GEO and one from ArrayExpress that met the requirements were recruited. The expression of miR-330-5p in HCC and adjacent liver control tissues was also obtained from TCGA miRNA-sequencing (Supplementary Figure 1) (36,37) . The violin charts were drawn using ggplot2 packages for each dataset to visualize the different expressions of miR-330-5p. Box plots were also inserted to re ect the average value and quartile. Student's t-test was used to evaluate the different levels of miR-330-5p between groups. The results of RT-qPCR, miRNA-microarray, and miRNA-sequencing were merged to show the overall expression level of miR-330-5p in HCC. The standard mean difference (SMD) and a 95% con dence interval (95%CI) were determined using meta R packages and STATA 12.0. A random effect model or xed-effect model was chosen according to the I 2 test and Chi-square test; high heterogeneity existed if I 2 >50% or P<0.05. Sensibility analysis was conducted to investigate the highly heterogeneous datasets. A funnel plot was also drawn to detect publication bias. The ROC curves were performed for each set, and the area under curves (AUCs) were calculated using SPSS 22.0 to access the predictive capacity of miR-330-5p for HCC. Next, sROC was performed for all datasets via STATA 12.0. The summarized AUC was calculated along with the sensitivity and speci city.

Analysis of the clinicopathological role of miR-330-5p in HCC
Information on related clinicopathological parameters, such as gender, age, tumor grading, tumor-regional lymph node-metastasis (TNM) staging, alcohol addiction, HBV infection, HCV infection, and albumin value was extracted to analyze the relationships between miR-330-5p and development of HCC based on miRNA-sequencing and in-house RT-qPCR data. The prognostic data of miR-330-5p was extracted from available datasets to perform an integrated prognostic analysis. The result was visualized by K-M curves, and the single hazard ratio (HR) and the summarized HR were evaluated.
2.2 Screening of overlapping target genes 2.2.1 Differently expressed genes (DEGs) of HCC To obtain the differentially expressed genes (DEGs) in HCC tissues for intersection with the predicted target genes of miR-330-5p and narrow the range of potential target genes of miR-330-5p, we obtained raw data from 63 high-throughput RNA-sequencing and microarrays (data not shown). Eventually, a total of 11 different detection platforms were included. The DEGs of RNA-sequencing data were calculated using Limma Voom and microarrays using Limma (logfc>1, adj P<0.05). After removing the batch effect of different platforms, 1,478 over-expressed genes appearing over eight times in 11 platforms were collected for subsequent use.

Prediction of miR-330-5p target genes
MiRWalk 2.0 was used for miRNA target prediction (38) . The target genes of miR-330-5p were predicted by 12 prediction software programs from MiRWalk. The target genes of miR-330-5p were intersected with the above 1,478 DEGs of HCC tissues. The overlapping genes were regarded as the most likely potential targets of miR-330-5p in HCC.

Signaling pathway enrichment analyses of miR-330-5p potential targets
The ClusterPro ler package was used to provide biofunctional enrichment data required for Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The GO items were selected at P<0.05 and FDR<0.05, and the KEGG items were selected at P<0.05. The top 10 signal pathways enriched by DEGs were listed according to the P value. The results of GO and KEGG analyses were exhibited in a bubble chart and chord chart via the "GO plot" package. The amount of enrichment genes in each term were represented by the size of the bubble, and the P value of each term was shown in gradual color. The genes enriched in the top three pathways were screened to analyze protein-protein interaction (PPI) networks. Protein interactions included direct physical interactions and indirect functional correlations; the interaction pairs in a PPI network were considered with a combined score >0.9 (39)(40)(41)(42)(43) .

Expression of miR-330-5p in HCC based on in-house RT-qPCR and other datasets
MiR-330-5p was markedly over-expressed in HCC tissues assessed by RT-qPCR, which was consistent with data shown by the high-throughput databases of different platforms, including Illumina Human v2 MicroRNA expression beadchip, Applied Biosystems Human MicroRNA Array A v2.0, Agilent-046064 Unrestricted_Human_miRNA_V19.0_Microarray, Affymetrix Multispecies miRNA-2 Array, and Agilent-070156 Human miRNA. Because of the uneven distribution of the results from a single platform, we integrated all the data and found that the SMD was 0.29 [0.05-0.54]. This implied that the expression level of miR-330-5p was predominantly higher in HCC than in non-tumor tissues under a random effect model with 1,315 cases of HCC and 880 control cases ( Figure 2A). Sensibility analysis was performed to discover the heterogenous chips ( Figure 2C). After removing the high-heterogeneity studies, the SMD rose to 0.47 and the trend of high miR-330-5p expression in HCC was more obvious ( Figure 3A). There was no publication bias in our study either before or after removing the high-heterogeneity studies.
The ROCs were used to judge the prospective diagnostic signi cance of miR-330-5p in each dataset ( Figure 4). To further validate whether miR-330-5p could be a predictor of HCC, 25 studies, including the in-house RT-qPCR, were summarized to calculate tp, fp, fn, and tn. The nal AUC of SROC was 0.76 [0.72-0.80], the sensitivity was 0.72, and the speci city was 0.69. This showed that miR-330-5p could be considered as a tool for discerning HCC from non-HCC liver tissues.

Clinicopathological role of miR-330-5p based on RNA-sequencing and RT-qPCR data
When looking at the clinical role of miR-330-5p in HCC based on a single cohort, we noted that there were close relations between the miR-330-5p expression and the deterioration of HCC. From the RNA-seq data, when the cases were younger and at more advanced tumor grading and clinical TNM stages, the miR-330-5p expression tended to increase ( Table 2). Similar up-regulated trends of miR-330-5p expression were con rmed in our own RT-qPCR data. When a tumor had metastasis or tumor thrombus in the blood vessels, miR-330-5p expression turned to be signi cantly higher. More importantly, several other highthroughput datasets could be leveraged to display the integrated role of miR-330-5p expression in tumor grading and TNM staging. The TNM stage was divided into advanced stage (stage III-IV) and early stage (stage I-II). Tumor grading was divided into the poorly-differentiated group (G3-G4) and the welldifferentiated group (G1-G2). The expression of miR-330-5p in 165 advanced-stage cases was higher than in 343 early-stage cases (SMD=0.27, 95%CI=0.07-0.47) ( Figure 8A). Similarly, the expression of miR-330-5p in 178 cases of poorly differentiated HCC tissues was higher than in 353 cases of welldifferentiated HCC (SMD=0.39, 95%CI=0.21-0.58) ( Figure 8B). The AUC of the miR-330-5p level for separate clinical stages was 0.88, the sensitivity was 0.67 (0.26-0.92), and the speci city was 0.94 (0.36-1.00). The AUC of miR-330-5p expression to differentiate tumor grades was 0.63, the sensitivity was 0.70 (0.57-0.80) and the speci city was 0.34 (0.16-0.59) ( Figure 9A-B).
3.1.3 Prognostic implications of miR-330-5p expression based on RNA-seq data After screening, since no prognostic data were available in other datasets, the RNA-seq data from TCGA was the only available source to assess the clinical role of miR-330-5p expression in HCC. Both OS and RFS data were extracted and used for prognostic analysis, and cutoff values of K-M curves were set to mean or quartile. The prognostic results were statistically signi cant in both the mean group and quartile group in the OS group; high expression of miR-330-5p was signi cantly related to unfavorable HCC prognosis, and the HR of the quartile group (HR=1.91) was higher than that of the mean group (HR=1.5, Figure 10A-B). However, in RFS, the expression of miR-330-5p was not signi cantly associated with HCC prognosis ( Figure 10C-D).

Targets of miR-330-5p in HCC
According to the result of the DEGs and miRWalk2.0, 499 genes overlapping the predicted miRNA target genes and the DEGs were achieved for GO and KEGG pathway analysis and miRNA-mRNA network construction ( Figure 11A). MiR-30-5p was over-expressed in HCC, so the low-expressed genes in the DEGs were used to intersect predicted target genes to narrow the possibilities. The target genes appeared at least ve times on every predictive platform, and DEGs appeared at least six times among all platforms.
The top 10 items in the biological processes, cell components, and molecular function categories were visualized in bubble charts (Table 4, Figure 11B). In biological processes, these overlapping genes are mainly involved in response to drugs, cholesterol metabolic processes, and positive regulation of gene expression. In cell components, these DEGs were enriched in the extracellular exosome, an integral component of plasma membrane, and external side of plasma membrane. In molecular function, these DEGs had the function of cytoskeletal protein binding, steroid hormone receptor activity, and transcription factor binding. The top 10 KEGG signal pathways were also listed (Table 5), and the chord diagrams were introduced ( Figure 12A). The results indicated that these DEGs were mainly enriched in alanine, aspartate, and glutamate metabolism, metabolic pathways, and cancer pathways. The three PPI networks were obtained by STRING ( Figure 12B-D). The hub genes in each PPI network were screened out via Cytoscape 3.6.1. The hub genes in the alanine, aspartate, and glutamate metabolism pathways were GLUD1, GOT1, and GLS2. Then, GLUD1, ALDH4A1, and GLS2 were considered as hub genes in the PPI network of metabolic pathways. In the network of cancer pathways, CDH1, CDKN2A, and IGF1 were screened as hub genes.

Discussion
The expression levels of miR-330-5p in HCC cells and tissues have been documented by several research groups. However, the results have been inconsistent. In this study, combining in-house real time RT-qPCR, miRNA-sequencing, miRNA-microarray, and in-silico integrated analysis, we noted that miR-330-5p levels were apparently upregulated in HCC tissues based on 1,095 cases from multiple research centers. More interestingly, we also pointed out that miR-330-5p may accelerate hepatocarcinogenesis via metabolic pathways, cancer pathways, and the signaling pathways of alanine, aspartate, and glutamate metabolism.
Contradictory miR-330-5p levels have been found in HCC cells and tissues by different groups. Downregulation of miR-330-5p was observed in HCC HepG2, HuH-7, Bel-7402, and SMMC-772 cells, compared to L-02 (24) . Two studies showed that ectopic miR-330-5p had a protective function to suppress the cell growth of HCC in vitro and in vivo. The restoration of miR-330-5p could also suppress the migration, invasion, and angiogenesis capacity of the HCC cells (24,29) . Hence, the above studies concluded that miR-330-5p may play a proactive role in the HCC process, as the increase of miR-330-5p could slow down or even block HCC development. However, higher expression levels of miR-330-5p have also been detected in HCC tissues and cells by three independent groups (26-28) . Among these three studies, 35 pairs of HCC and non-HCC cases from TCGA were detected by RNA-sequencing and miRNAmicroarray (26-28) . More importantly, two groups also veri ed the oncogenic role of miR-330-5p in HCC through in vitro and in vivo experiments (27,28) . Overexpression of miR-330-5p by transfection of miRNA mimics could assist the cell proliferation of HCC and the formation of HCC tumors, as well as promote cell in ltration ability (27,28) .
In the previous studies, The insu cient sample size in uenced the accuracy of the above incompatible results. To have a comprehensive view of the clinicopathological value of miR-330-5p levels in HCC, the current study combined various detecting methods, including in-house real time RT-qPCR, miRNAsequencing, miRNA-microarray, and SMD and sROC calculation. The RT-qPCR with the clinical samples showed a signi cantly higher level of miR-330-5p in 26 cases of HCC, with the relative expression of 7.51, 2.52 folds from the non-cancerous controls. To verify this nding, we mined the data from highthroughput datasets. Upregulation of miR-330-5p levels was consistently noted, as the SMD was 0.3 for 1,095 cases of HCC. With the integration of all available data, the nal SMD increased to 0.44, which was comparable to the reports of previous studies (26)(27)(28) . Hence, together with previous work, the current ndings con rm a marked increase in miR-330-5p expression levels in HCC tissues, compared to non-HCC liver tissues. This upregulation may play an oncogenic function in the tumorigenesis of HCC.
In previous studies, con icting results existed in the expression of miR-330-5p in HCC tissues. Beyond that, the association between miR-330-5p and the clinical parameters of HCC were contradictory. In the studies in which miR-330-5p was down-regulated in HCC, the expression value of miR-330-5p was negatively related to the progression parameters of HCC. The higher the TNM stage of HCC, the lower the expression value of miR-330-5p was (24,29) . However, the expression of miR-330-5p was positively associated with the TNM stage of HCC patients in the studies where miR-330-5p was up-regulated in HCC; in these studies, the higher expression of miR-330-5p was connected with poorer prognosis of HCC patients (26)(27)(28) . Although this contradiction existed in multiple studies, the researchers only used their own clinical data, including TNM stage and prognostic data; public data were not used in their studies. This lack of independent cohort veri cation was likely to be the cause of the contradictory results. In addition, Yu et al. reported that the expression of miR-330-5p in lower migrating HCC exosomes was greater than in higher ones (25) .
To obtain a comprehensive conclusion, clinicopathological parameters and prognostic data were collected from TCGA and GEO databases. In the TCGA data, the expression of miR-330-5p had no relationship with the TNM stage. However, there was a positive correlation between miR-330-5p levels and TNM stage in our PCR data. Because the TNM stage information existed in some GEO data, a metaanalysis was computed based on TCGA, GEO, and RT-qPCR data. The result showed that a higher expression of miR-330-5p was correlated with a more advanced TNM stage. Moreover, the AUC of miR-330-5p to differentiate advanced TNM stages from early stages showed that the speci city was 0.94, meaning the expression of miR-330-5p could be considered an indicator for tumor progression. Neither the included GEO data nor our PCR data had prognostic information, so the prognostic analysis of miR-330-5p in HCC was performed based on TCGA data. The results showed that miR-330-5p was closely related to patients' overall survival in HCC, and high expression of miR-330-5p led to poorer prognosis, but the expression of miR-330-5p did not affect patients' RFS. Previously, a relationship between miR-330-5p and RFS of HCC patients had been reported, but the literature used only its own data (27) and no public data. Hence, according to TCGA data and previous studies (27,28) , the overexpression of miR-330-5p could cause poorer overall survival of HCC patients.
The molecular function of miR-330-5p in HCC was still not clear, so we used various bioinformatics platforms to inquire into the putative molecular mechanism of miR-330-5p in HCC. First, the target genes of miR-330-5p were forecasted via miRWalk2.0, and the differential expression genes were calculated via 63 high-throughput RNA-sequencing and microarrays. All of the following operations were based on the intersection genes of the two gene sets. The result of KEGG pathway analysis revealed that miR-330-5p targeted genes might be involved in alanine, aspartate, and glutamate metabolism pathways, metabolic pathways, and pathways in cancer. In the top three pathways, there were two pathways associated with metabolism, particularly in glutamate metabolism. Numerous reports had indicated that the glutamate metabolic pathway was extremely important in the development of cancer. Glutamine can be broken down into glutamate, and glutamate can be further decomposed into α-ketoglutarate to promote the TCA cycle (44)(45)(46) . Thus, miR-330-5p may play an essential part in liver cancer; miR-330-5p may affect the occurrence and development of liver cancer by targeting related genes of glutamate metabolism to change the metabolic pathway. To further reveal this phenomenon, the hub genes (GLUD1 GOT1 GLS2) of alanine, aspartate, and glutamate metabolism pathways were screened. GLUD1, glutamate dehydrogenase 1, converted glutamate to alpha-ketoglutarate to promote the TCA cycle (47) . Glazer et al. found that the expression of GLUD1 in HCC tissues was lower than in non-tumorous liver tissues; combined with our study, we could conclude that targeting of GLUD1 by miR-330-5p led to low expression of GLUD1 in HCC (48) . Similarly, as a glutamate metabolism-related gene, GOT1 (glutamic-oxaloacetic transaminase 1), also called AST1, is an extremely important indicator of liver function and could also convert glutamate to alpha-ketoglutarate (49,50) . Nwosu et al. reported that the expression of GOT1 was down-regulated in poorly differentiated hepatocellular carcinoma cell lines, and the researchers found that GLUD1 was also down-regulated in HCC cell lines (51) . GLS2 (Glutaminase 2) is located on human chromosome 12, and it mainly encodes glutaminase (GA) which converts glutamine into glutamate in the liver (52) . The expression of GLS2 in HCC was evidentiarily lower than in normal liver cells, and the high expression of GLS2 could suppress the growth of HCC cells (47) . GLS2 could also be regarded as a molecular marker to assess the outcome of HCC patients (53) . These results show that all three hub genes could encode key enzymes in glutamine and glutamate metabolic pathways, and all three genes had low expression in HCC. Therefore, the highly expressed miR-330-5p was highly likely to target these three hub genes to form the regulatory axis of miR-330-5p-GLUD1/GOT1/GLS2.
There were some shortcomings in our study. First, although the three key genes (GLUD1, GOT1, GLS2) were highly likely to be the target genes of miR-330-5p, a dual-luciferase experimental con rmation was still de cient. Second, in vivo and in vitro work must still be conducted to assess the regulatory axis of miR-330-5p-GLUD1/GOT1/GLS2.

Conclusion
Contradictory results had been found in previous studies of miR-330-5p, so a series of integrated analyses were implemented to judge the facticity. The result showed a high-expressed level of miR-330-5p was obvious in the HCC tissue compared with non-tumor tissue, and the overexpression of miR-330-5p indicated poor progression and prognosis in HCC. In addition, GLUD1, GOT1 and GLS2 were highly likely to be targeted by miR-330-5p, and the progression and prognosis of HCC might be in uenced via the miR-330-5p-GLUD1/GOT1/GLS2 axis.

Declarations
Ethics approval and consent All the experimental protocols were approved by the Ethics Committee of the First A liated Hospital of Guangxi Medical University.

Consent
Both the clinicians and the patients signed the consent forms for the use of samples for the research.

Availability of data and materials
Our data is availible in the manuscript.

Competing interests
The authors declare that they have no competing interests.

Funding
The study was supported by   The correlation between clinicopathological variables and miR-330-5p expression in HCC from in house RT-qPCR. *P<0.05. Table 4 Figure 1 The           The expression value of miR-330-5p was divided into high expression group and low expression group based on the quartile, and the prognosis was compared between the two groups in RFS.

Figure 11
The overlapping genes and GO analysis. (A) The overlapping genes were screened between DEGs and the target genes of miR-330-5p. (B) Bubble plot was performed for GO analysis. Green bubbles were BP terms, red bubbles were CC terms and blue bubbles were MF terms. The size of bubble was associated with the gene counts.