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.

Introduction

Shiga toxins (Stxs) are bacterial cytotoxins mainly produced by enteric pathogens. It is a family of related toxins with two major groups, Stx1 and Stx2, encoded by genes belonging to the lambdoid prophages genome (Friedman and Court, 2001). The sources of Stxs includes Shigella dysenteriae, the Shiga toxigenic group of Escherichia coli (STEC) including serotypes O157:H7, O104:H4, and other enterohemorrhagic E. coli (EHEC) (Beutin, 2006; Spears et al., 2006) responsible for millions of cases of severe dysentery, food poisoning, gastroenteritis and bowel necrosis (Mims et al., 1993; Vogt and Dippold, 2005; Todar, 2007; World Health Organization, 2016), and this situation is worse in underdeveloped countries.

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” (Torgersen et al., 2010). 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 (Sandvig et al., 2010; 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 caspases- can 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 pro-inflammatory 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 THP1-monocytic 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.

Materials and Methods

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)¹. 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 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 statistic). To assess the quality of RNA in samples, AffyRNAdeg, summaryAffyRNAdeg, and plotAffyRNAdeg Bioconductor packages were used for degradation analysis (Affymetrix, 1999, 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 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)² and Cancer GeneticsWeb database (www.cancer-genetics.org)³ to filter the significant cancer-related genes. The roles of sorted genes in cancer were evaluated through reported experimental studies using PubMed⁴ 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., 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)⁵. 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 cancer-related, 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 over-represented regulatory motif of target matrices (Pavesi et al., 2004) in promoter regions of expressed genes using default parameters.

Results

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.

FIGURE 1

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₂(I₁) – log₂(I₂), A = 1/2 (log₂(I₁)+log₂(I₂)), where I₁ is the intensity of the array studied, and I₂ is the intensity of a “pseudo”-array that consists of the median across arrays.

FIGURE 2

Figure 2. Side-by-side plot produced by plotAffyRNAdeg representing 5′-3′ trend to assess the severity of degradation of RNA and significance level.

TABLE 1

Table 1. A summary statistic for each array in the batch, assessing the severity of RNA degradation and significance level.

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

TABLE 2

Table 2. List of 50-differentially expressed genes.

TABLE 3

Table 3. k-fold cross validation by bioconductor “boot” package using dispersion parameter of Gaussian family.

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.

TABLE 4

Table 4. The differentially expressed cancer-related genes curated from OMIM and Cancer Genetics. Web databases and their role has been referenced from PubMed.

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 7-cancer 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 down-regulated 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).

FIGURE 3

Figure 3. Cluster analysis of 7 cancer-related DEGs with Euclidean distance (Binning method: Quantile). Green corresponds to a small distance and Red to a large distance. Lines indicate the boundaries of the clusters in the level of the tree. Dotted lines indicate down or up-regulated genes.

TABLE 5

Table 5. Functional and GO* analysis of 7-differentially expressed genes.

FIGURE 4

Figure 4. Functional enrichment analysis of expressively cancer-related DEGs using FunRich annotation tool (A) Biological pathways analysis (B) clinical phenotypes for GEGs.

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.

FIGURE 5

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.

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.

FIGURE 6

Figure 6. The role of 7 DEGs in various cancers with Absolute Pearson correlation. Red color indicates the maximum involvement of a gene in various cancers and Blue color represents the least association. This association was curated using OMIM, Cancer Genomic Web, and PubMed databases.

De Novo Prediction of Regulatory Motifs

Cancer-related differentially expressed genes were used for over-representation 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 top-ranked matrix under both the Fisher-score (8.746e⁻⁰²) 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.

FIGURE 7

Figure 7. Over-representation of cancer-related DEGs using oPOSSUM with 80% Matrix score.

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 macrophage-like 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 down-regulated 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 Patatin-like 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 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 IL6-expression 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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2018.00380/full#supplementary-material

Supplementary Figure 1. Framework our study indicates the effects of Shiga toxins on cancer progressions.

Supplementary Figure 2. Transcription factor binding sites (TFBS) of cancer-related differentially expressed genes with their Jaspar core profiles. Top 10% of conserved regions (min. conservation 70%) are using oPOSSUM tool.

