A Three-Gene Signature for Predicting the Prognosis of Patients Treated with Transarterial Chemoembolization (TACE) and Identification of PD-184352 as a Potential Drug to Reverse Nonresponse to TACE

Background Transarterial chemoembolization (TACE) is a first-line treatment for patients with unresectable hepatocellular carcinoma (HCC). Owing to differences in its efficacy across individuals, determining the indicators of patient response to TACE and finding approaches to reversing nonresponse thereto are necessary. Methods Transcriptome data were obtained from the GSE104580 dataset, in which patients were marked as having TACE response or nonresponse. We identified differentially expressed genes (DEGs) and performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. We screened genes with a prognostic value for TACE in the HIF-1 signaling pathway by univariate regression analysis. By using least absolute shrinkage and selection operator (LASSO) Cox regression, we established a multigene signature in GSE14520, which we verified using a drug sensitivity test. The Connectivity Map (CMap) database was used to find potential drugs to reverse nonresponse to TACE. Results We constructed a prognostic signature consisting of three genes (erythropoietin (EPO), heme oxygenase 1 (HMOX1), and serine protease inhibitor 1 (SERPINE1)) that we validated by drug sensitivity test. After dividing patients treated with TACE into high- and low-risk groups based on this new signature, we showed that overall survival (OS) of the high-risk group was significantly lower than that of the low-risk group and that the risk score was an independent predictor of OS in patients treated with TACE. Based on our CMap findings, we speculated that PD-184352, an inhibitor of mitogen-activated protein kinase (MEK), had potential as a drug treatment to reverse nonresponse to TACE. We confirmed this speculation by using PD-184352 in a cell promotion experiment in a TACE environment. Conclusion We constructed a TACE-specific three-gene signature that could be used to predict HCC patients' responses to and prognosis after TACE treatment. PD-184352 might have potential as a drug to improve TACE efficacy.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common malignant tumors, with the fifth highest incidence rate and the third highest mortality rate [1]. More than two thirds of HCC patients cannot undergo surgical resection [2]. According to the Barcelona clinic liver cancer (BCLC) staging system, transarterial chemoembolization (TACE) is recommended as first-line treatment for midstage or advanced HCC [3,4]. Although TACE can improve outcomes in some patients, the objective remission rate after the first TACE treatment is less than 40% [5]. Therefore, it is particularly important to reveal the causes of TACE nonresponse and distinguish patients who cannot benefit from TACE before surgery.
Professor Wang [6] of the American Cancer Research Center believes that hypoxia is a potential drug resistance mechanism of TACE nonresponse. In TACE nonresponse tumors, hypoxia reprogramming may have been experienced before TACE treatment or it may lead to enhanced response to hypoxia caused by TACE. The HIF (hypoxiainducible factor) family is significantly upregulated in a hypoxic environment, abnormal activation of the HIF-1 signaling pathway plays an important role in tumor progression, HIF-dependent gene expression controls epithelial mesenchymal transformation (EMT), invasion, migration, and angiogenesis [7]. TACE causes the hypoxic microenvironment of the tumor, and researchers have shown that the expression of HIF-1α is significantly increased after TACE and is related to the poor prognosis of TACE [8,9], so the current research on TACE nonresponse mainly focuses on the abnormal activation of the HIF-1 signaling pathway. Our study also found that the HIF-1 signaling pathway is enriched in TACE nonresponse, which is why we decided to further study the HIF-1 signaling pathway.
Due to the heterogeneity of HCC, patient response to TACE is variable. Therefore, there is continuous interest in preoperative screening of those who could benefit from this procedure [10]. In recent years, scoring systems have been designed based on certain routine examination items to predict prognosis in TACE, such as the hepatoma arterial embolization (HAP) prognostic score, which applies simple laboratory and imaging indicators [11]. The HAP score has the advantages of ease of application and simplicity; however, because its calculation is based only on certain macroscopic data, it cannot reflect deeper indicators of prognosis after TACE treatment, which is a disadvantage in the era of precision medicine. Furthermore, the most current scoring systems are limited by their lack of specificity to TACE. Therefore, developing a method to assess whether patients would benefit from TACE before they are ready for such treatment could improve survival in certain patients.
It is important to clarify why some patients do not benefit significantly from TACE. In this study, using a Gene Expression Omnibus (GEO) dataset, we established a signature consisting of three genes (EPO, HMOX1, and SER-PINE1) to predict the prognosis in TACE. We divided cell lines into TACE response and nonresponse groups accord-ing to the gene signature and then evaluated the differences between the two groups to verify the signature's effectiveness. Using the Connectivity Map (CMap) database, we screened PD-184352, an inhibitor of mitogen-activated protein kinase (MEK), and proved that it could improve the sensitivity of HCC cells and the effects of TACE treatment in vitro. We believe that our initial exploration of a new intraoperative or postoperative administration pattern for TACE is worthy of further clinical trials.

