Identification of DEGs and influential genes perturbed by L-pampo (Step 1)Quantifying gene expressionOur main objective in this study is to determine whether L-pampo drives the differential responses of cellular networks in different cell types. Particularly, to investigate gene signatures associated with L-pampo in immune cells and cancer cells, we chose a differentiated human monocyte cell line, THP-1, that highly expresses TLR2/3 and induces various immune responses to L-pampo12,13. Additionally, we utilized two cancer cell lines, PC-3 and SW620, that can express TLR2/314,15 but manifest differential responses by L-pampo treatment. The cell death was induced to PC-3 in a dose-dependent manner, while L-pampo did not affect the cell death on SW620 (Figure S1a-b). Given that cancer cell death is an essential cellular process in the tumor microenvironment that leads to enhancing anti-tumor immune responses, we sought to examine transcriptomes related to cell death induced by L-pampo. We treated cells with L-pampo and isolated RNA samples at different time points (0, 3h, and 6h), and performed the RNA-seq to obtained gene expression profiles. As expected, the TLR2/3 are more highly expressed in immune cell (THP-1) compared to the cancer cells (Figure S1c). Also, we observed that transcriptomes of TLR2/3 were more highly expressed in PC-3 than SW620, suggesting that the signal transmission from the TLR receptors may depend on TLR2/3 expression in both immune cells and cancer cells (Figure S1c).Differentially expressed gene (DEG) and functional analysesFirst, we analyzed differentially expressed genes (DEGs) in three cell lines (THP-1, SW620, and PC-3) and three time points (0h, 3h, and 6h after treating L-pampo). The analyses revealed that 3,732 genes, 501 genes, and 1,415 genes upregulated or downregulated in THP-1, SW620 and PC-3 cell lines across different time points, respectively (Figure S2a).To understand the biological roles of the putative genes induced by L-pampo, we performed over-representation analysis for each DEG at THP-1, SW620, and PC-3 using KEGG pathway database16. Only 20 DEGs were identified among three cell lines and related to toll-like receptor (TLR) signaling pathway and cytokine-cytokine receptor interaction pathway. The 558 genes in the intersection between THP-1 and PC-3 were related with ribosome and oxidative phosphorylation (OXPHOS) (Figure S2b). The over-representation analysis showed that the common DEGs from 0 to 6h in THP-1 and PC-3 were related to JAK-STAT signaling and cytokine-cytokine receptor interaction besides ribosome and OXHOS, while no significant changes were observed in SW620 (Fig. 2).Figure 2KEGG pathway analysis for DEGs from THP-1, SW620, and PC-3. The columns represent cell lines and time changes, and the rows show the selected significant KEGG pathway terms from the cell lines. In each column, the length of each bar indicates the number of DEGs contained by the respective KEGG pathway. The bar color shows the significance of pathway enrichment on the -log10 scale for the adjusted p-value. The bold characters indicate the KEGG pathways of our interest perturbed by L-pampo in three cell lines.The gene signatures related to reactive oxygen species (ROS), NOD-like receptor signaling, IL-17 signaling, TGF-beta signaling, and lipid and atherosclerosis were also upregulated in PC-3 between 0 and 6h. Interestingly, RIG-1, a pathogen pattern recognition receptor, was also found to be significantly upregulated in SW620 and PC-3 compared to THP-1 (Figure S2c). On the other hand, the DEGs related with necroptosis which is a programmed inflammatory cell death17 were increased between 3 and 6h in PC-3 (Figure S2d). As an example, NLRP3, a well-known key regulator in necroptosis varied in the same way as necroptosis genes18.Influential gene analysis by network propagationDEG analysis is very powerful since it led us to a probably list of perturbed genes by L-pampo treatment across cell types and time points. However, those genes/proteins that play critical roles given a molecular stimulus (e.g., diseases) are often not caught by only DEG analysis. In many cases, DEGs may reflect significant transcriptomic changes as an effect of perturbation not a cause. On the other hand, conventional pathway enrichment analysis based on DEGs often yields a large number of pathways, making it challenging to organize the results and draw clear conclusions from fragmented evidence. Moreover, our DEG analysis only showed 20 DEGs common to the three cell lines in the 0h-6h comparison (Figure S2b), which account for just 0.5% of the total DEGs. This implies that hidden regulators may be responsible for the differences among the cell lines.To tackle this challenge, adopting gene–gene interaction networks is a potent strategy to augment DEG analysis and this approach has been used to identify functional genetic modules and even drug targets19. Fundamentally, using gene–gene interaction networks can allow us to search the vicinity of DEGs for hidden players that may cause or mediate important molecular signals. Therefore, we devised a network-based analysis and investigated the complex gene–gene interactions in immune cells and cancer cells to understand the L-pampo induced mechanism. To identify hidden players, we utilized network propagation that prioritizes genes by iteratively considering the influence of the seed genes (i.e., DEGs) on their neighbor genes in a protein–protein interaction (PPI) network19. We ranked the genes that may highly affect DEGs, naming the top-ranked genes as “influential genes” and this method as “influential gene analysis”.A PPI network is a collection of all known and curated interactions between genes/proteins. Analyzing the connections of DEGs and transcription factors on a PPI network can provide insights into key regulators at the transcriptome level10. However, all the interactions/correlations reported within a PPI network are not functionally meaningful to our specific problem solving, as a PPI network contains any known biological interactions accumulated to date.Thus, we also computed gene–gene correlations from our data and prune those edges (i.e., interactions reported in the PPI network) that did not have high correlation coefficients (Pearson correlation coefficient < 0.7)20 . As a result, we constructed our own gene–gene interaction network reflecting cellular specificity of our data (THP-1, SW620, and PC-3) and we then identified the top 200 influential genes using DEGs as seeds on the constructed network for the following steps.Instantiation of a template network for comparison (Step 2)Construction of a template networkAlthough consideration of temporal dynamics and cell line characteristics is essential for our study design, the DEG and influential gene analyses are constrained in multifactorial comparisons due to distinct baselines for each condition. To address this limitation, we built and used a universal gene–gene interaction network to compare transcriptomic profiles of all different cells and time points on the same baseline and referred to this as “template network”. We utilized DEGs, influential genes, and transcription factors to construct the template network so that it can reflect molecular signaling flows induced by L-pampo on the three cell lines. In this network, we particularly added transcription factors from CollecTRI database20 to further consider transcriptional regulations over DEGs and influential genes. As a result, 7,252 genes were selected and used as input for WGCNA’s hard thresholding. In contrast to WGCNA’s soft thresholding, hard thresholding applies a cutoff to remove all weak correlations, simplifying the network structure and making it less sensitive to outliers. When we chose an optimal cutoff of 0.9 based on the scale-free topological fit, we obtained the final template network with 6,585 nodes and 1,214,811 edges.Identification of perturbed modules with L-pampoNext, we constructed hierarchical clustering using WGCNA on 4,582 DEGs identified across all conditions, showing the results using a dendrogram (Figure S3a). WGCNA requires a power parameter for a soft threshold because it assumes that increasing the power can reduce noise in the gene–gene correlation matrix. We chose a power of 18, which achieves a scale-free topological fit by R2 larger than 0.9 (Figure S3b). To understand the relationship between gene modules and traits (i.e., cell type and time point), we performed module-trait relationship analysis, which informed the modules change in different directions depending on the cell type or time (Figure S3c).Different co-expression patterns over cell lines and time points guided us to select gene modules differentially perturbed by L-pampo. Among the ten identified modules, we finally chose five ‘perturbed modules’, named as A, B, C, D, and E, respectively. The correlation of Module A went up at 3h and stayed high until 6h only in THP-1. Module B had high correlation values in PC-3 compared to the other two cell lines. The correlation of Module C also increased at later time points in all cell lines, but it moves up faster in PC-3 than THP-1 (i.e., 3h vs. 6h). On the other hand, those of Modules D and E clearly decreased by time in THP-1 but not in the other two cancer cell lines.Next, we conducted a gene set enrichment analysis on the perturbed modules and found that there was an increase in the levels of PI3K-AKT signaling (Modules A and B) and oxidative phosphorylation (OXPHOS) (Modules B and C) (Figure S3d), respectively. To identify finer expression patterns within Modules A and B, we conducted a community detection analysis using the Walktrap algorithm to divide each module into submodules. This algorithm clusters densely connected gene communities given a network by comparing the neighborhood genes through simulated random walks21. We identified three submodules within Module A and named them as Modules A-1, A-2, and A-3. We repeated the same process for Module B, identifying two submodules and marking them as Modules B-1 and B-2.Upon a deeper analysis of the submodules, it was found that Module A-1 contained genes related to JAK-STAT, OXPHOS, and TLR signaling pathways. In contrast, Modules A-2 and A-3 were more associated with RAP-1 and PI3K-AKT pathways, respectively. JAK-STAT is closely linked to TLRs and involved in various cellular processes such as immune response, inflammation, or carcinogenesis22. Module B has a significantly large number of genes (676 genes), and most of these genes are highly expressed in PC-3. Modules B-1 and B-2 were related to PI3K-AKT and OXPHOS, whereas Module C was uniquely characterized by OXPHOS. These findings demonstrate that PI3K-AKT, JAK-STAT, and OXPHOS pathways can play essential roles in L-pampo-mediated response in THP-1 and PC-3.Oxidative phosphorylation (OXPHOS) and reactive oxygen species (ROS) as a key process in the L-pampo response by comparing three cellular networks (Step3)To identify the key molecular signals triggered by L-pampo, we used the TOPAS algorithm23 to analyze the shortest paths from the TLR pathway genes to the perturbed modules within the template network. We obtained six subnetworks, namely TLR-Module A, TLR-Module B-1, TLR-Module B-2, TLR-Module C, TLR-Module D, and TLR-Module E, that summarize the molecular interactions among the TLR pathway genes and the perturbed modules (Figure S4a). In addition to the qualitative assessment, we also attempted to quantify the differences of all subnetworks between cell lines. We used network-based Jensen-Shannon divergence (nJSD) to estimate the distance between probability distributions of the networks of interest24. Interestingly, our analysis revealed that TLR-Module A, TLR-Module B, and TLR-Module C subnetworks have longer distances between THP-1 and SW620 than THP-1 and PC-3. We also observed that these subnetworks showed larger differences between THP-1 and PC-3 compared to the other subnetworks (Fig. 3a).Figure 3Comparison of detected subnetworks. (a) Quantifying the relative differences of cell line pairs in terms of subnetwork between TLR genes and selected gene modules. The x-axis represents the cell line pairs and the y-axis the average of nJSD for the subnetworks. The average of nJSD was obtained by taking the mean of nJSDs at different time points. (b–d) Visualizing the expression changes in the subnetwork connected from TLR signaling genes to Module C (TLR-Module C) for the three cell lines: (b) THP-1, (c) SW620, (d) PC-3. Each circle represents a gene, and it has three time points of 0h, 3h, and 6h (counterclockwise). The node color represents the normalized TPM (Transcripts Per Million) using z-score across different cell lines and time points per gene. The color of the outer ring means the module ID (Module A: red, Module B: brown, Module C: magenta, Module D: green, and Module E: black) and the nodes marked by asterisks are TLR signaling genes.Especially, TLR-Module C subnetwork exhibited different expression patterns across the cell lines, which reflect the cell line specificity against L-pampo treatment (Fig. 3b-d). THP-1 showed the highest expression levels for TLR signaling genes that are related to pro-inflammatory cytokines, chemokines and activation markers compared to SW620 and PC-3 (Fig. 3b, 5 o’clock directed clusters with asterisk mark). SW620 displayed reduced expression levels in TLR-related genes within the subnetwork. However, some TLR signaling genes such as MAPK8, MAPK12, MAPK13, MAP2K6, TAB1, and AKT2 were relatively highly expressed in SW620 compared to THP-1 and PC-3, but their expression levels did not appear to be affected by L-pampo treatment (Fig. 3c, 1 o’clock directed cluster with asterisk mark). PC-3 showed increased gene expression levels in all genes within Module C and some Module B genes including MDH1, LSM3, and PSMA2 (magenta and brown color outer rings, respectively) compared to THP-1 and SW620 (Fig. 3d). Interestingly, MDH1 plays an essential role in the conversion of malate to oxaloacetate during the tricarboxylic acid (TCA) cycle in the OXPHOS process25.In particular, we observed that PIK3CD and IL-12A in Module B-1 played a crucial role in connecting the above major subnetworks such as TLR-Module A, TLR-Module B-2, and TLR-Module C (Figure S4a). TLR-Module A is enriched mostly by PI3K-AKT pathway genes, and PIK3CD and IL-12A acted as a bridge connecting Modules A-1 and A-3 (Figure S4a). These two genes are involved in TLR, PI3K-AKT, and JAK-STAT signaling pathways16,26. Genes related to JAK-STAT signaling, such as IL-6, CSF2, CSF3, and IL-7R in TLR-Module A-1 (Figure S4a), were connected to TLR-Module C through hub genes IL-6 and S100A6 (Fig. 3b-d). Interestingly, PIK3CD and IL-12A simultaneously worked as a bridge connecting not only TLR-Module B-1 and TLR-Module B-2 but also TLR-Module C and TLR-Module B-2 (Figure S4a). Therefore, TLR-Module C can be bridged with TLR-Module B through these two genes as a hub.Given that OXPHOS and ROS pathways were enriched in Modules A-1, B-2, and C (Figure S3d), we examined the expression changes of OXPHOS and ROS genes in each module over time. This showed that the average expression levels of OXPHOS and ROS genes in Modules A-1, B-2, and C notably increased in PC-3 (Fig. 4a). When collecting all the OXPHOS and ROS genes from the modules detected by WGCNA, genes belonging to Modules A-1, B-2, and C in PC-3 exhibited distinct expression patterns compared to THP-1 and SW620 (Figure S5, the genes marked by stars). More specifically, those genes in Module C showed substantially increased expression levels at the earlier time point (3h) in PC-3. However, their expression levels were low but continued to increase until the latest time point (6h) in THP-1 (Figure S5).Figure 4Reconstructing functional networks by combining different modules. (a) Average of RNA expression of OXPHOS and ROS genes of each module. Each point and error bar indicates mean ± standard deviation of the normalized-TPM using z-score across different cell lines and time points per gene. Line colors denote three different cell lines: THP-1 (red), SW620 (green), and PC-3 (blue). The number under each module title indicates the number of OXPHOS and ROS genes in the that module. (b) The reconstructed network of Modules A-1, B-2, and C with OXPHOS and ROS genes (PC-3 DEGs). Each node represents a gene, and it has three time points of 0h, 3h, and 6h (counterclockwise). The node color represents the normalized TPM using z-score across different cell lines and time points per gene. The color of the outer ring means the module ID (Module A: red, Module B: brown, and Module C: magenta). The nodes marked by asterisks are OXPHOS and ROS genes.To elucidate the OXPHOS and ROS process in differential L-pampo response, we conducted network analyses using TOPAS on the OXPHOS and ROS pathway genes with the combination of Modules A-1, B-2, and C (Fig. 4b). Our analysis revealed that genes related to OXPHOS and ROS, such as NADH dehydrogenase (NDUF) and ATP synthase, which form complexes of the mitochondrial respiratory chain27, were predominantly upregulated in Modules A-1, B-2 and C (Fig. 4b, asterisk marked). As shown in Fig. 4a, these genes had higher basal expression levels in PC-3 than other cell lines and further activation of these genes due to L-pampo may have led PC-3 to cell death (Figure S5).Collectively, our study revealed that L-pampo treatment triggered PI3K-AKT and JAK-STAT pathways via TLRs and likely led to cell death through OXPHOS and ROS pathways. We found that PIK3CD and IL-12 in Module B-1 connected PI3K-AKT to OXPHOS and ROS. Also, these two genes worked as a hub by linking Modules B-2 and C as well. Furthermore, OXPHOS and ROS also communicated with JAK-STAT pathway through IL-6 and S100A6 connecting Modules A and C (Fig. 5). These findings suggest that the OXPHOS and ROS process may play a critical role in distinguishing the differential responses in cancer cell death by L-pampo.Figure 5An illustrative diagram represents the inferred signaling transduction for L-pampo-induced PC-3 cell death L-pampo engages Toll-like receptors (TLRs) of cancer cells and triggers two primary signaling pathways: PI3K-AKT and JAK-STAT. Activation of these pathways appears to increase oxidative phosphorylation (OXPHOS) and reactive oxygen species (ROS) production through hub genes, which leads to L-pampo-induced cancer cell death, especially in PC-3. The color of the box means the module ID (Module A: red, Module B: brown, and Module C: magenta).
