Noninvasive Identification of Immune-Related Biomarkers in Hepatocellular Carcinoma

Primary liver carcinoma is one of the most common malignant tumors with a poor prognosis. This study aims to uncover the potential mechanisms and identify core biomarkers of hepatocellular carcinoma (HCC). The HCC gene expression profile GSE49515 was chosen to analyze the differentially expressed genes from purified RNA of peripheral blood mononuclear cells, including 10 HCC patients and 10 normal individuals. GO and KEGG pathway analysis and PPI network were used, and the enrichment of the core genes out of 15 hub genes was evaluated using GEPIA and GSEA, respectively. We employed flow cytometry to count mononuclear cells to verify our opinions. In this analysis, 344 DEGs were captured, including 188 upregulated genes and 156 downregulated genes; besides that, 15 hub genes were identified. GO analysis and KEGG analysis showed the DEGs were particularly involved in immune response, antigen processing and presentation, and infection and inflammation. The PPI network uncovered 2 modules were also mainly involved in immune response. In conclusion, our analysis disclosed the immune subversion was the major signature of HCC associated closely with JUN, VEGFA, TNFSF10, and TLR4, which could be novel noninvasive biomarkers in peripheral blood and targets for early diagnosis and therapy of HCC.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common malignancies, especially in the aged, which accounts for approximately 90% of all primary liver cancers severely threatening public health [1]. e mechanism of HCC is a complex process associated with the incremental accumulation of gene mutation, giving rise to abnormal immune subversion, cell cycle, and angiogenesis [2][3][4]. As for immune subversion, effector immune cells could execute immune control of HCC, which efficiently decrease malignant transformed cells. However, progression of HCC clearly certifies failure of tumor immune control suggesting inhibition of anticancer immune responses [5]. Especially, tumor-related mononuclear cells collaborate within an inflammatory network, which result in the immune privilege in the tumor environment [6]. erefore, immunosuppressive mononuclear cells are equivalent to heterogeneous cell lines, including lymphocytes and monocytes cooperating by direct cell contact, secretion of cytokines, or production of extracellular matrix, which lead to the suppression of the immune response in the tumor milieu [7].
Currently, imageological examination and pathological biopsy are the conventional diagnostic methods of HCC [8]. However, imaging displays poor specificity, and pathological biopsy is an invasive method which may result in iatrogenic injury [9]. erefore, serum biomarkers are routinely used for tumor diagnostic. For example, alpha-fetoprotein (AFP) has been widely used in clinical practice [10]. Although many studies have reported the accuracy of AFP for HCC, solely AFP still has some false-positive or false-negative rate [11]. Hence, the identification of specific and sensitive biomarkers is necessary in order to achieve accurate diagnosis and treatment of HCC as early as possible, especially noninvasive biomarkers.
High-throughput gene microarray is increasingly being widely used, which can analyze cancer and noncancer samples indicating us tumor-related genes at multiple levels from molecular diagnosis and pathological classification to therapeutic evaluation and prognosis prediction, as well as drug sensitivity and neoplasm recurrence [12][13][14]. However, the use of microarrys in clinical application is restricted by countless number of genes identified by gene profiling, lack of both repeatability and independent verification, and requirement for complex statistical analyses. Moreover, most of the microarrys are based on the genes in tissues which are difficult to detect except by invasive methods [15]. erefore, in order to put these expression profiles into clinical applications as soon as possible, it is necessary to identify an appropriate amount of serum genes and develop a suitable way that can be done by routine assay.
In this study, we downloaded the HCC gene expression profile GSE49515 in the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/), an online public collection database for microarray data and used GEO2R online software to compare gene expression profiles of tumor cells with normal liver cells to identify differentially expressed genes (DEGs). en, we constructed the protein-protein interaction (PPI) network of the DEGs and selected 15 hub genes according to a high degree of connectivity. Following this, we analyzed gene ontology (GO) and pathway enrichment including the biological process (BP), molecular function (MF), cellular component (CC), and KEGG pathway of the DEGs. Moreover, we performed two modules and confirmed their enriched pathways. e core genes out of the 15 hub ones were found, and the interactions between any of them were detected with the help of GEPIA. After the analysis of the core status and biological function of any hub gene, we performed flow cytometry (FCM) to count mononuclear cells, confirming the findings.