Patient Cohorts and Cell Lines.
We screened genes in the HIF-1 signaling pathway from the GSE104580 cohort. This dataset contains data from 147 patients treated with TACE for whom TACE was the primary treatment. Of these 147 patients, 81 responded well (TACE response) and 66 responded poorly (TACE nonresponse) [12]. Ribonucleic acid (RNA) extracted from tumor biopsies of HCC patients to be treated with TACE was analyzed using Affymetrix gene arrays (Affymetrix, Santa Clara, CA, USA).
From the GSE14520 development cohort, we selected 74 patients with liver cancer who had received adjuvant TACE after liver resection and 30 who had received TACE after recurrence. In addition, we included 85 patients who had undergone liver resection alone. The details of this cohort were reported by Fako et al. [6].
Transcriptomic data of liver cancer cell lines were obtained from the Cancer Cell Line Encyclopedia (CCLE). We divided cell lines into responder-like and nonresponder-like groups by the calculating risk score based on the transcriptomic data described above.

Establishment of a Potential Prognostic HIF-1 Signaling
Pathway-Related Gene Signature. Gene expression data were normalized by formula log 2 using the limma package [13] in R software v4.1.0 (R Foundation for Statistical Computing, Vienna, Austria) to calculate differentially expressed genes (DEGs) between the TACE response and nonresponse groups in GSE104580 under the following conditions: Benjamini-Hochberg adjust P < 0:05, jlog 2 fold change ðFCÞj > 0:5, and P < 0:05. We then performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the dysregulated genes in the TACE nonresponse group using the clusterprofiler package in R [14]. Seventeen genes enriched in the HIF-1 signaling pathway were selected and then analyzed by univariate regression in GSE14520. Finally, we screened eight survival-related genes: aldolase, fructose-bisphosphate C (ALDOC); egl-9 family hypoxia-inducible factor 3 (EGLN3); enolase 2 (ENO2); EPO; hexokinase 2 (HK2); HMOX1; 6phosphofructo-2-kinase/fructose-2,6-biphosphatase (PFKFB3); and SERPINE1. Overall survival (OS) and recurrence-free survival (RFS) were calculated using the Kaplan-Meier (K-M) method; the survival package was used for statistical analysis and survminer package was used for visualization. With the glmnet package, we performed least absolute shrinkage and selection operator (LASSO) regression and established a prognostic signature to compress some regression coefficients, preserving the benefits  Journal of Oncology of subset shrinkage [15,16]. Subsequently, we established the following formula to calculate the risk score: where n represents the number of genes in the signature and β represents the coefficient of genes obtained from LASSO regression.  [17]. We input significantly upregulated genes of the high-risk group in 3 Journal of Oncology GSE14520 into the L1000 platform of CMap to find compounds that might cause reversal reactions. Subsequently, we took compounds with connectivity scores of less than −98, which could cause the opposite reaction of the Hep G2 cell line, as drugs that could potentially reverse nonresponse to TACE.    Journal of Oncology  1  2  2  2  2  2  3  3  3  3  3  3  3  3  3  3  3  4  5  5  5  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  7  7  7  7  7  7  7  7  7 Partial likelihood deviance 2.9. Statistical Analysis. All statistical tests were two tailed, and we considered P < 0:05 statistically significant. Univariate and multivariate Cox regression models were applied to confirm characteristics that were valuable in predicting the prognosis of patients treated with TACE. We used an independent t-test to assess differences between the two groups. A correlation matrix was constructed with ggplot2.

