{'fraction_non_alphanumeric': 0.036069123951946856, 'fraction_numerical': 0.0207635240924189, 'mean_word_length': 5.277395058326444, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characterization of heterogeneous redox responses in hepatocellular carcinoma patients using network analysis

Background Redox metabolism is often considered a potential target for cancer treatment, but a systematic examination of redox responses in hepatocellular carcinoma (HCC) is missing. Methods Here, we employed systems biology and biological network analyses to reveal key roles of genes associated with redox metabolism in HCC by integrating multi-omics data. Findings We found that several redox genes, including 25 novel potential prognostic genes, are significantly co-expressed with liver-specific genes and genes associated with immunity and inflammation. Based on an integrative analysis, we found that HCC tumors display antagonistic behaviors in redox responses. The two HCC groups are associated with altered fatty acid, amino acid, drug and hormone metabolism, differentiation, proliferation, and NADPH-independent vs -dependent antioxidant defenses. Redox behavior varies with known tumor subtypes and progression, affecting patient survival. These antagonistic responses are also displayed at the protein and metabolite level and were validated in several independent cohorts. We finally showed the differential redox behavior using mice transcriptomics in HCC and noncancerous tissues and associated with hypoxic features of the two redox gene groups. Interpretation Our integrative approaches highlighted mechanistic differences among tumors and allowed the identification of a survival signature and several potential therapeutic targets for the treatment of HCC.


Introduction
Hepatocellular carcinoma (HCC) is a prevalent form of primary liver cancer and represents the second leading cause of worldwide cancer mortality [1,2]. Its prevalence is predicted to increase in the next two decades, and the prognosis is usually poor, with most diagnosed cases resulting in death [3]. Due to the high inter-and intra-tumor heterogeneity, it is challenging to establish a robust HCC classification and develop targeted therapies [4]. Hence, there is an urgent need to stratify patients to reveal the underlying molecular mechanisms of the disease that may be used in the development of targeted and effective treatment strategies.
Reactive oxygen species (ROS) such as hydrogen peroxide and superoxide are central players in redox metabolism and have important cellular functions in which they control signaling, cell cycle progression, proliferation, inflammation and immune responses [5]. However, the deregulation of redox and ROS metabolism may be toxic for cells, disrupting pathways and promoting mutagenesis and apoptosis [6,7], which are important features in cancer [8]. Due to their pivotal cellular functions, recent efforts have sought to devise ROS-based cancer treatments [9]. However, effective redox-based HCC therapeutics are currently lacking due to our incomplete understanding of the underlying redox mechanisms of the disease in different patient groups, and the extensive crosstalk between redox and central metabolism [5]. Systems level analyses have previously been employed to unravel cancer heterogeneity, to identify novel prognostic genes, therapeutic targets, and drug candidates and to enable patient stratification [10][11][12][13][14], including for HCC [15][16][17].
Here we used systems biology approaches (Fig. S1A) to reveal the heterogeneous redox responses in HCC combining transcriptomics (RNA-seq and an independent set of microarray data), proteomics and metabolomics data. We first identified several patterns of coexpression among genes involved in redox metabolism, liver-specific genes, and genes associated with immunity and inflammation that have prognostic properties based on integrative network analyses. Second, we identified several associations between the genes involved in redox and central metabolism. Next, we identified specific functional patterns between important ROS scavengers and how they are coupled to the central metabolism of the tumors. Tumor stratification based on the expression of redox genes and analysis of clinical data allowed the discovery of potential markers for effective treatment strategies of HCC patients and novel survival signatures based on the expression of key redox genes. Finally, we validated our results by performing experiments using HCC mouse models and analyzing independent human HCC tumors.

Identification of redox genes
We identified 132 genes directly involved in redox metabolism (Dataset 1). This list included genes encoding for proteins 1. with antioxidant and/or ROS-dependent signaling activities (e.g., superoxide dismutases, catalase, peroxidases, peroxiredoxins, thioredoxins); 2. involved in reactions coupled with antioxidant defense or ROS-mediated signaling transduction (e.g., thioredoxin and glutathione reductases); 3. involved in producing compounds with relevant redox roles (e.g., glutathione synthetase, glutamate-cysteine ligase); 4. enzymes involved in reducing equivalent biosynthesis from pentose phosphate pathway (e.g., glucose-and hexose-6-phosphate dehydrogenases), folate metabolism (e.g., tetrahydrofolate synthase, bifunctional methylenetetrahydrofolate dehydrogenase), malate metabolism (malic enzymes), and other processes (e.g., NAD(P)H dehydrogenase, NAD(P) transhydrogenase, aldehyde and glutamate dehydrogenases); 5. transcription factors (e.g., hypoxia-inducible factor 1-alpha, nuclear factor erythroid 2-related factor 2) that coordinate responses to oxidative stress. This list was used to perform the initial hierarchical clustering, and co-expression analysis with all the genes in the genome. Due to the high co-expression between redox genes and those involved in central metabolic pathways, we expanded the list to 174 redox genes (Dataset 1), including those indirectly involved in redox metabolism with known importance for cancer and/or that control or are regulated by the above pathways (e.g., glycolytic enzymes).

Data selection and gene expression analysis
RNA-seq gene expression data for primary liver hepatocellular carcinoma in 360 subjects and 50 noncancerous samples were retrieved from NCI's Genome Data Commons [18] as raw counts or Fragments Per Kilobase of transcript per Million mapped reads (FPKM). Only one sample was considered per patient, and samples belonging to control or normal types were excluded. Patient survival, expressed in Living Days, was identified by considering age at diagnosis for deceased patients, or days to last follow up for those alive. Unless otherwise stated, all analyses were performed based on FPKM data, and genes displaying a median (FPKM) b 1 across patients were considered as not detected.
In the cancer progression analyses, we considered grade and not stage classification of tumors because the former is a histological classification of the tumor and the latter takes into account other features that are not related to the primary tumor. In total, 55 Grade 1, 175 Grade 2, 118 Grade 3 and 12 Grade 4 patients were considered and divided according to early (G1), intermediate (G2) and late (G3 + G4) cancer.
For a smaller subset of the tumor samples (186), more detailed clinical information was available [13] and used to assess relationships with redox behavior. Liver-specific and enriched genes, and group-enriched genes, were attained from ref. 25 . Genes involved in the immune system (adaptive and innate immunity, and cytokine signaling in immune system) and inflammation (CLEC7A/inflammasome pathway, DEx/H-box helicases activate type I IFN and inflammatory cytokines production, inflammasomes) were identified using REACTOME [64].

Hierarchical clustering, differential expression and co-expression analysis
Hierarchical clustering was performed in R considering Ward.D2 and the Euclidean distance. Whenever clustering was performed based on gene expression matrices, genes were Z-normalized before clustering. Differential expression analyses were performed based on raw counts through R using the DESeq2 package [65] using default methods and considering genes with absolute log fold changes higher than 1 by setting the results function with arguments altHypothesis = "greaterAbs" and lfcThreshold = 1 and the DEseq function with argument betaPrior = FALSE. An FDR of 5% was set during the differential expression. Although the filtering based on the log fold change and FDR may be performed a posteriori, specifying these arguments upon performing the differential expression analysis is recommended for independent gene filtering and significance of the Wald testing [65].

Research in context
Evidence before this study Imbalances in redox metabolism are prevalent features in cancer, influencing tumorigenesis and proliferation. Because of this, several recent efforts have sought to target redox and antioxidant metabolism in cancer. However, a systematic examination and characterization of redox metabolism has never been performed in hepatocellular carcinoma, hindering our understanding of redox behavior and its relationship with metabolism, signaling, or patient clinical data.

Added value of this study
We performed a systems level multi-omics analysis of redox metabolism and identified two major opposing clusters of genes and their associations with central metabolism, differential metabolic responses, immunity and inflammation. Tumors stratified according to redox gene clusters display substantial progression and survival differences and subtype-specific potential therapeutic targets. Experiments with mouse models validated the antagonistic behavior and highlighted differential hypoxic associations for redox genes. We finally identified and validated a survival signature derived based on genes associated with redox metabolism.

Implications of all the available evidence
Our observations highlighted, for the first time, the antagonistic relationships between key genes in redox metabolism. Furthermore, while redox genes are heterogeneously expressed, HCC tumors may be stratified into two groups with distinct phenotypes, potential therapeutic targets, and associated patient survival.
Differential expression analysis was performed by partitioning subjects according to survival (in days), or between tumors displaying high vs low expression levels of redox genes. Because gene expression expressed in raw counts could not be compared across different samples, we sorted subjects according to the FPKM expression of the target gene and then identified the two groups that displayed highest and lowest gene expression, each with 50 individuals. Genes displaying median FPKM b 1 across individuals were ignored.
Spearman rank correlations were computed for gene vs gene expression across tumors, and Q values were used to find co-expressed genes. Genes with expression that was not statistically correlated (Q N 0.01) with at least one gene in the genome were disregarded.

Gene set enrichment analysis and reporter metabolite analysis
GSEA was performed through PIANO [66] using the whole-genome log 2 -fold changes and adjusted P values attained from DESeq2 as gene-level statistics, with geneSetStat = reporter, and nPerm = 1000. Gene Ontology biological processes were downloaded from MSigDB [67]. Gene Ontology processes were considered as enriched with an FDR of 5% and with clear direction (i.e., non-mixed). Additionally, we ignored Gene Ontology processes related to tissues/organs of different embryonic origin (e.g., brain, bone, hair). Reporter metabolites [28] were identified through PIANO using the same input as above through the function runGSA with default arguments and based on fold changes. Gene set collections were determined by assigning ensemble ids to each reaction's metabolites using the iHCC model [15]. Only those metabolites with clear direction were considered. The concept of reporter metabolites and their application in the determination of gene-set statistics and associated P values was detailed by Refs. [28,66].
After selection of liver-specific genes, as well as genes involved in the immune system and inflammation co-expressed with redox genes, gene set enrichment analysis was performed with BINGO in Cytoscape [68], and the significantly enriched biological processes (Q b 0.05) were reported.

Tumor stratification and comparison
Tumor stratification was carried out through Consensus Clustering [29,31] for those genes displaying a median FPKM N 1 across all samples after row-normalization of gene expression. This unsupervised method permits the determination of an optimum (i.e., stable) number of nonoverlapping clusters. Briefly, the data were resampled 1000 times by considering 80% sample (i.e., tumor) and feature (genes) resampling to achieve robust clustering. Resampled data were then transformed into a similarity matrix, the consensus matrix. Agglomerative hierarchical clustering was performed to stratify tumors using the consensus matrix based on Pearson correlation distances through the R package ConsensusClusterPlus [31]. The optimum cluster number was determined by testing 2 to 10 clusters, and based on CDF and Δ(K), the area that increased in the cumulative distribution function increased as the number of clusters.
2.6. Genome-scale metabolic modeling and essentiality analysis RNA expression data were integrated into the iHCC models to generate personalized genome-scale metabolic models using tINIT [15] and RAVEN [69]. The following thresholds for gene levels were considered: no expression (FPKM b 1), low expression (1 ≤ FPKM b 10), medium expression (10 ≤ FPKM b 50), and high expression (FPKM ≥ 50). Upper and lower bounds of the exchange reactions considered were based on Ref. [70].
To model HCC progression or the differences between tumor clusters, we used MADE [34] and TIGER [35] and used as input the genespecific fold changes and Q values. MADE integrates expression data and a metabolic model without the prior definition of activity thresholds. MADE utilizes an optimization-based approach to generate a sequence of expression states reflecting the most significant changes in the series of gene expression measurements, simultaneously solving the Flux Balance Analysis problem for each condition, leading to a sequence of functioning models. To implement MADE, we used TIGER, which uses a custom MATLAB class CMPI to solve mathematical programming problems. The Common Mathematical Programming Interface defines a consistent structure for mixed-integer linear programming and mixed-integer quadratic programming problems, providing independence from the underlying mixed-integer linear programming solver software, and MADE is formulated as a single mixedinteger linear programming problem. CPLEX solver was also employed. In this study, genes were considered multilevel instead of in a binary state during the implementation of MADE. Using MADE and TIGER, we reconstructed grade-specific genome-scale metabolic models on the basis of their corresponding gene expression data. We integrated the expression data achieved from DESeq2 [65] and considered the upper and lower bounds for exchange reactions in GEMs based on experimental liver data [70].
Here we considered three processes as important cellular objectives by considering the following fake reactions as objective functions for maximization: biomass production, NADPH production and H 2 O 2 scavenging (see below). These objective functions maximize the fluxes of NADPH oxidation and H 2 O 2 production, such that the optimization method redirects the cell's resources towards reducing NADP +

Survival analysis and signature identification
Kaplan-Meier survival analysis was performed to compare two groups of stratified tumors and based on the 10-year survival data for all 360 patients in R through the package survival [71].
To determine the best survival signatures, we generated all possible signatures consisting of 3 or 4 genes by taking the average expression of the selected genes in tumors from the hALDH2 cluster. Then, all 360 patients were stratified based on the correlation between their expression of selected genes and the signature using a cutoff that was maximally selected similarly to our previous work [27]. The minimum group size considered was 20% of the total sample size. The signature that presented the lowest log-rank P value after multiple hypothesis correction among all the investigated gene combinations was then selected as the best one. KM analysis had been performed for single gene predictions [27], and was corrected for multiple hypothesis testing considering all redox genes.

Microarray and metabolomic analyses
Three separate microarray cohorts were used to validate the coexpression and gene expression patterns observed using transcriptomics. The microarray for 61 HCC samples was downloaded from GEO with accession number GSE76297 (Ref. [14]) for the Affymetrix Human Transcriptome Array 2.0 platform. Row data were normalized using the Robust Multi-array Average followed by quantile normalization. Log 2 -transformed data was then analyzed in R, with probes annotated through the hta20sttranscriptcluster.db package. Normalization of the second cohort of 91 HCC samples with accession GSE1898 was carried out as previously indicated [32,33]. The final microarray cohort of 221 hepatocellular carcinoma patients was retrieved from GSE14520 and tumor normalized data was used for all analyses [42,43].
The metabolomics data of the 127 most variable metabolites between HCC and intrahepatic cholangiocarcinoma, after ignoring unidentified metabolites, were analyzed for the 54 HCC samples with both gene expression and metabolomics data [14]. Correlations between gene vs gene, and gene vs metabolite were computed across the 61 and 54 HCC samples, respectively.
2.9. mRNA expression of redox genes in murine HCC samples Total RNA was isolated from tumor and noncancerous liver tissues from six sets of DEN-treated (HCC) CB1R +/+ mice and CB1R −/− mice [46]. Next, 100 ng/sample of rRNA-depleted RNA was treated with RNase III to generate 100-to 200-nt fragments, which were pooled and processed for RNA sequencing. All data were normalized based on housekeeping genes through the Program CLC Genomics Workbench (version 5.1, CLC Bio). Absolute numbers were extracted from the reads, and the data were adjusted to non-HCC biopsies for each gene.

Statistical analyses
Statistical approaches were chosen according to the diverse data types in this study and indicated throughout. Hypothesis testing was performed by considering the null hypothesis of the absence of an association between the compared variables. The association with clinical features was tested according to the nature of the data: continuous vs continuous (Spearman rank correlation test); continuous vs categorical (Kruskal-Wallis test); categorical vs categorical (Fisher's exact test for 2 × 2 tables, or Chi-squared test). Q values were computed from the P values (false discovery rate with Benjamini-Hochberg procedure, family-wise error rate 0.05).

Data availability
The GEMs used in this work are found in SBML format (Computer Code 1). The authors declare that data supporting the findings of this study are available within the paper and its supplementary information files, or otherwise clearly point to the online resources where they were obtained.

Redox genes are co-expressed with liver-enriched, immune and inflammation genes, and differentiation markers
To clarify the role of genes involved in redox metabolism, we identified a panel of 132 protein-coding genes directly involved in this process through a literature search (Methods, Dataset 1). Specifically, we selected genes involved in i) ROS production, scavenging, and metabolism; ii) the metabolism of reducing equivalents and iii) genes that are essential for the response to oxidative stress using a generic human Genome-scale Metabolic Model (GEM) [15] and other publicly available biochemical reaction databases. We retrieved the gene expression data of 360 HCC tumors together with patient metadata from NIH's TCGA [18] and found that HCC tumors display highly heterogeneous expression of genes that are directly involved in redox metabolism (Fig. 1A).
Gene co-expression analysis is a useful approach to reveal functional relationships [19,20]. We calculated pairwise Spearman rank correlations between these 132 genes and all other genes expressed in HCC tumors [21]. This unsupervised clustering approach identified genes showing coordinated co-expression and are likely co-regulated or functionally associated. The 132 genes were co-expressed with each other and with many other genes in the genome (Dataset 2), including genes involved in carbohydrate, lipid, nucleotide and amino acid metabolism. For instance, aldolases (ALDOA and ALDOB), pyruvate kinases (PKM, PKLR) and phosphofructokinases (PFKP) were highly coexpressed with these 132 genes (highest Spearman's ρ 0.71, Q b 10 −50 ). Importantly, some enzymes involved in central metabolism (e.g., pyruvate kinases, isocitrate dehydrogenases, glutaminases) are not only important players in cancer but are also indirectly involved in redox metabolism [22]. Aiming to explore the relationships between genes directly or indirectly involved in redox metabolism, we expanded the list to 174 genes by considering those involved in the metabolism of carbohydrates, lipids or other compounds that are closely associated with redox metabolism (Dataset 1). This 174 highly heterogeneous gene set (Fig. 1A, Dataset 1) is hereafter referred to as "redox genes".
Among the significantly correlated redox genes, we identified two major clusters with positively co-expressed genes in the same cluster but negatively co-expressed with those of the opposite cluster (Q b 0.01, Fig. 1B, Dataset 3). The first cluster of the top 10 correlated genes ( Fig. 1B inset) comprised genes involved in the pentose phosphate pathway (PPP) (G6PD, RPIA), glycolysis (PKM, ALDOA), differentiation (TXNDC9), response to hypoxia (HIF1A), and metabolism of folate (MTHFD1L), glutathione (GLRX3), purines (ATIC), and amino acids (GLS) and is hereafter referred as the G6PD cluster. The second cluster also comprised genes involved in glycolysis and the PPP (ALDOB), metabolism of folate (MTHFD1, MTHFS, SHMT1), amino acids (ALDH5A1, ALDH6A1, ALDH7A1) or other compounds (ALDH2, ALDH8A1), and responses to oxidative stress (CAT) and is hereafter referred to as the ALDH2 cluster. Glucose 6 phosphate dehydrogenase (G6PD) controls the flux into the PPP and provides reducing power and ribose phosphates to maintain the redox balance and biosynthesis of nucleotides and lipids. These reducing equivalents may be used towards peroxidase-catalyzed antioxidant defense, which is linked to apoptosis, angiogenesis, and the efficacy of anti-cancer therapy and is identified as a promising target in cancer therapy [23]. ALDH2 is involved in alcohol metabolism and redox homeostasis, and it has recently been associated with HCC progression and a poor prognosis in mice [24].
Interestingly, genes in the ALDH2 and G6PD clusters were coexpressed (absolute Spearman ρ N 0.5, Q b 10 −20 ) with 1) 14 genes involved in inflammation; 2) 696 genes involved in immune mechanisms, including innate and adaptive immunity; and 3) 60 liver-enhanced and 89 liver-enriched genes [25]. Many of these were co-expressed with N5 genes in the ALDH2 and G6PD clusters but negatively and positively with respect to each cluster ( Fig. 1C inset). For instance, IL1, TNFα and NFκB, key genes involved in inflammation, were negatively coexpressed with the ALDH2 cluster (Q b 10 −4 ) but positively coexpressed with the G6PD cluster (Q b 0.001). In turn, liver-specific genes were negatively co-expressed with genes in the G6PD cluster but positively co-expressed with the ALDH2 cluster, in which we identified the liver-specific gene ALDOB [25]. This finding suggests that the G6PD cluster is associated with a poorly differentiated tumor phenotype, as opposed to the ALDH2 cluster. For instance, the G6PD cluster gene MTHFD2 has been previously identified as a marker of poor differentiation in other cancers [26]. Indeed, vimentin, a marker of poor differentiation and mesenchymal transition, was positively co-expressed with 9 genes in the G6PD cluster and negatively co-expressed with 7 genes in the ALDH2 cluster (absolute Spearman's ρ N 0.2, Q b 10 −5 , Dataset 4). Alpha-fetoprotein, also a marker of undifferentiated tumor cells, was positively co-expressed with 7 genes in the G6PD cluster and negatively co-expressed with 9 genes in the ALDH2 cluster (absolute Spearman's ρ N 0.2, Q b 10 −4 ). Similar observations were obtained with respect to other markers of differentiation such as cytokeratins (e.g., KRT8, KRT19).
Gene set enrichment analysis (GSEA) for liver-specific or immune genes (Fig. 1D, Dataset 5) indicated that the liver-specific genes were involved in lipid, steroid, small molecule and amino acid metabolism; small molecule metabolism, including hormones, cofactors and ketones; and oxidation/reduction processes. In turn, the immune system genes were involved in cell cycle and division; regulation of transport, metabolism and cell death; ubiquitination, regulation of protein metabolism and phosphorylation; proteasomal-mediated protein catabolism, cytoskeleton and organelle organization. Finally, we observed that several genes were co-expressed with 10 or more genes in both clusters: liver-specific genes, including ABAT, CDO1, CYP4F3, F11, SEC14L2, SPDYC, and TTC36, were positively and negatively co-expressed with genes in the ALDH2 and G6PD clusters, respectively. Genes associated with the adaptive immune system, including PPP2R1A, UBE2L3, VASP, and YWHAZ, and innate immune system, including ARPC2 -4, CFL1, GYG1, IKBKE, LIMK1, RHOG, and STK11IP, were negatively co-expressed with genes in the ALDH2 cluster but positively co-expressed with those in the G6PD cluster.

