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

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). Supplementary Information The online version contains supplementary material available at 10.1186/s12885-022-09357-y.

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.

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 hyperactive Akt stimulation.

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 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 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).

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.
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, upand 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 NCB-P2AS2, NBP and NAG5 and ATP2A1 downregulated in the treatment of PIK-75 alone while they were upregulated all the other Huh7 treatments.
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. 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, 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.

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.

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 Fig. 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 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.

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 timezero 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 Nan-oDrop and A260/A280, A260/A230 ratios were checked for RNA quality and purity. Total RNAs were provided to BGI Tech (https:// en. genom ics. 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.
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 plugin 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.