Footnotes

1. ^Gene Expression Omnibus: http://www.ncbi.nlm.nih.gov/geo/. Accessed 2016-05.25.

2. ^NCBI: http://www.ncbi.nlm.nih.gov/omim. Accessed 2015-08.28.

3. ^Cancer Genetics Web: http://www.cancer-genetics.org/index.htm. Accessed at Thursday Sep 10, 2016.

4. ^PubMed: http://www.ncbi.nlm.nih.gov/pubmed. Accessed 2015-08-28.

5. ^National Cancer Database: http://www.facs.org. Accessed 2016-01.01.

References

Affymetrix (1999). Affymetrix Microarray Suite User Guide, 4th Edn. Santa Clara, CA: Affymetrix.

Affymetrix (2001). Affymetrix Microarray Suite User Guide, 5th Edn. Santa Clara, CA: Affymetrix.

Ali, M., Yopp, A., Gopal, P., Beg, M. S., Zhu, H., Lee, W., et al. (2016). A variant in PNPLA3 associated with fibrosis progression but not hepatocellular carcinoma in patients with hepatitis C virus infection. Clin. Gastroenterol. Hepatol. 14, 295–300. doi: 10.1016/j.cgh.2015.08.018

CrossRef Full Text | Google Scholar

Andreoli, M., Persico, M., Kumar, A., Orteca, N., Kumar, V., Pepe, A., et al. (2014). Identification of the first inhibitor of the GBP1:PIM1 interaction. Implications for the development of a new class of anticancer agents against paclitaxel resistant cancer cells. J. Med. Chem. 57, 7916–7932. doi: 10.1021/jm5009902

PubMed Abstract | CrossRef Full Text | Google Scholar

Barooei, R., Mahmoudian, R. A., Abbaszadegan, M. R., Mansouri, A., and Gholamin, M. (2015). Evaluation of thymic stromal lymphopoietin (TSLP) and its correlation with lymphatic metastasis in human gastric cancer. Med. Oncol. 32:217. doi: 10.1007/s12032-015-0653-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B Stat. Methodol. 57, 289–300.

Google Scholar

Beutin, L. (2006). Emerging enterohaemorrhagic Escherichia coli, causes and effects of the rise of a human pathogen. J. Vet. Med. B Infect. Dis. Vet. Public Health 53, 299–305. doi: 10.1111/j.1439-0450.2006.00968.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolstad, B. M., Irizarry, R. A., Astrand, M., and Speed, T. P. (2003). A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19, 185–193. doi: 10.1093/bioinformatics/19.2.185

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandelli, J. R., Griener, T. P., Laing, A., Mulvey, G., and Armstrong, G. D. (2015). The effects of shiga Toxin 1, 2 and their subunits on cytokine and chemokine expression by human macrophage-like THP-1 cells. Toxins 7, 4054–4066. doi: 10.3390/toxins7104054

PubMed Abstract | CrossRef Full Text | Google Scholar

Campioni, M., Ambrogi, V., Pompeo, E., Citro, G., Castelli, M., Spugnini, E. P., et al. (2008). Identification of genes down-regulated during lung cancer progression: a cDNA array study. J. Experiment. Clinic. Cancer Res. 27:38. doi: 10.1186/1756-9966-27-38

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, A. W., Bhatti, M., DuPont, H. L., Nataro, J. P., Carlin, L. G., and Okhuysen, P. C. (2017). Clinical features and molecular epidemiology of diarrheagenic Escherichia coli pathotypes identified by fecal gastrointestinal multiplex nucleic acid amplification in patients with cancer and diarrhea. Diagn. Microbiol. Infect. Dis. 89, 235–240. doi: 10.1016/j.diagmicrobio.2017.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J. Y., Mamidipalli, S., and Huan, T. (2009). HAPPI: an online database of comprehensive human annotated and predicted protein interactions. BMC Genomics 10:S16. doi: 10.1186/1471-2164-10-S1-S16

PubMed Abstract | CrossRef Full Text | Google Scholar

Cline, M. S., Smoot, M., Cerami, E., Kuchinsky, A., Landys, N., Workman, C., et al. (2007). Integration of biological networks and gene expression data using cytoscape. Nat. Protoc. 2, 2366–2382. doi: 10.1038/nprot.2007.324