Several redox genes are associated with survival outcome in HCC patients
To reveal the association between the expression of redox genes and patient survival, we selected 50 subjects with the lowest and highest survival based on patient metadata and performed a differential expression analysis between the two groups (Methods, Dataset 6). We observed that 11 and 14 redox genes (Table 1) were significantly (Q b 0.05) up-and downregulated between patients in the low vs high survival groups, respectively. For instance, TXNRD1 and GSR were upregulated in the low-survival group and found to encode NADPdependent reductases involved in peroxidase-dependent H 2 O 2 defense. G6PD, which was upregulated in the low-survival group, encodes glucose-6-phosphate dehydrogenase, one of the main NADPH sources in cancer cells. In addition, many redox genes have been identified as potential prognostic genes (Table 1) [27]. Here, ALDH2 was identified as a potential prognostic gene since its downregulation in the low survival group is in agreement with its identification as a favorable prognostic gene in liver cancer [27]. These observations indicated that several genes involved in redox metabolism were significantly correlated with survival outcome in HCC patients.
GSEA between low-and high-survival patients (Q b 0.05, Dataset 6) indicated that subjects with low survival tended to over represent genes involved in cell cycle regulation and related processes (e.g., chromatid/chromosome segregation and regulation, DNA replication, initiation and integrity checkpoint, spindle organization and recombination, regulation of cell division). In turn, these patients underrepresented processes such as metabolic and energy regulation (e.g., electron transport chain, aerobic respiration, oxidative phosphorylation, metabolism of amino acids, fatty acid β oxidation and catabolism, lipid oxidation). Only two upregulated redox genes (GART and HIF1A) were associated with overrepresented processes (Fig. 2, Dataset 7). However, several downregulated genes were associated with underrepresented processes, including energy and reducing equivalent metabolism (ALDH5A1, NNT), folate and amino acid metabolism (MTHFS, SHMT1), and detoxification (GSTO1, GSTK1), among other metabolic processes (Dataset 7). These observations suggested that the biological differences displayed between subjects with low vs high survival involved several redox genes associated with metabolic processes that are important in the occurrence of HCC.
3.3. Genes in ALDH2 and G6PD clusters are associated with opposing functional responses, prognostic genes, isoform-encoding genes, and hypoxia behavior Based on our integrative analysis, we observed that HCC tumors displayed different gene expression patterns of redox genes in each cluster, which might translate into functional metabolic differences. To understand the key role of each redox gene in cancer, we classified HCC tumors based on expression and compared the global gene expression profiling of the 50 tumors with the highest and lowest expression of individual genes in each cluster (see Methods). Differential expression with respect to majority of the genes in the ALDH2 and G6PD clusters indicated similar expression patterns between redox genes that were distinct from those in the opposing cluster (Appendix Fig. S1B), in agreement with the above-described positive intra-cluster and negative inter-cluster correlations (Fig. 1B).
GSEA revealed underlying metabolic changes between tumors with high expression of genes in each cluster (Fig. 3). We found that tumors displaying high expression levels of genes in the ALDH2 cluster, or low expression of those in G6PD cluster, were enriched in several processes related to drug metabolism, lipid oxidation and metabolism, amino acid metabolism and biosynthesis, and carbohydrate metabolism. In turn, Fig. 1. Highly heterogeneously expressed redox genes comprise two major clusters of co-expressed redox genes, highly co-expressed with genes that are liver-specific, involved in immune system or inflammation. A. Hierarchical clustering of 174 genes involved in redox metabolism across 360 HCC subjects showed highly heterogeneous redox responses. Genes (rows, mean FPKM N 1) were Z-normalized and clustered. B. Two redox gene clusters were highly positively co-expressed with genes within the clusters (red) but negatively co-expressed with genes of the opposing cluster (blue). Only statistically significant (Q b 0.01) correlations are presented, ordered according to first principal component. The inset displays the 10 genes in each cluster with the highest correlation coefficients. C. Genes in G6PD and ALDH2 clusters were strongly co-expressed with liver-specific genes and genes involved in inflammation (Q b 10 −20 , absolute Spearman ρ N 0.5). Inset: genes that were simultaneously co-expressed positively and negatively with N5 genes in ALDH2 and G6PD clusters and at least one gene in each cluster. Red edges indicate positive co-expression, blue edges indicate negative co-expression. Inflammasome categories: 1. Inflammasome; 2. CLEC7A/Inflammasome pathway; 3. DEx/H-box helicases that activate type I IFN and inflammatory cytokine production. D. Gene set enrichment analysis performed for liver-specific and immune system genes co-expressed with N5 genes in the ALDH2 and G6PD clusters (C inset). All colored nodes are statistically enriched functional terms (Q b 0.05). Table 1 Favorable and unfavorable survival and prognostic genes in HCC involved in redox metabolism. Significantly altered redox genes between 50 high vs 50 bottom survival subjects are displayed, where fold changes N 1 indicate an upregulated gene among the poor survival group, and prognostic genes in liver cancer [27] are shown. Favorable and unfavorable prognostic genes are those respectively up-and downregulated in the high survival group. CAT was identified as a prognostic gene in liver cancer in agreement with its expression variation between survival groups, despite being marginally not statistically differentially expressed.

