Being a global health challenge, the incidence of liver cancer is continuously growing worldwide. It is the third most common cause of cancer-related deaths with 758,725 deaths and ranks sixth in terms of incidence rate with 866,136 new cases in 2022 [1]. Hepatocellular carcinoma (HCC) is the most common form of liver cancer, accounting for nearly 90 % of all cases [2]. Due to late diagnosis of the disease, the median survival of patients with advanced HCC is only approximately 6 to 20 months [3]. In the United States, the current 5-year survival is only 10 % [4]. While the treatment of many other types of tumors has progressed, the five-year survival rate for advanced HCC has not been significantly improved during the last decades [5].
The treatment for HCC is determined based on tumor stage and the expected benefits of major interventions, following the Barcelona Clinic Liver Cancer (BCLC) staging system. Generally, patients with early-stage HCC tumors are recommended for resection, transplantation and local ablation, whereas those with intermediate stages are the first candidates for transarterial chemoembolization (TACE) [6]. For those with advanced disease, patients first receive systemic therapies. To preclude relapse, adjuvant therapies are applied, but remain unmet medical need since randomized controlled trials yielded negative results so far. In fact, as many as 70 % of these patients develop tumor recurrence after five years [7–9].
Systemic therapies for advanced HCC utilize either immune-checkpoint inhibitors (ICIs), tyrosine kinase inhibitors (TKIs), or monoclonal antibodies. A substantial improvement in the development of systemic therapies recently emerged, with studies reporting an increase in overall survival and in the quality of life of patients [10]. While other regimens also showed survival benefits, including regorafenib, cabozantinib, and ramucirumab, the TKIs sorafenib and lenvatinib remain as the most effective single-drug therapies in the first-line setting [11,12].
Lenvatinib has become an important player in HCC therapeutics in 2018, overcoming the decades-long status of sorafenib as a favored treatment. Lenvatinib, an oral small molecule multi-kinase inhibitor similar to sorafenib, was approved for patients with unresectable HCC. It exerts its anticancer properties by inhibiting vascular endothelial growth factor receptor 1-3 (VEGFR1-3), fibroblast growth factor receptors 1-4 (FGFR 1-4), platelet-derived growth factor receptor α (PDGFR-α), RET, and KIT protooncogenes [13,14]. It inhibits angiogenesis, cell proliferation, and tumor growth. In the landmark phase III non-inferiority REFLECT trial comparing lenvatinib and sorafenib, lenvatinib showed non-inferiority in overall survival (OS) with duration of 13.6 vs. 12.3 months (hazard ratio, 0.92; 95 % CI, 0.79-1.06) [12,15]. Significant improvements on the median progression-free survival (PFS), median time-to-progression, and overall response rate (ORR) were seen in patients that received lenvatinib. Based on the promising results of this trial, the Food and Drug Administration (FDA) approved lenvatinib as first-line treatment of patients with unresectable HCC in 2018.
But despite lenvatinib being a key player in HCC therapy, its efficacy still remains poor with the overall median OS of ∼1 year and the ORR of approximately 40 % [12,16]. This poor prognosis is mainly attributed to drug resistance, which lowers the clinical benefits of lenvatinib. Over 60 % of HCC patients develop resistance to lenvatinib within 1 year, decreasing its clinical utility; and in fact only a fraction of patients obtains long-term benefits [17]. Therefore, a better understanding of the molecular mechanism of lenvatinib resistance could offer new insights to overcome resistance and increase its clinical benefits. As such, this study examined public datasets to mine novel gene players that drive lenvatinib resistance. Differentially expressed genes (DEGs) were processed for functional analysis and the most promising candidates were validated in in vitro cell models.
2Materials and Methods2.1Dataset search and detection of differentially expressed genes (DEGs)Four public datasets (GSE223201, GSE186191, GSE214324, and GSE211850) were retrieved from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (Table 1). These datasets compare lenvatinib-resistant cells to control parental cells, as well as patients that received lenvatinib as neoadjuvant treatment after surgical resection. DEGs were extracted using GEO2R, where significance level cut-off was set at 0.05 and was adjusted by Benjamini & Hochberg (False Discovery Rate, FDR). The DEGs were then filtered by selecting those genes with a fold change (FC) of ≤1.5 which were considered as upregulated genes. Those genes with FC of ≥0.5 were considered downregulated genes.
2.2Combination of DEGsVenn diagrams were drawn to determine overlapping genes in the four GEO datasets using Bioinformatics & Evolutionary Genomics (https://bioinformatics.psb.ugent.be/webtools/Venn/). To avoid missing critical genes, DEGs that overlapped in at least two datasets were selected for further analysis.
2.3Functional analysis of DEGsFor enrichment analysis, the two set of DEGs (FC≥1.5 and ≤0.5) were subjected to Gene Ontology (GO) and pathway enrichment analyses to explore the biological functions and pathways associated to lenvatinib resistance in hepatocellular carcinoma. GO enrichment analyses were performed using ShinyGO 0.80 (http://bioinformatics.sdstate.edu/go/), categorizing DEGs into Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Pathway Enrichment analyses were conducted using KEGG (Kyoto Encyclopedia of Genes and Genomes) and WikiPathways databases to identify significantly enriched pathways. Both GO terms and pathways with a p-value of <0.05 were considered statistically significant.
2.4Survival analysis of HCC patientsThe prognostic relevance of the candidate gene drivers of lenvatinib resistance on overall survival was evaluated using The Human Protein Atlas (https://www.proteinatlas.org/) that uses the Cancer Genome Atlas for Liver Hepatocellular Carcinoma (TCGA LIHC) cohort (https://portal.gdc.cancer.gov/) [18]. Patients were classified into two groups (high and low expression) based on the best expression cut-off, which refers to the FPKM value that yields maximal difference with regards to survival at the lowest log-rank P-value.
2.5Cell cultureHuman liver cell line Huh7 and Hep3B were cultured in DMEM (high glucose) and MEM media, respectively, supplemented with 10 % fetal bovine serum (FBS), 1 % penicillin-streptomycin, and 1 % L-glutamine. Cells were maintained at 37°C in a humidified incubator with 5 % CO2. For lenvatinib treatment, 2.0 × 105 cells were seeded in a 6-well plate. At 24 hours after seeding, cells were treated with lenvatinib (HY-10981, MedChemExpress, USA) at 1 µM and 5 µM doses for 48 hours. All experiments were performed in three biological replicates.
2.6Reverse transcription - quantitative polymerase chain reaction (RT-qPCR)Total RNA was extracted from cells using TriFast II (EMR517100, EuroClone, Italy) and was quantified using a Nanodrop 2000c Spectrophotometer (Thermo Scientific, Waltham, MA, USA). cDNA were obtained by RNA reverse transcription using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Waltham, MA, USA) following the manufacturer’s instructions. Quantitative PCR was conducted in CFX 9600 real-time PCR system using SYBR Green Supermix protocol (Bio-Rad Laboratories, Hercules, CA, USA) and relative mRNA expression was calculated using CFX Maestro Software 2.0 (Bio-Rad Laboratories) normalizing to the expressions of 18S and β-actin. The primer sequences are listed in Table 2.
Primer sequences for RT-qPCR.
Statistical significance was calculated using GraphPad Prism version 8.0 (GraphPad Software, Inc., La Jolla, CA, USA). The in vitro data were obtained from at least three biological replicates and were expressed as mean ± SD. Statistical significance was set to p-value <0.05 and reported as indicated here: * p < 0.033, ** p < 0.002, and *** p < 0.001.
2.8Ethical StatementsThe study did not involve human or animal specimens, and therefore, ethical approval was not applicable.
3Results3.1Screening of differentially expressed genes (DEGs) associated with lenvatinib resistanceIn this study, four datasets were analyzed to acquire a novel gene signature that might be associated with lenvatinib resistance in HCC (Table 1). GSE223201 reported expression profiles of five patients that underwent surgical resection after neoadjuvant lenvatinib treatment and 10 matched patients. The matched patients were based on age, sex, tumor size, and etiology. GSE186191 compared the transcriptional profile of parental (Huh7 and Hep3B cells) against lenvatinib-resistant cells which were exposed to lenvatinib for 6 months. GSE214324 measured gene expression profiles of MHCC-97H cells treated with lenvatinib at different concentrations of 40, 80, and 100 µM. GSE211850 compared the gene expression profiles of parental Huh7 cells and lenvatinib-resistant cells which received lenvatinib for 10 months.
A total of 110, 8844, 9346, 166, and 8438 DEGs were acquired from the four datasets (Fig. 1A, Table 3). Of these genes, DEGs with a FC ≥1.5 were considered upregulated genes, whereas those with FC ≤0.5 were considered downregulated genes. These genes in each dataset were selected for further analysis. Venn diagrams were drawn for the set of genes with FC ≥1.5 and FC ≤0.5 to identify common genes in all the four datasets (Fig. 1B). KIF26B was upregulated in 4 of 5 datasets, whereas TNFSF10, PRSS22, NYNRIN, AGPAT4, UGT2B11, and C6 were upregulated in 3 datasets. NPC1L1 and EBNA1BP2 were downregulated in 4 of the 5 datasets.
Differential Gene Expression in four public datasets. A. Volcano plots illustrating the differential expression of genes in four datasets (GSE223201, GSE186191, GSE214324, and GSE118850). The x-axis represents the log2 fold change in gene expression between the two conditions, while the y-axis denotes the -log10 p-value, which indicates the statistical significance of the observed differences. B. Venn diagrams depicting common differentially expressed genes (DEGs) across the four datasets with fold change of ≥1.5. C. Venn diagrams depicting common DEGs across the four datasets with fold change of ≤0.5.
Number of differentially expressed genes in the four datasets.
| Code | Dataset | DEGs | FC ≥1.5 | FC ≤0.5 |
|---|---|---|---|---|
| S1 | GSE223201 | 110 | 37 | 68 |
| S2 Huh7 | GSE186191 Huh7 | 8844 | 515 | 6643 |
| S2 Hep3B | GSE186191 Hep3B | 9346 | 897 | 6346 |
| S3 | GSE214324 | 166 | 24 | 86 |
| S4 | GSE211850 | 8438 | 643 | 5804 |
DEGs, differentially expressed genes; FC, fold change.
GO Enrichment Analysis was performed using ShinyGO 0.80 (http://bioinformatics.sdstate.edu/go/) and significant enrichment in several biological processes, cellular components, and molecular functions were found in the two gene sets in the four datasets. The upregulated genes (FC ≥1.5) were enriched in BP terms such as system development, regulation of biological quality, and cell differentiation. In CC terms, extracellular matrix, external encapsulating structure, and intrinsic component of plasma membrane were the most enriched. For MF terms, signaling receptor binding, calcium ion binding, and lipid binding were the most enriched (Fig. 2).
Gene Ontology (GO) Enrichment Analysis of upregulated (FC≥1.5) and downregulated genes (FC≤0.5). The graphs are categorized into Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The x-axis shows the -log10 of the p-value, indicating the significance of the enrichment for each GO term. The color of the bars corresponds to fold enrichment, while the size of the circle represents the number of genes.
For downregulated genes (FC ≤0.5), RNA processing, organonitrogen compound biosynthesis, and organelle organization were most enriched in BP. For CC, nuclear lumen, nucleoplasm, and nuclear protein-containing complex were the most represented while for MF, RNA binding, nucleic acid binding, and enzyme binding were the most significantly enriched (Fig. 2).
Notably, pathway analysis using KEGG and WikiPathways identified steroid hormone biosynthesis, retinol metabolism, and chemical carcinogenesis-DNA adducts, metapathway biotransformation, and glucuronidation were highly represented in upregulated genes, suggesting the potential activation of these pathways in lenvatinib resistance (Fig. 3). For downregulated genes, metabolic pathways and VEGFA-VEGFR2 signaling were the most represented, indicating that the dysregulation of these pathways may be the driving force of lenvatinib resistance in HCC (Fig. 3).
KEGG and WikiPathways Enrichment Analysis of Differentially Expressed Genes. The bar plots depict the significantly enriched pathways identified through KEGG and WikiPathways analysis for upregulated genes (fold change ≥1.5) and downregulated genes (fold change ≤0.5) in lenvatinib resistance. The KEGG and WikiPathways enrichment analysis of upregulated (in pink) and downregulated genes (in blue) shows the top pathways, with the x-axis representing the -log10 of the p-value, indicating the significance of enrichment, and the y-axis listing the enriched pathways.
The enrichment of certain GO terms and pathways underscores the relevance of these biological processes and pathways in the context of lenvatinib resistance in HCC. These insights provide a foundation for further exploration of the molecular mechanisms driving the HCC tumors to develop resistance.
3.3Expression of DEGs in the cancer genome atlas (TCGA)To obtain a more intuitive insight into the several target genes that contributes in lenvatinib resistance, DEGs were analysed to identify dysregulated genes present in the different datasets. Those genes which were differentially expressed in at least two datasets were considered, which narrowed down the candidates to 24 upregulated and 47 downregulated genes (Supplementary Fig. S1 and S2). Of these dysregulated genes, 12 genes were consistently dysregulated in liver HCC from TCGA LIHC and GTEx cohorts. SECTM1, SEZ6L2, FBLN7, IFI6, UGT2B22, and AKR1C1 were upregulated both in lenvatinib resistance and in HCC tumor tissues, indicating that the expressions of these genes were enhanced during lenvatinib resistance (Fig. 4A). More so, APOA4, ATOH8, CPEB3, IGFBP3, NPC1L1, and PALM3 were downregulated in lenvatinib resistance and in HCC tumor (Fig. 4B). Of these genes, five genes emerged to be good candidates for further study because they were not previously studied in HCC. These candidate genes were SEZ6L2, SECTM1, FBLN7, IFI6, and NPC1L1.
Gene expressions of lenvatinib-associated genes in TCGA Liver Hepatocellular Carcinoma (TCGA LIHC) cohort. Normalized expressions of candidate genes in hepatocellular carcinoma (HCC) tumor tissues (n=369) vs normal liver tissues (n=160) from TCGA and GTEx datasets, accessed through GEPIA2 database (http://gepia2.cancer-pku.cn/#index). A. Candidate genes with consistent high expressions in TCGA LIHC cohort. B. Candidate genes with consistent low expressions in TCGA LIHC cohort. Statistical significance *p-value<0.05.
To acquire more in depth information about these five genes in lenvatinib resistance, their TPM normalized expression across the datasets were acquired (Fig. 5). SEZ6L2, SECTM1, FBLN7, and IFI6 were significantly upregulated in HCC tissue and in lenvatinib-treated patient and cells. The low expression of NPC1L1 was seen both in HCC tumor tissues and in lenvatinib-resistant cells.
Expressions of five candidate genes in public datasets. TPM (transcript per million) normalized expression of SEZ6L2, SECTM1, IFI6, NPC1L1, and FBLN7 in public datasets. LENV, lenvatinib; LENV-R, lenvatinib-resistant; Student t-test, *p < 0.033, **p < 0.02, ***p < 0.01. Data was retrieved from GSE223201, GSE186191, GSE214324, and GSE211850 from the Gene Expression Omnibus.
As in vitro validation, we treated Huh7 and Hep3B cells with lenvatinib at two doses for 48 hours. Our own results show increased expressions of SECTM1 and IFI6 in Huh7 cells, which is consistent with the public datasets (Fig. 6). In Hep3B cells, all resistance genes increased after lenvatinib treatment. The raw data for Figs. 5 and 6 are available in Supplementary Tables T1 and T2, respectively.
Effect of lenvatinib in the expression of candidate genes in HCC cell lines. Huh7 and Hep3B cells were treated with low (1 µM) and high (5 µM) dose of lenvatinib for 48 hours. RNA expression levels of target genes were normalized to 18S and β-actin. Statistical analysis: one-way ANOVA with post-hoc test, *p < 0.033, **p < 0.02, ***p < 0.01 against DMSO-treated cells. Data is presented from three biological replicates.
The survival curves were obtained from The Human Protein Atlas that utilizes TCGA data (Fig. 7). For overall survival, poor prognosis was seen in patients with high expressions of SECTM1, SEZ6L2, IFI6, and NPC1L1, whereas patients with low FBLN7 expression had poorer prognosis, even though it was not significant (Fig. 7). In lenvatinib resistant HCC cells, SECTM1, SEZ6L2, and IFI6 were also upregulated in HCC tumor, indicating that the enhanced expression drives tumors to develop resistance, possibly resulting to even poorer prognosis. NPC1L1 expression is lower in HCC tumor tissues, but patients with lower NPC1L1 expression tend to have better 5-year prognosis (58 % versus 45 %). More so, NPC1L1 expression tends to be lower in lenvatinib-resistant patient and HCC cells.
Survival probability based on the five candidate genes of liver cancer patients. The Kaplan-Meier overall survival (OS) curves illustrate the survival probability of liver cancer patients stratified by high (pink line) and low (blue line) expression levels of the five candidate genes. The OS curves were obtained from The Human Protein Atlas (https://www.proteinatlas.org/).
Collectively, these results suggest that these genes could both play a pivotal role in HCC development and progression by modulating certain signaling pathways, affecting the overall survival probability. However, their molecular mechanisms should further be evaluated using in vitro and in vivo assays, to fully discriminate their functions in lenvatinib resistance.
4DiscussionPrimary resistance to the initial drug treatment is highly associated with genetic heterogeneity of HCC tumors [19]. This acquired resistance remains a major challenge for patients with HCC since it lowers the overall clinical utility of the drug. Over a certain of period, drug resistance can either manifest or be enhanced after exposure to the same pharmacological treatment [20]. In fact, majority of HCC patients do not obtain long-term benefit from systemic therapy due to acquired drug resistance [21]. The mechanism of drug resistance is complex and can be due to a number of events, including alterations in cell signaling pathways, dysregulations of apoptosis and other cell death programs, the heterogeneous nature of the tumor microenvironment, cancer stem cell renewal, tumor hypoxia, alterations in drug metabolism, drug efflux and uptake, DNA repair, and epigenetic control [21]. Overall, chemoresistant tumor cells possess a survival advantage and thus create a tumor whose genomic makeup confers resistance to drug therapy [22].
Lenvatinib, an oral multi-kinase inhibitor, was approved for patients with unresectable HCC. However, HCC tumors could develop lenvatinib resistance, limiting its long-term efficacy. This could be due to dysregulations of certain signalling pathways, mutations in key genes, epithelial-mesenchymal transition, or alteration of tumor microenvironment [23–25].
In this study, we found a number of genes which were differentially expressed in lenvatinib-resistant patients and cells. Our analyses have found that certain GO terms and KEGG and WikiPathways were more represented in upregulated and downregulated genes in lenvatinib-resistance group. For instance, VEGFA-VEGFR2 signalling was evidently enriched in downregulated genes, which was not surprising since lenvatinib inhibits the kinase activities of VEGF receptors. After applying several filtering conditions on the datasets, 12 genes were significantly altered in HCC tumor tissues based on TCGA and GTEx datasets, whose expression patterns were similar in lenvatinib resistance. Of these genes, five were selected because their roles in lenvatinib resistance have not been fully elucidated until now. These genes were SEZ6L2, SECTM1, FBLN7, IFI6, and NPC1L1.
SEZ6L2 (seizure related 6 homolog like 2) is a type 1 transmembrane protein that is primarily expressed in the brain and is affiliated with the seizure-related gene 6 (SEZ6) family. The overexpression of SEZ6L2 can drive tumor progression in CRC patients [26] and is correlated to tumor size, TNM stage, tumor number, and poor prognosis in HCC [27]. The high expression of SEZ6L2 has been found to be an independent prognostic indicator of thyroid carcinoma [28].
SECTM1 (secreted and transmembrane 1) is highly expressed in various normal cells, including epithelial cells, dendritic cells, and neutrophils [29], but also in many tumors, including melanoma, breast, prostate, and myeloid leukemias [30]. SECTM1 is proposed as a ligand for CD7 and promotes T cell proliferation [31]. In glioblastoma, SECTM1 promotes cell proliferation, migration, and invasion by promoting EMT through TGFβ1/Smad signalling pathways [32].
FBLN7 (fibulin 7) is a member of fibulin protein family, a group of cell-secreted glycoproteins that functions as a cell adhesion molecule and interacts with other extracellular matrix proteins as well as cell receptors [33]. It is overexpressed in glioblastoma tissue among astrocytic tumors and binds to angiopoietin-1 through interaction between the N-terminal portions of FBN7 and angiopoietin-1. This binding inhibits the Ang1-Tie2 interaction and the subsequent phosphorylation of the Tie2 receptor. This suggests that FBLN7 may contribute to the aberrant vessel formation via modulation of the angiopoietin-1/angiopoietin-2-Tie2 axis in glioblastoma [34]. Its close relative fibulin-5 (FBLN5) can inhibit HCC cell migration and invasion by downregulating matrix metalloproteinase-7 (MPP7) expression [35].
IFI6 (interferon alpha inducible protein 6) plays important roles in various biological processes including antiviral response, apoptosis regulation, and tumor suppression [36,37]. It has been extensively studied in various cancer types including lung adenocarcinoma, myeloma, pancreatic, breast, gastric, prostate, and ovarian cancers [38–42]. In colorectal cancer, downregulation of IFI6 could reverse oxaliplatin resistance by activating ROS-induced p38MAPK signalling pathway [43].
NPC1L1 (Niemann-Pick C1-like 1) is a protein that is essential for intestinal cholesterol absorption and plays vital roles in dietary cholesterol absorption and biliary cholesterol resorption. In the small intestine and liver, NPC1L1 primarily acts as a sterol transporter protein that controls lipid homeostasis [44]. In colorectal cancer, high NPC1L1 expression is associated with disease development and poor prognosis, and has value as a standalone prognostic factor [45]. In pancreatic ductal adenocarcinoma, NPC1L1 inhibition using ezetimibe significantly decreases cell viability and tumor volume [46]. In HCC, NPC1L1 is not highly expressed in HCC liver tissues as in the peritumoral liver tissues, both at protein and mRNA levels [47].
5ConclusionsPrimary resistance remains a major challenge for HCC patients. Here, we highlighted new genes (SECTM1, SEZ6L2, FBLN7, IFI6, and NPC1L1) that were involved in lenvatinib resistance in patient and in vitro models. Alteration in the various signaling pathways should be explored since their dysregulation could drive lenvatinib resistance. We would recommend further investigations to elucidate the specific mechanisms of the candidate genes, especially SECTM1 and IFI6, in lenvatinib resistance and to describe the functional outcomes of these genes in terms of tumor growth, therapy response, and prognostic relevance.
FundingCD is funded by a fellowship of the Department of Science and Technology and the Philippine Council for Health Research and Development (DOST-PCHRD).
Data availability statementThe data that support the findings of this study are available on reasonable request from the corresponding author.
Author contributionsConceptualization: CD and CS. Data Curation: CD. Formal Analysis: CD. Investigation: CD. Project Administration: CS and CT. Supervision: CS and CT. Writing-original draft: CD. Writing-review & editing: CD, CS and CT
CD would like to thank Department of Science and Technology-Philippine Council for Health Research and Development (DOST-PCHRD). CS is supported by the Fondazione Veronesi, Milan, Italy.