Results
3.1. The HIF-1 Signaling Pathway Was Enriched in the TACE Nonresponse Group. By analyzing differences in gene expression between 81 TACE nonresponders and 66 TACE responders in GSE104580, we identified 735 genes that were upregulated in nonresponders and 749 that were upregulated in responders (jlog 2 FCj > 0:5, P < 0:05). Based on these genes, we drew a heat induction map (Figure 1(a)) and a volcano map (Figure 1(b)). Then, we performed KEGG enrichment analysis on all 1484 genes and found that the HIF-1 signaling pathway was significantly enriched in TACE nonresponders (Figure 1(c)) but not in TACE responders (Figure 1(d)). Since cutting off the blood supply is an important way for TACE to kill cancer cells, the hypoxia microenvironment is one of the most important features after TACE. The hypoxia microenvironment is an important factor in the activation of the HIF-1 signaling pathway. More importantly, many studies have proved that HIF-1α is closely related to poor prognosis in TACE, so we chose the HIF-1 signaling pathway as our next research focus.

3.3.
Construction of the Prognostic Signature. Next, we performed LASSO regression analysis to screen the appropriate genes from the abovementioned eight in order to determine their prognostic characteristics (Figures 3(a) and 3(b)). The following formula was used to calculate the risk score: Risk score = SERPINE1 × 0:252860249076791 + HMOX1 × 0:0830953460565714 We divided patients into high-risk (n = 52) and low-risk (n = 52) groups based on the median risk score (Figure 3(c)). K-M curve analysis showed that the high-risk group had a poorer prognosis than the low-risk group (Figure 3 9 Journal of Oncology reached 0.826, 0.794, and 0.681 at 1, 3, and 5 years, respectively, indicating that the prognostic signature showed excellent predictive ability (Figure 3(e)). Then, we used the same formula to calculate the risk scores of samples in the GSE104580 dataset in order to detect any differences between TACE nonresponders and responders. As shown in Figure 3(f), the risk score in the TACE nonresponse group was higher than that in the TACE response group (P < 0:0001; 95% confidence interval (CI), 0.5971-1.235).

Our Signature Was Specific to TACE.
To verify whether this signature was specific to patients treated with TACE rather than applicable to all HCC patients, we divided patients treated with postrecurrence TACE, adjuvant TACE, and resection only into high-and low-risk groups by median risk score in GSE14520. One hundred four patients were treated with TACE, 74 with adjuvant TACE, and 30 with postrecurrence TACE. K-M curve analysis showed that the high-risk group had a poorer OS than the low-risk group among patients treated with postrecurrence TACE (Figure 4(a)) and adjuvant TACE (Figure 4(b)). However, we saw no difference in patients treated with resection only (Figure 4(c)). Then, we analyzed recurrence-free survival (RFS), except for patients treated with postrecurrence TACE, because the recurrence time is measured from resection to recurrence. As expected, the high-risk group had poorer RFS than the low-risk group among patients treated with TACE (Figure 4(d)) and adjuvant TACE (Figure 4(e)). There was still no difference in patients treated with resection only (Figure 4(f)). These results showed that our signature was specific to TACE.

Verification of the Signature's Distinguishing Ability in
Different Clinical Subgroups. Next, we verified the signature's distinguishing ability in different clinical subgroups. To this end, we divided patients treated with TACE into subgroups according to their clinical characteristics (age ≤ 50 and >50; tumor size ≤ 5 and >5 cm; AFP > 300 ng/mL; BCLC stage A; cirrhosis; TNM stages I, II, and III). K-M curves for OS after such division showed that significant differences remained in all subgroups except for TNM III, suggesting that this signature also had a good ability to distinguish between different clinical subgroups ( Figure 5(a)-5(k)). Because we enrolled too few patients with no cirrhosis, BCLC stage B, or BCLC stage C, we did not analyze OS for these patients.  3.6. Independent Prognostic Value of the Signature. To illustrate the clinical application value of this signature, we compared the risk scores of various clinical subgroups in GSE14520. That of patients with recurrence within 24 months was higher than that of patients with recurrence after 24 months (Figures 6(a) and 6(b)), which was consistent with our clinical observations. The risk score of patients with cirrhosis of the liver was higher than that of patients without (Figures 6(c) and 6(d)), and patients with tumor size > 5 cm had a higher risk score than those with tumor size ≤ 5 cm (Figures 6(e) and 6(f)), which suggested that TACE offered better prognosis in patients with tumor size ≤ 5 cm and those without cirrhotic livers. In addition, we found remarkable differences among different TNM grades: risk score increased along with the TNM grade, indicating that the signature could effectively distinguish different TNM grades of HCC and therefore might be a potential biomarker of this cancer (Figures 6(g) and 6(h)). We performed univariate Cox analyses to explore the predictive value of a series of clinical metrics and risk score in GSE14520. As shown in Table 2, a larger main-tumor size (>5 cm), higher TNM stage (III), higher BCLC stage (stages B and C), recurrence months (≤ 24), and risk score are significant risk factors. We then performed multivariate Cox analysis including these four variables and found that the risk score remained an independent predictor (hazard ratio   13 Journal of Oncology the relative mRNA expression of EPO, HMOX1, and SER-PINE1 in Hep G2 is much higher than that in Hep 3B. To verify the effectiveness of the signature, we tested the IC 50 of lobaplatin, a platinum chemotherapeutic drug commonly used in TACE, under both normoxia and hypoxia to mimic an environment caused by TACE. We showed that the IC 50 of lobaplatin in Hep 3B was 1.861 (range, 1.153-2.57) μg/mL under normoxia (Figure 7(a)) and 1.643 (range, 1.017-2.268) μg/mL under hypoxia (Figure 7(b)), while the IC 50 of lobaplatin in Hep G2 was 3.575 (range, 2.718-4.431) μg/ mL under normoxia (Figure 7(c)) and 3.87 (range, 2.999-4.741) μg/mL under hypoxia (Figure 7(d)). The nonresponder-like cell line Hep G2 had a higher IC 50 than the responder-like cell line Hep 3B under both conditions, and the IC 50 of lobaplatin on HepG2 increased slightly under hypoxia. Taken together, these results showed that our risk score could distinguish between TACE responders and nonresponders.
3.8. PD-184352 Could Be Used to Reverse Nonresponse to TACE. To find a way to break nonresponse to TACE, we determined DEGs between the high-and low-risk groups. We found that 114 genes were significantly upregulated in the high-risk group (Supplementary file 1, jlog 2 FCj > 0:5, P < 0:05). We input them into CMap to identify compounds that could potentially reverse nonresponse to TACE. We paid particular attention to drugs that might elicit the opposite response in the human hepatoma cell line Hep G2, because in this study, we used Hep G2 as the TACE nonresponse-like cell line. A total of 69 chemicals with connectivity scores of less than −98 were detected (Supplementary file 2). We found that MEK inhibitors ranked first, indicating that they might reverse nonresponse to TACE (Figure 8(a)). Interestingly, MEK is just upstream of the HIF-1 signaling pathway (Figure 8(b)). Among MEK inhibitors, PD-184352 ranked second and it is easily available; therefore, we identified PD-184352 as a potential reverser of nonresponse to TACE. Subsequently, we detected the IC 50 of PD-184352 on Hep G2 in a hypoxic environment, which was about 0.43 μM (Figure 8(c)).
Finally, we used lobaplatin and 1% O 2 to simulate TACE in vitro and added PD-184352 to observe its effect on the efficacy of TACE. As shown in Figure 8(d), one-way analysis of variance (ANOVA) showed that cell proliferation of the responder-like cell line Hep G2 was significantly reduced by using either lobaplatin alone (P = 0:03) or PD-184352 alone (P < 0:0001). When PD-184352 and lobaplatin were combined, the inhibitory effect on cell proliferation was further enhanced (P < 0:0001).

Discussion
With the increasing emphasis on personalized treatment and the lower cost of genome sequencing, the era of precision medicine has arrived. Specifying treatment options and selecting drugs according to genome sequencing results have become a standardized treatment model. In this study, we established a prognostic signature specific to TACE using three genes according to the expression pattern of DEGs between TACE nonresponders and responders, which we verified in cell experiments. More importantly, we selected PD-184352, an inhibitor of MEK, as a drug that could potentially reverse nonresponse to TACE.
The difference between TACE and traditional treatment is that TACE embolizes the blood supply artery when chemotherapeutic drugs are administered. Therefore, the current research has focused on hypoxia-and angiogenesisrelated pathways so as to explore differences in TACE response among individuals. Some studies have found that levels of hypoxia-related biomarkers before TACE, such as vascular endothelial growth factor (VEGF) and HIF-1α, are negatively correlated with survival and prognosis [6,20,21]. HIF-1α has been shown to be a key molecule in nonresponse to TACE. In 2018, Liu et al. [22] found that knocking down HIF-1α enhanced the efficacy of transarterial embolization. Therefore, the HIF-1 signaling pathway, to which HIF-1α belongs, is worthy of in-depth research. In this study, we found that this pathway was significantly enriched in TACE nonresponders compared with TACE responders. Then, using univariate Cox regression analysis, we found that eight genes were related to OS probability in GSE14520. The final prognostic signature was composed of three genes: EPO, HMOX1, and SERPINE1. The EPO protein encoded by the EPO gene, mainly expressed in the liver and kidneys, is a hypoxia-reactive protein. A hypoxic environment increases the activity of HIF-1α, which binds to cis-acting deoxyribonucleic acid (DNA) hypoxia-response elements (HREs) to activate EPO transcription [23]. Much evidence has shown that high expression of EPO is associated with poor prognosis in HCC and breast cancer [23,24], but its role in nonresponse to TACE is still unknown. In view of the important role of EPO in hypoxia, we think that it deserves further studies. Heme oxygenase (HO) is the rate-limiting enzyme in heme metabolism [25]. HMOX1 is a downstream target of the transcription factor HIF-1αlimiting enzyme in heme metabolism [26]. The role of HMOX1 in cancer remains controversial, and further exploration is urgently needed. SERPINE1, also known as plasminogen activator inhibitor-1 (PAI-1), has long been considered an important factor in poor prognosis of cancer [27,28]. In 2015, Divella et al. [29] directly proved that high expression of SERPINE1 is related to poor prognosis in TACE. The functions of these genes, especially their effects on TACE nonresponse, also need further research.
In this study, we obtained the transcriptome sequencing data of liver cancer cell lines in the CCLE database. We treated Hep G2 as a TACE nonresponse-like cell line and Hep 3B as a TACE response-like cell line. Then, we used lobaplatin and 1% O 2 culture to simulate a TACE environment. We found that Hep G2 cells were more tolerant of the toxic effects of lobaplatin than Hep 3B cells in both hypoxic and normoxic environments. More interestingly, the IC 50 of Hep G2 was higher in the hypoxic than in the normoxic environment.
We queried CMap build 1.0 based on an L1000 assay using the genes that were significantly upregulated in the GSE14520 high-risk group and screened the MEK inhibitor PD-184352 on the L1000 platform as a drug that could 14 Journal of Oncology reverse nonresponse to TACE. MEK, also known as MAP2K, regulates HIF-1α through the MEK/ERK/HIF-1α axis and is an important regulatory factor of the HIF-1 signaling pathway [30,31]. PD-184352 could simultaneously inhibit the activities of MEK1 and MEK2. Previous studies have demonstrated that PD-184352 shows antitumor activity in thyroid cancer and lung adenocarcinoma [32,33], but its role in HCC has not been studied. We demonstrated that the combination of lobaplatin and PD-184352 inhibited cell line proliferation more than either treatment alone under hypoxic conditions; therefore, this combination could become a novel intraoperative or postoperative TACE medication strategy. Finally, because we chose only lobaplatin to mimic TACE in our in vitro model in this study, we cannot rule out that the use of different chemotherapeutic agents combined with MEK inhibition might have different effects.
Our future work will be devoted to explaining the principle underlying clinical application of MEK inhibitors to TACE and providing a theoretical basis for such application.

Conclusion
We constructed a TACE-specific three-gene signature that could be used to predict HCC patients' responses to and prognosis after TACE treatment. PD-184352 might have potential as a drug to improve TACE efficacy.

Data Availability
The transcriptome in this research is available in GEO database (https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE14520 and GSE104580. The cell line transcriptome in this research is available in CCLE database (https://sites.broadinstitute.org/ccle).