Identification of DEGs and Hub Genes.
A comparison of 10 HCC samples with 10 normal samples in our study was performed by employing the GEO2R online analysis tool based on P value < 0.05 and logFC ≤ − 2 or logFC ≥ 2 criteria. A total of 344 DEGs were picked up after analyzing GSE49515, 188 of which were upregulated while 156 were downregulated (Figure 1(b)). e expression levels of the top 50 DEGs were displayed in a heat map to visualize the results (Figure 1(a)).

GO Function-and KEGG Pathway-Enrichment Analysis.
To gain a more extensive and in-depth knowledge of those selected DEGs, we use DAVID to analyze significantly enriched GO function and KEGG pathways. After inputting all of the DEGs to DAVID online analysis tool, we obtained the GO analysis of these upregulated DEGs and downregulated DEGs. e results showed that these DEGs were mainly enriched in biological processes (BP), including apoptotic process, immune response, and inflammatory response, among which were positive regulation of NF-κB transcription factor activity and cell-cell signaling for downregulation, positive regulation of angiogenesis, negative regulation of cell proliferation, positive regulation of cell proliferation, mitotic spindle organization, and neutrophil chemotaxis for upregulation. Concerning molecular function (MF), the downregulated DEGs were particularly related to receptor binding, iron ion binding, haptoglobin binding, oxygen transporter activity, and peroxidase activity, while the upregulated DEGs were mainly implicated with nucleotide binding and ubiquitin protein ligase binding. Besides, GO cell component (CC) analysis indicated that the downregulated DEGs were mainly enriched in cytosol, extracellular exosome, Golgi membrane, blood microparticle, and nuclear chromosome (telomeric region) and the upregulated DEGs were principally enriched in nucleus, nucleoplasm, platelet alpha granule and extracellular space (Table 1).
Afterwards, we analyzed the most significantly enriched KEGG pathway of the upregulated and downregulated DEGs, which is shown in Table 2. e downregulated DEGs were involved in measles, influenza A, rheumatoid arthritis, antigen processing and presentation, and legionellosis, while the upregulated DEGs were involved in bladder cancer, rheumatoid arthritis, malaria, herpes simplex infection, and osteoclast differentiation. e scatter plots in Supplement 1 (A, B, and C) show a GO and KEGG pathway-enrichment plot of HCC.

Hub Genes and Module Screening from PPI Network.
Besides, 15 hub genes from Cytoscape software were identified in accordance with a high degree of connectivity selected (Table 3). We built the PPI network of the top 15 hub genes via the information of the STRING protein query from public databases (Figure 2(a)). e top 15 hub genes with a higher degree of connectivity are as follows : JUN, IL8,  VEGFA, TLR4, IFNG, TNFSF10, EHHADH, ATF3, FUS,  DUSP1, HSPA1A, CUL1, FPR2, POLR2H, and RHOB. Based on the GO function and KEGG pathway analysis of the profile, we uncovered JUN, VEGFA, TNFSF10, and TLR4 that were enriched in immune response-related pathway.
Moreover, we used MCODE plug-in to detect the highest modules in the PPI network. We chose the top 2 modules, and GO function-and KEGG pathway-enrichment analysis indicated that Module 1 was related to immune response, apoptotic process, and epithelial cell migration, while Module 2 was associated with a signaling pathway and cellular response to various substances (Figures 2(b) and 2(c), Table 4).

e Kaplan-Meier Plotter and Expression Level of Hub
Genes.
ese hub genes in PBMC which can be the biomarkers for early diagnosis may not be easy to detect in liver tissue. en, we employed DAVID to analyze the correlation of 15 hub genes. We found JUN, IFNG, VEGFA, TLR4, and TNFSF10 are the 5 high-degree-of-connectivity genes which are fully associated with immune response, inflammatory response, and HIF-1 signaling pathway (Table 5). en, in order to confirm the most relevant hub genes, we used correlation analysis in GEPIA, and we detected JUN and VEGFA, JUN and ATF3, and JUN and RHOB are distinctly correlated (P value � 0, R � 0.47; P value � 0, R � 0.71; P value � 0, R � 0.69), which means JUN may be the core gene of HCC (Figures 4(e)-4(g)).