PubMed Abstract | CrossRef Full Text | Google Scholar

DeRisi, J., Penland, L., Brown, P. O., Bittner, M. L., Meltzer, P. S., Ray, M., et al. (1996). Use of a cDNA microarray to analyse gene expression patterns in human cancer. Nat. Genet. 14, 457–460. doi: 10.1038/ng1296-457

PubMed Abstract | CrossRef Full Text | Google Scholar

DesRochers, T. M., Kimmerling, E. P., Jandhyala, D. M., El-Jouni, W., Zhou, J., Thorpe, C. M., et al. (2015). Effects of shiga toxin type 2 on a bioengineered three-dimensional model of human renal tissue. Infect. Immun. 83, 28–38. doi: 10.1128/IAI.02143-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, H., Zhu, G., Tamada, K., and Chen, L. (1999). B7-H1, a third member of the B7 family, co-stimulates T-cell proliferation and interleukin-10 secretion. Nat. Med. 5, 1365–1369. doi: 10.1038/70932

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisen, M. B., Spellman, P. T., Brown, P. O., and Botstein, D. (1998). Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. U.S.A. 95, 14863–14868. doi: 10.1073/pnas.95.25.14863

CrossRef Full Text

Elmore, S. (2007). Apoptosis: a review of programmed cell death. Toxicol. Pathol. 35, 495–516. doi: 10.1080/01926230701320337

PubMed Abstract | CrossRef Full Text | Google Scholar

Fontanini, G., Campani, D., Roncella, M., Cecchetti, D., Calvo, S., Toniolo, A., et al. (1999). Expression of interleukin 6 (IL-6) correlates with estrogen receptor in human breast carcinoma. British J. Cancer. 80, 579–584. doi: 10.1038/sj.bjc.6690394

CrossRef Full Text | Google Scholar

Fraser, M. E., Chernaia, M. M., Kozlov, Y. V., and James, M. N. G. (1994). Crystal structure of the holotoxin from Shigella dysenteriae at 2.5 Å resolution. Nat. Struct. Biol. 1, 59–64. doi: 10.1038/nsb0194-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, D. I., and Court, D. L. (2001). Bacteriophage lambda: alive and well and still doing its thing. Curr. Opin. Microbiol. 4, 201–207. doi: 10.1016/S1369-5274(00)00189-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujita, A., Sato, J. R., de Rodrigues Lde, O., Ferreira, C. E., and Sogayar, M. C. (2006). Evaluating different methods of microarray data normalization. B Bioinformatics 7:469. doi: 10.1186/1471-2105-7-469

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, L. M., van den Hoogen, C., van Haaften, W. C., and Tesh, V. L. (2005). Chemokine expression in the monocytic cell line THP-1 in response to purified Shiga toxin 1 and/or lipopolysaccharides. Infect. Immun. 73, 403–412. doi: 10.1128/IAI.73.1.403-412.2005

CrossRef Full Text | Google Scholar

Hattori, T., Watanabe-Takahashi, M., Shiina, I., Ohashi, Y., Dan, S., Nishikawa, K., et al. (2016). M-COPA, a novel Golgi system disruptor, suppresses apoptosis induced by Shiga toxin. Genes Cells 21, 901–906. doi: 10.1111/gtc.12386

PubMed Abstract | CrossRef Full Text | Google Scholar

Hengeveld, P. J., and Kersten, M. J. (2015). B-cell activating factor in the pathophysiology of multiple myeloma: a target for therapy? Blood Cancer J. 5:e282. doi: 10.1038/bcj.2015.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirano, T., Yasukawa, K., Harada, H., Taga, T., Watanabe, Y., Matsuda, T., et al. (1986). Complementary DNA for a novel human interleukin (BSF-2) that induces B lymphocytes to produce immunoglobulin. Nature 324, 73–76. doi: 10.1038/324073a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho Sui, S. J., Mortimer, J. R., Arenillas, D. J., Brumm, J., Walsh, C. J., Kennedy, B. P., et al. (2005). oPOSSUM: identification of over-represented transcription factor binding sites in co-expressed genes. Nuc. Acids Res. 33, 3154–3164. doi: 10.1093/nar/gki624

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang da, W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, T., and Pines, J. (1994). Cyclins and cancer II: cyclin D and CDK inhibitors come of age. Cell 79, 573–582. doi: 10.1016/0092-8674(94)90543-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ibsen, M. S., Gad, H. H., Andersen, L. L., Hornung, V., Julkunen, I., Sarkar, S. N., et al. (2015). Structural and functional analysis reveals that human OASL binds dsRNA to enhance RIG-I signaling. Nuc. Acids Res. 43, 5236–5248. doi: 10.1093/nar/gkv389