Module
Gene tumors with low expression levels of genes in the ALDH2 cluster, or high expression of those in the G6PD cluster, were enriched in several processes related to development, differentiation, and morphogenesis and, to a smaller extent, fatty acid biosynthesis (Dataset 8). No functional terms related to NADPH biosynthesis were downregulated with ALDH2 expression and upregulated with GPX and PRDX expression. Tumors enriched in genes co-expressed with those in the ALDH2 and G6PD clusters (Fig. 4) tended to display similarly antagonistic responses (Appendix Fig. S2). Similar observations were also identified in terms of CAT and G6PD, two important genes for defense against H 2 O 2 (Doc S1, Appendix Fig. S3).
To reveal the detailed metabolic differences between tumors with high expression of genes in each cluster, we performed reporter metabolite analysis [28] using a generic human HCC GEM [15]. We found significant differences between reporter metabolites, with the most significant transcriptional changes observed in their associated enzymes, considering the full cellular metabolic network [28]. We identified 1999 reporter metabolites (out of 5134) that displayed the opposite behavior with respect to genes in ALDH2 and G6PD clusters (Q b 0.05, Dataset 9). Reporter metabolites identified in our study involved in redox metabolism included H 2 O 2 , NAD(P)H, NAD(P) + , glutathione, CoA, and O 2 • − (Dataset 9), as well as several metabolites involved in lipid metabolism (e.g., hydroperoxyeicosatetraenoic acids). These trends were consistently observed across cellular compartments, where 1566 reporter metabolites displayed similar antagonistic behaviors to those described above. For instance (Fig. 3B), genes involved in reactions involving NAD(P)H and NAD(P) + were consistently upregulated in all compartments, with the exception of the nucleus, in subjects displaying high expression of genes in the ALDH2 cluster, or low expression of those in the G6PD cluster. Likewise, genes involved in H 2 O 2 -associated reactions were upregulated and downregulated in subjects displaying high ALDH2-and G6PD-co-expressed genes, respectively.
The above-described antagonistic behavior (Fig. 1) extended to the genes co-expressed with those in the two clusters, which tended to be negatively co-expressed with those of the opposing cluster (Fig. 4A,  Appendix Fig. S4). Importantly, the ALDH2 and G6PD clusters were co-expressed with the genes encoding for enzyme isoforms or enzymes that catalyze alternative reactions in the same pathway (Fig. 4B, Q b 10 −3 ). For instance, among the genes in the G6PD cluster or their co-expressed genes (Fig. 4 right), we found that PKM, MTHFD2, MTHFD1L, and ALDOA, as opposed to genes in the ALDH2 cluster or their first neighbors PKLR, MTHFD1, and ALDOB (Fig. 4  left). This opposite behavior was not observed for GLRX3 vs GLRX5, IDH3ABG vs IDH1, or ACSS1 vs ACSS2, as previously observed [16]. Transcription factors and genes involved in the response to oxidative stress, TP53, HIF1A, and MTH1, were co-expressed with genes in the G6PD cluster, whereas FOXO1 was co-expressed with genes in the ALDH2 cluster.
The observation that HIF1A was part of the G6PD cluster suggested that HCC tumors and their redox behavior responded to hypoxia. We observed that PKM and HIF1A were co-expressed, consistent with the response of PKM to hypoxia [16]. Interestingly, we found that HIF1A expression reflected the antagonistic responses of ALDH2 and G6PD clusters observed herein (Fig. 4C, Q b 0.02), in which the only genes that showed no significant difference with respect to HIF1A expression were MTHFD1, MTHFS and CAT (Q N 0.05).
Additionally, we identified 7 unfavorable prognostic genes and favorable prognostic genes among the G6PD cluster and its coexpressed genes, and 9 favorable prognostic genes among the ALDH2 cluster and its co-expressed genes. Functional enrichment analysis with respect to tumors with high expression of each gene showed similar observations to those described above (Appendix Fig. S2).