Gene Set Enrichment Analysis.
In order to make the further function of the hub gene clear, GSEA was used to map into GO analysis and KEGG pathways database. Under  the cutoff criteria nominal P value < 0.05, │enrichment score (ES)│ > 0.6, and gene size ≥ 100, six functional gene sets were enriched in total, which were particularly centralized on pathway associated with immune response and inflammatory response. Six pathways were "inflammatory response to antigenic stimulus," "regulation of T cell migration," "antigen processing and presentation via MHC class IB," "negative regulation of B cell activation," "negative regulation of NF-kB signaling," and "positive regulation of interleukin-1 production" ( Figure 5).

Identification of Biomarkers.
To confirm the results we have stated above, we test the mRNA level of four key genes we predicted (JUN, TLR4, VEGFA, and TNFSF10). We found the mRNA level of JUN, TLR4, VEGFA, and TNFSF10 in the PBMC of HCC patients was significantly upregulated in the PBMC of HCC patients (P < 0.05), of which VEGFA increased obviously (P < 0.01) (Figures 4(a)-4(d)).

Immune Subversion in HCC.
According to our prediction, the negative regulation of immune cells in HCC patients' peripheral blood significantly occurred detected by FCM (Figures 6(a)-6(h)). In detail, compared with the healthy subjects, the levels of T lymphocytes in peripheral blood, both helper T cells and cytotxic T cells, were significantly lower in HCC patients, which is mainly the immune mechanism of tumor patients (P < 0.05) ( Figure 6(k)). Meanwhile, the B lymphocytes as well as NK cells also decreased in HCC patients, especially NK cells (P < 0.05) (Figures 6(i) and 6(j)). Taken together, our experiment of FCM demonstrated that T cell migration and B cell activation in adaptive immune and NK cells inhibition in innate immune were important mechanisms and could do further research in future.