PubMed Abstract | CrossRef Full Text | Google Scholar

Keepers, T. R., Psotka, M. A., Gross, L. K., and Obrig, T. G. (2006). A murine model of HUS: Shiga toxin with lipopolysaccharide mimics the renal damage and physiologic response of human disease. J. Am. Soc. Nephrol. 17, 3404–3414. doi: 10.1681/ASN.2006050419

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. S., Kim, M. H., and Tesh, V. L. (2013). Shiga toxins expressed by human pathogenic bacteria induce immune responses in host cells. J. Microbiol. 51, 724–730. doi: 10.1007/s12275-013-3429-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Leyva-Illades, D., Cherla, R. P., Galindo, C. L., Chopra, A. K., and Tesh, V. L. (2010). Global transcriptional response of macrophage-like THP-1 cells to shiga toxin type 1. Infect. Immun. 78, 2454–2465. doi: 10.1128/IAI.01341-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Long, T., Lu, Y., Ouyang, Q., and Tang, C. (2004). The yeast cell-cycle network is robustly designed. Proc. Natl. Acad. Sci. U.S.A. 101, 4781–4786. doi: 10.1073/pnas.0305937101

PubMed Abstract | CrossRef Full Text | Google Scholar

Lingwood, C. A. (2003). “Shiga toxin receptor glycolipid binding. Pathology and utility,” in Methods in Molecular Medicine E. coli Shiga Toxin Methods and Protocols, Vol. 73, eds D. Philpott and F. Ebel (Totowa, NJ: Humana Press, Inc.), 165–186.

Google Scholar

Lingwood, C. A., Binnington, B., Manis, A., and Branch, D. R. (2010). Globotriaosyl ceramide receptor function—where membrane structure and pathology intersect. FEBS Lett. 584, 1879–1886. doi: 10.1016/j.febslet.2009.11.089

PubMed Abstract | CrossRef Full Text | Google Scholar

Lukyanenko, V., Malyukova, I., Hubbard, A., Delannoy, M., Boedeker, E., Zhu, C., et al. (2011). Enterohemorrhagic Escherichia coli infection stimulates Shiga toxin-1 macropinocytosis and transcytosis across intestinal epithelial cells. Am. J. Physiol. Cell Physiol. 301, C1140–C1149. doi: 10.1152/ajpcell.00036.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Malyukova, I., Murray, K. F., Zhu, C., Boedeker, E., Kane, A., Patterson, K., et al. (2009). Macropinocytosis in Shiga toxin 1 uptake by human intestinal epithelial cells and transcellular transcytosis. Am. J. Physiol. Gastrointest. Liver Physiol. 296, G78–G92. doi: 10.1152/ajpgi.90347.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Manda, K. R., Tripathi, P., His, A. C., Ning, J., Ruzinova, M. B., Liapis, H., et al. (2015). NFATc1 promotes prostate tumorigenesis and overcomes PTEN loss-induced senescence. Oncogene 35, 3282–3292. doi: 10.1038/onc.2015.389

PubMed Abstract | CrossRef Full Text | Google Scholar

Melton-Celsa, A. R. (2014). Shiga toxin (Stx) classification, structure, and function. Microbiol. Spectr. 2, 1–20. doi: 10.1128/microbiolspec.EHEC-0024-2013

CrossRef Full Text | Google Scholar

Mims, C. A., Playfair, J. H. L., Roitt, I. M., Williams, W. (1993). Medical Microbiology, 1st Edn. Mosby, MO: Cambridge University Press.

