Simulation Study of cDNA Dataset to Investigate Possible Association of Differentially Expressed Genes of Human THP1-Monocytic Cells in Cancer Progression Affected by Bacterial Shiga Toxins

Shiga toxin (Stxs) is a family of structurally and functionally related bacterial cytotoxins produced by Shigella dysenteriae serotype 1 and shigatoxigenic group of Escherichia coli that cause shigellosis and hemorrhagic colitis, respectively. Until recently, it has been thought that Stxs only inhibits the protein synthesis and induces expression to a limited number of genes in host cells, but recent data showed that Stxs can trigger several signaling pathways in mammalian cells and activate cell cycle and apoptosis. To explore the changes in gene expression induced by Stxs that have been shown in other systems to correlate with cancer progression, we performed the simulated analysis of cDNA dataset and found differentially expressed genes (DEGs) of human THP1-monocytic cells treated with Stxs. In this study, the entire data (treated and untreated replicates) was analyzed by statistical algorithms implemented in Bioconductor packages. The output data was validated by the k-fold cross technique using generalized linear Gaussian models. A total of 50 DEGs were identified. 7 genes including TSLP, IL6, GBP1, CD274, TNFSF13B, OASL, and PNPLA3 were considerably (<0.00005) related to cancer proliferation. The functional enrichment analysis showed 6 down-regulated and 1 up-regulated genes. Among these DEGs, IL6 was associated with several cancers, especially with leukemia, lymphoma, lungs, liver and breast cancers. The predicted regulatory motifs of these genes include conserved RELA, STATI, IRFI, NF-kappaB, PEND, HLF, REL, CEBPA, DI_2, and NFKB1 transcription factor binding sites (TFBS) involved in the complex biological functions. Thus, our findings suggest that Stxs has the potential as a valuable tool for better understanding of treatment strategies for several cancers.