Discussion
In recent decades, the morbidity and mortality of HCC have been increasing worldwide. Although the early diagnosis and treatment have developed a lot recently, the overall survival rates of HCC is still poor [16]. erefore, the sensitive and specific biomarkers for HCC are urgently necessary. High-throughput studies can develop the thorough exploration of the vital mechanisms which lead to HCC. In our study, we identified DEGs between 10 HCC samples and 10 normal samples from the GEO database of GSE49515. In order to increase the statistical power of these DEGs, we defined that the absolute value of the logarithm (base 2) fold change (logFC) greater than 2 and a total of 344 DEGs were captured, including 188 upregulated genes and 156 downregulated genes. In order to have an in-depth detection of these DEGs, we employed GO function, KEGG pathway, PPI network, and connectivity analysis of these DEGs, via which we found that HCC-related genes and pathways have great importance in cancer initiation and progression.
ere were plenty of mechanisms uncovered to contribute to the development of HCC, but the predominant mechanism implicated with tumorigenesis is still controversial, thereby causing difficulties to the diagnosis and treatment of HCC. GO term enrichment and PPI analysis in our study disclosed that downregulated DEGs were mainly associated with immune response. Our experimental results detecting the mononuclear cells in peripheral blood mononuclear cells of HCC patients and healthy individuals further confirmed the important role of immune subversion.
Immune response is well recognized to play a vital part in the initiation and progression of carcinogenesis because the development of HCC obviously records failure of tumor immune control which stands for immune subversion by the tumor environment [2]. To our knowledge, protective immune surveillance of tumor is mainly conducted by tumordirected NK cells and lymphocytes, which can effectively identify and eradicate malignant cells [6]. On the one hand, for the reason that innate immune cells are capable of eliminating malignant cells and pertaining to the first-line defense to restrain tumor initiation and progression [17], they act as an essential player in the immunological surveillance. On the other hand, adaptive immune cells including B lymphocytes and T lymphocytes develop along with innate immune responses, which finally target tumor-associated antigens (TAAs) [18]. To be more specific, once the immune cells in peripheral blood are unable to generate or be recruited, specific immunity and nonspecific immunity function will degrade, eventually causing the rapid growth of carcinoma tissue. e quantity of immune cells is associated with many immune-related genes which may be differentially expressed during tumor manifestation and progression.
NK cells (CD3− /CD56+) reside in the liver and account for about 30% to 50% of the hepatic lymphocytes in humans [17]. erefore, NK cells constitute a corresponding effector cell population of innate immune cells contributing to tumor surveillance within the liver. In our data from FCM, we identified that NK cells in HCC patients' peripheral blood were significantly lower than those in the blood from normal people. Additionally, our GO analysis indicated that the downregulated DEGs were mainly related to immune response including immune cells manifestation and recruitment Journal of Oncology 5 and release of cytokines. We found that there were 4 genes (JUN, VEGFA, TNFSF10, and TLR4) closely related to the downregulation of immune response which occurs in HCC. Mgrditchian et al. [19] found that the phosphorylation of JUN could induce NK cell infiltration into the tumor bed by inducing the transcription of CCL5, which eventually results in targeted autophagy. However, NK cells produce vascular endothelial growth factor A (VEGFA) in tumor tissues, which may enhance the formation of tumor by the way of angiogenesis. e presence of VEGFA-secreting NK cells is associated with a poor prognosis and depends on the type and stage of the tumor [20]. In our GO analysis, we detected that  Our result suggested that innate immune cells, especially the function of NK cells, are downregulated via some hub genes in HCC.
Adaptive immunity including humoral immunity and cellular immunity is highly specific elimination of transformed cells, which protects from tumor manifestation and progression. B lymphocytes, CD4+ T lymphocytes, and CD8+ T lymphocytes are vital members in adaptive immunity [25]. Equally, we employed FCM to detect the B lymphocytes and two T lymphocyte subsets which verified the reduction of adaptive immune function in HCC. We predicted the mechanisms of B lymphocytes and T lymphocyte declining may also be associated with the four hub genes (JUN, VEGFA, TNFSF10, and TLR4) via GO and KEGG analysis. c-JUN NH 2 -terminal kinase (JNK) signaling pathway was implicated in various T cell functions. e JNKs are synergistically activated by stimulation of the TCR with antibodies to its CD3 component and the CD28 auxiliary receptor, which was related to Tcell activation [26]. Patterson et al. [27] discovered that Igα non-ITAM tyrosine 204 promoted T-independent  B cell proliferation and differentiation via phosphorylation of JUN. In addition, tumor can produce VEGFA, which increases expression of inhibitory immune checkpoints mediating Tcell exhaustion on intratumoral CD8+ T cells [28]. e previous report identified that VEGFA directly triggers regulatory T cells (Treg) proliferation resulting in tumors escape, which can become a therapeutic goal in the future [29]. TNFSF10 influences T cells function, which cannot only enhance the maximal suppressive function of Treg cells leading to the antitumor effect but also be involved in T cellmediated killing of cancer cells [30,31]. Furthermore, Fang et al. demonstrated that TLR4 appears to induce antitumor T cell response by activating dendritic cells [32]. Similarly, there is research confirming that TLR4 can promote tumorspecific cytotoxic T cell responses [33]. Taken together, from our perspective based on the GO and KEGG analyses and FCM, the regulation of immune-related genes including JUN, VEGFA, TNFSF10, and TLR4 can influence the immune response and be involved in the occurrence and progression of HCC. Drugs targeting immune response of peripheral mononuclear cells may be the potential candidates for therapy of HCC. e PPI network could form a visible framework for a better understanding of the function of the proteome [34]. From the enriched pathways of top 2 modules, we discovered that the interactions among the proteins in HCC are particularly associated with pathways of immune response and cytokine-cytokine receptor interaction. It emphasizes that immune-related gene interaction can regulate the tumor-associated immune surveillance. erefore, we utilized DAVID to uncover the correlation of 15 hub genes. Coincidentally, the 4 genes, JUN, VEGFA, TNFSF10, and TLR4, have the high degree of connectivity, especially JUN and VEGFA. We predicted that JUN can regulate mononuclear cells to release VEGFA, which may promote tumor angiogenesis, which was proved as the reason for tumor initiation and progression. Furthermore, we analyzed the gene set enrichment to make sure of the primary function of the 4 hub genes. We found these 4 hub genes were inextricably linked with the various processes of immune and inflammatory responses like the regulation of T and B cell migration and antigen processing via MHC class I B. Hence, monitoring the immune process of tumor immunity including immune-related genes and cytokines is of great importance for the diagnosis and treatment of HCC.
Moreover, HCC patients always used to suffer from hepatitis infected by hepatitis B virus (HBV) and hepatitis C virus (HCV). Patients who suffer cirrhosis or HCC from chronic virus hepatitis account for 90% [35]. It is consistent with our GO and KEGG analyses that the inflammation is an important process of HCC. JUN, VEGFA, TNFSF10, and TLR4 are connectivity genes involved in the inflammatory response. For instance, there was a report that stated JUN can act as an intermediate in antiviral immunity and TLR4 recognizes bacterial ligands to constrain the ability of antiviral immunity which combine bacterial and viral infections [36,37]. As a result, we hypothesized that JUN, VEGFA, TNFSF10, and TLR4 in PBMC can also regulate virus infection of HCC patients. Detecting these genes can early protect patients from HCC occurrence and development as well as treat them early.
In conclusion, we provide a comprehensive and novel analysis of gene expression profiles to recognize DEGs which may play a core part in the development and prognosis in patients with HCC. Genes implicated with immune and inflammatory response were significantly altered in HCC patients. In order to acquire more precise correlation results, we plan to carry out subsequent authentication experiment later to prove these predictive results. Taken together, we found that JUN, VEGFA, TNFSF10, and TLR4 in PBMC play core roles in the immune response of HCC. Hence, these 4 genes may serve as potential serum biomarkers combining with AFP and targets of immunotherapy. We expect this analysis method will offer accurate and valuable information for future study on the molecular mechanisms of HCC and supply evidences for the detection of new diagnosis biomarkers and therapeutic strategies.