Muhammad, S. A., Ahmed, S., Ali, A., Huang, H., Wu, X., Yang, X. F., et al. (2014). Prioritizing drug targets in Clostridium botulinum with a computational systems biology approach. Genomics 104, 24–35. doi: 10.1016/j.ygeno.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Muhammad, S. A., Fatima, N., Syed, N.-I-H., Wu, X., Yang, X. F., and Chen, J. Y. (2015). MicroRNA expression profiling of human respiratory epithelium affected by invasive Candida infection. PLoS ONE 10:e0136454. doi: 10.1371/journal.pone.0136454

PubMed Abstract | CrossRef Full Text | Google Scholar

Nam, D., and Kim, S. Y. (2008). Gene-set approach for expression pattern analysis. Brief. Bioinformatics 9, 189–197. doi: 10.1093/bib/bbn001

PubMed Abstract | CrossRef Full Text | Google Scholar

Naz, A., Awan, F. M., Obaid, A., Muhammad, S. A., Paracha, R. Z., Ahmad, J., et al. (2015). Identification of putative vaccine candidates against Helicobacter pylori exploiting exoproteome and secretome: a reverse vaccinology based approach. Infec. Genet. Evol. 32, 280–291. doi: 10.1016/j.meegid.2015.03.027

PubMed Abstract | CrossRef Full Text | Google Scholar

NCBI (2015). NCBI: http://www.ncbi.nlm.nih.gov/gene/3624. Accessed August 29, 2015.

Obenchain, V., Lawrence, M., Carey, V., Gogarten, S., Shannon, P., and Morgan, M. (2014). Variant annotation: a bioconductor package for exploration and annotation of genetic variants. Bioinformatics 30, 2076–2078. doi: 10.1093/bioinformatics/btu168

CrossRef Full Text | Google Scholar

Pathan, M., Keerthikumar, S., Ang, C. S., Gangoda, L., Quek, C. M. J., Williamson, N. J., et al. (2015). FunRich: a standalone tool for functional enrichment analysis. Proteomics 15, 2597–2601. doi: 10.1002/pmic.201400515

PubMed Abstract | CrossRef Full Text

Pavesi, G., Mereghetti, P., Mauri, G., and Pesole, G. (2004). Weeder web: discovery of transcription factor binding sites in a set of sequences from co-regulated genes. Nucleic Acids Res. 32, W199–W203. doi: 10.1093/nar/gkh465

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramegowda, B., Samuel, J. E., and Tesh, V. L. (1999). Interaction of Shiga toxins with human brain microvascular endothelial cells: cytokines as sensitizing agents. J. Infect. Dis. 180, 1205–1213. doi: 10.1086/314982

PubMed Abstract | CrossRef Full Text | Google Scholar

Ripley, B. (2010). Package ‘Boot’. Available online at: CRAN: http://cran.rproject.org/web/packages/boot/index.html (Accessed September 10, 2015).

Ryoo, H. D., and Bergmann, A. (2012). The role of apoptosis-induced proliferation for regeneration and cancer. Cold Spring Harb. Perspect. Biol. 4:a008797. doi: 10.1101/cshperspect.a008797

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandvig, K., Bergan, J., Dyve, A., Skotland, T., and Torgersen, M. L. (2010). Endocytosis and retrograde transport of Shiga toxin. Toxicon 56, 1181–1185. doi: 10.1016/j.toxicon.2009.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandvig, K., Grimmer, S., Lauvrak, S. U., Torgersen, M. L., Skretting, G., van Deurs, B., et al. (2002). Pathways followed by ricin and Shiga toxin into cells. Histochem. Cell Biol. 117, 131–141. doi: 10.1007/s00418-001-0346-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Schweppe, C. H., Bielaszewska, M., Pohlentz, G., Friedrich, A. W., Büntemeyer, H., Schmidt, M. A., et al. (2008). Glycosphingolipids in vascular endothelial cells: relationship of heterogeneity in Gb3Cer/CD77 receptor expression with differential Shiga toxin 1 cytotoxicity. Glycoconj. J. 25, 291–304. doi: 10.1007/s10719-007-9091-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Seymour, G. (1993). Predictive Inference. New York, NY: Chapman and Hall.

