Comparative RNA-Sequencing Analysis Reveals High Complexity and Heterogeneity of Transcriptomic and Immune Profiles in Hepatocellular Carcinoma Tumors of Viral (HBV, HCV) and Non-Viral Etiology

Background and Objectives: Hepatocellular carcinoma (HCC), the most common type of primary liver cancer, is the leading cause of cancer-related mortality. It arises and progresses against fibrotic or cirrhotic backgrounds mainly due to infection with hepatitis viruses B (HBV) or C (HCV) or non-viral causes that lead to chronic inflammation and genomic changes. A better understanding of molecular and immune mechanisms in HCC subtypes is needed. Materials and Methods: To identify transcriptional changes in primary HCC tumors with or without hepatitis viral etiology, we analyzed the transcriptomes of 24 patients by next-generation sequencing. Results: We identified common and unique differentially expressed genes for each etiological tumor group and analyzed the expression of SLC, ATP binding cassette, cytochrome 450, cancer testis, and heat shock protein genes. Metascape functional enrichment analysis showed mainly upregulated cell-cycle pathways in HBV and HCV and upregulated cell response to stress in non-viral infection. GeneWalk analysis identified regulator, hub, and moonlighting genes and highlighted CCNB1, ACTN2, BRCA1, IGF1, CDK1, AURKA, AURKB, and TOP2A in the HCV group and HSF1, HSPA1A, HSP90AA1, HSPB1, HSPA5, PTK2, and AURKB in the group without viral infection as hub genes. Immune infiltrate analysis showed that T cell, cytotoxic, and natural killer cell markers were significantly more highly expressed in HCV than in non-viral tumors. Genes associated with monocyte activation had the highest expression levels in HBV, while high expression of genes involved in primary adaptive immune response and complement receptor activity characterized tumors without viral infection. Conclusions: Our comprehensive study underlines the high degree of complexity of immune profiles in the analyzed groups, which adds to the heterogeneous HCC genomic landscape. The biomarkers identified in each HCC group might serve as therapeutic targets.