Screen Genes of Differential
Expression. e analysis of DEGs between HCC samples and normal liver tissue samples was performed by using GEO2R (https://www.ncbi. nlm.nih.gov/geo/geo2r/), an online analysis method for GEO database on the basis of R language. We regulated DEGs as differentially expressed with logFC < − 2 (upregulated genes) or logFC > 2 (downregulated genes), according to the criteria described in [38,39]. e adjusted P value < 0.05 was treated as statistically significant in order to reduce the false-positive rate.
ereafter, 344 DEGs were picked up, containing 188 upregulated genes and 156 downregulated genes. We utilized visual hierarchical cluster analysis to get the heatmap and volcano plot of two groups by ImageGP (http://www.ehbio.com/ImageGP/index.php/ Home/Imdex/index.html) after the correlative raw data of TXT files were downloaded.

Gene Ontology and KEGG Pathway Analysis of DEGs.
Gene ontology (GO) analysis is a common framework which can annotate genes and gene products including functions of cellular components, biological pathways, and molecular function [40]. Kyoto Encyclopedia of Genes and Genomes (KEGG) contains a set of genomes and biological pathways related to disease and drugs online database, which essentially is a resource for systematic understanding of biological system and certain high-level genome functional information [41]. e Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.ncifcrf.gov) is an online bioinformatics database. It has covered many biological data and relevant analysis tools, and then provided tools for the biological function annotation information for plenty of genes or proteins [42]. P < 0.05 was considered as the cutoff criterion with significant difference. We could visualize the key biological processes, molecular functions, and cellular components and pathways of DEGs by using DAVID online database. And the scatter plot was performed by ImageGP according to the results of GO and KEGG pathway.

PPI Network and Module Analysis.
Search Tool for the Retrieval of Interacting Genes (STRING) is an online tool for the assessment and integration of the protein-protein interaction (PPI) information, containing physical and functional associations. It covered 9,643,763 proteins from 2031 organisms in STRING version 10.0 [43]. We drew DEGs using STRING to evaluate the interactional associations among them, thereby utilizing the Cytoscape software to build a PPI network. In the meantime, we set maximum number of interactors � 0, confidence score ≥0.4 as the cutoff criterion. Moreover, the Molecular Complex Detection (MCODE) app was utilized to screen modules of the PPI network in Cytoscape in line with degree cutoff � 2, kcore � 2, node score cutoff � 0.2, and max. depth � 100. And we chose the top 15 genes with a high degree of connectivity as hub genes. e pathway analysis of genes in every module was worked out according to DAVID. en, 15 hub genes were also mapped into STRING with maximum number of interactors ≤5 and confidence score ≥0.4. GO and KEGG pathway analyses were also used to explain the potential information.

Comparison of the Hub Genes Expression Level.
GEPIA (http://gepia.cancer-pku.cn/index.html) is a newly developed interactive web server which analyzes the RNAsequencing expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx projects, employing a standard processing pipeline. GEPIA provides customizable functions such as tumor and normal differential expression analysis, profiling according to pathological stage or cancer types, patient survival analysis, similar gene detection, correlation analysis, and dimensionality reduction analysis [44]. In our study, we mainly used the correlation to reveal the relevance out of any two hub genes in HCC and normal people's peripheral blood. Moreover, two suspicious genes that demonstrated a good manner in the scatter diagram were chosen.

Survival Analysis of Hub Genes.
We used GEPIA database to analyze the relapse-free and overall survival information related to the hub genes. e hazard ratio (HR) with 95% confidence intervals and log rank P value were computed and showed on the plot. P < 0.05 was statistically significant.

Gene Set Enrichment Analysis (GSEA).
We divided 20 HCC samples from GSE49515 into two groups (high and low) on the basis of expression level of several key genes, and median expression value was regarded as the cutoff point. In order to explain the potential function of these key genes, GSEA (http://software.broadinstitute.org/gsea/index.jsp) was operated between the two groups. Annotated genes were selected as the reference gene sets. Nominal P value < 0.05, gene size ≥ 100, and│enrichment score (ES)│ > 0.5 were regarded as the cutoff criteria.

Identification of Biomarkers.
Based on the information in the individual MCODE modules, the node with the highest score was selected as the hub gene in GSE49515. According to our analysis of KEGG and PPI, we confirm four key genes (JUN, TLR4, VEGFA, and TNFSF10) which can be the noninvasive biomarkers in early HCC. To verify the vital role of the four key genes in HCC, the mRNA level of these four genes were measured in the PBMC of 20 HCC patients and 20 healthy individuals. Primers used for PCR are listed as follows: JUN-F: ACAGAGCATGACCCTGAACCT; JUN-R: TGTGCCCGTTGCTGGAC; TLR4-F: ATCAAG-GACCAGAGGC; TLR4-R: CACTGAGGACCGACAC; VEGFA-F: GACGGACAGACAGACAGACACC; VEGFA-R: GAAGCGAGAACAGCCCAGA; TNFSF10-F: AGTGG-CATTGCTTGTTT; TNFSF10-R: GAGCTGACGGAGTTGC.

Reidentification of Immune Subversion in HCC.
Meantime, in order to confirm the key role of immune subversion in HCC, we screened peripheral blood mononuclear cells (PBMCs) using flow cytometry (FCM), which were obtain from 20 patients with HCC and 20 healthy individuals (ethical approval for the study was obtained from the Zhongnan Hospital of Wuhan University (no. 2019016)). ereafter, we detected the level of some typical immune cells in PBMCs of HCC and healthy people via the molecules on the immune cell surface using BD FACSCanto II. e samples we detected were from anticoagulated blood with BD Multi TEST IMK Kit (Catalog No. 340503). e total level of T-lymphocytes was screened through CD3 molecules which includes helper lymphocyte T subsets detected by CD3 and CD4 molecules and inhibitory T lymphocyte subsets detected by CD3 and CD8 molecules. B lymphocytes and natural killer (NK) cells were determined by CD19, CD16, and CD56 molecules, respectively. is study was carried out according to legal requirements and supported by the Ethics Committee of the Zhongnan Hospital of Wuhan University. e HCC patients we chose are stage 0 and A of the Barcelona Clinic Liver Cancer (BCLC) stage as well as T1 and T2 of the TNM stage. It means the HCC patients are in the early stage of HCC. Besides, all patients did not take any invasive therapies, like tumor excision, interventional therapy, and hepatic artery embolism. e clinical characteristics of the participants we choose are summarized in Table 6.

Statistical Analysis.
All values were reported as means ± SD. Statistical significance was analyzed by SPSS 19.0 software. Differences were considered significant when P < 0.05.

Data Availability
e figure and table data used to support the findings of this study are included within the article.