Tumors stratified based on redox genes display substantial metabolic differences
The above observations indicated potentially important and opposite responses in HCC subjects differentially expressing redox genes. We then questioned whether HCC tumors could be stratified based on genes in ALDH2 and G6PD clusters. Hierarchical clustering of tumors has some pitfalls that have been overcome by recent methods that achieve robust stratification even considering small gene subsets [29,30], and does not permit a clear separation of tumors (Fig. 1A). Unsupervised clustering methods permit the identification of unknown groups in a sample based on intrinsic features without external information [31]. Here, we employed consensus clustering for partitioning patients into stable groups through repeated subsampling and clustering [29,31]. Briefly, considering genes in the ALDH2 and G6PD clusters (Fig. 1B inset), we determined an optimum number of clusters of 2 ( Fig. 5A) by determining the highest area increase in the cumulative distribution function as the number of clusters increased.
We repeated these analyses using gene expression data in independent human HCC cohorts [14,32,33] and identified 2 tumor clusters (Appendix Fig. S5). The two clusters (Fig. 5A) comprised 190 and 170 patients with high-expression of ALDH2-and G6PD-co-expressed genes, respectively, and hereafter are referred to as the hALDH2 and hG6PD groups. Consensus clustering based on the full set of redox genes yielded similar clusters. Importantly, among the top differentially expressed genes between the two clusters, we identified ALDOB, PKM, G6PD and MTHFS (Appendix Fig. S6, Dataset 10).