Introduction
Liver cancer is the sixth most commonly diagnosed cancer and the third most common cause of cancer deaths worldwide [1] (GLOBOCAN 2020 report (https://gco.iarc.fr/today (accessed in May 2021)); Table 1). E October 2013). Ribosomal depletion was performed with 1 µg of total RNA using Ribo-Zero Gold before a heat-fragmentation step that targeted generating libraries with insert sizes ranging from 120 to 200 bp; thus, the remaining RNA was purified, fragmented, and primed for cDNA synthesis. The purified and fragmented RNA was then used to generate cDNA using SuperScript II Reverse Transcriptase (Invitrogen, catalog no. 18064) and random primers. The synthesized cDNA was then transformed into double-stranded DNA incorporating dUTP in place of dTTP to prevent subsequent amplification of the second strand and therefore improve the library's strand specificity. Libraries were subjected to 15 cycles of PCR after 3 adenylation and adaptor ligation steps to generate selectively enriched RNA-Seq libraries suitable for sequencing.
The RNA-seq libraries were evaluated prior to sequencing using an Agilent BioAnalyzer 2100 System and an Agilent DNA 1000 kit (Agilent, part no. 5067-1504) to determine the quality and distribution of DNA fragments. Next, the final library concentration was determined using the Qubit dsDNA HS assay (Thermo Fisher Scientific).
The RNA-seq libraries were standardized to 10 nM, denatured with 0.2 N NaOH, then diluted to 20 pM for downstream sequencing. Sequencing of denatured libraries was carried out in accordance with the manufacturer's standard (Illumina, document no. 15048776 v04, May 2018) using a NextSeq500 platform and NextSeq 500 High Output Kit v2 (150 cycles; up to 400M reads) kits (Illumina, San Diego, CA, USA, catalog no. FC-404-2002).

Analysis of Sequencing Data-Identification of Differentially Expressed Genes (DEGs)
The analysis of sequencing data was performed in collaboration with Illumina by quantifying gene expression against the human genome (version GRCh37 (hg19)) in the The CummeRbund program was used for initial data exploration, analysis, and visualization [23].

Functional Enrichment Analysis
Functional enrichment analysis was performed using METASCAPE [24]. For each gene list, pathway and process enrichment analyses were carried out using: the KEGG Pathway, GO Biological Processes, Reactome Gene Sets, Canonical Pathways, Cell Type Signatures, CORUM, TRRUST, DisGeNET, PaGenBase, Transcription Factor Targets, WikiPathways, PANTHER Pathway, and COVID. All genes in the genome were used as the enrichment background. Terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor >1.5 (the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities. More specifically, p-values were calculated based on cumulative hypergeometric distributions, and q-values were calculated using the Benjamini-Hochberg procedure to account for multiple testing. Kappa scores were used as the similarity metrics when performing hierarchical clustering of the enriched terms, and sub-trees with a similarity of >0.3 were considered as clusters. The most statistically significant term within a cluster was chosen to represent the cluster [24]. Relevant gene functional analysis was performed using GeneWalk, PMID: 33526072 [25]. Genes of interest (hub genes or moonlighting genes) were selected using the following thresholds: global_padj < 0.1; ncon_gene ≥ 50; and ncon_go ≥ 50 [25].

Protein-Protein Interaction (PPI) Network Analysis of DEGs
Interactions of 3 group of DEGs were shown in the STRING online database (http:// string-db.org (accessed on 5/6 December 2021)) [26]. Network nodes represented proteins and edges represented protein-protein associations.

Immune Infiltrate Analysis
The immune infiltrates in HCC samples were investigated using Immunome [27]-a compendium of immune cell markers preferentially expressed in the majority of immune subtypes infiltrating tumors.
The raw count data for tumors and normal samples were analyzed. The quality control revealed that there was no batch effect on the samples. Genes with fewer than 10 counts were removed. Matrices with raw counts were log2-normalized using the TMM method "edgeR". Differential expression analysis was performed with Limma-Voom. Enrichment analyses were performed with Cytoscape Apps [28], ClueGO [29], and CluePedia [30].

Validation of Target mRNA Levels Using Quantitative Real-Time PCR
We carried out validation by quantitative reverse transcriptase (RT) real-time PCR for a selected group of genes (BIRC5 and SLC22A1 in the HCV group; CLEC1B in the HBV group; FGFR4, HSF1, RNF187, HSP90AB1, and HSPB1 in the group without viral infection; and HGF, COLEC10, and CYP17A as genes common to all groups) to validate our next-generation sequencing (NGS) results. These genes were selected based on their up-or downregulation in our HCC samples.
A High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA) was used to synthesize 2000 ng cDNA by incubation as follows: 25 • C for 10 min, 37 • C for 120 min, 85 • C for 5 min, and 4 • C for 5 min. The amplification steps were performed using SYBR Green PCR Master Mix (Applied Biosystems, Thermo Fisher Scientific) with the following thermocycler protocol: 95 • C for 10 min + (95 • C for 15 s; 60 • C for 1 min) for 40 cycles. The ABI PRISM 7300 Detection System (Applied Biosystems, Thermo Fisher Scientific) was used to analyze the relative expression all target genes normalized to b-actin. The expression levels of all target genes were related as fold changes 2 −∆∆Ct . The primers were designed and synthesized by Kaneka Eurogentec S.A. Liège (Supplementary Table S2).

Gene Expression Validation
We verified our gene expression results by qPCR by comparing them with curated databases, such as HCCDB [32], UALCAN [33], CTdatabase [34], and TCGA, and data curated from published research articles.

Differentially Expressed Gene (DEG) Analysis
Although initially the number of differentially expressed genes (DEGs) was very high after setting the filtering threshold of the q-value to 0.05, the number of DEGs significantly dropped.
A first observation was that in each group the gene expression in tumoral tissues was significantly different from that in non-tumoral tissues.
In addition, the total number of DEGs varied between the three groups. Table 2 shows the number of up-and downregulated genes in the HBV, HCV, and non-viral (non-B, non-C) groups.
The highest numbers of upregulated and downregulated genes were identified in the Total HCV group (with 465 upregulated and 226 downregulated) and the lowest number of DEGs was identified in the HBV group (120 upregulated and 102 downregulated). 3.2. Identification of "Common" and "Unique" Genes Further comparative analysis revealed that there were DEGs present in all three groups or just in two out of three, which we called "common"/overlapped (though they presented variable Log2 ratios/fold changes between tumor groups), as follows: 26 upregulated and 17 suppressed genes were found to be common to all three groups (Table 3, common/overlapped HBV, HCV, and non-B, non-C). Figure 1 shows the GO functional enrichment by STRING.
The HBV group had 36 upregulated and 23 downregulated genes that were common to/overlapped with the HCV groups (Supplementary Table S3, data common to HBV and HCV) and another 14 upregulated and 9 downregulated genes common to the non-B, non-C group (Supplementary Table S4, data common to HBV and non-B, non-C).
The HCV group had 59 upregulated and 37 downregulated genes common to/overlapped with the non-B, non-C group (Supplementary Table S5, data common to HCV and non-B, non-C). The HBV, HCV, and non-viral gene lists are shown in Supplementary Tables S6-S8. In addition, the total number of DEGs varied between the three groups. Table 2 shows the number of up-and downregulated genes in the HBV, HCV, and non-viral (non-B, non-C) groups.
The highest numbers of upregulated and downregulated genes were identified in the Total HCV group (with 465 upregulated and 226 downregulated) and the lowest number of DEGs was identified in the HBV group (120 upregulated and 102 downregulated).  Upregulated  Downregulated  Total HBV  120  102  Total HCV  465  226  Total non-B, non-C  441  187  Unique HBV  39  50  Unique HCV  331  145 Unique non-B, non-C 332 121

Identification of "Common" and "Unique" Genes
Further comparative analysis revealed that there were DEGs present in all three groups or just in two out of three, which we called "common"/overlapped (though they presented variable Log2 ratios/fold changes between tumor groups), as follows: 26 upregulated and 17 suppressed genes were found to be common to all three groups ( Table 3, common/overlapped HBV, HCV, and non-B, non-C). Figure 1 shows the GO functional enrichment by STRING.
The HBV group had 36 upregulated and 23 downregulated genes that were common to/overlapped with the HCV groups (Supplementary Table S3, data common to HBV and HCV) and another 14 upregulated and 9 downregulated genes common to the non-B, non-C group (Supplementary Table S4, data common to HBV and non-B, non-C).
The HCV group had 59 upregulated and 37 downregulated genes common to/overlapped with the non-B, non-C group (Supplementary Table S5, data common to HCV and non-B, non-C). The HBV, HCV, and non-viral gene lists are shown in Supplementary Tables S6-S8.  Besides the "common"/overlapped genes, in every group we identified also genes that we called "unique" because they were differentially expressed in only one of the etiological groups (Table 3).  The overlaps between differentially expressed genes (DEGs-upregulated and downregulated) among tumor types are shown in Figure 2A,B as Venn diagrams.
For the visualization of gene expression across the samples from the RNA-Seq results, heatmaps and volcano plots were generated (for each tumor group) ( Figure 3A

Validation by RT-PCR
Validation of the gene expression levels obtained by RNA-seq specific to the HCV group (BIRC5 and SLC22A1), the HBV group (CLEC1B), and the group without infection (FGFR4, HSF1, RNF187, HSP90AB1, and HSPB1) was performed using the real-time PCR technique. In addition, common genes among all groups, such as HGF, COLEC10, and CYP17A, were validated. The results indicated that the gene expression data obtained by RNA-seq were consistent with the expression data determined by real-time PCR ( Figure 4A-D). For the visualization of gene expression across the samples from the RNA-Seq results, heatmaps and volcano plots were generated (for each tumor group) ( Figure 3A-F).   group (BIRC5 and SLC22A1), the HBV group (CLEC1B), and the group without infection (FGFR4, HSF1, RNF187, HSP90AB1, and HSPB1) was performed using the real-time PCR technique. In addition, common genes among all groups, such as HGF, COLEC10, and CYP17A, were validated. The results indicated that the gene expression data obtained by RNA-seq were consistent with the expression data determined by real-time PCR ( Figure  4A-D).

Differential Gene Expression of Solute Carrier Transporters (SLC), ATP Binding Cassette, and Cytochrome Genes
In our comparative study, we identified consistent differential expression (fold change) levels for various SLC, ATP binding cassette, and Cytochrome 450 genes in our three tumor groups (Tables 4 and 5).
The liver is the basic organ of drug and xenobiotic metabolism, transport, and excretion, and it expresses a variety of enzymes and transporters involved in these processes.
The expression of genes encoding these enzymes can be influenced by liver pathologies, such as viral infection, alcoholic liver disease, primary sclerosis, cholangitis, nonalcoholic fatty liver disease, and hepatocellular carcinoma (HCC) [35]. In the processes associated with carcinogenesis, some of the principal classes of macromolecules (carbohydrates, proteins, lipids, and nucleic acids, etc.) could be modified. Consequently, some

Differential Gene Expression of Solute Carrier Transporters (SLC), ATP Binding Cassette, and Cytochrome Genes
In our comparative study, we identified consistent differential expression (fold change) levels for various SLC, ATP binding cassette, and Cytochrome 450 genes in our three tumor groups (Tables 4 and 5).
The liver is the basic organ of drug and xenobiotic metabolism, transport, and excretion, and it expresses a variety of enzymes and transporters involved in these processes.
The expression of genes encoding these enzymes can be influenced by liver pathologies, such as viral infection, alcoholic liver disease, primary sclerosis, cholangitis, non-alcoholic fatty liver disease, and hepatocellular carcinoma (HCC) [35]. In the processes associated with carcinogenesis, some of the principal classes of macromolecules (carbohydrates, proteins, lipids, and nucleic acids, etc.) could be modified. Consequently, some genes involved in transport and metabolism might be differentially expressed in tumoral tissues. For example, in our study, the expression of many solute carriers (SLCs) and cytochromes was modified, with different fold changes presented in the three groups of tumors.
Solute carriers (SLCs) represent a major and important class of cellular transporters. They could become valuable targets in cancer therapeutic strategies (e.g., by blocking or activating them) [36,37].  We present here some of the differentially expressed SLCs and cytochromes identified in our groups of tumors compared with data in the literature.

•
SLC44A5 is an intermediate-affinity choline transporter; high expression of SLC44A5 demonstrates its important role in the development and progression of HCC [38]. • SLC26A6 belongs to an anion transporter family [39]. SLC26A6 expression is an independent prognostic factor for HCC, and its upregulation is correlated with poor prognosis [40]. • SLC38A4 transporter is found predominantly in the liver and transports both cationic and neutral amino acids. Low expression of SLC38A4 is associated with poor prognosis of HCC patients [41]. • SLC22A1 codes for one of the three organic cation transporters, OCT1, an integral transmembrane protein involved in metabolic processes and detoxification. OCT1/SLCCA1 transports a wide range of substances, such as catecholamines, toxins, and anticancer drugs, and is of pharmaceutical interest [42]. In HCC, the expression of OCT1 is significantly reduced and associated with tumor progression and worse patient survival [43]. • CYP17A1 codes for an enzyme involved in the synthesis of steroid hormones, mineralocorticoids, and glucocorticoids [44]. CYP17A1 is significantly increased in human HCC. CYP17A1, as well as CYP19A1, is targeted by inhibitors in cancer treatments. • CYP39A1 studies revealed that total bile acid, total bilirubin, and direct bilirubin were significantly increased in patients with low CYP39A1, and survival analysis of HCC patients indicated that lower CYP39A1 expression was associated with poorer overall survival. The downregulation of CYP39A1 is associated with HCC carcinogenesis, tumor differentiation, and poor overall survival. CYP39A1 may serve as a tumor suppressor gene and a novel biomarker for HCC patients [45]. • CYP2C9 codes for one of the most important drug metabolizing enzymes in humans. Substrates for CYP2C9 include fluoxetine, losartan, phenytoin, tolbutamide, torsemide, S-warfarin, numerous NSAIDs, etc. [46]. In the TCGA database, low expression of CYP2C8, CYP2C9, and CYP2C19 in tumor tissue was associated with short median survival [47]. CYP2C9 could be used as a new biomarker for diagnosis. • CYP2C8 plays an important role in oxidative metabolism; the enzyme metabolizes certain chemicals that contain steroids, arachidonic acids, and retinoids and the anionic parts of some drugs.
Polymorphism in this gene is associated with variable ability to metabolize drugs. CYP2C19 influences metabolism (particularly the detoxification of carcinogens) as a tumor suppressor [48]. CYP2C19 is downregulated in HCC [44] and, consequently, detoxification processes are lower and exposure to carcinogens is higher. As a result, carcinogenesis and proliferation easily occur, leading to aggressive manifestations and poor prognosis in HCC [48].

•
CYP3A4 is an important mono-oxygenase that metabolizes xenobiotics (drugs, toxins, etc.) to eliminate them from the body. The enzyme is predominantly found in the liver but also in the intestines. Dowregulation of CYP3A4 in HCC is associated with poor prognosis. It may be a novel biomarker for HCC [49]. In our study, we identified CYP3A4 as being dowregulated only in the HCV tumor group.
In our study, the expression of cytochromes was variable between groups, and the most affected was the group of HCV-related tumors.
• ABC transporters are mostly exporters; they transport a large variety of molecules using the energy generated by hydrolysis of ATP against their electrochemical gradient. They regulate cellular levels of lipids, ions, xenobiotics, and other small molecules. Studies have revealed that the members of the ABCA subfamily are significantly involved in membrane lipid trafficking (ABCA1, A3, A5, and A9 are detected in almost all tissues) and in cholesterol homeostasis and they have been associated with some inherited diseases [50]. In hepatocellular cancer, ABCA8 and ABCA9 are downregulated, and HCC patients had significantly shorter survival times [51].
Overexpression of many ABC transporters mediates multidrug resistance (MDR) in cancer. In hepatocellular carcinoma, MDR is mediated by ABCB1, ABCB5, ABCC1, ABCC2, and ABCG2 [52]. Beyond the augmentation of the capacity to efflux various therapeutic cytotoxic drugs, the ABC dysregulated genes are being increasingly associated with cancer development and evolution processes (angiogenesis, apoptosis, proliferation, invasion, metastasis, etc.) [50]. ABCB5 has been reported to be overexpressed in HCC and as being associated with chemoresistance, cancer stemness properties, and poor recurrence-free survival [53].
A recent study revealed the upregulation of ABCF1 in drug-selected chemoresistant HCC cells. ABCF1 is a hepatic oncofetal protein that modulates migration, epithelialmesenchymal transition (EMT), and cancer stemness properties and is considered a novel potential therapeutic target for HCC treatment [54].

Differential Expression of Cancer-Testis-Specific Genes
A series of CT genes are upregulated in hepatocellular carcinoma, involved in cellcycle regulation, cancer progression, and signaling pathways and can serve as biomarkers for diagnostics, prognostics, and treatment [55,56].
In our study, we identified a series of dysregulated cancer testis DEGs (Table 6).

•
The Melanoma Antigen Gene (MAGE) family was reported to participate in the progression of multiple cancers in humans, including HCC [57]. • ACTL8 is a member of the sugar kinase/heat shock protein 70/actin superfamily. It is upregulated and contributes to invasion and metastasis in many cancers [58][59][60].
ACTL8 could be involved in epithelial cell differentiation and may be a potential prognostic marker and novel therapeutic target. • ATAD2 (ATPase family AAA domain-containing 2) participates in carcinogenic processes. ATAD2 is overexpressed in various human malignancies, including HCC; it is a potential proliferation marker for liver regeneration and is a poor prognostic marker for hepatocellular carcinoma after curative resection [61,62].

Differential Expression of Heat Shock Proteins and Heat Shock Factors
The heat shock proteins (HSPs) or stress proteins are cellular constitutive, ubiquitous, highly evolutionarily conserved molecules that have cytoprotective properties; they are the main basic elements of the cellular proteoprotection system. HSPs have multiple functions in physiological conditions: chaperoning functions, disaggregation and possibly even refolding of damaged proteins, and, more important, protection of newly synthesized proteins and help with their folding (in an ATP-dependent manner) into functional forms. HSP expression might be induced by many stress factors, e.g., thermic stress (heat shock), chem-ical stress (heavy metals), oxidative stress (free radicals), denatured proteins, antibiotics, immunosuppressive drugs, hypoxia, nutritional inadequacy, and pathological conditions (inflammation, fever, viral and bacterial infections, carcinogenesis) [63][64][65]. Deregulation of stress gene expression is associated with various human diseases, including malignancies. Heat shock proteins (HSPs) are found to be overexpressed in tumor cells, where they protect oncogenic proteins. Stress induction of HSPs plays a crucial role in tumorigenesis, metastasis, and therapeutic resistance [66].
Our analysis identified significant dysregulation of stress proteins and heat shock factors, thus confirming perturbed metabolic homeostasis in HCC (Table 7).
• HSF4 (Heat Shock Transcription Factor 4) HSF4 is a member of the heat shock transcription factor family and is expressed in human tissues. Dysregulation of HSF4 expression might induce carcinogenesis. HSF4 was found to be upregulated in HCC tissues and, more important, elevated in primary HCC tissues derived from recurrent patients; consequently, HSF4 was considered an independent poor-prognosis predictor after resection [67,68]. Our analysis identified HSF4 as being upregulated only in the HCV tumor group.

Functional Enrichment Analysis of DEGs
Here, we present the most relevant GO processes and pathways enriched with DEGs as found in our study for the HBV (Figure 5A,B), HCV ( Figure 5C,D), and non-B, non-C groups respectively ( Figure 5E,F).

Identification of Regulator Genes, Hub Genes, and Moonlighting Genes
Regulator genes (or regulatory genes) are genes that regulate the expression of one or more structural genes. Their function is to ensure that gene products, such as enzymes, structural proteins, and RNA molecules, are synthesized when they are needed and in the proper amounts. Regulator genes also allows cells to react quickly to changes in their environments [69].
Moonlighting proteins comprise a class of multifunctional proteins which singly perform multiple physiologically relevant biochemical or biophysical functions that are not due to gene fusion, multiple RNA splice variants, or pleiotropic effects (e.g., soluble enzymes that also bind to DNA or RNA to regulate translation or transcription) [70,71].

Functional Enrichment Analysis of DEGs
Here, we present the most relevant GO processes and pathways enriched with DEGs as found in our study for the HBV ( Figure 5A Numerous studies have demonstrated that individual proteins can moonlight, meaning that they can have multiple functions based on their cellular or developmental contexts.
Moonlighting may be particularly relevant in the context of human disease, especially in cancer [72].
In our study, the DEG analysis using GeneWalk showed 145 regulator genes in the HCV-related tumor group and 1 moonlighting gene (ECT2). In the non-viral-infected group of tumors, 106 regulator genes (graphic representations in Figure 6A,B) and 5 moonlighting genes (ACTN4, FLNA, NOTCH1, TOP2A, and PDGFRA) were detected ( Figure 7A,B). The regulator genes were identified as those with a wide connectivity to other input genes and high fractions of relevant GO annotations. The list of regulator genes is available in Supplementary Table S9 (HCV group) and Supplementary Table S10 (non-viral group). Moonlighting genes were identified as those with many GO annotations of which only a small fraction are relevant [25]. No significant regulatory or moonlighting genes were identified in our HBV-related tumor group.
From the regulator gene list, we further selected the most significant hub genes using the following thresholds: global_padj < 0.1; ncon_gene ≥ 50; and ncon_go ≥ 50. The hub genes and their connectivity degrees are presented in Table 8 and Supplementary Figure S4, respectively, for the non-viral group and in Table 9 and Supplementary Figure S3, respectively, for the HCV group.
We further describe the functional role and involvement of some identified hub genes in normal and pathological conditions.

Non-Viral Group HUB Genes and Proteins
• HSF1 (Heat Shock Transcription Factor 1). In vertebrates, the prototype of heat shock transcription factor is HSF1, which mediates the induction of heat shock gene expression in response to environmental stress [73].
As a mitotic regulator, HSF1 is a major contributor to cancer morbidity. It allows a series of cell-level tumorigenic processes (deregulation of cell-cycle progression, increased cell survival, etc.) and modulates tumor-level tumorigenic features (invasion, angiogenesis, and metastasis) [74]. HSF1 participates in the initiation, development, and progression of various cancers, including hepatocellular carcinoma. HSF1 exhibits high expression in HCC and in other malignancies [75][76][77][78].
• HSPB1/HSP27 is a stress-inducible chaperone which belongs to the small heat shock protein family [79]. HSPB1 has multiple functions and regulates many cellular processes, such as cytoskeleton organization, maintenance of cellular proteostasis, inhibition of apoptosis, modulation of autophagy induction of resistance to anticancer drugs, etc. [80][81][82]. Numerous studies have revealed that HSPB1 promotes tumorigenesis [66,83,84] and is dysregulated in different malignancies. HSPB1 is upregulated in HCC, and it was identified as a hub gene [85]. The overexpression of HSPB1 was associated with a worse prognosis in HCC patients and it was considered a possible target of immunotherapy in HCC [86]. • HSP90AA1: The heat shock protein 90 (HSP90) family perform a large number of cellular regulatory functions in normal and pathological processes. In vertebrates, the two major paralog isoforms are HSP90AA1 and HSPAB1 [87].
moonlighting genes (ACTN4, FLNA, NOTCH1, TOP2A, and PDGFRA) were detected (Figure 7A,B). The regulator genes were identified as those with a wide connectivity to other input genes and high fractions of relevant GO annotations. The list of regulator genes is available in Supplementary Table S9 (HCV group) and Supplementary Table S10 (nonviral group). Moonlighting genes were identified as those with many GO annotations of which only a small fraction are relevant [25]. No significant regulatory or moonlighting genes were identified in our HBV-related tumor group.  Hsp90 is an essential element for malignant transformation and progression as a cancer supporter that assists and interacts with oncogenic proteins [88]. HSP90 acts as an important regulator of autophagy that leads to inhibited apoptosis and increased drug resistance [89]. HSP90AA1 suppression results in increased sensitivity to chemotherapy [90].
Hepatocellular HSP90 is positively involved in HCC development by increasing liver cancer cell invasion, inhibiting cancer stem cells, apoptosis, etc. [90]. It is considered a potential biomarker for detection/screening, prognostics, and supervision of human hepatocarcinogenesis [91] and a valuable target in cancer therapy [92]. From the regulator gene list, we further selected the most significant hub genes using the following thresholds: global_padj < 0.1; ncon_gene ≥ 50; and ncon_go ≥ 50. The hub genes and their connectivity degrees are presented in Table 8 and Supplementary Figure  S4, respectively, for the non-viral group and in Table 9 and Supplementary Figure S3, respectively, for the HCV group.      • HSPA1A gene codes for molecular chaperons proteins, belonging to the HSP70 Heat Shock Protein Family A. HSPA1A is a major stress-induced member, having crucial roles in protein homeostasis and cell survival [93].
In oncogenesis, HSPs play an essential, facilitating role through the accumulation of overexpressed and mutated oncogenes through their cytoprotective functions (inhibition of apoptosis, as well as HSP27) [81,94] and have multiple implications for the hallmarks of cancer [95].
• HSPA5 (GRP78, BiP) is a chaperone protein constitutively expressed in the endoplasmatic reticulum (ER); GRP78 maintains normal ER functions and is the principal regulator of cellular response to ER stress [101,102]. A series of studies have demonstrated that GRP78/HSPA5 is anti-apoptotic and has a critical cytoprotective role in oncogenesis (protects tumor cells from ER stress) [103].
It was demonstrated that GRP78 is a novel obligatory component of autophagy in mammalian cells [103][104][105]. GRP78 is involved in tumor proliferation, survival, tumor angiogenesis, metastasis, and drug resistance, and overexpression of GRP78 was observed in the progression of many human cancers, including hepatocellular carcinoma [105][106][107].
• ACTB (Beta-actin) is a highly conserved cytoskeleton structural protein generally upregulated and involved in the development and metastasis of various cancers, including HCC [108,109].
is an important member of the glucose metabolism enzyme family. Glucose metabolism dysfunction is one of the most important characteristics of cancers. High expression of ALDOA is associated with the initiation and progression of many cancers. ALDOA contributes to moonlighting functions; under hypoxia, ALDO A regulates cell proliferation, invasion, and apoptosis, being an essential driver in HCC [110,111]. • GAPDH (Glyceraldehyde-3-phosphate dehydrogenase) is an essential regulator of glycolysis overexpressed in numerous cancers, including HCC, and enabling tumor progression. GAPDH is functionally active in the nucleus, cytoplasm, and plasma membrane and also carries out numerous, non-glycolytic "moonlighting" functions. Glycolytic enzymes have gained increasing attention as potential anticancer therapeutic targets. • CTTN (Cortactin) is an important actin-binding and assembly protein involved in cytoskeletal regulation. It is found at sites of dynamic actin assembly, in cellular protrusions, such as invadopodia, and is associated with cell motility and invasion. CTTN enhances cell migration, invasion, and tumor cell metastasis and is overexpressed in many cancers, including HCC [112,113].

HCV Group HUB Genes and Proteins
• ACTN2-Alpha actinin is an actin-binding cytoskeletal protein. The alpha actinin isoform, which is concentrated in the cytoplasm, is thought to be involved in metastatic processes.
Studies have revealed that ACTN2 overexpression in HCC stimulates invasion abilities by enhancing cellular motility, demonstrating a pro-metastatic role in tumorigenesis [114]. ACTN2 was also cited as an HCC hub gene in other studies [115].
• ANXA 2-Annexin A2 belongs to a protein family (annexins) whose members bind anionic phospholipids in a calcium-dependent manner and have the ability to aggregate membranes.
ANXA 2 is overexpressed in many human cancers, including hepatocellular carcinoma (HCC), and has multiple regulatory roles and is correlated with proliferation, cell migration, adhesion, angiogenesis, apoptosis, etc. [116].
Upregulated ANXA2 in HCC plays an important role in tumor immune escape and is proposed as a target in cancer treatment [117,118].
• AURKA (Aurora Kinase A) is a serine/threonine kinase that plays essential roles in regulating cell division during mitosis. Abnormal activity of AURKA promotes tumorigenic progression [119] and is highly expressed in various cancers, including HCC [120]. It might be a reliable predictor of early-stage HCC, a crucial biomarker for HCC development, and a reliable target for cancer therapy [121][122][123]. • AURKB (Aurora Kinase B) is a serine/threonine fundamental kinase (as is AURKA) involved in the regulation of cell mitosis, especially in chromosomal segregation [124]. The dysregulation of aurora kinase genes has been reported in many cancers. The expression of AURKB was found to be higher in HCC than in a control and was consistently correlated with patient tumor stage [125]. • BRCA1 (BRCA1 DNA Repair Associated) encodes a nuclear phosphoprotein that plays an important role in the correct repair of damaged DNA and maintaining genomic stability. BRCA1 is overexpressed in many type of cancers, including HCC, where its expression correlates with immune cell infiltration [125,126]. • CCNB1 (Cyclin B1) belongs to the cyclin family. Eukaryotic cell-cycle progression is regulated by cyclin-dependent kinases (Cdks) and their regulatory cyclin subunits. Cyclin/Cdk complexes activate transcription, enable DNA replication, and catalyze mitosis [127,128]. Overexpression of CCNB1 can promote proliferation in human HCC cells and was identified as a hub gene in HCC in others studies as well [129]. • CDK1 (Cyclin Dependent Kinase 1) is a member of the serine/threonine protein kinase family and has a crucial role in cell proliferation initiating mitosis [127,128]. Overexpression of CDK1 has been observed in different type cancers, including hepatocellular carcinoma [130,131], where it is correlated with poor OS. Moreover, expression levels of CDK1, CCNB1, and CCNB2 were positively correlated with infiltrating levels of CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells in HCC [132,133].
Other studies also identified CDK1 and CCNB1 as Hub genes for HCC [129]. • CYP3A4 (Cytochrome P450 Family 3 Subfamily A Member 4) codes for an important mono-oxygenase, metabolizing xenobiotics (drugs, toxins, etc.) to eliminate them from the body [134]. The enzyme is predominantly found in the liver but also in the intestines. Dowregulation of Cyp3A4 in HCC was associated with poor prognosis [135]. • CYP1A2 (Cytochrome P450 family 1 Subfamily A Member 2) is the major hepatic isoform of the human CYP1A subfamily. It is involved in the clearance mechanisms for important drugs (tizanidine, theophylline, clozapine, caffeine, etc.) and participates in the biotransformation processes of different procarcinogens. CYP1A2 is markedly decreased in primary HCC tumors and is an independent predictor for post-surgical recurrence in early-stage HCC patients [136,137]. • EGF (Epidermal Growth Factor) is a growth factor secreted by tumors and inflammatory cells in the tumor microenvironment. EGF binds to a transmembrane glycoprotein, its receptor EGFR (epidermal growth factor receptor), and activates/triggers regulatory signal transduction pathways of proliferation, differentiation, survival, and migration.
Overexpression of EGF was reported in many human cancers, including HCC. EGF is highly expressed in HCC and facilitates DNA synthesis, regeneration, tumor growth and progression, and promotes metastasis [138,139].
• TERT (Telomerase reverse transcriptase) is the catalytic subunit of telomerase. In early stages of cancer, because of the increased cell proliferation, telomeres are shortened, but with tumor progression telomerase is reactivated and the capacity for infinite cell division (immortalization) is gained [140]. TERT upregulation is a critical event in hepatocarcinogenesis. It has been shown that TERT expression increases in hepatocyte cultures after overexpression of HCV core protein as compared to normal human liver and uninfected cells [141]. The HCV core protein is a transcriptional activator of a number of host genes [142] and it has been suggested that it interferes with telomerase expression and might be essential for malignancy [141]. • TOP2A (DNA Topoisomerase II Alpha) is a gene that encodes a nuclear enzyme which catalyzes the transient breaking and rejoining of two strands of duplex DNA (which allows the strands to pass through one another). It is involved in processes that occur during DNA transcription and replication, such as relief of torsional stress, chromosome condensation, and chromatid separation [143]. TOP2A is involved in many cancers, being a prognostic biomarker and potential therapeutic target for bladder cancer, lung adenocarcinoma, prostate cancer, colon cancer, breast cancer, and HCC [144,145].

Analysis of Tumor Immune Infiltrate
Based on the transcriptomic differences shown between etiological tumor groups analyzed in our study, we further examined the tumor immune infiltrate transcriptomic profiles.
In general, in HCC the immune landscape of tumoral tissue is significantly altered compared to the normal one; consequently, the evaluation of immune infiltration patterns contributes to the establishment of new HCC immunotherapy strategies in personalized medicine.
The immune infiltrate in HCC samples was investigated using Immunome [27]-a compendium of immune cell markers preferentially expressed in the majority of immune subtypes infiltrating tumors.
We first compared the tumoral tissues with the normal samples, regardless of HBV and HVC infection, and obtained genes significantly differentially expressed (adjusted p-value < 0.005) with a fold change greater than 2.5. Among these genes, markers preferentially expressed in Th2 cells and eosinophils had significantly higher expression in tumors compared to normal samples (Supplementary Figure S5A). In contrast, in normal tissue markers of B cells, Th1, TFH, iDC, NK cells, mast cells, and macrophages had significantly higher expression than in tumors. Similar results were observed in the TCGA cohort (n = 421 patients; Supplementary Figure S5B).
In the next step, we analyzed the Immunomes from tumors with HBV, HCV, and tumors without viral infection. Immune gene data were extracted and clustered ( Figure 8A, Supplementary Table S11) and the biological roles of genes with the highest expression in HBV (Cluster 3), HCV (Cluster 2), or tumors without viral infection (Cluster 1) were investigated with ClueGO [29] and CluePedia [30]. An over-representation of T-cell-, cytotoxic-, and natural-killer-cell-related GO terms was observed for Cluster 2 genes ( Figure 8B, Supplementary Table S11). Many of these genes are known to be involved in protein-protein interactions leading to the activation or inhibition of expression or to immune cell activation, while others are chemokine-receptor binding pairs (Supplementary Figure S5C). Cluster 3 genes were associated with GO terms involved in primary adaptive immune response and complement receptor activity, while Cluster 1 genes were associated with monocyte activation involved in immune response and other metabolism-related terms.  Fusion was applied. Network shows GO terms after multiple testing correction. The size of the nodes shows the significance of the terms. Nodes are colored based on the proportions of associated genes from Cluster 1, Cluster 2, or Cluster 3. Equal proportions of genes from the three clusters are shown in gray. Bar charts showing the expression levels of (C) T cells (D) CD8 T cells (E) cytotoxic cells, (F) immune checkpoint markers, and (G) cytokine markers in HBV (red), HCV (blue), and non-B, non-C (green) tumors. Differential expression analysis was performed with Limma-Voom. Significance levels are shown as: * p < 0.05, ** p < 0.01, *** p < 0.005.
As previously reported [146], a similar immune profile was observed for HBV and HCV tumors ( Figure 8C). Many immune genes had significantly higher expression in these tumors than in tumors without viral infection, as was illustrated for markers of T cells, CD8 T cells, and cytotoxic cells ( Figure 8C-E).
The immune checkpoint inhibitors PDCD1 and CTLA4 were also part of Cluster 2 and had significantly higher expression in HCV compared to non-B, non-C tumors ( Figure 8F). A similar trend was observed for HBV tumors.
The expression of cytokines, such as CXCL9, CXCL10, CXCL11 and CXCL13, that attract specific immune cells at the tumor site was significantly higher in HBV and HCV tumors than in non-B, non-C tumors ( Figure 8G).

Discussion
Our study analyzed three original NGS whole-transcriptome datasets and revealed consistent differential gene expression between non-tumoral and tumoral tissues, including 222 DEGs (120 upregulated and 102 downregulated) in HBV-related tumors, 691 DEGs (465 upregulated and 226 downregulated) in HCV-related tumors, and 628 DEGs (441 upregulated and 187 downregulated) in non-viral-infected tumors. In the HBV group, a smaller number of DEGs were consistently identified. We identified common (overlapped) DEGs (present in all three etiological groups or in two of three) and unique DEGs (present only in one of the three groups) in all three analyzed groups.
Further analysis showed variable fold change values for common DEGs between the tumor groups.
In the HCV-related group, we identified a higher number of dysregulated DEGs as SLC (solute carrier), cytochrome p450, cancer testis, oncogenes, tumor suppressor genes, etc.
SLC and cytochrome (CYP) genes are involved in drug absorption, distribution, metabolism, and excretion (ADME). Their correlated actions control liver drug metabolism and clearance and, consequently, the efficacy of therapy. In pathological conditions (such as viral infection, alcohol abuse, HCC, etc.), the expression and activity of these genes are modified.
The solute carriers (SLCs) are important cellular carriers, having consistent roles in different metabolic processes and tumorigenesis, and may become cellular targets for new therapeutic agents [36,37]. The analysis of SLC DEGs identified a series of dysregulated genes in the three tumor groups.
Our results are in accordance with data reported in the literature concerning the upregulation/overexpression of SLC genes in HCC, including SLC44A5 and SLC26A6 [38,40], and the downregulated expression of SLC38A4 and SLC22A1 [41,43]. A recent study evidenced SLC26A6 as a novel oncogene in HCC [147].
Cancer testis (CT) genes are restrictedly expressed in normal tissues except for the testis and are aberrantly expressed in tumor tissues and often trigger humoral or cellular antitumor responses in patients with cancer. Expressed testis-derived antigens might be immunogenic because they have never been encountered by the immune system. Consequently, cancer testis (CT) genes have been indicated as a group of potential targets for TAA-specific immunotherapy [148].
A series of studies have demonstrated that CT antigens are regulators of cancer hallmarks, e.g., sustaining proliferative signaling, resisting cell death, and evading growth suppressors [149].
The MAGE and GAGE families are groups of tumor-specific antigens that can be used as molecular markers for early diagnosis but are also appropriate targets for vaccine-based cancer immunotherapy in human HCC [56,57].
Recent studies have revealed the involvement of members of the MAGE family in stress response pathways. CT MAGE-B2 is normally expressed in the testis but is highly expressed in tumors. It was reported that MAGE-B2 enhances cellular stress threshold by suppressing stress granule (SG) assembly and promoting cellular stress tolerance. This allows cancer cells to continue to grow even in stressful conditions [150].
Our analysis identified high expression of MAGEB2, MAGEC3, and MAGEB17 only in the non-viral-infected tumor group, in which we observed the overexpression of genes involved in stress response. The overexpression of other CT genes, including MAGEA1, PAGE4, PAGE5, GAGE2A, BAGE, BAGE3, BAGE4, BAGE5, ACTL8, TPTE, HSPB9, ATAD and DSCR8, was shown in the three groups. Our results for the CT gene expression analysis are in accordance with data reported in the literature.
The next step in the analysis releveled common pathways and GO terms upregulated in the HCV and HBV groups, such as cell-cycle, mitotic, and PLK1 pathways. These results support the idea that hepatitis viruses deregulate the cell cycle in infected cells to promote an environment that can sustain viral replication. The persistence of viral infection and ineffective T cell responses maintain an inflammatory state in the liver that probably culminates in the onset of fibrosis, cirrhosis, and HCC.
In the non-viral-infected group, functional enrichment analysis revealed a significant upregulation of cellular stress response.
In HCV-related tumors, the downregulated genes were mainly enriched in "carcinogenesis-DNA adducts KEEG pathway" but also in "monocarboxylic acid metabolism" and "steroid metabolic processes" (demonstrating a reprogramed metabolism).
Considering that the liver is the site of major metabolic processes, it is not surprising that we identified altered metabolic pathways in our samples, as these may be highjacked by tumors in order to help them survive and proliferate, especially against a background of viral infection. It has been shown that HCV infection may have a direct effect on lipid metabolism and that cholesterol metabolism is decreased in HCV-infected HCC through the downregulation of genes involved in cholesterol synthesis, absorption and transport, and bile acid synthesis [151].
Immunome analysis revealed the immune expression patterns of the three patient groups that were analyzed. For example, the HCV tumors showed significant upregulation of T cell genes, as well as of genes with cytotoxic properties, markers of CD8 T cells, and cytotoxic cells, compared with tumors without viral infection. The HBV tumors showed a similar trend. High expression of such immune markers was previously reported to be associated with prolonged survival in many cancer types, such as colorectal cancer [27,[152][153][154]. A similar expression pattern was observed for cytokines that modulate intratumoral immune infiltrate and also for immune checkpoint CTLA4 and PD-1 (PDCD1).
CD96, among T cell genes, is considered a novel immune checkpoint receptor target [155], and higher expression was associated with a poorer clinical outcome [156], while low LCK expression has been identified as a potential prognostic biomarker for immunotherapy in HCC [157]; however, the study did not segregate patients based on viral infection. In our study, we found LCK to be over-expressed in the HCV group but to have a lower expression in the non-B, non-C (non-viral) group.
We would like to acknowledge one limitation of the study, namely, the small number of patients. However, we wanted to give a close representation of real-world data by including the most common viral infections that represent risk factors for HCC and further compare these with data from non-viral-infected patients.
Further validation with bigger cohorts is needed, but we were able to determine both differences as well as commonalities between the three different HCC groups analyzed in our study. The HBV group was characterized by the smallest number of genes. Metascape analyses showed the upregulation of the cell-cycle pathway, similar to the HCV group.
Additionally, the HCV group showed overexpression of a series of immune cells (including CTLA-4 and PDCD1) and downregulation of Cyp and SLC. The third group of patients (non-B, non-C/non-viral) showed the presence of upregulated Heat Shock Factor 1 and a series of HSPs as hub genes (upregulation of response to stress pathway), possibly pointing to HSF1/HSPs as modulators of apoptosis/autophagy in this group.

Conclusions
Our comparative RNA-seq analysis of liver cancer tumors revealed heterogeneity among HBV, HBV, and non-B, non-C (non-viral) tumors with regard to transcriptomic and immune profiles and contributes to a better understanding of the pathogenesis and progression mechanism of HCC.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/medicina58121803/s1, Figure S1: Histology (hematoxylin-eosin staining) for the HCC patients (non-viral, HBV and HCV etiology), Figure S2: Tuxedo pipeline, Figure S3: HCV HUB and moonlighting genes enrichment, Figure S4: nonBnonC HUB and moonlighting genes enrichment, Figure S5: Immunome analysis in tumor tissues versus non-tumoral tissues and in TCGA cohort; Table S1: Patients' features,  Table S6: List of up and down-regulated genes in HBV related HCC, Table S7: List of up and down-regulated genes in HCV related HCC, Table S8: List of up and down-regulated genes in non-viral HCC, Table S9: List of regulatory genes in HCV related group, Table S10: List of regulatory genes in non-viral group, Table S11: Clustering of immune gene data by ClueGO.