Role of the antineoplastic drug bleomycin based on toxicogenomic-DNA damage inducing (TGx-DDI) genomic biomarkers data: A meta-analysis

Objectives: Accurately identifying the cellular, biomolecular, and toxicological functions of anticancer drugs help to decipher the potential risk of genotoxicity and other side effects. Here, we examined bleomycin for cellular, molecular and toxicological mechanisms using next-generation knowledge discovery (NGKD) tools. Methods: This study was conducted at the Faculty of Applied Medical Sciences, King Abdulaziz University (KAU), Jeddah, Saudi Arabia in October 2022. We first analyzed the raw Toxicogenomic and DNA damage-inducing (TGx-DDI) gene expression data from Gene Expression Omnibus (GEO) (GSE196373) of TK6 cells treated with 10 µM bleomycin and TK6 cells treated with DMSO for four hours using the GEO2R tool based on the Linear Models for Microarray Analysis (limma) R packages to derive the differentially expressed genes (DEGs). Then, iPathwayGuide was used to determine differentially regulated signaling pathways, biological processes, cellular, molecular functions and upstream regulators (genes and miRNAs). Results: Bleomycin differently regulates the p53 pathway, transcriptional dysregulation in cancer, FOXO pathway, viral carcinogenesis, and cancer pathways. The biological processes such as p53 class mediator signaling, intrinsic apoptotic signaling, DNA damage response, and DNA damage-induced intrinsic apoptotic signaling and molecular functions like ubiquitin protein transferase and p53 binding were differentially regulated by bleomycin. iPathwayGuide analysis showed that the p53 and its regulatory gene and microRNA networks induced by bleomycin. Conclusion: Analysis of TGx-DDI data of bleomycin using NGKD tools provided information about toxicogenomics and other mechanisms. Integration of all “omics” based approaches is crucial for the development of translatable biomarkers for evaluating anticancer drugs for safety and efficacy.


INTRODUCTION
Bleomycin is a glycopeptide broad-spectrum antibiotic originally isolated from the bacterium Streptomyces verticillus. 1 Bleomycin is a cytotoxic chemotherapy drug approved by the Food and Drug Administration (FDA) to treat various malignancies such as Hodgkin's disease or Hodgkin's lymphoma, non-Hodgkin's lymphoma, head and neck cancer, ovarian carcinoma, squamous cell carcinoma of the vulva, testes and penis either alone or in combination with other drugs both as curative and palliative treatment. 2 The reactive oxygen species (ROS) increase carbohydrate oxidation, lipid peroxidation and affect or influence the synthesis and breakdown of prostaglandins. 2 Bleomycin interferes with the growth and proliferation of cancer cells and helps to inhibit cancer. 2,3 Also, bleomycin is useful in treating cancerinduced pleural effusions and growth. However, like most antineoplastic agents' bleomycin can also affect normal cells leading to other unwanted or undesirable side effects. Bleomycin use may cause decreased appetite, vomiting, nausea, skin darkening, hair loss, but most importantly it can induce severe pulmonary toxicity in cancer patients. 3,4 The transcriptome profiles of the immortalized lymphoblastoid cells (TK6 cells) in vitro after treatment with different drugs or chemicals were obtained to decipher their genotoxicity. 5 Unlike most cancer cells which are routinely used for basic or translational research, the TK6 cells do not carry many genetic changes or mutations. 6 The TGx-DDI biomarkers were originally determined from transcriptome profiles of TK6 cells treated with chemical agents or drugs with well-categorized modes of action. 7 Here, we assessed the genotoxicity and other adverse molecular effects of bleomycin in TK6 cells using the NanoString TGx-DDI nCounter® assay data derived from Gene Expression Omnibus (GEO) based on next-generation knowledge discovery (NGKD) approaches.