Simulation of the metabolic differences between HCC tumor types
To reveal the metabolic differences between HCC subpopulations, we generated two HCC GEMs, one per patient cluster, using the MADE [34] and TIGER [35] algorithms by incorporating gene-specific fold changes and Q values. Considering tumor growth as an objective function, we observed that models for hG6PD and hALDH2 tumors (Fig. 5B top and bottom) displayed major differences in terms of H 2 O 2 consumption. NADPH-independent and -dependent systems, i.e., reactions catalyzed by catalase but not by peroxidases, displayed fluxes in hALDH2 but not in hG6PD tumors (compare NADPH-independent reaction 4 with NADPH-dependent reaction 6). NADPH-production by glucose-6phosphate dehydrogenase mirrored this trend, displaying low flux in hALDH2 but high flux in hG6PD (Dataset 11, Appendix Fig. S7). Additionally, mitochondrial superoxide dismutase-catalyzed H 2 O 2 production (Fig. 5B, reaction 7) displayed flux and promoted H 2 O 2 permeation to the cytoplasm (reaction 8) in hALDH2. This was not observed in hG6PD, where cytoplasmic superoxide dismutase produces H 2 O 2 that enters mitochondria and is therein consumed by peroxidases (reaction 6). Finally, reactions involved in fatty acid β-oxidation (reactions 1-3), and purine metabolism (reactions 13 and 19) and that produce H 2 O 2 displayed fluxes in hALDH2 but not in hG6PD. NADPH metabolism through the folate pathway displayed fluxes in hALDH2 but not in hG6PD (reactions 15, 25, 28 and 33, Appendix Fig. S7). In turn, reactions involved in mitochondrial but not cytoplasmic reduction of oxidized glutathione (reaction 27) displayed fluxes in hALDH2 but not in hG6PD, reflecting the increased peroxidase activity in the latter group of tumors. Similar observations were observed through patient-specific GEMs, reflecting the major redox patterns of the two stratified groups (Dataset 12). Reporter metabolite analysis performed between hALDH2 vs hG6PD qualitatively reproduces (Dataset 12) the observations found for gene-specific reporter metabolites (Fig. 3C).
To validate these observations, we repeated these analyses using gene expression data from an independent cohort [32,33]. Overall, these cohorts displayed metabolic features that were qualitatively similar to those observed above, including fluxes in glutathione-reductase driven H 2 O 2 scavenging and glucose-6-phosphate dehydrogenase, and fatty acid biosynthesis but not β-oxidation in hG6PD, and not in hALDH2 (Dataset 13).
Essentiality analysis performed on each of the genes in the ALDH2 and G6PD clusters indicated that they were essential for hALDH2 and hG6PD tumors, respectively, but not the converse (with the only exception being GLRX3, which is not in the model).

