Skip to main content
Download PDF

Context dependent isoform specific PI3K inhibition confers drug resistance in hepatocellular carcinoma cells

BMC Cancer volume 22, Article number: 320 (2022) Cite this article

Abstract

Background

Targeted therapies for Primary liver cancer (HCC) is limited to the multi-kinase inhibitors, and not fully effective due to the resistance to these agents because of the heterogeneous molecular nature of HCC developed during chronic liver disease stages and cirrhosis. Although combinatorial therapy can increase the efficiency of targeted therapies through synergistic activities, isoform specific effects of the inhibitors are usually ignored. This study concentrated on PI3K/Akt/mTOR pathway and the differential combinatory bioactivities of isoform specific PI3K-α inhibitor (PIK-75) or PI3K-β inhibitor (TGX-221) with Sorafenib dependent on PTEN context.

Methods

The bioactivities of inhibitors on PTEN adequate Huh7 and deficient Mahlavu cells were investigated with real time cell growth, cell cycle and cell migration assays. Differentially expressed genes from RNA-Seq were identified by edgeR tool. Systems level network analysis of treatment specific pathways were performed with Prize Collecting Steiner Tree (PCST) on human interactome and enriched networks were visualized with Cytoscape platform.

Results

Our data from combinatory treatment of Sorafenib and PIK-75 and TGX-221 showed opposite effects; while PIK-75 displays synergistic effects on Huh7 cells leading to apoptotic cell death, Sorafenib with TGX-221 display antagonistic effects and significantly promotes cell growth in PTEN deficient Mahlavu cells. Signaling pathways were reconstructed and analyzed in-depth from RNA-Seq data to understand mechanism of differential synergistic or antagonistic effects of PI3K-α (PIK-75) and PI3K-β (TGX-221) inhibitors with Sorafenib. PCST allowed as to identify AOX1 and AGER as targets in PI3K/Akt/mTOR pathway for this combinatory effect. The siRNA knockdown of AOX1 and AGER significantly reduced cell proliferation in HCC cells.

Conclusions