Data Source:
The TGx-DDI gene expression data was obtained from GEO [(GSE196373 (super series) and GSE196369 (sub series)]. 8 As specified in these GEO datasets, the nCounter® experiments were conducted following the instructions outlined in the NanoString nCounter custom TGx-DDI Assay Panel (GEO platform GPL31920). 8 Ethical Statement: This study was conducted using publicly available gene expression data from the Gene Expression Omnibus (GEO). 8 It was conducted at the Faculty of Applied Medical Sciences, King Abdulaziz University (KAU), Jeddah, Saudi Arabia in October 2022. It was exempt from Institutional Review Board (IRB) approval because it did not involve animal models or human subjects. 9 GEO2R Analysis of NanoString nCounter® Data: We analyzed the data sets of human lymphoblastoid cell line TK6 treated with 10 µM bleomycin (treatment) and TK6 cells treated with DMSO (control) for four hours using the GEO2R analysis tool. Data were obtained using the TGx-DDI nCounter® assay using high quality, purified RNA. DMSO served as a control while bleomycin (BL; 10 µM) was used as a treatment group in the nCounter® experiments performed as previously reported. 8,10 The GEO data was compared using GEO2R based on the GEOquery and Linear Models for Microarray Analysis (limma) R packages obtained from the Bioconductor project, which provides R-based tools for highthroughput analysis of genomic data. The DEGs were identified using the limma R package using a log2 fold change (Log2FC) cut-off of 2 and p value cut off of 0.05. The DEGs were further analyzed using iPathwayGuide for NGKD analyzes as previously described. 9 Next Generation Knowledge Discovery using iPathwayGuide: The DEGs were subjected to iPathwayGuide analysis as previously described. 9 The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used for the pathway analysis in iPathwayGuide. 11 In short, the iPathwayGuide evaluates differentially regulated pathways based on the impact analysis method as described previously. 12 The overrepresentation of DEGs was used in the impact analysis method in a given pathway and the expression changes calculated used for the perturbation of that pathway topology. 12 The probability scores of over-representation analysis (pORA) and total pathway accumulation (pAcc) were applied to compute the scores and combined into a specific pathway-specific p-value using Fisher's method. 13 This p-value was then corrected for multiple comparisons using the false discovery rate (FDR) 14 and Bonferroni corrections. 15 In the iPathwayGuide Gene Ontology (GO) analysis for biological processes, molecular function, and cellular components, the number of DEGs annotated with each GO term was compared to the random expected number of DEGs observed in this study. 16 iPathwayGuide applies an over-representation approach (ORA) to determine the statistical significance of observing at least the  13 The p-value is estimated using the hypergeometric distribution as defined for pORA of pathways. This p-value is corrected for multiple comparisons using FDR and Bonferroni. 14,15 Upstream Regulator Analysis: The enrichment of DEGs identified in the present study and their network of regulatory interactions obtained from Advaita proprietary knowledge base was used to predict the upstream gene regulators in the iPathwayGuide analysis. The prediction of active miRNAs is based on the enrichment of differently downregulated target genes of the miRNAs. miRNAs mostly inhibit their targets and therefore the iPathwayGuide method calculates the ratio between the number of differentially down-regulated genes or targets and all DEGs or targets and compares it to the ratio of all downwardly expressed targets to all targets. Overall, iPathwayGuide assesses the likelihood of observing at least the number of differentially downregulated target genes for a given miRNA purely by probability. This p-value is computed using the hypergeometric distribution as described for pORA. 13

RESULTS
In the present study, DEGs identified by GEO2R analysis (Table-I) were used to decipher biological knowledge and associations from a total of 54,683 genes available in the Advaita Knowledge Base (AKB).
These data were associated with pathways from the KEGG database (Release 100.0+/11-12, November 21), 11 gene ontologies from the Gene Ontology Consortium database analyzed (2021-Nov4), 16 18 About 66 signaling pathways were found to be significantly affected. It was also found that 616 gene ontology terms (GO) and 117 gene upstream regulators were significantly enriched before correction for multiple comparisons. The five major differentially regulated pathways by bleomycin include the p53 pathway, transcriptional dysregulation in cancer, the forkhead box O (FOXO) pathway, viral carcinogenesis, and pathways in cancer (Table-II). The p53 pathway (KEGG: 04115), triggered by a range of stress signals (Fig.1), and the FoxO pathway (Fig.2), in which a family of transcription factors regulates genes involved in cell cycle, apoptosis and resistance to oxidation involves stress and longevity, glucose metabolism, etc. The five major biological processes differentially regulated by bleomycin are signaling through a p53 class mediator, intrinsic apoptotic signaling through a p53 class mediator, DNA damage response, signaling through a p53 class mediator  and intrinsic apoptotic signaling in response to DNA damage by a p53 class mediator (Table-III). Approximately 11 genes were differentially regulated by bleomycin in signal transduction through the p53 class mediator. We identified six DEGs in the intrinsic pathway of apoptosis in response to DNA damage by a p53-class mediator (GO:0042771). About seven DEGs were involved in the DNA damage response, signal transduction through the p53 class mediator (GO:0030330). Molecular functions that were significantly differentially regulated by bleomycin included ubiquitin protein transferase activity (GO:0004842) and p53 binding (GO:0002039) ( Table-IV). There were three DEGs that regulated ubiquitin protein transferase activity and three DEGs involved in p53 binding in bleomycin treated TK6 cells. The iPathwayGuide analysis predicted that the major upstream regulatory gene induced by bleomycin (Fig.3A) is p53 and its corresponding regulatory gene networks ( Fig.3B and Fig.3C) and microRNAs such as hsa-miR-218-5p, hsa-miR-200a-3p, hsa-miR-141-3p, hsa-miR-323a-3p and hsa-miR-202-5p (Fig.3D).

DISCUSSION
In the present study, we used the NGKD approach to delineate the molecular mechanism of genotoxicity and other adverse molecular effects of bleomycin in TK6 cells using the NanoString TGx-DDI nCounter® assay data derived from GEO. 10 Bleomycin differentially regulated genes involved in the DNA damage response and signal transduction through the p53 class mediator. Here,  bleomycin induced the stress response-associated p53 gene in TK6 cells. The signals of stress response such as DNA damage, oxidative stress and activated oncogenes induce the p53 gene. The tumor suppressor protein (TP53) encoded by p53 has three domains, namely the transcription activation domain, the DNA binding domain, and the oligomerization domain. TP53 is a transcriptional regulator for a number of genes involved in cell cycle arrest, cellular senescence, DNA repair, and apoptosis. Importantly, mutations in p53 are typically associated with a variety of cancers. 19 Bleomycin activates the intrinsic apoptotic signaling pathway in response to DNA damage by the p53 class mediator, where an array of intracellular molecular signals is conveyed to trigger apoptosis mediated cell death. 19 Bleomycin-induced FOXO family of transcription factors are critical for cell cycle regulation, resistance to oxidative stress, programmed cell death, etc. 20 FOXOmediated regulation involves phosphorylation by the serine-threonine kinase Akt/protein kinase B (Akt/ PKB) and the downstream phosphatidylinositol 3-kinase (PI3K) in response to various stimuli. 19 Moreover, FOXO phosphorylation at three conserved residues causes the export of FOXO proteins from the nucleus to the cytoplasm, reducing gene-induced FOXO. On the other hand, stress-activated c-Jun N-terminal kinase (JNK) and energy-sensitive AMP-activated protein kinase (AMPK) phosphorylates and activate FOXO proteins that are dependent on oxidative and nutritional stress stimuli. Apart from PKB, JNK and AMPK, FOXOs are regulated by a number of post-translational modifications such as acetylation, methylation and phosphorylation. 20 Cyclindependent kinase inhibitor 1A (CDKN1A) encodes a potent cyclin-dependent kinase inhibitor that inhibits cyclin-cyclin-dependent kinase2 (CDK2) and regulates the cell cycle in G1 phase. 21 CDKN1A expression is tightly regulated by TP53 and leads to cell cycle arrest in G1 phase in response to bleomycin. This CDKN1A protein was cleaved by CASP3-like caspases, which led to activation of CDK2, triggering apoptosis after caspase activation. 21 Growth inhibition and DNA damage-inducible alpha (GADD45A) belongs to a group of genes induced by stressful conditions of growth inhibition and treatment with DNA-damaging agents. The protein encoded by GADD45A responds to environmental stress by mediating the activation of p38/JNK signaling via the MTK1/ MEKK4 kinase. 22 DNA damage-induced transcription of this gene is mediated by both p53-dependent and -independent mechanisms. MDM2 is transcriptionally regulated by p53 and is a proto-oncogene encoding a nuclear localized E3 ubiquitin ligase that targets tumor suppressor proteins such as p53 for proteasomal degradation. Overexpression or amplification of the MDM2 locus is found in a number of cancers. 19 The antitumor effect of the upstream miRNA hsa-miR-218-5p has been demonstrated in numerous studies. 23,24 The upstream miRNA has-miR-181c-5p has also shown antitumor potential in numerous studies. 25 Our study sheds light on the rapid identification of toxicological functions and mechanism of action(s) of drugs or drug candidates based on TGx-DDI gene expression data using NGKD tools to unravel the safety and efficiency and thus help in the anti-cancer drug discovery and development which is of pivotal importance and requirement for successfully treating patients suffering from cancer. Furthermore, the integration of all "omics" based methodologies is essential for the development of translatable biomarkers for evaluating anticancer drugs for safety and efficacy.
Limitations: We conducted our study using publicly available NanoString TGx-DDI nCounter® datasets to evaluate the genotoxicity and other adverse molecular effects of bleomycin in TK6 cells based on NGKD approaches. However, in the future, integration of miRNA and metabolome profiling and other "omics" based approaches such as genomics, proteomics, glycomics, etc. may be critical for the development of translatable biomarkers for antineoplastic drug safety and efficacy. In addition, in vitro and in silico characterization of known and novel compounds may provide clues that aid in the development and discovery of much safer drugs in the near future.

CONCLUSIONS
The NGKD analysis of bleomycin TGx-DDI gene expression data revealed that a variety of genotoxic stress pathways, genes, and miRNAs were differentially regulated in TK6 cells. Our study demonstrated that NGKD tools can be used to unravel the molecular mechanism of genotoxicity and other adverse molecular effects of bleomycin. This provided important clues for the identification of additional biomarkers for the critical evaluation of safety and efficacy of druggable candidates. The identification of toxicological functions of drug candidates based on TGx-DDI gene expression data using NGKD tools is essential for evaluating the potential risk of genotoxicity and helps in the process of developing safe and effective drugs.