Redox responses vary with tumor progression and subtype and identify survival signatures across HCC patients
The above observations indicated that HCC tumors display substantial metabolic differences, leading us to question how these differences relate to clinical properties of the patients. For instance, survival analysis indicated that hALDH2 patients displayed significantly higher survival than hG6PD patients (median 5-year survival of 0.57 and 0.39 respectively, log-rank test P = .0017, Fig. 5C). We also identified significantly altered redox gene expression between patients with different tumor grades. Cancer grades identify the dedifferentiation stage of HCC with respect to normal cells, and they include well (Grade 1), moderately (Grade 2), poorly (Grade 3) differentiated and undifferentiated (Grade 4) tumors. Here we identified differentially expressed genes throughout cancer progression between early (G1), intermediate (G2) and late (G3 + G4) HCC. We observed that few redox genes were significantly (Q b 0.05) different between early and intermediate cancer, but many of them were significantly differentially expressed between intermediate and late cancer (Fig. 5D, Dataset 14). Notably, the differences followed the shifts in gene expression observed between the opposing G6PD and ALDH2 clusters. Namely, expression of genes in the G6PD cluster was upregulated throughout cancer progression, whereas expression of genes in the ALDH2 cluster decreased with cancer progression. Metabolic modeling and accounting grade-specific gene expression for early and late HCC showed trends similar to those observed in hALDH2 and hG6PD tumors above (Fig. 5B). These observations indicated that through HCC progression, and particularly in more advanced cancer, there was a switch from NADPH-independent to NADPH-dependent H 2 O 2 defense supported by glucose-6-phosphate dehydrogenase activity.
Detailed clinical information, including the tumor subtype, viral infection and inflammation status, is available for 186 patients [13]. Analysis of clinical data (Fig. 5E) indicated that the proportions of tumors in hALDH2 and hG6PD groups were not independent of clinical traits (Dataset 15) such as TP53 mutation incidence, histological grade and stage (Q b 0.05), tumor growth pattern (Q b 0.005), and molecular data including mRNA and microRNA expression and protein levels (Q b 0.01). The proportion of tumors in the two clusters also varied with the tumor subtype: NCI proliferation (NCIP) [32], cholangiocarcinomalike (CCL) [36] and hepatoblastoma 16 gene (HB16) [37] signatures (Q b 10 −5 ); Hoshida [38] and 65-gene recurrence risk score (RS65) [39] signatures (Q b 10 −4 ); Seoul National University recurrence (SNUR) [40], and Hippo pathway [41] signatures (Q b 0.05). In contrast, the proportion of tumors in hALDH2 and hG6PD did not vary (Q N 0.2) with Fig. 3. Opposite functional responses are observed in HCC tumors differentially expressing redox genes. A. Using metabolomics and gene expression data [14], we identified several metabolites with abundances that positively (red) or negatively (blue) correlated with the expression of genes in the ALDH2 and G6PD clusters (Q b 0.05, absolute Spearman's ρ N 0.4). Genes with expression that was not quantified by Ref. 14 are not indicated. B. GSEA for subjects expressing high and low redox genes in the ALDH2 (blue) and G6PD (orange) clusters. Only processes that were significant (Q b 0.05) in N4 columns are displayed (see Methods). Log 2 -normalized test statistics indicate the directionality of the processes (positive are upregulated, negative are downregulated, null are not statistically significant). C. Reporter metabolites were determined, and the minimum Q was computed for those with associated genes that were consistently up-or downregulated across compartments (Q b 0.05). Q values were normalized by computing −1 × log 10 (minimum Q), or log 10 (minimum Q), respectively, for up-or downregulated changes. race, gender, alcohol liver disease or consumption, viral infection (hepatitis B and/or C), inflammation (including intra-tumoral), cirrhosis, fibrosis, tumor ploidy, purity, or DNA methylation. At the gene expression level (Fig. 5E, Dataset 15), we observed that expression of most genes in the ALDH2 and G6PD clusters (≥17 genes, Q b 10 −4 ) varied with the tumor subtype (CCL, HB16, Hoshida, NCIP, SNUR, Hippo), RS65, and mRNA expression. The expression of many genes (≥11) also varied with miRNA expression and protein levels.
Negative correlations were identified between RS65 index score and the expression of genes in the ALDH2 cluster (ALDH2 and MTHFD1, ρ b −0.3, Q b 0.0005), but positive correlation were observed with those in the G6PD cluster (PKM, MTHFD1L, G6PD, ATIC, GLRX3, ALDOA, RPIA and Fig. 4. HCC tumors antagonistically co-express genes in ALDH2 and G6PD clusters, associated with opposite hypoxic environments. A. Genes in the two clusters (yellow) and their first neighbors (gray) are indicated. Redox genes that were positively co-expressed with genes in each cluster are displayed (absolute Spearman's ρ N 0.3, Q b 10 −7 ). Green-and purplecircled nodes respectively indicate favorable and unfavorable prognostic genes ( Table 1). Genes that were positively co-expressed with those in both clusters are not presented (e.g., NFE2L2, NNT, XCT, TXNRD1). Refer to Dataset 3 for the full list of correlations. B. Genes encoding enzyme isoforms or enzymes that catalyze alternative reactions in the same pathway were alternatively expressed. Boxplots display the expression of genes (FPKM) as a function of HCC samples displaying high (H) and low (L) expression of each gene (50 samples per group), and the respective Q values (Mann-Whitney U test). No significant difference (Q N 0.05) was observed for the other isoforms, with the exception of IDH1 as a function of IDH3A (Q b 0.005), for which high expression of the latter gene was associated with high expression of Idh1. Similar observations were obtained through differential expression analysis. C. Gene expression as a function of HIF1A expression. Note that among all genes in the ALDH2 and G6PD clusters, only the expression of MTHFD1, MTHFS and CAT did not significantly (Q b 0.05) differ between samples with high vs low HIF1A expression. For all boxplots, whiskers extended up to 1.5 times the interquartile range. TXNDC9, ρ N 0.3, Q b 0.0006). Only the expression of genes in the G6PD cluster (G6PD, RPIA, PKM, ALDOA, ATIC, and MTHFD1L) strongly varied with the TP53 mutation incidence (Q b 10 −4 ), which was also observed to a lower extent for genes in the ALDH2 cluster (ALDH2 and ALDH5A1, Q b 0.01).
The above observations indicated that HCC tumors were successfully stratified into different groups with distinct redox behavior and survival (Fig. 5) and that several redox genes were potentially favorable and unfavorable prognostic genes (Table 1). We sought to determine a signature that would permit a more efficient prediction of patient survival by using a subset of genes in the redox clusters. We examined all 3 by 3 and 4 by 4 combinations of genes in the redox clusters and, through Kaplan-Meier analysis, identified PKM, G6PD, MTHFS, and HIF1A as achieving the best survival signature where the group size comprised made up at least 20% of the total sample (Fig. 6A, log-rank P b 10 −16 , Q b 10 −16 , Dataset 16). PKM, MTHFS and G6PD also presented the most significant 3 by 3 combination (log-rank P b 10 −15 , Q b 10 −12 ) and were substantially better than any single gene prediction (logrank P b 10 −12 , Q b 10 −9 , Appendix Fig. S8, Dataset 16, analysis based on ref. [27]). This 4-gene signature included genes from pentose phosphate pathway and folate metabolism, suggesting that the survival of HCC patients was associated with these pathways.
We performed Kaplan-Meier analysis in independent cohorts [42,43] and validated the prognostic value of this gene signature (Fig. 6A, log-rank P b .01, Q ≈ 0.023). Taken together, these observations indicated that tumors displayed substantial differences in redox behavior that were related to the clinical and survival differences of their respective patients.

Validation of the redox response in independent human HCC cohorts
The above relationships were validated using several independent cohorts at gene expression protein and metabolite levels. First, the antagonistic behavior observed among co-expressed redox genes (Fig. 1) was also displayed at the gene expression level in 61 HCC tumors (Appendix Fig. S9A, Ref. [14]) and 91 HCC tumors [42,43] (Appendix Fig. S2A), and at the proteomic level in 19 HCC tumors [44] (Dataset 17).
Second, the above antagonistic behavior was observed in an independent microarray and metabolomics tumor dataset [14] (Fig. 3A). The authors used multi-omics data, including gene expression and metabolomic data that quantified 718 metabolites using LC/MS and GC/MS, to identify and identified molecular subtypes and driver genes shared by HCC and intrahepatic cholangiocarcinoma. The integrative approach employed by the authors permitted identifying allowed the identification of potential biomarkers and shared metabolic features in common in the two liver diseases despite their intrinsic heterogeneity. Using gene expression and metabolomics data from this cohort, we observed that the abundance of glycolytic/gluconeogenic compounds (e.g. glucose and oxaloacetate), pyrimidines (e.g. thymine) and other compounds was directly or inversely proportional to the expression of genes in the ALDH2 and G6PD clusters. Notably, we observed that the reporter metabolite analysis between the clusters was consistent with the observed gene-metabolite associations. For instance, we observed that 75% of the metabolites simultaneously found by the reporter metabolite analysis and gene-metabolite association analysis (Fig. 3) showed concordant behavior between the two independent analyses on the two datasets, i.e., positive correlations with ALDH2-cluster genes and upregulated reporter metabolites in hALDH2 tumors, and the same with respect to G6PD cluster/hG6PD. Analyses performed using the additional cohort supported these observations (Dataset 13). The metabolites identified included glucose, oxaloacetate, and cholate derivatives (positively correlated with ALDH2 cluster genes), and guanosine, thymidine and deoxyinosine (positively correlated with G6PD cluster genes). While the reporter metabolite analysis alone did not permit inferences of metabolite abundance, they pointed towards metabolic hotspots with likely increases or decreases in metabolic fluxes [28]. The observed metabolite-gene associations suggested that aside from representing these flux-alteration hotspots, the abundances of those metabolites simultaneously identified in the 2 analyses were associated with redox genes and thus might represent potential biomarkers in ALDH2 or G6PD clusters.
Third, we observed that patient stratification according to genes in ALDH2 and G6PD clusters was associated with differences in survival (Fig. 5C) and known HCC subtypes (e.g., NCIP, Fig. 5D), which were reproduced using an independent cohort (Appendix Fig. S9B) [32,33].
In addition to these observations, we have previously shown that similar antagonistic patterns are observed for other genes associated with the G6PD cluster (ACSS1, PKM) and ALDH2 cluster (ACSS2) [16,45]. Finally, we have identified the expression of G6PD, HIF1A, PFKP, CAT, and MTHFD1 using the antibodies HPA000834, HPA001275, HPA018257, HPA055838, and HPA050052, respectively [25,27], in human HCC and noncancerous samples (Fig. 6B). The G6PD-clustered genes HIF1A, G6PD and PFKP show moderate/high protein expression in HCC samples but low expression in noncancerous samples, whereas the ALDH2-clustered genes CAT and MTHFD1 show medium protein expression in HCC but high expression in noncancerous tissues.

Validation of redox responses in independent mouse HCC models
Retrospective analysis of murine and human samples [16,46] showed agreement between hypoxic behavior and progression for the 2 redox gene clusters. G6PD-clustered genes were associated with hypoxic behavior, a less favorable prognosis or accumulated throughout progression when compared with ALDH2-clustered genes. Our observations (Fig. 6C, D) indicated that most genes associated with the G6PD cluster displayed more marked expression increases between HCC vs noncancerous samples (e.g., Hif1a, G6pd, Pfkp, Gart, Me1) than those in the ALDH2 cluster (e.g., Aldh2, Cat, Shmt1, Mthfd1) in N-nitrosodiethylamine (DEN)-induced HCC in wild-type (wt) mice. These changes were more pronounced when comparing expression levels of HCC tumors in CB 1 R +/+ vs CB 1 R −/− mice since the wild-type mice were associated with hypoxic tumor responses when compared with the mutant littermates (Fig. 6C, D). Aldob presented as an exception to this behavior, possibly because no other aldolase was detected in either HCC or noncancerous samples. Overall, these observations indicate that the antagonistic redox Fig. 5. HCC tumors stratified based on redox genes display substantial differences in the survival and metabolism of H 2 O 2 , purines and fatty acids. A. Tumors were stratified based on the expression of genes in the ALDH2 and G6PD clusters through consensus clustering into two clusters that respectively displayed high expression of genes in ALDH2 and G6PD clusters (cyan and blue horizontal bars, respectively, hALDH2 and hG6PD). Clustering displayed the likelihood that a tumor was found in that cluster. Similar partitioning was achieved if all 174 redox genes were considered, reinforcing the robustness of the clustering. Note that N94% of the samples were consistently clustered in the two clusters (likelihood N 0.9). B. Metabolic differences displayed by the two groups of tumors as determined through metabolic modeling. GEMs for hG6PD (top) and hALDH2 (bottom) tumors were integrated with the gene expression data. Cell compartments are indicated by colored boxes. Numbers indicate enzyme-catalyzed reactions or transport reactions, as indicated in Dataset 11. Thin arrows indicate reactions with low or null fluxes (B: 6, 10, 11; C: 1-5, 7, 13, 19). Bidirectional reactions (e.g., transport reactions) indicate the estimated direction based on the flux value. Common metabolites (e.g., H 2 O and H + ) are not displayed. C. Kaplan-Meier survival plot for hALDH2 (cyan) and hG6PD (blue) clusters indicates substantially lower survival for patients in the latter cluster (respective median 5-year survival of 0.57 and 0.39, P ≈ 0.0018). D. Differential expression analysis showing the opposing behavior between ALDH2-co-expressed and G6PD-coexpressed genes. Cumulative log 2 (fold changes) were computed for the significantly differentially expressed redox genes (Q b 0.05, DESeq2) between late (G3 + G4) vs intermediate (G2) vs early (G1) HCC and noncancerous samples. Genes in the ALDH2 (cyan) and G6PD (blue) clusters (Fig. 1B, inset) behavior observed in human HCC is also displayed in a mouse model of HCC.

Discussion
Redox metabolism regulates cell proliferation and aberrant redox metabolism, and it is now recognized as one of the hallmarks of cancer [22,47]. Antioxidants scavenge ROS that may otherwise have toxic effects, such as promoting oxidative DNA damage and tumorigenesis. This phenomenon has led many researchers to propose the use of antioxidants as therapeutic agents against cancer [48,49], but several clinical trials have failed to demonstrate the beneficial effects of antioxidants and instead showed that they may increase cancer risk and promote progression and metastasis [5]. Due to the highly heterogeneous responses of HCC, including in terms of global antioxidant gene expression [5], the commonalities and relationships between redox metabolism, patient survival, and global metabolism had never been systematically investigated. Here we extended our previous pancancer analysis [27] to systematically analyze heterogeneous redox behavior across HCC tumors and used systems biology approaches to analyze transcriptomics, proteomics and metabolomics data. We considered genes directly and indirectly involved in redox metabolism, i.e., redox genes, and found that they were co-expressed with liverspecific genes and genes involved in the immune system and inflammation, all of which are involved in important processes associated with cancer, such as cell cycle progression, apoptosis, amino acid and lipid metabolism. The close relationships between redox metabolism, inflammation and immune defense are well known in other contexts [7,50], and the immune system is targeted for the treatment of HCC [51]. Our observations emphasized the necessity of jointly targeting immunity, inflammation and redox metabolism for the effective treatment of HCC.
Importantly, our unsupervised approach highlights opposite behaviors among redox genes through two antagonistic and highly coexpressed groups of genes co-expressed with ALDH2 or G6PD. ALDH2 is a gene with important roles in redox homeostasis associated with survival and lower malignancy [24]. In turn, the role of G6PD in cancer is widely known [22]. Their co-expressed genes have known relationships in tumorigenesis. For instance, ALDOA is co-expressed with G6PD and has been previously associated with the response to hypoxia and with HIF1A [52,53], similarly to our observations in HCC. Tumors with high ALDH2-co-expressed genes displayed high fatty acid and lipid oxidation, catabolism, and lower lipid biosynthesis, whereas those with high G6PD-co-expressed genes (hG6PD) showed the opposite behavior. Our previous analyses indicated that the metabolic alterations in HCC promote NADPH biosynthesis through the PPP, lipid biosynthesis and decreased fatty acid oxidation [16]. Here we observed that HCC tumors differentially displayed the above metabolic changes and were closely related to redox behavior (Fig. 7). The correlation between the expression of genes in the ALDH2 cluster with alterations in lipid, amino acid and carbohydrate metabolism additionally suggested that tumors with high expression of those genes might be metabolically more active. In turn, tumors with high expression of G6PD-clustered genes, including HIF1A, were associated with several markers of dedifferentiation, poor prognosis and low survival, consistent with the higher aggressiveness of those tumors.
HCC tumors displayed antagonistic redox responses at the gene expression, protein and metabolite level. For instance, hG6PD tumors displayed lower levels of glycolytic metabolites (including oxaloacetate and glucose) and higher levels of nucleoside and nucleobases (thymine, deoxyinosine and orotidine) and amino acid intermediates (acetylglycine, phenylalanine, γ-glutamyl tryptophan). Although other important metabolites, such as NADP and ROS, were not available in the metabolomics data [14], our modeling analyses indicated that genes involved in the metabolism of NAD(P)H, H 2 O 2 , and, to a smaller extent, O 2 •-, glutathione, and CoA, were significantly upregulated in hALDH2 tumors but downregulated in hG6PD tumors. This finding suggested that tumors responded differently to reducing equivalents and ROS and implied that targeting redox metabolism might lead to different outcomes depending on whether the tumors display the hALDH2 or hG6PD subtypes. Using metabolic models and essentiality analysis for hALDH2 and hG6PD tumors, we observed that the former group was not sensitive to inhibition of genes in G6PD cluster but displayed no growth if ALDH2-clustered genes ALDH2 are inhibited, and vice versa.
It is generally considered that cancer cells rely on high antioxidant activity, together with promoted NADPH biosynthesis that is used towards peroxidase-catalyzed defense, enabling tumors to proliferate under sustained oxidative stress [54]. HCC displayed lower glutathione content, catalase and glutathione peroxidase activities in comparison with normal cells [55,56]. Our analyses showed that hALDH2 tumors displayed substantial catalase-mediated H 2 O 2 scavenging, as opposed to hG6PD tumors that instead rely on glutathione/glutathioneperoxidase-mediated H 2 O 2 scavenging. Additionally, we observed that mitochondria were a significant source of H 2 O 2 in hALDH2 tumors, but not in hG6PD tumors, which displayed high mitochondrial glutathione-peroxidase activity. These observations highlighted important NADPH-independent and -dependent differences in ROS scavenging in HCC. Furthermore, they indicated that targeting antioxidants and redox metabolism for HCC treatment is highly dependent on the cell phenotype, i.e., hALDH2 or hG6PD, and thus requires subgrouptailored drug choices.
Major metabolic and functional alterations that are hallmarks of malignant transformation and often linked to different cancer stages are associated with differential redox metabolism [8]. Expression and functional switches are associated with genes in the ALDH2 and G6PD clusters and thus hALDH2 and hG6PD phenotypes. Importantly, a switch between the two phenotypes was observed from early to late HCC, associated with two main features. First, genes in the ALDH2 and G6PD clusters tended to be respectively co-expressed with favorable and unfavorable prognostic genes. Second, liver-specific genes and markers of differentiation loss respectively were co-expressed positively and negatively with genes in the ALDH2 and G6PD clusters. These two features suggested a malignant association with the expression of genes in the G6PD cluster. They were consistent with the low survival of patients with hG6PD tumors and the observed inverse relationship between the expression of ALDOB, a liver-specific gene in the ALDH2 cluster [57], and ALDOA, a G6PD cluster gene repressed in adult liver and expressed in other tissues [25]. Furthermore, a decrease in serum aldolase B and malignant liver cancer has been previously observed [58], whereas a positive relationship between ALDOA and malignancy and progression has been observed in other cancers [53,59].
Importantly, genes in ALDH2 and G6PD clusters were co-expressed with genes encoding alternative enzymes but that catalyze the same or alternative reactions, sometimes exhibiting the same flux directions. This suggests a compensatory behavior in the metabolic switching of the hALDH2 to hG6PD phenotype. This phenomenon was observed for the expression of MTHFD1, ALDOB and PKLR, as well as all parts of the ALDH2 cluster or their co-expressed genes, which were replaced by MTHFD2, MTHFD1L, ALDOA and PKM, and all parts of the G6PD cluster or co-expressed genes. Such switching Fig. 6. Predominant expression of G6PD-and ALDH2-clustered genes in mouse and human HCC and noncancerous cells. A. Kaplan-Meier plots for the 4-gene signature (G6PD, HIF1A, PKM, MTHFS) using RNA-seq data for 360 patients from TCGA (left), and from an independent microarray cohort with 222 patients [42,43]. B. Immunohistochemistry staining of G6PD, HIF1A, PFKP, CAT, and MTHFD1 through antibodies HPA000834, HPA001275, HPA018257, HPA055838, and HPA050052, respectively [25,27]. Protein expression is displayed in blue, with counterstaining in blue. C. Gene expression profile in murine samples revealed substantially higher expression of G6PD cluster genes in HCC than noncancer samples and higher expression in CB1R +/+ (hypoxia-associated) compared with CB1R −/− samples, compared with the expression of genes in the ALDH2 cluster. mRNA expression profile measured in six pooled wild-type HCC tumor and noncancerous samples using RNA-seq in CB1R +/+ and CB1R −/− six murine samples. For all boxplots, whiskers extended up to 1.5 times the interquartile range. D. The distribution of differences for genes in the two clusters in HCC and noncancerous (NC) murine samples. Fig. 7. HCC tumors display antagonistic metabolic behavior associated with redox responses. Genes of ALDH2 and G6PD clusters (blue and red nodes) promoted several reactions and processes, as summarized at the bottom of the figure. Reducing equivalents (NADPH and NADH) are highlighted (black nodes). Dashed lines indicate regulatory processes. Genes in the G6PD and ALDH2 clusters (black-circled red and blue nodes) and some of their co-expressed genes (Fig. 4). The reaction catalyzed by H6PD occurred in the endoplasmic reticulum. Note the differences in ROS metabolism in tumors with high expression of ALDH2-clustered genes (blue) vs G6PD-clustered genes (red). Prdx5 was expressed by hALDH2 and hG6PD tumors, but mitochondrial ROS consumption by this gene was only active in hG6PD. Expression profiles of G6PD, HIF1A, PKM and MTHFS, genes in the two redox clusters, identified a signature that more effectively predicted patient survival than predictions from a single gene. has been previously observed in acetate metabolism in HCC [16], and in glycolysis in rat hepatomas [60].
Tumor stratification based on the expression of genes in the ALDH2 and G6PD clusters, where ALDOB is the best stratifying gene, has shown not only differences at the functional and genomic levels but also at the clinical level, including TP53 mutation and known tumor subtypes [36]. Through patient metadata, we also identified several markers of favorable and unfavorable survival among genes involved in redox metabolism. The 4-gene signature herein introduced (PKM, G6PD, MTHFS and HIF1A) improved the commonly used single-gene survival markers and showed great potential for clinical use to predict the prognosis of HCC patients. These 4 genes and their individual association with redox metabolism and tumorigenesis are known in other cancers [22,[61][62][63]. However, here we showed for the first time the relationship of their joint behavior with metabolic and functional responses in HCC and translation to survival differences.
The current study also revealed gene expression behaviors of ALDH2 and G6PD clusters in a DEN-induced HCC murine model. In human HCC, high expression of the hG6PD cluster resulted in low survival and was dominant in late-stage HCC. In mice, G6pd cluster genes were more robustly increased in HCC tumors than Aldh2 clusters in wild-type mice. A similar polarization trend between humans and mice with a dominance of the G6PD cluster increased the translational value of DEN-induced HCC murine models for studying the preclinical efficacy of pharmacological agents in HCC as metabolism mimicking the low-survival stratified cluster. Furthermore, this study also exemplified an effect of deletion of a druggable gene, Cnr1, on the regulation of G6PD cluster genes in DEN-induced HCC in mice. Overactivation of CB 1 R has been shown to promote HCC initiation and progression [46]. Deletion of the Cnr1 gene significantly attenuated the activation of G6PD cluster genes in HCC. This finding supports previous observations of the potential involvement of CB 1 R in the stratification of HCC patients [16] and also reinstates the therapeutic potential of peripheral CB 1 R antagonism in HCC [46]. In conclusion, our systems level analyses of multi-omics data [72] highlighted substantial differences in redox behavior in HCC tumors, relationships between redox genes, and genes that are liver specific, involved in the immune system, or inflammation. Different redox responses in subgroups of HCC tumors and their relationships with clinical outcome were also identified. We finally identified several novel markers of survival and prognosis and showed that they may be targeted for the treatment of HCC in specific patient groups, depending on the redox response.