Simultaneously constructed and analyzed differentially expressed cellular networks presented in this study, revealed distinct consequences of isoform specific PI3K inhibition in PTEN adequate and deficient liver cancer cells. We demonstrated the importance of context dependent and isoform specific PI3K/Akt/mTOR signaling inhibition in drug resistance during combination therapies. (https://github.com/cansyl/Isoform-spesific-PI3K-inhibitor-analysis).

Peer Review reports

Background

According to WHO-Global cancer observatory (GCO) that one-fifth of men and one-sixth of women will be diagnosed with cancer throughout their lives and one-eighth of men and one-eleventh of women will die of it worldwide before the age of 75 years. Hepatocellular cancer (HCC) which constitutes the 75% of Primary liver cancers is the 5th most common and the 3rd most lethal cancer in the world [1, 2]. While the death rates from other cancers are decreasing due to advances in diagnosis and therapeutics, the incidence and the mortality of HCC follow an increasing trend due to high rate of obesity associated liver diseases [3, 4].

Development of HCC is multi-factorial and complex biological process, where the chronic liver disease is initiated due to hepatic injury, followed by continuous inflammation and cell death, which in turn leads to the regeneration of hepatocytes and the increased rate of mutations along with genomic instability [5]. The increased number of proliferating cells evokes the activation of several cell signaling pathways involved in liver regeneration, such as growth factor signaling, cell differentiation, angiogenesis and cell survival. Stimulation of these pathways is mostly associated with tyrosine kinases which are usually the members of PI3K/Akt/mTOR cell signaling [6]. Studies show that the heterogeneous nature of HCC is mainly caused by the variations of mutations and alterations in expression levels of these key proteins [7].

Currently there is no effective therapy for patients suffering from HCC, the survival rates is only 7% for 5 years [1]. There are two FDA approved small molecule drug treatments for HCC; Sorafenib (Nexavar, BAY43–9006) and Regorafenib (Bayer, BAY73–4506), are receptor tyrosine kinase inhibitors targeting Raf, VEGFR and PDGFR kinases. They inhibit tumor cell proliferation and angiogenesis while promoting apoptosis. However, in most of the cases they are not capable of eliminating the cancer cells primarily because of the heterogeneous nature of HCC [8, 9]. Moreover, the signaling pathways involved in proliferation, growth, angiogenesis and metastasis are redundant, compensating each other through some key molecular regulations. Which makes them with superfluous functions due to the potential cross-talks between them, which could be another reason for the ineffectiveness of these two multi-kinase inhibitors [6].

The constitutive activation of PI3K/Akt/mTOR signaling pathway is frequently observed in liver cancer due to inactivating mutations or loss of heterozygosity in a tumor suppressor protein Phosphatase and tensin homolog (PTEN). PTEN dysfunction is observed in nearly 50% of the HCC cases and correlated with poor prognosis, drug resistance and low patient survival [10,11,12]. PTEN prevents the Akt activation by dephosphorylating PIP3, or mutations activating PIK3CA gene, or damage in the negative-feedback loop from mTOR signaling pathway in various epithelial cancers including HCC [13,14,15,16].

The influence of isoform diversity on responses to drugs with respect to large number of GPCR receptors has been demonstrated at systems level recently [17]. Furthermore, there are resent studies on the association of isoform specific differential involvement of AKT in the pathophysiology and therapeutic responses of cancer cells [18,19,20]. Here in this study, we focused on the response of HCC cells to isoform specific PI3K inhibitors. PI3Ks are grouped into three classes based on their structures [21, 22] but two of Class I members of PI3Ks have heterodimeric class IA p110-α (p110) and class IB p110-β (p85) regulatory subunits are well studied enzymes in cancer. PIK3CA gene encoded PI3K isoform p110-α, is activated through receptor tyrosine kinases (RTKs) and Ras oncogene. In cancer, signaling though PI3K predominantly depends on alpha isoform regulating cellular growth, metabolism and angiogenesis. The other PI3K isoform encoded by PIK3CB, p110-β is regulated mostly by G protein-coupled receptors (GPCRs) and has critical functions in inflammatory cells [23, 24].

In this study, we demonstrated that context (PTEN function) dependent isoform specific PI3K inhibition confers drug resistance by their antagonistic and synergistic effects with Sorafenib on HCC cells at network level in and studies focusing on the discovery of agents against HCC aim to identify target proteins that escape from regulatory signaling mechanisms of the cell.

Results

Molecular and cellular characterization Huh7 and Mahlavu cells in the presence of small molecule isoform specific PI3K inhibitors

Well-differentiated Huh7 cell line with adequate PTEN and poorly-differentiated PTEN deficient Mahlavu cells were selected to exploit throughout this study. The expression levels and the phosphorylation status of key proteins in PI3K/Akt/mTOR and RAF/MEK/ERK signaling pathways were reported by our group, and in correlation with their PTEN status, Mahlavu cells display hyper-activated cell survival proteins [25]. Initially Sorafenib, LY294002, PI3K inhibitor p110α subunit specific (PIK-75) and PI3K inhibitor p110β subunit specific (TGX-221) were analyzed for their cytotoxic bioactivity and their effect on cell cycle progression on Huh7 and Mahlavu cells (Fig. 1A). G1, S and G2/M cell cycle phases were analyzed separately to calculate viable cell distributions among them (Fig. 1B). Sub-G1 percentage demonstrating apoptotic cells were also calculated. Cell cycle distribution remained stable for both cell lines and all inhibitor treatments. In both cell lines, Sorafenib and PIK-75 treatments showed stimulation of apoptosis through increase in sub-G1 population. In Huh7, Sorafenib seems to be more active while PIK-75 functioned more in Mahlavu cells which was more aggressive than Huh7 cell line by PTEN-loss based hyper-active Akt stimulation.

Fig. 1
figure 1

Characterization of HCC cells in the presence of small molecules inhibitors. Real time cell growth analysis of Huh7 and Mahlavu cells with increasing concentrations (40 μM, 20 μM, 10 μM, 5 μM, 2.5 μM) of Sorafenib, PI3K inhibitor LY294002, PI3Ki-β inhibitor (TGX-22) and PI3Ki-α (1 μM, 0.5 μM, 0.25 μM, 0.125 μM, 0.0625 μM) PI3Ki-α (PIK-75) along with DMSO vehicle control (Control is black and increasing drug concentrations is given in grey level, highest concentration is being the darkest) (A). Cell cycle analysis with flow cytometry. Sub-G1 population represents apoptotic cells (B). Wound healing assay for 24 and 48 h for cell migration. (C). 10 μM of Sorafenib, LY294002 and PI3Ki-β (TGX-221) and 0.1 μM of PI3Ki-α (PIK-75) were used for cell cycle and migration assays

Full size image

Migration analysis of the inhibitors

In order to analyze the effects of selected inhibitors on cell migration, wound-healing assay was performed. The percentages of wound closures after 48 h of initial scratch were calculated for Huh7 and Mahlavu. We observed that Sorafenib and PIK-75 reduced migration significantly (p < 0.001) in both Huh7 and Mahlavu (Fig. 1C).

Synergistic cytotoxicity analysis

Since none of the treatments alone was fully effective in inhibiting growth and stimulate apoptosis, we addressed the value of co-treatments of Sorafenib with PIK-75 and TGX-221 through real-time cell growth analysis (Fig. 1C). A synergistic effect of Sorafenib and PIK-75 treatments was observed on growth of both cell lines. TGX-221 combinatory treatment with Sorafenib also resulted in synergistic growth inhibition on Huh7 cell line. On the other hand, TGX-221 displayed a growth inhibition of Mahlavu, TGX-221 co-treatment with Sorafenib resulted in an antagonistic effect and stimulated cellular growth. Furthermore, Sorafenib and PIK-75 treatment had more drastic effect on Mahlavu compared to Huh7. Therefore, these findings indicated that in PTEN deficient Mahlavu cells, constitutive activation of PI3K/Akt signaling mainly depends on p110-α (Fig. 2).

Fig. 2
figure 2

Real-time cell growth analysis. Human liver cancer cells Huh7 (A, B) and Mahlavu (MV) (C, D) were treated with the Sorafenib, PI3Ki-α and PI3Ki-β alone or in combination with increasing concentrations as indicated. Cell index measurements were obtained by RT-CES software. DMSO was used as negative control A B. 72 h of the percent growth inhibition values were used to calculate drug interactions with The SynergyFinder web application. Positive delta score reflects synergistic and negative score reflects antagonistic drug interactions. Experiments were performed in triplicate

Full size image

Network level analysis of isoform specific combinatory effects of PI3Ks

In order to describe the molecular events in the differential response of PTEN adequate (Huh7) and deficient (Mahlavu) cells displaying differential PI3K/Akt/mTOR pathway activities toward PI3K-α inhibitor (PIK-75) and PI3K-β inhibitor (TGX-221) alone or in combination with Sorafenib, we performed RNA sequencing experiments. Further, network-based data analysis using systems biology approaches which is represented in Fig. 3A is applied.

Fig. 3
figure 3

Systems Biology flowchart for RNA-seq Data analysis. Systems level methodology flowchart for differential PIK3K/Akt/mTOR pathway activities in Huh7 and Mahlavu calls treated with Sorafenib, PI3K-α inhibitor PIK-75 and PI3K-β inhibitor TGX-221 alone or in combination (A). Differentially Expressed Genes (DEG) Table summarize the abbreviations of samples as the treatments to HCC cells and differentially expressed gene (DEG) numbers. DEG filtration for A and B as follows; Huh7 cells:logFC ≥2.0, ≤ − 2.0 and p ≤ 0.01, Mahlavu cells: logFC ≥1.5, ≤ − 1.5 and p ≤ 0.01 (B). Pearson correlations for gene expressions in Huh7 and Mahlavu, no filtration (C). Dendrogram analysis on logFC for top 50 DEGs, red and blue color represented for up- and downregulated genes (D)

Full size image

Initially, we identified differentially expressed genes (DEGs) by using empirical analysis of digital gene expression data in R (edgeR tool) [26]. While PIK-75 (ALPHA) treated Huh7 cells had more upregulated genes, combination with Sorafenib reversed the number of genes in favor of gene downregulation. The response of Mahlavu to the same combination treatment (SALPHA) is significantly different with highly increased number of up and downregulated genes when compared to Sorafenib or PIK-75 alone. In both cell lines, TGX-221 treatment had minor action yet in PTEN deficient Mahlavu, single treatment of Sorafenib had a neutral effect. Interestingly, TGX-221 and Sorafenib combination resulted in higher number of downregulated genes in Mahlavu compared to the number of upregulated genes in Huh7 (Fig. 3B).

Using Pearson correlation analysis, we demonstrated the similarities in the overall gene expressions using the corresponding logFC values. A significant correlation was observed between Sorafenib alone treatment and its combinatory treatment with TGX-221 (0.92) in Huh7 and (0.63) in Mahlavu indicating ineffectiveness of single TGX-221 treatments in both cells. The similarity between PIK-75 and its combinatory treatment with Sorafenib in Mahlavu (0.74) is also high (Fig. 3C), which may be an evidence of the underfilling action of Sorafenib alone treatment in PTEN deficient cells.

50 most commonly differentially regulated genes ranked by the sum of absolute logFC values were represented through a dendrogram in Fig. 3D. In Huh7, up- and down-regulated genes were well clustered. DUSP5, PCNA, VCAN, GADD45B and DUSP8 genes in Huh7 were the shared mostly. In Mahlavu, DEGs were not well separated like Huh7, some of the genes like EGR1, LINC00641, MIR6723, FOSB and ACTA2 were found to vary in different treatments. ESRG, CYP1B1-AS1, LDLR and TPTEP1 genes were the most common DEGs in Mahlavu.

Gene enrichment analysis of differential expression patterns in HCC cells

Considering the high correlation between specific treatments, in order to investigate the expression patterns between different inhibitory treatments in HCC cell lines, Huh7 and Mahlavu genes were clustered separately using their corresponding logFC values.

Heatmap analysis of DEGs revealed expression pattern of HCC cells and the gene enrichment analysis to these patterns (gene clusters) exploit the functional consequences (Fig. 4). In Huh7 cell line, expressions of single Sorafenib and its combined treatment with TGX-221 were highly correlated, which was also visualized in heatmap analysis. For all treatments in HCC, a positive regulation of extracellular matrix organization and developmental processes were observed, while regulation of cell proliferation and actin filament bundle assembly ontologies were more active in single PIK-75 treatment. PIK-75 and Sorafenib combined treatment resulted in downregulation of genes enriched in negative regulation of biosynthetic processes and cell fate commitment ontologies. Likewise, cholesterol metabolic process was downregulated for TGX-221 and its Sorafenib combinatory treatment. We also identified a group of genes involved in apoptosis stimulation process. The genes are 2 histone family proteins, 1 long intergenic non-translating RNA, uncharacterized proteins FAM184B and NCBP2AS2, NBP and NAG5 and ATP2A1 downregulated in the treatment of PIK-75 alone while they were upregulated all the other Huh7 treatments.

Fig. 4
figure 4

Gene Expression Patterns. Heatmaps of gene expressions illustrated as dendrograms separately for Huh7 and Mahlavu cells lines. We removed single PI3K-β inhibitor treatments for both cell lines considering its ineffectiveness. Sample sets for Huh7 and Mahlavu were separately joined, and united sets included 11,033 and 11,615 genes in total before filtration. Gene enrichment analysis was performed using BiNGO (FDR ≤ 0.05) and significant gene ontologies were selected according to the context. Hence, dendrogram analysis were performed on 581 genes for Huh7 and 583 genes for Mahlavu cells. For more detailed analysis and to view interactive dendrogram please see CanSyL github repository. Clusters were generated by heatmaply and colored; 8 for Huh7 and 6 for Mahlavu. Clusters not showing any significant enrichment were excluded. Up- and downregulated gene expression levels are colored as red and blue respectively, the intensity of the color indicates how strong the logFC value is. ALPHA; PIK-75, SALPHA; PIK-75 and Sorafenib, SBETA; TGX-221 and Sorafenib, SOR; Sorafenib treatments

Full size image

Immune response was downregulated more significantly in Mahlavu cells in treatments combined with Sorafenib treatment. Cation binding was enhanced for PIK-75 and its Sorafenib combination. Cholesterol metabolic processes, angiogenesis and vascular endothelial growth factor receptor 2 binding were downregulated for all treatments. A group of genes were upregulated in both single Sorafenib and combinatory TGX-221 treatments while they were downregulated in single PIK-75 and combinatory PIK-75 treatments. Because of this opposite effect, we anticipated the relation of these genes with antagonistic action of combined therapy of TGX-221 and Sorafenib. In this group, most of the genes were mitochondrial pseudo-genes. It is known that mitochondrial dysfunctions are mostly associated with apoptotic resistance and metabolism of tumor cells and one of HCC hallmarks points out the mitochondrial mutations in cancer development. Furthermore, downregulation of enzymes mediating oxidative phosphorylation for TGX-221 and Sorafenib treatment with respect to PIK-75 treatments confers the previous antagonistic nature in Mahlavu cells [27].

Network based interpretation with omics integrator

A traditional way of RNA-seq analysis is to use only DEG sets for gene enrichment analysis, which generally restricts the capture of complete cellular events. However, application of a conventional method to connect DEGs in a network though their known protein-protein interactions can reveal intersecting/hidden regulation patterns. Using the Omics Integrator tool, we adapted Prize Collecting Steiner Tree (PCST) algorithm to associate DEGs by adding intermediate genes (or Steiner nodes) aiming the construction of the most optimal gene-to-gene network. As reference network, we converted protein nodes in STRING human protein interaction network to gene nodes. Using Steiner nodes together with DEG sets introduced more specific gene ontologies. Distribution and relations of PCST enriched gene sets were presented in 5 sets Edwards–Venn diagrams for each treatment (Fig. 5A and D). NHBE was found to be differentially expressed in PIK-75 (ALPHA) and TGX-221 (BETA), and Sorafenib-PIK-75 (SALPHA). LAMP3, SIRPG and CD83 genes were differentially expressed in TGX-221 (BETA), and Sorafenib-PIK-75 (SALPHA). FOSB was differentially expressed in TGX-221 (BETA), Sorafenib-PIK-75 (SALPHA), and Sorafenib-TGX221 (SBETA). RICTOR and STK32A were differentially expressed in Sorafenib-TGX-221 (SBETA) and Sorafenib (SOR) treatments.

Fig. 5
figure 5

Network based interpretation of DEGs. Venn diagram scheme of Huh7 network nodes (A). Dendogram of GO enrichments for Huh7 (B). Network representation of PI3Ki-β and Sorafenib treated Huh7 cells (C). Dendogram of GO enrichments for Mahlavu (D). Venn diagram scheme of Mahlavu network nodes (E). Network representation of PI3Ki-β and Sorafenib treated Mahlavu cells (F). ALPHA; PI3Ki-α inhibitor, SALPHA; PI3Ki-α inhibitor and SOR, SBETA; PI3Ki-β inhibitor and SOR, SOR; Sorafenib treatments

Full size image

Although, DEG lists were found to be highly correlated, network based functional analysis revealed that kinase inhibitor treatments in both cell lines resulted in different biological processes. A better comparison of the networks was provided through a functional encolouring and sizing of the nodes and a systematic usage of the network centrality measures for clustering. PCST predicted optimal gene-to-gene networks were imported into Cytoscape, and gene logFC values were used to color the nodes to represent up- and downregulated branches. We arranged the sizes of the nodes according to their betweenness centrality to better organize hub genes. PCST optimal input parameters (final network statistics) were summarized in Supplementary Table 3. In Fig. 5B and E, network nodes including both DEGs and Steiners were compared. As a result of PCST, hidden expression patterns were identified. Since input DEG numbers for TGX-221 treated Huh7 and Mahlavu cells and Sorafenib treated Mahlavu cell were low, their networks were smaller. Sorafenib and TGX-221 combined treatments in Huh7 and single PI3K-α inhibitor (PIK-75) and combined treatments in Mahlavu were separated from the other treatments in the dendrograms.

We optimized the networks by limiting the number of trees to one. By minimizing their overall degrees to avoid hairballs we had more than one central hub nodes in the network generating more branches for the analysis (Fig. 5C and F). In order to understand the relatedness between the gene nodes in the networks, the functions of the branches should be exploited. Yet, the significance of biological enrichment analysis highly depend on the input size. The power of statistical analysis is low for large DEG sets. Eventually, we clustered networks using gene nodes’ betweenness centralities by Glay algorithm creating branches and applied BiNGO for each cluster/branch to get their enriched gene ontologies. Finally, we selected the significant gene ontologies for clusters, and connected them through the network. The ultimate network visualizations allowed us to analyze overall effect of up- or downregulation of genes and provided a comprehensive space for network comparisons through clusters. Other network representations can be reached in Supplementary Figs. 1–10 and associated Cytoscape files are in the referred CanSyL github repository.

Selection and validation of genes from optimal networks based on centrality metrics

In order to reveal candidate genes to be target for drug studies, optimal PCST networks were analyzed with their centrality metrics. Since hub nodes in the optimal networks were mostly the well-studied genes, we decided to eliminate them to find novel targets in branches. Although Omics Integrator scales the optimal networks avoiding hub node bias, we also filtered out the nodes that have betweenness centrality values higher than 0.001 after filtration of random nodes (frequency ≥ 0.01).

Then, we used centrality properties of optimal networks which were calculated by Networkx python library. Each network was filtered by degree, eigenvector and betweenness centralities that were higher than 0.001 allowing us both not to select the nodes at the end of the branches. The remaining nodes were sorted by inhibitor treatments and eigenvector centrality, and at most six genes were selected for each treatment (Fig. 6A). In the final sets, we have come up with 20 genes for each cell line (Fig. 6B). For Huh7 cells: CDC27, CCDC80, AARS2, ACSBG2 and CITED2 genes in PIK-75 inhibitor treatment, RIMKLA in TGX-221 treatment, CEBPB, DNAJC10, DLK1, EDEM1, ATP6V1D and DUSP8 genes in PIK-75 and Sorafenib combined treatment, LIN7C gene in TGX-221 and Sorafenib combined treatment, EXOC7, FEZ1, GAB2, HOXA10, BIRC7 and ANKRD28 genes in Sorafenib inhibitor treatment and for Mahlavu cells: ATP1B1, CACNA1H, CAPNS1, CCT7, ATG9A and BOLA2B genes in PIK-75 treatments, CGA and TNFRSF4 genes in TGX-221 inhibitor treatment, ALMS1, AOX1, BCL3, CD276, ANKRD1 and ASIC1 genes in PIK-75 and Sorafenib combined treatment, HMGCS1, GDF15, AGER, FABP1, ACOT12, CRHR1 genes in TGX-221 and Sorafenib combined treatment were prioritized for further investigations. For PTEN deficient Mahlavu cells along with single Sorafenib treatment, our prioritization strategy found no significant genes. Yet, it is interesting to find some of the targets from Steiner nodes (white boxes) since they cannot be exploited using classical differential expression analysis.

Fig. 6
figure 6

Proposed drug target genes. Prioritized nodes for Huh7 and Mahlavu was ranked by betweenness centrality values of randomized networks for each inhibitor treatment (A). Expressions of the genes in the cell lines. Prioritized treatments were pointed for corresponding drug treatment target (B). Relative expression profile of Mahlavu (C) and Huh7 (D) cells for AOX1 and AGER genes determined by qRT-PCR. Expression values were normalized with RPL19. Experiment was performed as triplicates and for statistical analysis, unpaired t-test with Welch correction was performed. *p < 0.05; **p < 0.01. ALPHA; PI3Ki-α inhibitor, HBETA; PI3Ki-β inhibitor, SALPHA; PI3Ki-α inhibitor and SOR, SBETA; PI3Ki-β inhibitor and SOR, SOR; Sorafenib treated cells

Full size image

In vitro validation of selected target genes AOX1 and AGER

We selected AOX1 and AGER genes to be validated by mRNA expression though qPCR experiments (Fig. 6C and D). Mahlavu and Huh7 cells were treated with Sorafenib or its combinations with PI3K-α and β inhibitors. AGER was selected because it is a pure Steiner node and not found in our DEG list. Whereas AOX1 was in both Steiner node and part of the DEG list. Hence, our qPCR results correlate and validated our network analysis results.

Furthermore, we performed AOX1 and AGER gene knockdown experiments on HCC cells to investigate the effects of the loss of these genes on cell proliferation. Knockdown with specific siRNAs (Fig. 7A) resulted in significant effect on cell proliferation (Fig. 7B). Hence, real-time cell proliferation analysis has shown that silencing AGER and AOX1 significantly inhibited growth of these cells with respect to negative control siRNA treatments. Overall, results from our in vitro experiments have supported and validated the systems level PCST network analysis.

Fig. 7
figure 7

Effects of knockdown of AGER and AOX1 genes on proliferation of HCC cell lines. Knockdown of AGER and AOX1 genes in Mahlavu and Huh7 cells using 25 nM siRNA validated with q-RT-PCR experiments. Experiments were performed in triplicates and results were represented as means ± SEM. One-Way ANOVA was performed for statistical analysis (A). Cell growth analysis of Mahlavu and Huh7 cells for which siRNA treatments targeting AGER and AOX1 genes were done and cell index values reflecting cell growth were monitored for 72 h. Knockdown of both genes resulted in significant drop in cell proliferation in both cells (B). Experiments were performed in triplicates, and results were represented as means ± SEM. Two-Way ANOVA was performed for statistical analysis. *p < 0.05, **p < 0.01, ***p < 0.001

Full size image

Discussion

Recently the influence of isoform diversity on responses to drugs with respect to large number of GPCR receptors is demonstrated at systems level [17]. Current studies focusing on the discovery of agents against HCC aim to identify target proteins that escape from regulatory signaling mechanisms of the cell. Conventionally, these studies concentrate on a single gene or locus which result in a comprehensive investigation of a new tumor driver gene, yet, in most cases single driver gene analyses are inadequate to solve the complex network of cancer pathogenesis as the interaction of several signaling pathways and their interlaced connections including HCC which represents high rate of tumor heterogeneity.

PI3K/Akt/mTOR signaling is involved in cellular growth, proliferation, and cell cycle progression in various cancer cells [6]. Frequent mutations and loss of function in PTEN tumor suppressor gene leads to constitutive activation of Akt protein hence activation of PI3K/Akt/mTOR cell survival pathway. Sorafenib, which is the most studied therapeutic agent for HCC, that targets Raf/MEK/ERK cascade. This is compensated by PI3K/Akt signaling activation in favor of cell survival in cancer cells which is one of the reasons for limited effectiveness of this drugs [8]. Sorafenib treatment frequently results in resistance to the treatment within nearly 6 months due to release of pro-inflammatory cytokines and chemokines in tumor microenvironment which promotes cancer stemness, tumor proliferation, and angiogenesis [28, 29]. Hence signaling pathways activated by these molecules in favor of Sorafenib resistance must be targeted in order to alter resistance toward Sorafenib.

There are several drugs in clinical trials targeting PI3Ks inhibiting tumor progression [22]. PI3K inhibitors are classified as pan-PI3K inhibitors and isoform specific inhibitors. Current molecular and clinical trial studies focus on the effectiveness of these inhibitors as well as the mechanisms of resistance to PI3K inhibition. In this study, we investigated the molecular alterations in combinatory treatment of isoform PI3kinase inhibitors targeting PI3K/Akt/mTOR pathway with Raf/MEK/ERK signaling inhibitor Sorafenib either on PTEN adequate (Huh7) or deficient (Mahlavu) HCC cells (Fig. 1). PI3K inhibitors are usually combined with mTOR inhibitors in order to increase the effectiveness of the treatment [22, 30]. However, there are limited studies which targets alternate cell signaling survival pathways (i.e. Raf/MEK/ERK vs PI3K/Akt/mTOR) with the aim of revealing the genes involved in synergistic or antagonistic resistance to inhibitors.

We showed synergistic inhibition of cell growth in both cell lines treated with PI3K-α inhibitor PIK-75 and Sorafenib. The synergistic cytotoxicity was more effective in PTEN-adequate Huh7 cells. While combinatory PI3K-β inhibitor TGX-221 treatment also synergistically inhibited cellular growth in Huh7, we observed a strong antagonistic effect in Mahlavu cells indicating the importance of isoform specific actions of Phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit isoform kinase inhibitors. We then investigated the molecular mechanisms involved in this differential phenotypic response to isoform specific inhibition by comprehensive network analysis of RNA-seq data upon treatment with drugs in combination with Sorafenib. Both inhibitors resulted with up regulation of key pathways in inflammation and immune response like; BCL3/NF-κB, cell proliferation (CDH2 and CCD1), Jun kinase and osmotic stress in Huh7 cells. Moreover, in TGX-221 and Sorafenib treatment, genes having role on regulation of programmed cell death and apoptotic mitochondrial changes were downregulated. It is known that interactions of GTPases with PI3K are isoform specific, while Ras cannot bind to p110-β, RAC1 and CDC42 proteins can activate p110-β [31]. Our network analysis showed that PI3K-α inhibitor (PIK-75) and its combined treatment with Sorafenib also showed a correlation in Huh7 cells. They both resulted with negative regulation of Erk1/Erk2 signaling and activation of MAPKK activity. Combinatory treatment with Sorafenib, mostly resulted with downregulation of hub proteins JUN, INSIG1, MDM2 and SOX9 associated with cancer cells. Hence, in PTEN adequate Huh7 cells, PIK-75 and Sorafenib treatment could decrease cell proliferation and decrease Sorafenib dependent immune response. Our data indicates that targeting PI3Ki-α isoform in an inhibited MAPK pathway background with Sorafenib would be a better therapeutic approach in both PTEN deficient and adequate Hepatocellular cancer cells.

On the other hand, combination of PI3K-β inhibitor (TGX-221) and Sorafenib in Mahlavu cells showed a strong antagonistic action, which probably depends on PI3Ki-α isoform activity. Our data with PTEN deficient Mahlavu cells demonstrated that constitutive PI3K/Akt pathway activation makes these cells more resistant due to PI3K-α isoform activation since the inhibition of PI3K-β with its specific inhibitor TGX-221 makes these cells resistant to Sorafenib. When TGX-221 and Sorafenib was combined, MAPK and nuclear factor kappa B (NF-κB) signaling upregulated and increasing activity in response to immune stress and inflammatory injuries genes were enriched (Figs. 4 and 5). In Mahlavu, TGX-221 combined with Sorafenib treatment shows a decreased level of Bcl-3 responsible for antagonistic action.

In this study we also performed system level network analysis in order to identify genes involved in isoform specific actions of PI3K inhibitors using DEG genes from our RNA-seq transcriptome data on STRING human protein interaction networks. Our prioritization strategy using topological features of the optimized networks identified hub and Steiner nodes representing genes involved in differential synergistic or antagonistic effects of isoform specific PI3K-α (PIK-75) or PI3K-β (TGX-221) inhibitors combined with Sorafenib. Many of these (DLK1, GAB2, BOLA2B, AOX1 and AGER) were closely related to cell proliferation and tumor progression, and associated with poor prognosis in HCC [32,33,34]. Steiner nodes prioritized in Mahlavu cells treated with PI3K-α (PIK-75) and Sorafenib identified Aldehyde oxidase 1 (AOX1), and Mahlavu cells treated with PI3K-β (TGX-221) and Sorafenib identified Advanced glycosylation end product-specific receptor (AGER) genes. Differential expression of these genes were validated by qPCR experiments as shown in Fig. 6C and D. Furthermore, siRNA knockdown of these two genes negatively affected cell proliferation significantly (Fig. 7). Both genes are associated with glucose metabolism and generation of reactive oxygen species and are involved in proinflammatory actions in liver carcinogenesis [35, 36]. AOX1 is considered one of the key biomarkers in HCC and abnormal expression of AOX1 is correlated with the poor prognosis [37]. AGER, also is shown as one of the main responsible factors in tumorigenesis of HCC cells in the presence of high glucose for diabetes [36].

Conclusions

Combination of targeted drugs to inhibit alternative compensatory pathways holds great promise for effective treatment of cancer including HCC. As we clearly demonstrated and validated both in silico and in vitro, in this study system level analysis of cellular networks in response to combination treatments and the investigation of the regulation signaling pathways are of necessity, because such treatments may result in an opposite of the desired effect. The importance of context dependent (PTEN status) PI3K/Akt/mTOR signaling inhibition must be taken into consideration during the use of isoform specific or pan-PI3K inhibitors in combination therapies with Sorafenib with respect to resistance in HCC cells.

Material and methods

Cell lines and kinase inhibitors

Mahlavu and Huh7, HCC cell lines were cultured in DMEM medium, supplemented with 10% fetal bovine serum (FBS), 1% penicillin/streptomycin (P/S) and 1% non-essential amino acids (NEA) and incubated in humidified 37 °C incubator with 5% CO2. Mahlavu and Huh7 cell lines were treated with the inhibitors which are listed in Supplementary Table 1 (Sorafenib (Nexavar) was from Bayer Healthcare Pharmaceuticals, Inc., NJ USA, Inhibitors PIK-75 (cat#528116), TGX-221 (cat#528113), LY294002 (cat#440202) were from Calbiochem).

Cytotoxicity and cell cycle

Huh7 (2000cell/well) and Mahlavu (1000 cell/well) cell lines were seeded into 96-well plates in 150 μl of medium/well. The next day, cells were treated with the drugs (Sorafenib (Nexavar) Bayer Healthcare Pharmaceuticals, Inc., NJ USA, or PIK-75 (PI3K-α inh.), TGX-221 (PI3K-β inh.), LY294002 (PanPI3Ki) and control (DMSO) in triplicates. After 72 h, the media was discarded, and the wells were washed with PBS, and 50 μl of 10% cold TCA (Merck, Germany) was added for fixation and incubated with TCA at + 4 °C in dark for 1 h. Then cells were washed with ddH2O for 4 times, and the plates were air-dried at room temperature. Finally, 50 μl of 0.4% sulphorhodamine B (SRB) (Sigma Aldrich) solution in 1% acetic acid was applied to each well, and the plates were incubated for 10 min in dark at room temperature. Excess dye was washed off with 1% acetic acid (4 or 5 washes). Finally, 200 μl of 10 mM cold Tris-Base was applied to each well to solubilize SRB. Then, the absorbance values were measured at 515 nm and were analyzed to determine the effect of each drug on cell proliferation compared to control [38]. Sorafenib, PIK-75 and TGX-221 were used to treat HCC cells in concentrations which described in respective figure legend and the cells were incubated for 96 h. Cell viability and DNA content calculations in flow cytometric cell cycle analysis was performed using propidium iodide.

Real-time cell electronic sensing (RT-CES) system for cell growth and cytotoxicity analysis and synergy analysis

50 μl Huh7 (2000 cell/well) and Mahlavu (1000 cell/well) cells were seeded into 96-well plates in 100 μl of medium/well. The next day, cells were treated with the drugs and CI (Cell Index) values were taken every 10 min for 4 h to get the fast drug response and then every 30 min to obtain the long-term drug response. Impedance measurements displayed as CI values reflect cell growth. The SynergyFinder Zero Interaction Potency (ZIP) model is used for the evaluation of the combined effect of PIK-75, TGX-221, and their combinations with Sorafenib (Sorafenib (Nexavar) Bayer Healthcare Pharmaceuticals, Inc., NJ USA). ZIP model defines the effect of combining two compounds by comparing the change in the dose-response curves between individual drugs and their combinations [39]. For monitoring the effects of siRNA treatment for AOX1 or AGER genes, Mahlavu and Huh7 cells were seeded onto 96-well E-Plate in triplicates. After overnight incubation, siRNA treatments were done as described previously. Cell index (CI) were recorded every 30 min for total of 96 h. Data was normalized using time-zero CI values (when siRNA treatment was performed). For statistical analysis, Two-Way ANOVA was performed using GraphPad Prism 8.

RNA extraction and sequencing

Total RNA was isolated with NucleoSpin RNA II Kit (Macherey-Nagel) according to the manufacturer’s protocol (MN, Duren, Germany) with small modifications such as 30 min of DNA digestion instead of 15 min and 2-step elution with 20 μl water instead of one elution with 60 μl. RNA concentration was measured with NanoDrop and A260/A280, A260/A230 ratios were checked for RNA quality and purity. Total RNAs were provided to BGI Tech (https://en.genomics.cn/) for sequencing. RIN values are acquired Agilent Bioanalyzer system and they were above 0.8 for all samples. Details of RNA-seq experiment and data can be found at PRJNA556552.

Wound healing

In the wound-healing assay, a wound was made in the middle of a confluent cell monolayer and the migration of cells to this area was assessed by taking photos at different time points and calculating the wound closure with respect to the initial wound width. Sorafenib, and TGX-221 were used at a concentration of 10 μM, except PIK-75 inhibitor, which was used at 0.1 μM concentration. Photos of the wounds were taken after 24 h and 48 h. The sizes of the wounds were calculated at all time points. At least 12 different wound distances were noted for each condition at each time point and the averages were used for analysis to construct the graphs.

Quantitative RT-PCR (qRT-PCR)

Mahlavu (100.000 cells/dish) and Huh7 (250.000 cells/dish) were seeded into 10 cm culture dishes. After 24 h, cells were treated with the inhibitors and incubated for 48 h and then collected for RNA isolation. RNA-purification kit (Qiagen, cat#74106) and cDNA synthesis (ThermoFischer, cat#K1621) kit were used according to manufacturer’s protocol. Total RNA amounts were measured with Nanodrop One (ThermoFisher). qPCR was initiated with 50 ng cDNA and performed with FastStart Essential DNA Green Master (Roche, cat#6402712001) via Roche LightCycler 96 Instrument, according to manufacturer’s protocol optimized for this instrument. Primer sequences are: AOX1-f: 5′-ggggtgttccgtgtttttcg-3′, AOX1- r: 5′-caggttcatctctcggaatcattt-3′, AGER-f: 5′-agcatcagcatcatcgaacca-3′, AGER-r: 5′-gcctttgccacaagatgacc-3′ and RPL19-f: 5′-gctctttcctttcgctgctg-3′, RPL19-r: 5′-ggatctgctgacgggagttg-3′. All reactions were performed in triplicates. The Ct (cycle threshold) values were normalized against RPL19 reference gene [40]. To determine the relative expression of target genes in inhibitor treated cells to that of DMSO treated cells, the ∆ ∆Ct method was used. Results were analyzed with GraphPad Prism 9.0.

siRNA transfection

Mahlavu (13.000 cells/well) and Huh7 (30.000 cells/well) were seeded onto 24-well plates in Penicillin-Streptomycin free DMEM; supplemented with 10% FBS, 1x L-Glutamine, and 1x non-essential amino acid solution. After overnight incubation, media was removed, and siRNA targeting AOX1 (ON-TARGETplus Human AOX1 (316) siRNA, SMARTpool cat#L-008291-00-0005) or AGER mRNA (ON-TARGETplus Human AGER (177) siRNA, SMARTpool, cat#L-003625-00-0005) were prepared in 1x siRNA Buffer (cat#B-002000-UB-100) and administered in 25 nM concentrations in the presence of FBS and Penicillin-Streptomycin free DMEM according to manufacturer’s protocol. As negative controls, 25 nM non-targeting pool siRNA (cat#D-001810-10-05) (NC-siRNA) and DharmaFECT 4 Transfection Reagent (cat#T-2004-02) treated groups (UNT) were used. After 12 h of incubation, treatment media was replaced with DMEM including 10% FBS, 1x Penicillin-Streptomycin, 1x L-Glutamine and 1x non-essential amino acid and incubated for the following 48 h. Cell pellets were collected via trypsinization, frozen in liquid nitrogen, and stored at − 800°. RNA isolation, cDNA conversion and q-RT-PCR experiments for AOX1 and AGER genes were performed as described previously. For statistical analysis, One Way ANOVA was performed using GraphPad Prism 8.

Bioinformatics methods

RNA-Seq analysis

RNA reads were processed by Illumina Hiseq 2000 (SE50). 12 FASTQ files (PR-JNA556552), were first analyzed through a well-known quality assessment tool; FASTQC [41]. Then, without any trimming, single-end reads were aligned to the reference human genome (GRCh38/hg38) using a split read aligner algorithm TopHat V2.1.0 [42]. TopHat itself features an ultrafast mapper Bowtie v2.2.6 algorithm [43]. After that, aligned reads were quantified by HTSeq-count v0.6.1 [44] for given human gene split regions (GRCh38 v84) to count how many transcripts map to each gene, which generates a gene level count matrix.

FASTQC analysis of this study results are included into the referred CanSyL github repository and other alignment metrics summarized in Supplementary Table 2.

Differential expression analysis for sequence count data

EdgeR [26], from Bioconductor package, is a widely used method for differential expression analysis. We used gene level count matrices of 12 RNA-seq treatment sets as input of EdgeR. DMSO treated Huh7 and Mahlavu cells were used as negative controls. EdgeR constructs a negative binomial model using the RNA count data. In our experimental design, there was no biological replicates of the samples to inherit the in-sample variation. EdgeR solves no-replication problem by suggesting a different dispersion calculation method to estimate variation within each sample compared to housekeeping genes. A set of housekeeping genes in Hepatocellular carcinoma was well characterized in Ersahin T. et al. [45]. We used these housekeeping genes to estimate biological coefficient of variation (BCV) value manually.

Before EdgeR analysis, genes with less than 5 readings were filtered out using counts per million constraint (cpm ≤ 5). A biological model was constructed by taking BCV as 0.045. Differential analysis performed using exactTest function of EdgeR package. Then, we limited logFC (log2 of fold change) to −/+ 2 for both cell lines. However, with this limitation, Mahlavu had no significant numbers of DEGs for downstream analysis. On the other hand, Huh7 resulted greater number discounting for single PI3Ki-β treatment. Therefore, for Mahlavu cell lines, a less stringent logFC value (− 1.5/+ 1.5) was used for further analysis. Finally, we selected the top DEGs according to following filters; p-value ≤0.01, FDR ≤ 0.01, and logFC ranges (− 2/+ 2) for Huh7 cell and (− 1.5/+ 1.5) for Mahlavu cell. Gene annotations were obtained using org. Hs.eg.db R package [46] from Bioconductor.

Dendrogram analysis

Heatmap representation is one of the most popular graphical methods for visualization of bigdata providing color encoding cells that represent numbers. Heatmaply [47] is a very powerful way of investigating clusters in a high dimensional data since final heatmap result is visualized as interactive graph offering inspection over the cells making possible for zoom-in. In our study, all dendrograms were visualized through heatmaply using default hclust clustering by using Euclidian as distance measure.

Gene ontology (GO) analysis

Given a set of genes on the network, Cytoscape plug-in BiNGO tool [48] maps functional terms to enriched genes to output GO terms and their statistical features. In order to have a better understanding of the processes that selected genes having role on, statistically over-represented GO terms were characterized using BiNGO in our analysis. We have used a very stringent Benjamini&Hochberg False Discovery Rate (FDR ≤ 0.005) to filter out non-significant GO-terms. GO sets containing redundant and electronically annotated terms generated noise for functional comparisons. In order to avoid those suspicious GO terms, we have only used GOs with experimentally validated codes (EXP, IDA, IEP, IGI, IMP, and IPI). The codes were matched though Gene Ontology Annotation (GOA) database.

Network construction and optimization of DEGs using PCST approach

PCST (Prize Collecting Steiner Tree) [49, 50] aims to identify sub-networks from an interaction network given a set of weighted genes. By using PCST, we have extracted the biologically meaningful interactions between the DEGs from human protein-protein interaction data. We used Omics Integrator software to implement PCST algorithm. We used Forest module in our analysis to determine multiple sub-pathways in the human interactome. PCST algorithm finds an optimal tree, including the terminal nodes (from DEG lists in our case) with prizes travelling through the interactome nodes which have costs of edges only if they are included. The task is to find the shortest paths between the prize nodes avoiding the costs on the edges. The algorithm minimizes the cost of all edges by passing through as many prize nodes as possible. In order to construct meaningful trees using DEGs, forest parameters must be fine-tuned. The size and degree of the forests are expected to vary as the number of genes in the input files changes. Forest parameters depend highly on the distribution of prizes and numbers of the nodes. The best combinations of parameters for each DEG set were explored using forest-tuner [51] which is PCST algorithm parameter tuner for ω, β and μ parameters. This script was used to find the best arrangements of the parameters to be used in Forest module for each treatment. We had searched the parameters in the following ranges: ω (1–10.0 or 5–15), β (1–15.0), μ (0.01–0.05). Here, ω parameter tunes the number of trees in the network, β parameter increases the number of prices entering the tree and μ is another parameter that arranges the dominance of hub proteins in the network. Among all of the possible solutions, we have selected the combination which generates a network with minimum mean degree. Optimal PCST parameters are summarized in Supplementary Table 3.

The interactome set given to Forest module was derived from STRING protein- protein interaction database v10 [52]. In STRING, network edges were scored according to a confidence score (range of 0 to 1) determined through an algorithm by the database. The confidence score gets higher as it gets more experimental proofs basically. In our analysis, we used the interactions only with high confidence proofs (at least 0.7). Omics Integrator performs -log10 conversion to the confidence score, so the cost negatively correlates to the confidence score.

Randomization tests

In order to test the significance of the nodes appearing in the optimal nodes, each PCST network was subjected to randomization tests using forest module (−randomTerminals 100). The tests were performed using random set of terminals with respect to keeping node numbers, and original interactome set with same edge weights and optimization penalties. The probability that a node randomly to be connects in the network was expressed by its frequency of randomness in the network. Therefore, less frequent nodes would be the most specific ones to the network. Throughout the analysis, we had used nodes that appeared only once in the random networks.

Network centrality

Centrality measures are the indicators of most valuable vertices in the graph for network analysis and they are often used to identify influential nodes of the network providing a ranking which identifies the important nodes in the network. We had used degree, eigenvector and betweenness centralities in order to estimate network topology. Networkx python library [53] was used to calculate centrality measures.

Effective visualization and clustering of the networks

Omics Integrator Forest module generates networks in .sif format which is compatible format for Cytoscape visualizations. Cytoscape includes many add-on for biological network analysis, therefore we both analyzed and visualized our networks on this tool. For our study, after .sif file was imported into the Cytoscape yFiles layout algorithms was implemented and hierarchical layout was selected for visualizations. In order to cluster the networks, Glay algorithm [54] using edge betweenness centrality was implemented. For proper annotation of the clusters AutoAnnotate plug-in [55] was used. Application of this strategy resulted in the most connected patterns in the networks. Then, in order to better understand what processes of these patterns have role on, statistically over-represented Gene Ontology (GO) terms were characterized using BiNGO for each cluster in networks. Only experimentally proven GOs (FDR ≤ 0.005) were used. Selected GO terms were also imported into the network using their gene maps. The final visualizations represented all gene relationships, up- and downregulated genes and internal Steiner nodes. Also, highly connected groups and GO annotations were provided for an easy and efficient way to compare networks with each other.

Prioritization of nodes in PCST generated networks

Based on the network topology, we developed a prioritization strategy for further investigation as drug targets. A node is treatment specific only if it occurs in the branches on random networks while present in more central areas on the optimal networks. In order to accomplish these nodes, we used the least frequent nodes (0.01) resulting from randomization test. Here, hub nodes of optimal networks were selected through using degree, eigenvector and betweenness centralities greater than 0.001. From these nodes, we eliminated the predominant nodes in the random network using degree centrality of random networks smaller than 0.001. Finally, the top 20 nodes for each treatment were selected and represented.

Availability of data and materials

FASTQ reads of the RNA-seq experiments are available in the NCBI-SRA repository [https://www.ncbi.nlm.nih.gov/sra/PRJNA556552]. The datasets generated and analysed during the current study are available in the CanSyL lab github repository [https://github.com/cansyl/Isoform-spesific-PI3K-inhibitor-analysis].

Abbreviations

HCC:

Hepatocellular carcinoma

DEGs:

Differentially expressed genes

GO:

Gene Ontology

FDR:

False Discovery Rate

logFC:

log2 of fold change

References

  1. Bray F, Ferlay J, Soerjomataram I. Global Cancer Statistics 2018: GLOBOCAN Estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68:394–424. https://doi.org/10.3322/caac.21492.

    Article  PubMed  Google Scholar 

  2. Perz JF, Armstrong GL, Farrington LA, Hutin YJF, Bell BP. The contributions of hepatitis B virus and hepatitis C virus infections to cirrhosis and primary liver cancer worldwide. J Hepatol. 2006;45(4):529–38. https://doi.org/10.1016/j.jhep.2006.05.013.

    Article  PubMed  Google Scholar 

  3. Sun B, Karin M. Review obesity, inflammation , and liver cancer. J Hepatol. 2012;56(3):704–13. https://doi.org/10.1016/j.jhep.2011.09.020.

    Article  CAS  PubMed  Google Scholar 

  4. Aleksandrova K, Stelmach-Mardas M, Schlesinger S. Obesity and liver cancer. Recent Results Cancer Res. Fortschritte der Krebsforschung. Progres dans les recherches sur le cancer. 2016;208:177–98. https://doi.org/10.1007/978-3-319-42542-9_10.

    Article  CAS  PubMed  Google Scholar 

  5. Farazi PA, DePinho RA. Hepatocellular carcinoma pathogenesis: from genes to environment. Nat Rev Cancer. 2006;6(9):674–87. https://doi.org/10.1038/nrc1934.

    Article  CAS  PubMed  Google Scholar 

  6. Ersahin T, Tuncbag N, Cetin-Atalay R. The PI3K/AKT/mTOR interactive pathway. Mol Bio Syst. 2015;11(7):1946–54. https://doi.org/10.1039/c5mb00101c.

    Article  CAS  Google Scholar 

  7. Moeini A, Cornell’a H, Villanueva A. Emerging signaling pathways in hepatocellular carcinoma. Liver Cancer. 2012;1(2):83–93. https://doi.org/10.1159/000342405.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc J-F, et al. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med. 2008;359(4):378–90. https://doi.org/10.1056/NEJMoa0708857.

    Article  CAS  PubMed  Google Scholar 

  9. Villanueva A, Llovet J. Targeted therapies for hepatocellular carcinoma. Gastroenterology. 2013;140(5):1410–26. https://doi.org/10.1053/j.gastro.2011.03.006.

    Article  CAS  Google Scholar 

  10. Targeted Network CGAR. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017;169(7):1327–134123. https://doi.org/10.1016/j.cell.2017.05.046.

    Article  CAS  Google Scholar 

  11. Dimri M, Satyanarayana A. Molecular signaling pathways and therapeutic targets in hepatocellular carcinoma. Cancers. 2020;12(2):491. https://doi.org/10.3390/cancers12020491.

    Article  CAS  PubMed Central  Google Scholar 

  12. Marin JJG, Macias RIR, Monte MJ, Romero MR, Asensio M, Sanchez-Martin A, et al. Molecular bases of drug resistance in hepatocellular carcinoma. Cancers. 2020;12(6):1–26. https://doi.org/10.3390/cancers12061663.

    Article  CAS  Google Scholar 

  13. Ward A, Shukla K, Balwierz A, Soons Z, König R, Sahin O, et al. MicroRNA-519a is a novel oncomir conferring tamoxifen resistance by targeting a network of tumour-suppressor genes in ER + breast cancer. J Pathol. 2014;233(June):368–79. https://doi.org/10.1002/path.4363.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Bae J-J, Rho J-W, Lee T-J, Yun S-S, Kim H-J, Choi J-H, et al. Loss of heterozygosity on chromosome 10q23 and mutation of the phosphatase and tensin homolog deleted from chromosome 10 tumor suppressor gene in Korean hepatocellular carcinoma patients. Oncol Rep. 2007;18(4):1007–13. https://doi.org/10.3893/or.18.4.1007.

    Article  CAS  PubMed  Google Scholar 

  15. Buontempo F, Ersahin T, Missiroli S, Senturk S, Etro D, Ozturk M, et al. Inhibition of Akt signaling in hepatoma cells induces apoptotic cell death independent of Akt activation status. Investig New Drugs. 2011;29(6):1303–13. https://doi.org/10.1007/s10637-010-9486-3.

    Article  CAS  Google Scholar 

  16. Engelman JA. Targeting PI3K signalling in cancer: opportunities, challenges and limitations. Nat Rev Cancer. 2009;9(8):550–62. https://doi.org/10.1038/nrc2664.

    Article  CAS  PubMed  Google Scholar 

  17. Marti-Solano M, Crilly SE, Malinverni D, Munk C, Harris M, Pearce A, et al. Combinatorial expression of GPCR isoforms affects signalling and drug responses. Nature. 2020;587(7835):650–6. https://doi.org/10.1038/s41586-020-2888-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wang J, Zhao W, Guo H, Fang Y, Stockman SE, Bai S, et al. AKT isoform-specific expression and activation across cancer lineages. BMC Cancer. 2018;18(1):742. https://doi.org/10.1186/s12885-018-4654-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Wadhwa B, Paddar M, Khan S, Mir SA, Clarke P, Grabowska AM, et al. Akt isoforms have discrete expression in triple negative breast cancers and roles in cisplatin sensitivity. Oncotarget. 2020;11(45):4178–94. https://doi.org/10.18632/oncotarget.27746.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Riggio M, Perrone MC, Polo ML, Rodriguez MJ, May M, Abba M, et al. AKT1 and AKT2 isoforms play distinct roles during breast cancer progression through the regulation of specific downstream proteins. Sci Rep. 2017;7(1):44244. https://doi.org/10.1038/srep44244.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Vanhaesebroeck B, Guillermet-Guibert J, Graupera M, Bilanges B. The emerging mechanisms of isoform-specific PI3K signalling. Nat Rev Mol Cell Biol. 2010;11(5):329–41. https://doi.org/10.1038/nrm2882.

    Article  CAS  PubMed  Google Scholar 

  22. Yang J, Nie J, Ma X, Wei Y, Peng Y, Wei X. Targeting PI3K in cancer: mechanisms and advances in clinical trials. Mol Cancer. 2019;18(1):26. https://doi.org/10.1186/s12943-019-0954-x.

    Article  PubMed  PubMed Central  Google Scholar 

  23. He Q, Johnston J, Zeitlinger J, City K, City K. PI3K p110α and p110β have differential effects on Akt activation and protection against oxidative stress-induced apoptosis in myoblasts. Cell Death Differ. 2015;33(4):395–401. https://doi.org/10.1038/nbt.3121.ChIP-nexus.

    Article  CAS  Google Scholar 

  24. Liu P, Cheng H, Roberts TM, Zhao JJ. Targeting the phosphoinositide 3-kinase pathway in cancer. Nat Rev Drug Discov. 2009;8(8):627–44. https://doi.org/10.1038/nrd2926.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Durmaz I, Guven EB, Ersahin T, Ozturk M, Calis I, Cetin-Atalay R. Liver cancer cells are sensitive to Lanatoside C induced cell death independent of their PTEN status. Phytomedicine. 2016;23(1):42–51. https://doi.org/10.1016/j.phymed.2015.11.012.

    Article  CAS  PubMed  Google Scholar 

  26. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009;26(1):139–40. https://doi.org/10.1093/bioinformatics/btp616.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Weilin W, Qiang S, Zehui W, Dongkai Z, Jianfeng W, Haiyang X, et al. Mitochondrial dysfunction-related genes in hepatocellular carcinoma Weilin. Front Biosci. 2013;18:1141–9. https://doi.org/10.2741/4169.

    Article  Google Scholar 

  28. Sethi N, Kang Y. Unravelling the complexity of metastasis — molecular understanding and targeted therapies. Nat Rev Cancer. 2011;11:735. https://doi.org/10.1038/nrc3125.

    Article  CAS  PubMed  Google Scholar 

  29. Kahraman DC, Kahraman T, Cetin-Atalay R. Targeting pi3k/akt/mtor pathway identifies differential expression and functional role of il8 in liver cancer stem cell enrichment. Mol Cancer Ther. 2019;18(11):2146–57. https://doi.org/10.1158/1535-7163.MCT-19-0004.

    Article  CAS  PubMed  Google Scholar 

  30. Calero R, Morchon E, Martinez-Argudo I, Serrano R. Synergistic anti-tumor effect of 17AAG with the PI3K/mTOR inhibitor NVP-BEZ235 on human melanoma. Cancer Lett. 2017;406:1–11. https://doi.org/10.1016/j.canlet.2017.07.021.

    Article  CAS  PubMed  Google Scholar 

  31. Pan Y-F, Tan Y-X, Wang M, Zhang J, Zhang B, Yang C, et al. Signal regulatory protein alpha is associated with tumor-polarized macrophages phenotype switch and plays a pivotal role in tumor progression. Hepatology (Baltimore, Md.). 2013;58(2):680–91. https://doi.org/10.1002/hep.26391.

    Article  CAS  Google Scholar 

  32. Xu X, Liu R-F, Zhang X, Huang L-Y, Chen F, Fei Q-L, et al. DLK1 as a potential target against cancer stem/progenitor cells of hepatocellular carcinoma. Mol Cancer Ther. 2012;11(3):629–38. https://doi.org/10.1158/1535-7163.MCT-11-0531.

    Article  CAS  PubMed  Google Scholar 

  33. Chen Y, Liu Q, Wu M, Li M, Ding H, Shan X, et al. GAB2 promotes cell proliferation by activating the ERK signaling pathway in hepatocellular carcinoma. Tumor Biol. 2016;37(9):11763–73. https://doi.org/10.1007/s13277-016-5019-9.

    Article  CAS  Google Scholar 

  34. Luo J, Wang D, Zhang S, Kuan H, Haijun W, Li J, et al. BolA family member 2 enhances cell proliferation and predicts a poor prognosis in hepatocellular carcinoma with tumor hemorrhage. J Cancer. 2019;4(6):1–32. https://doi.org/10.7150/jca.31829.

    Article  CAS  Google Scholar 

  35. Weigert J, Neumeier M, Bauer S, Mages W, Schnitzbauer AA, Obed A, et al. Small-interference RNA-mediated knock-down of aldehyde oxidase 1 in 3T3-L1 cells impairs adipogenesis and adiponectin release. FEBS Lett. 2008;582(19):2965–72. https://doi.org/10.1016/j.febslet.2008.07.034.

    Article  CAS  PubMed  Google Scholar 

  36. Qiao Y, Zhang X, Zhang Y, Wang Y, Xu Y, Liu X, et al. High glucose stimulates tumorigenesis in hepatocellular carcinoma cells through AGER-dependent O-GlcNAcylation of c-Jun. Diabetes. 2016;65(3):619–32. https://doi.org/10.2337/db15-1057.

    Article  CAS  PubMed  Google Scholar 

  37. Ma X, Zhou L, Zheng S. Transcriptome analysis revealed key prognostic genes and microRNAs in hepatocellular carcinoma. PeerJ. 2020;8:8930. https://doi.org/10.7717/peerj.8930.

    Article  Google Scholar 

  38. Onen-Bayram FE, Durmaz I, Scherman D, Herscovici J, Cetin-Atalay R. A novel thiazolidine compound induces caspase-9 dependent apoptosis in cancer cells. Bioorg Med Chem. 2012;20(17):5094–102. https://doi.org/10.1016/j.bmc.2012.07.016.

    Article  CAS  PubMed  Google Scholar 

  39. Yadav B, Wennerberg K, Aittokallio T, Tang J. Searching for drug synergy in complex dose-response landscapes using an interaction potency model. Comput Struct Biotechnol J. 2015;13:504–13. https://doi.org/10.1016/j.csbj.2015.09.001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Tomizawa M, Shinozaki F, Motoyoshi Y, Sugiyama T, Yamamoto S, Ishige N. 2-Deoxyglucose and sorafenib synergistically suppress the proliferation and motility of hepatocellular carcinoma cells. Oncol Lett. 2017;13(2):800–4. https://doi.org/10.3892/ol.2016.5510.

    Article  PubMed  Google Scholar 

  41. Andrews, S.: FastQC 0.11.4: A quality control tool for high throughput sequence data (2010). Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

    Google Scholar 

  42. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11. https://doi.org/10.1093/bioinformatics/btp120.9605103.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357. https://doi.org/10.1038/nmeth.1923.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Anders S, Pyl PT, Huber W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9. https://doi.org/10.1093/bioinformatics/btu638.

    Article  CAS  PubMed  Google Scholar 

  45. Ersahin T, Carkacioglu L, Can T, Konu O, Atalay V, Cetin-Atalay R. Identification of novel reference genes based on MeSH categories. Plos One. 2014;9(3). https://doi.org/10.1371/journal.pone.0093341.

  46. Marc Carlson: org. Hs.eg.db: Genome wide annotation for Human. R package version 3.4.0 (2016).

    Google Scholar 

  47. Galili T, O’Callaghan A, Sidi J, Sievert C. heatmaply: an R package for creating interactive cluster heatmaps for online publishing. Bioinformatics (Oxford, England). 2018;34(9):1600–2. https://doi.org/10.1093/bioinformatics/btx657.

    Article  CAS  Google Scholar 

  48. Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics. 2005;21(16):3448–9. https://doi.org/10.1093/bioinformatics/bti551.

    Article  CAS  PubMed  Google Scholar 

  49. Tuncbag N, Braunstein A, Pagnani A, Huang S-SC, Chayes J, Borgs C, et al. Simultaneous reconstruction of multiple signaling pathways via the prize-collecting steiner forest problem. J Comput Biol. 2013;20(2):124–36. https://doi.org/10.1089/cmb.2012.0092 PMID: 23383998.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Tuncbag N, Milani P, Pokorny JL, Johnson H, Sio TT, Dalin S, et al. Network modeling identifies patient-specific pathways in glioblastoma. Sci Rep. 2016;6:1–12. https://doi.org/10.1038/srep28668.

    Article  CAS  Google Scholar 

  51. Budak G. forest-tuner: GitHub Repository; 2018. Cited 2019 May 9. https://github.com/gungorbudak/forest-tuner

    Google Scholar 

  52. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(D1):447–52. https://doi.org/10.1093/nar/gku1003.

    Article  CAS  Google Scholar 

  53. Hagberg AA, Schult DA, Swart PJ. Exploring network structure, dynamics, and function using NetworkX. In: Proceedings of the 7th Python in Science Conference (SciPy 2008) (SciPy); 2008. p. 11–5.

    Google Scholar 

  54. Su G, Kuchinsky A, Morris JH, States DJ, Meng F. GLay: community structure analysis of biological networks. Bioinformatics. 2010;26(24):3135–7. https://doi.org/10.1093/bioinformatics/btq596.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Kucera M, Isserlin R, Arkhangorodsky A, Bader GD. AutoAnnotate: a Cytoscape app for summarizing networks with semantic annotations. F1000Res. 2016;5:1717. https://doi.org/10.12688/f1000research.9090.1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We are grateful to our colleagues part of CanSyL laboratory for their comments on the study and manuscript.

Funding

This work was supported by The Scientific and Technological Research Council of Turkey Grant #110S388 and The Turkish Ministry of Development project (#KanSil 2016 K121540). KN is supported by TUBITAK (2211) scholarship.

Author information

Authors and Affiliations

  1. Cancer System Biology Laboratory, CanSyL, Graduate School of Informatics, Middle East Technical University, 06800, Ankara, Turkey

    Kubra Narci, Deniz Cansen Kahraman, Altay Koyas, Tulin Ersahin, Nurcan Tuncbag & Rengul Cetin Atalay

  2. Present Address: Section of Pulmonary and Critical Care Medicine, the University of Chicago, Chicago, IL, 60637, USA

    Rengul Cetin Atalay

Authors
  1. Kubra Narci
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Deniz Cansen Kahraman
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Altay Koyas
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Tulin Ersahin
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Nurcan Tuncbag
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rengul Cetin Atalay
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

KN analyzed the data and wrote the manuscript. DCK, TE and AK performed the experimental analysis. DCK wrote the manuscript. NT provided Omics Integrator tool and helped through data analysis. RCA designed the experiments and wrote the manuscript. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Rengul Cetin Atalay.

Ethics declarations

Competing interest

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Supporting ResultsIC50 values, RNA-seq quality, optimal PCST networks and clustering and enrichment analysis of PCST generated networks are listed and figured in this additional file.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Narci, K., Kahraman, D.C., Koyas, A. et al. Context dependent isoform specific PI3K inhibition confers drug resistance in hepatocellular carcinoma cells. BMC Cancer 22, 320 (2022). https://doi.org/10.1186/s12885-022-09357-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-022-09357-y

Keywords

BMC Cancer

ISSN: 1471-2407

Contact us