Spears, K. J., Roe, A. J., and Gally, D. L. (2006). A comparison of enteropathogenic and enterohaemorragic E. coli pathogenesis. FEMS Microbiol. Lett. 255, 187–202. doi: 10.1111/j.1574-6968.2006.00119.x

CrossRef Full Text | Google Scholar

Stricklett, P. K., Hughes, A. K., Ergonul, Z., and Kohan, D. E. (2002). Molecular basis for up-regulation by inflammatory cytokines of Shiga toxin 1 cytotoxicity and globotriaosylceramide expression. J. Infect. Dis. 186, 976–982. doi: 10.1086/344053

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Kuhn, M., Simonovic, M., Roth, A., Minguez, P., et al. (2011). The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 39, D561–D568. doi: 10.1093/nar/gkq973

PubMed Abstract | CrossRef Full Text | Google Scholar

Tesh, V. L. (2010). Induction of apoptosis by Shiga toxins. Future Microbiol. 5, 431–453. doi: 10.2217/fmb.10.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorpe, C. M., Smith, W. E., Hurley, B. P., and Acheson, D. W. (2001). Shiga toxins induce, superinduce, and stabilize a variety of C-X-C chemokine mRNAs in intestinal epithelial cells, resulting in increased chemokine expression. Infect. Immun. 69, 6140–6147. doi: 10.1128/IAI.69.10.6140-6147.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Todar, K. (2007). Pathogenic E. coli. Online Textbook of Bacteriology. University of Wisconsin-Madison Department of Bacteriology. Cornell University, Ithaca, New York, NY.

Torgersen, M. L., Engedal, N., Bergan, J., and Sandvig, K. (2010). The intracellular journey of shiga toxins. Open Toxinol. J. 3, 3–12. doi: 10.2174/1875414701003020003

CrossRef Full Text | Google Scholar

Troyanskaya, O., Cantor, M., Sherlock, G., Brown, P., Hastie, T., Tibshirani, R., et al. (2001). Missing value estimation methods for DNA microarrays. Bioinformatics 17, 520–525. doi: 10.1093/bioinformatics/17.6.520

PubMed Abstract | CrossRef Full Text | Google Scholar

Tusher, V. G., Tibshirani, R., and Chu, G. (2001). Significance analysis of microarrays applied to the ionizing radiation response. Proc. Natl. Acad. Sci. U.S.A. 98, 5116–5121. doi: 10.1073/pnas.091062498

CrossRef Full Text

Uchida, H., Kiyokawa, N., Taguchi, T., Horie, H., Fujimoto, J., and Takeda, T. (1999). Shiga toxins induce apoptosis in pulmonary epithelium-derived cells. J. Infec. Diseas. 180, 1902–1911. doi: 10.1086/315131

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Kar, N. C., Monnens, L. A., Karmali, M. A., and van Hinsbergh, V. W. (1992). Tumor necrosis factor and interleukin-1 induce expression of the verocytotoxin receptor globotriaosylceramide on human endothelial cells: implications for the pathogenesis of the hemolytic-uremic syndrome. Blood 80, 2755–2764.

PubMed Abstract | Google Scholar

Vogt, R. L., and Dippold, L. (2005). Escherichia coli O157:H7 outbreak associated with consumption of ground beef, June–July 2002. Public Health Rep. 120, 174–178. doi: 10.1177/003335490512000211

PubMed Abstract | CrossRef Full Text | Google Scholar

World Health Organization (2016). Diarrheal Diseases: Shigellosis. Initiative for Vaccine Research (IVR). New York, NY: World Health Organization (Accessed August 4, 2016).

Yoon, D., Yi, S. G., Kim, J. H., and Park, T. (2004). Two-stage normalization using background intensities in cDNA microarray data. BMC Bioinformatics 5:97. doi: 10.1186/1471-2105-5-97

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Boone, T., Delaney, J., Hawkins, N., Kelley, M., Ramakrishnan, M., et al. (2000). APRIL and TALL-I and receptors BCMA and TACI: system for regulating humoral immunity. Nat. Immunol. 1, 252–256. doi: 10.1038/79802

PubMed Abstract | CrossRef Full Text | Google Scholar