Original Research ARTICLE
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
- 1Institute of Biopharmaceutical Informatics and Technologies, Wenzhou Medical University, Wenzhou, China
- 2Wenzhou Medical University 1st Affiliated Hospital, Wenzhou, China
- 3Institute of Molecular Biology and Biotechnology, Bahauddin Zakariya University, Multan, Pakistan
- 4Department of Computer and Information Science, Purdue University Indianapolis, Indianapolis, IN, United States
- 5Institute for Systems Biology, Seattle, WA, United States
- 6Department of Microbiology and Immunology, Indiana University School of Medicine, Indianapolis, IN, United States
- 7Informatics Institute, School of Medicine, The University of Alabama at Birmingham, Birmingham, AL, United States
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 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)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 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)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 PubMed4 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)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 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.
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. 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 = log2(I1) – log2(I2), A = 1/2 (log2(I1)+log2(I2)), where I1 is the intensity of the array studied, and I2 is the intensity of a “pseudo”-array that consists of the median across arrays.
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. 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 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. 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. 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.
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. 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. 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−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.
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.
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.
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.
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.
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
Ripley, B. (2010). Package ‘Boot’. Available online at: CRAN: http://cran.rproject.org/web/packages/boot/index.html (Accessed September 10, 2015).
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
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
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
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
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
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
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
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
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
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.
Keywords: cDNA dataset, Shiga toxin, differential expression, cancers, enrichment analysis, TFBS
Citation: Muhammad SA, Guo J, Nguyen TM, Wu X, Bai B, Yang XF and Chen JY (2018) 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. Front. Microbiol. 9:380. doi: 10.3389/fmicb.2018.00380
Received: 03 January 2018; Accepted: 20 February 2018;
Published: 13 March 2018.
Edited by:Andrea Masotti, Bambino Gesù Ospedale Pediatrico (IRCCS), Italy
Reviewed by:Grzegorz Wegrzyn, University of Gdansk, Poland
Vernon L. Tesh, Texas A&M University College of Medicine, United States
Copyright © 2018 Muhammad, Guo, Nguyen, Wu, Bai, Yang and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.