Shiga toxin (Stxs) is a family of structurally and functionally related bacterial cytotoxins produced by Shigella dysenteriae serotype 1 and shigatoxigenic group of Escherichia coli that cause shigellosis and hemorrhagic colitis, respectively. Until recently, it has been thought that Stxs only inhibits the protein synthesis and induces expression to a limited number of genes in host cells, but recent data showed that Stxs can trigger several signaling pathways in mammalian cells and activate cell cycle and apoptosis. To explore the changes in gene expression induced by Stxs that have been shown in other systems to correlate with cancer progression, we performed the simulated analysis of cDNA dataset and found differentially expressed genes (DEGs) of human THP1-monocytic cells treated with Stxs. In this study, the entire data (treated and untreated replicates) was analyzed by statistical algorithms implemented in Bioconductor packages. The output data was validated by the k-fold cross technique using generalized linear Gaussian models. A total of 50 DEGs were identified. 7 genes including TSLP, IL6, GBP1, CD274, TNFSF13B, OASL, and PNPLA3 were considerably (<0.00005) related to cancer proliferation. The functional enrichment analysis showed 6 down-regulated and 1 up-regulated genes. Among these DEGs, IL6 was associated with several cancers, especially with leukemia, lymphoma, lungs, liver and breast cancers. The predicted regulatory motifs of these genes include conserved RELA, STATI, IRFI, NF-kappaB, PEND, HLF, REL, CEBPA, DI_2, and NFKB1 transcription factor binding sites (TFBS) involved in the complex biological functions. Thus, our findings suggest that Stxs has the potential as a valuable tool for better understanding of treatment strategies for several cancers.
The Stxs are capable to cross epithelial, endothelial, leukocytic, lymphoid and other neuronal cells through transcytotic or paracellular mechanisms (Malyukova et al., 2009), may then associate with blood monocytes, macrophages and neutrophils to circulate the bloodstream and bind to toxin-binding glycosphingolipid Gb3 receptor by clathrin-dependent or clathrin-independent mechanisms (Lingwood, 2003;Schweppe et al., 2008). After internalization, the toxins undergo retrograde intracellular flow to reach the endoplasmic reticulum (Sandvig et al., 2002). Importantly, the toxin receptor, Gb3, has a limited expression in normal tissues but is overexpressed in several types of cancer and besides apoptosis, may cause "apoptosis induced proliferation" . Stxs are cytotoxic proteins that enter the host cell via macropinosome and function as an N-glycosidase, cleave a specific adenine nucleobase from the 28S RNA of the 60S subunit of the ribosome, thereby halting host cell protein synthesis Lukyanenko et al., 2011). This activity of the toxin resides in the A subunit and the pentamer of similar B subunit intercedes toxin binding to the Gb3 (Fraser et al., 1994;Lingwood et al., 2010;Melton-Celsa, 2014). These toxins can also activate host cell's signaling pathways and trigger apoptosis in many cell types. They induce apoptosis of epithelial, endothelial, leukocytic, lymphoid and neuronal cells (Tesh, 2010). After activation, the interleukin-1 (IL-1) and tumor necrosis factor alpha (TNF-α) sensitize endothelial cells to the action of Stx in vitro by increasing Gb3 expression (van de Kar et al., 1992;Stricklett et al., 2002). Upon sensitization, the innate immune response stimulated by Stxs may contribute to the development of vascular lesions (Ramegowda et al., 1999). It has been reported that when human monocytic THP-1 cells are treated with Stxs, they secrete tumor necrosis factor alpha (TNF-α)-, interleukin-1 (IL-1), and interleukin-6 (IL-6), which further alter the expression of these cytokines and chemokines (Leyva-Illades et al., 2010). Transcriptional regulation involves prolonged activation of stress-associated protein kinases JNK (c-Jun N-terminal kinases) and p38, extracellular signal-activated kinase1/2 (ERK1/2), and activation of transcription factors NF-kappaB (NF-κB) (nuclear factor) (Thorpe et al., 2001;Harrison et al., 2005). Recent evidence indicates that signaling through MAPK pathways and proapoptotic proteins -mostly caspasescan induce proliferation of neighboring surviving cells to replace apoptotic and dying cells. This "apoptosis induced proliferation" is critical for tissues regeneration and cancer progression (Tesh, 2010;Ryoo and Bergmann, 2012).
The expression profiling of treated and untreated THP-1 cells verifies that Stxs are responsible for genetic alterations (Leyva-Illades et al., 2010). Studies show that Stxs is associated with an elevated secretion of cytokines and other chemical mediators responsible for numerous diseases including cancer (DesRochers et al., 2015;Hattori et al., 2016). E. coli producing Stxs was isolated from cancer and diarrheagenic individuals (Chao et al., 2017). It has been observed that Stxs significantly caused increased expression of GRO, G-CSF, IL-1β, IL-8, and TNFα in THP-1 like cells with an increased level of proinflammatory cytokines and chemokines causing several diseases (Brandelli et al., 2015). The exact mechanism of influence of gene expression by Stxs as well as the functions of these influenced genes in respect to pathophysiology and molecular biology of cancer development have not been completely understood (Tesh, 2010). In this regard, cDNA microarray technology is a valuable tool to discover differential gene expression. The accumulated data from gene expression studies in public repositories provides an opportunity to construct pooled gene expression data sets from a larger number of individuals. In this study, we performed cDNA differential analysis of the publicly available dataset to discover the DEGs and to further explore the transcriptional responses of human THP1monocytic cells affected by bacterial Stxs. We uncovered the aberrantly expressed genes which are significantly linked to cancer. We also report the analysis of the pathophysiological pathways and predicted potential regulatory motif of these genes.

Accession of cDNA Dataset
The goal of this study was to identify aberrantly expressed cancer-related genes in the human THP1-monocytic cells lines affected by bacterial Shiga toxins (Supplementary Figure 1). We acquired the cDNA microarray dataset (NCBI, 2015) GSE19315 (Leyva-Illades et al., 2010) from Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE19315) 1 . This dataset examined the transcriptional changes in human THP-1 cells to Stxs toxins. It was performed by taking three untreated (controls) and six treated replicates (samples) followed by the real-time studies. The Stxs changed the genetic expression which is associated with the increased production of cytokines and other proinflammatory mediators. The GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Inc., Santa Clara, CA, 95051, USA, Technology: in situ oligonucleotide) platform was used, and the Annotation information (hgu133plus2) of probes was used to detect the gene expression.

Normalization and Identification of Differentially Expressed Genes (DEGs)
Pheno-data files of this dataset were prepared in identifiable format and missing values imputed (Troyanskaya et al., 2001). The Bioconductor packages including ArrayQuality Metrics on FIGURE 1 | Normalization and analysis of array quality metrics (A) Box-plot representing summaries of the signal intensity distributions of the arrays. Each box corresponds to one array. Outlier detection was performed by computing the Kolmogorov-Smirnov statistic Ka between each array's distribution and the distribution of the pooled data (B) Heatmap of the distances between arrays. The color scale is chosen to cover the range of distances encountered in the dataset (C) Histogram representing expression after normalization (D) MA plots. M and A are defined as: M = log 2 (I 1 ) -log 2 (I 2 ), A = 1/2 (log 2 (I 1 )+log 2 (I 2 )), where I 1 is the intensity of the array studied, and I 2 is the intensity of a "pseudo"-array that consists of the median across arrays.  R version 3.1.3 were used to execute the preprocessing steps of normalization and quality control (Bolstad et al., 2003;Fujita et al., 2006;Obenchain et al., 2014). Robust Multi-array Analysis (RMA) was used to adjust the background and normalization (Yoon et al., 2004) for perfect matches (PM) and mismatches (MM) in order to get a summary of intensities (summary To assess the quality of RNA in samples, AffyRNAdeg, summaryAffyRNAdeg, and plotAffyRNAdeg Bioconductor packages were used for degradation analysis (Affymetrix, 1999(Affymetrix, , 2001. After normalization, we performed the relative study and identified differentially expressed genes by pairwise comparison from genomic experiments (Tusher et al., 2001) and multiple testing corrections were completed by Benjamini-Hochberg method (Benjamini and Hochberg, 1995). The Bioconductor package LIMMA, a modified statistic that is proportional to the statistic with sample variance offsets, was used to shortlist the DEGs. By this package, duplicate spots and quality weights were Signif. codes: 0 "***" 0.001 "**" 0.01 "*" 0.05 "." 0. measured. The moderated statistics were calculated; genes were ranked with respect to the resulting scores and p-values. A false discovery rate (FDR) less than 0.05, p ≤ 0.05, Average Expression Level (AEL) ≥40% and an absolute log fold change (logFC) greater than 1 were set as the significant cutoffs.

K-Fold Cross Validation
We used the k-fold technique (Seymour, 1993) to validate the shortlisted differentially expressed genes using the Bioconductor "boot" package (Ripley, 2010). Bootstrapping can be helpful to correct some of the bias associated with the analysis. The generalized linear models (Gaussian) were applied and the "cv.glm" function was used to assess the k-fold cross validation for these cases. This function, firstly, predicts the error of the raw cross-validation followed by the adjusted cross validation estimation. The Gaussian function was trailed by the complete Leave-One-Out-Cross-Validation (LOOCV) procedure. The LOOCV method is instinctively named as one case is left out as the testing set and the rest of the data are used as the training set (Ripley, 2010).

Curation of Cancer-Related Genes
From the shortlist of DEGs, cancer-associated genes were curated from OMIM (Online Mendelian Inheritance in Man) database (www.ncbi.nlm.nih.gov/omim) 2 and Cancer GeneticsWeb database (www.cancer-genetics.org) 3 to filter the significant cancer-related genes. The roles of sorted genes in cancer were  evaluated through reported experimental studies using PubMed 4 database.

Cluster and Functional Enrichment Analysis of DEGs
We performed the cluster analysis (Eisen et al., 1998) based on expression values in each sample to verify the variations in gene expression levels between untreated and treated replicates. The functional Annotation and pathway enrichment analysis of cancer-related genes help us to reveal biological functions (Nam and Kim, 2008;Muhammad et al., 2015), and it was performed using the web-based DAVID (Database for Annotation Visualization and Integrated Discovery) (Huang da et al., 2009) and FunRich Annotation tools (Pathan et al., 2015).

Construction of Protein-Protein Interaction (PPI) Network
Proteins usually interact with each other to perform biological functions (Li et al., 2004;Muhammad et al., 2014;Naz et al., 4 PubMed: http://www.ncbi.nlm.nih.gov/pubmed. Accessed 2015-08-28. 2015). Therefore, the interacting partners of cancer-related genes were assessed with a high confidence score (0.999) from the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) (Szklarczyk et al., 2011) and HAPPI (Human Annotated and Predicted Protein Interaction) databases (Chen et al., 2009), and PPI-network were constructed using Cytoscape software (Cline et al., 2007) version 3.2.1. STRING and HAPPI connect major databases and predict interactions based on experiments, text mining and sequence homology. The role and association of these genes (source and target) in cancer were accessed from OMIM (www.omim.org), Cancer Genetics Web (www.cancerindex.org) and National Cancer Database (www. facs.org) 5 . The relationships between source proteins and target proteins (interacting proteins) were finally calculated.

De Novo Prediction of Regulatory Motifs
To better understand the complex mechanisms of cancerrelated, differentially expressed genes controlling important biological functions at transcriptional and post-transcriptional levels, oPOSSUM web tool (Ho Sui et al., 2005) was used to identify the transcription factor (TF) binding sites and overrepresented regulatory motif of target matrices (Pavesi et al., 2004) in promoter regions of expressed genes using default parameters.

Microarray Analysis, Normalization, and RNA Degradation Plots
The microarray data consists of 9-samples with 54675 probe ids from a study of the effect of bacterial Shiga toxins on human THP1 macrophage/monocytic cell lines (Leyva-Illades et al., 2010). The AffyBatch object comprises the size of the array 1,164 × 1,164 features with 54,675 affyids. Background correction and normalization by quantile normalization were completed to avoid systematic variation between samples. The normalized probe-level data represents expression levels of genes of entire DNA chip (Figure 1) involved in experimental condition and individual probes in a probe set were ordered by location relative to the 5 ′ end of the targeted RNA molecule. The 3 ′ /5 ′ intensity gradient has been shown to depend on the degree of competitive binding of specific and of non-specific targets to a particular probe. Poor RNA quality is linked with a reduced amount of RNA quantity hybridized to the array paralleled by a declined total signal level. Increasing degrees of saturation reduce the 3 ′ /5 ′ intensity gradient. The probe set contained the target gene at a position near the 3 ′ end of the corresponding transcripts. A side-by-side plot was produced by the function plotAffyRNAdeg (Figure 2), and the function summary of AffyRNAdeg produced a single summary statistic for each array in the batch (Table 1), presenting an assessment of the severity of RNA-degradation and significance level. We focused to sort out important factors in the cDNA dataset such as sequences biases, RNA degradation and quality of RNA to ensure the reliability of the dataset to identify transcriptional variations of the original samples. We standardized the sample handling procedures through normalization and assessed the optimal RNA reliability threshold using statistical and algorithm discrimination measures.

Identifying Differentially Expressed Genes (DEGs)
Our analyses revealed a total of 50 DEGs from Shiga effected human THP1-monocytic cells between 03 untreated and 06 treated to replicate samples, including 01 up-regulated (2%) and 49 down-regulated (98%) DEGs ( Table 2). For generalized linear models, the "cv.glm" function estimated the k-fold cross-validation prediction error. The dispersion parameter for Gaussian is 0.03268, which indicates the confidence models. The similar delta value of 0.032715 with 10 K-folds estimation is obtained when we used the LOOCV method (during raw cross-validation and then during adjusted cross-validation). The significant codes (0.1, 0.01, 0.001, and 0.05) with minimum deviance residuals are assessed (Table 3).

Identifying Cancer-Related Genes
Among differentially expressed genes, 7-cancer related genes were obtained including TSLP, IL6, GBP1, CD274, TNFSF13B, OASL, and pnpla3 after mapping (p < 0.00005) with OMIM and Cancer Genetics Web database. The references to the role of these genes in cancer have been shown in Table 4 using PubMed database.

Cluster and Functional Enrichment Analysis of DEGs
Clustering analysis has recognized to be helpful to understand gene function, gene regulation, cellular processes, and subtypes of cells. The cluster analysis of significantly expressed 7cancer related genes was studied with Euclidean distance (Figure 3). The genetic expression of Shiga effected cell samples is distinguished from the untreated replicates, indicating that obvious differences existed between the two groups (treated and untreated). The analysis showed that 6 downregulated and 1 up-regulated genes were significantly enriched. The total number of enriched terms was calculated with their significant cutoff parameters ( Table 5). The immune response (p < 1.68E-05), regulation of lymphocyte's proliferation (p < 0.001017), T-cell activation (p < 0.002007), cytokine activity (p < 0.004485), leukocyte homeostasis (p < 0.021097) and patatin-like phospholipase domain containing-3 were expressively enriched among dysregulated genes. These gene products are related to immune response system and signaling pathways, signaling of chemical mediators and regulation of a cell defense system ( Figure 4A). The dysregulation of these genes has been found associated with Kaposi sarcoma, cerebral malformation, familial polyposis, and other clinical phenotypes ( Figure 4B).

Protein-Protein Interaction Network
In PPI network, a total of 215 nodes was retrieved from STRING and HAPPI databases. In this network, IL6_Human as a source protein is connected with a number of target proteins (interactors) (28%) that have been reported in several cancers followed by the PD1L1_HUMAN (7%) (Figure 5).
Overexpression of PLPL3_HUMAN is associated with 4% of interactors that are responsible for various types of cancers. It has been shown that IL6_HUMAN is involved in the progression of several cancers, extremely responsible for leukemia, breast cancer and lymphoma (Figure 6). PD1L1_HUMAN has also been found the connection in a number of cancer cases, including lung cancer, leukemia, Hodgkin lymphoma, and head/neck cancers. The observed least association of OASL_HUMAN is with lymphoma and skin cancer cases.

De Novo Prediction of Regulatory Motifs
Cancer-related differentially expressed genes were used for overrepresentation and de novo analysis. Following transcription factors (TF) class were found in these DEGs: REL, Stat, TRP-CLUSTER, and bZIP. They are related to six target matrices in the Jaspar core profile. Opossum-tool found RELA as a topranked matrix under both the Fisher-score (8.746e −02 ) and the z-score (18.14) followed by the STAT1 which satisfied the thresholds for significant over-representation (Figure 7). This shows that genes containing conserved RELA, STATI, IRFI, NF-kappaB, PEND, HLF, REL, CEBPA, DI_2, and NFKB1 binding sites (Supplementary Figure 2) are indeed over-represented in cancer-related DEGs and this over-representation can be recuperated computationally. The rank of parameter settings for the analysis of oPOSSUM-tool could reflect properties how a specific TF controls its targets.

DISCUSSION
Monocytes/macrophages are key cellular components of the human immune system. Upon exposure to Shiga toxins produced by Shigella dysenteriae serotype 1 and enterohaemorrhagic Escherichia coli, these cells secrete chemical mediators, which initiate in?ammation and activate of immune system signals (Harrison et al., 2005). It has been known that Stx alters the genetic expression and several signaling pathways triggered by the toxin modify the expression of several genes in mammalian cells. The comparative cDNA studies showed that Stx toxins modulated the transcriptional response of human macrophagelike THP-1 cells (Leyva-Illades et al., 2010) associated with a number of abnormalities. In another study, it was observed that Stxs dysregulated the expression of Gb3/CD77 of lung carcinoma cell lines (Uchida et al., 1999). Although it is known that these toxins alter the genetic expression and affect the cellular division and apoptosis pathway (Elmore, 2007), the underlying mechanisms and its association with cancer development remain unclear. The purpose of the current study was to investigate differentially expressed genes of human THP1 cells in response to Stxs and the association of these genes with cancer development.
In this study, we found 50 DEGs, and 49 of which were downregulated and only a single gene being up-regulated by Stxs. Among these, 7-genes were mapped with cancer development that fell into the following functional categories: Immune response, regulation of lymphocyte's proliferation, cytokine activity and patatin-like phospholipase domain containing-3. These DEGs are linked with important biological pathways that are related to immune response signaling pathways, signaling of chemical mediators and regulation of cell defense system. Dysregulation of these genes can cause Kaposi sarcoma, familial polyposis, and other clinical phenotypes. These findings appear different from earlier microarray analyses performed by Keepers et al. (2006), which revealed Stxs-encoded proteins with functions related to transcriptional regulation, cell proliferation, and cell cycle regulation. Cell proliferation, cell morphogenesis, immune response signaling, anti-apoptosis, and regulation of cell death are strongly associated with the down-regulated genes, whereas cell cycle process, glycerolipid metabolism, and Patatinlike phospholipase activity are enriched in the up-regulated genes. These factors are known to cancer progression (Hunter and Pines, 1994;DeRisi et al., 1996).
Stxs have been reported to up-regulate the expression of cytokines and chemokines from human primary blood monocytes and to transform monocytic cell lines in vitro, and to affect cytokine expression via several mechanisms, including the activation of MAPK cascades leading to increased cytokine gene transcriptional activity and improved translation initiation efficiency (Harrison et al., 2005). However, in our analysis, we found that Stxs down-regulated IL6 gene, which is a potent inducer of the acute phase response and plays an important role in the final differentiation of B-cells into FIGURE 5 | Protein-Protein Interaction network of cancer-related differential expressed genes. Each node represents the protein, while lines indicate the interaction edges. The association of these source proteins and their interactors (target proteins) in cancers are shown in a bar graph. Blue and red color indicate down and up-regulated genes, respectively.
Ig-secreting cells responsible for lymphocyte and monocyte differentiation (Hirano et al., 1986). The downregulation of IL-6 and other DEGs is associated with highly malignant mammary carcinomas (Fontanini et al., 1999). Some genes such as MMP11 (extra-cellular matrix proteins), XRCC1 (DNA repair), VEGF (regulator of angiogenesis), Cyclin D1 (cell cycle regulators) and tumor-suppressor genes (Semaphorin 3B, WNT-5A) have been found as down-regulated genes during lung cancer progression (Campioni et al., 2008). Recently it has studied that Stxs elicit a proinflammatory response, including the elevated production of TNF-α, IL-1β, IL-8, IL-6, and growth-related oncogene α (Groα) both in vitro and in vivo (Lee et al., 2013). Furthermore, our analysis reveals that toxin may affect IL6expression via a number of mechanisms, including JAK-STAT cascade, cytosolic DNA-sensing, TNF signaling pathway and transcriptional misregulation in cancer cells. IL6 is associated with the progression of multiple cancers, including leukemia, lymphoma, breast cancer, kidney, and lungs cancer. TNFSF13B encodes Tumor Necrosis Factor (Ligand) Superfamily, Member 13b, also referred to as B-Cell-Activating Factor. It is important in the proliferation and differentiation of B-cells that can act as a transcription factor for its own gene in association with NF-κB p50 (Yu et al., 2000). Similarly, CD274 is another gene that is involved in the co-stimulatory signal, essential for T-cell proliferation and production of IL10 and IFNG. The dysregulation of this gene inhibits T-cell proliferation and cytokine production causing a number of cancer including lung, leukemia, Hodgkin lymphoma, and head/neck cancers (Dong et al., 1999). We found cytokines and chemokines, including IL-8, CSF2, GRO-1, GRO-2, and GRO-3, are affected by toxin exposure (Leyva-Illades et al., 2010). In addition, expression of TSLP, IL6, GBP1, CD274, TNFSF13B, OASL, and pnpla3 involved in the immune system, cell cycle, T-cell activation and lymphocyte's proliferation were also altered in THP1 cells by toxin treatment. Dysregulation of these genes cause increased apoptosis, cellular division and cancerous activities (Dong et al., 1999;Andreoli et al., 2014;Barooei et al., 2015;Hengeveld and Kersten, 2015;Ibsen et al., 2015;Manda et al., 2015;Ali et al., 2016). We performed motif enrichment to predict the transcription factor binding sites and over-represented regions of these cancer-related genes to better understand the complex biological mechanism. Several motifs, including RELA, STATI, IRFI, NF-kappaB, PEND, HLF, REL, CEBPA, DI_2, and NFKB1 were found in the gene set for factors that play critical roles in modeling required for immune response system, cellular cycles, chemical mediation and cell growth.
It should be noted that the publicly accessible microarray dataset studied here were subjected to rigorous statistical algorithms. Thus, there is a high level of confidence in the list of differentially expressed genes altered by Stxs. The present study did not evaluate an independent prognostic factor but suggest that differentially expressed genes should be evaluated as  potential biomarkers for further study to uncover their possible effects in the diagnosis, treatment, and prognosis of cancers.

CONCLUSION
This study identified differentially expressed genes in human THP1-monocytic cells upon treatment by bacterial Shiga toxins, and many of these genes are involved in the pathophysiological carcinomas. Stxs not only inhibited the protein synthesis of target cells but even induce the genetic expression and activate the multiple signaling pathways. Thus, cells treated with Stxs may lead to progression of multiple cancers. These findings may provide a valuable framework for developing diagnostic biomarkers and treatment strategies for cancer. Further molecular and real-time studies are warranted to investigate the role of the identified genetic elements in cancer development.

AUTHOR CONTRIBUTIONS
SM and JG performed the research and wrote the paper; TN and XW analyzed the data; BB and XY collected data and proofread the article; JC supervised and helped perform the research and commented on the draft manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.