A Comprehensive Transcriptomic Analysis of Arsenic-Induced Bladder Carcinogenesis

Abstract

Arsenic (sodium arsenite: NaAsO2) is a potent carcinogen and a known risk factor for the onset of bladder carcinogenesis. The molecular mechanisms that govern arsenic-induced bladder carcinogenesis remain unclear. We used a physiological concentration of NaAsO2 (250 nM: 33 µg/L) for the malignant transformation of normal bladder epithelial cells (TRT-HU1), exposed for over 12 months. The increased proliferation and colony-forming abilities of arsenic-exposed cells were seen after arsenic exposure from 4 months onwards. Differential gene expression (DEG) analysis revealed that a total of 1558 and 1943 (padj < 0.05) genes were deregulated in 6-month and 12-month arsenic-exposed TRT-HU1 cells. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that cell proliferation and survival pathways, such as the MAPK, PI3K/AKT, and Hippo signaling pathways, were significantly altered. Pathway analysis revealed that the enrichment of stem cell activators such as ALDH1A1, HNF1b, MAL, NR1H4, and CDH1 (p < 0.001) was significantly induced during the transformation compared to respective vehicle controls. Further, these results were validated by qPCR analysis, which corroborated the transcriptomic analysis. Overall, the results suggested that stem cell activators may play a significant role in facilitating the arsenic-exposed cells to gain a survival advantage, enabling the healthy epithelial cells to reprogram into a cancer stem cell phenotype, leading to malignant transformation.

Keywords: KEGG and GO enrichment; arsenic; bladder cancer; differential gene expression; stem cells.

Figure 1

Determination of cell viability, colony formation ability of transforming cells, and identification of differentially expressed genes from RNA seq datasets. (a) The cell viability was measured by the intake of trypan blue by dead cells in control and transforming cells (0, 2, 4, 6, 8, 10, and 12 months). (b) Colony-forming assays on 0, 2, 4, 6, 8, 10, and 12 months transforming TRT-HU1 cells were performed. Colonies were manually counted after being stained with crystal violet. All experiments were performed in triplicate. One-way ANOVA with multiple comparison tests was used to calculate the statistical significance between different experimental groups *, p < 0.05; and ****, p < 0.0001, # Not significant (c) PCA plot analysis of differentially expressed genes in three datasets. The volcano plot analysis of differentially expressed genes between (d) 6-month against 0-month, (e) 12-month against 0-month, and (f) 12-month against 6-month is plotted on the X axis, and false discovery rate (FDR) significance is plotted on the Y axis (-log10 scale). The grey dots represent no significant change, red dots represent logFC of >2 and FDR < 0.05, and blue dots represent logFC < −2 and FDR < 0.05.

Figure 2

Gene set enrichment analysis (GSEA) of 6-month and 12-month As-exposed cells compared to control with positive enrichment score. GSEA curves of the list of differentially expressed genes in 6-month vs. 0-month and 12-month vs. 0-month. Gene set enrichment at the top of the ranked list is indicated by a positive enrichment score (ES).

Figure 3

KEGG pathway and GO term analysis of differentially expressed genes based on RNA-seq data for 6-month and 12-month As-exposed cells compared to control cells. Each bubble’s color and size correspond to the amount of differentially expressed mRNAs enriched in GO or KEGG pathway, respectively. The cut-off for choosing to GO and KEGG keywords was p < 0.05.

Figure 4

Expression validation of selected genes (RNA-seq datasets) using qRT-PCR. Heatmap analysis showing the expression level of top 10 upregulated genes (a) in RNA-seq datasets. The expression for upregulated genes (b) was determined by quantitative RT-PCR in successive stages of transforming cells (0, 6, and 12 months). Expression levels were normalized to β-actin. Error bars indicate the standard deviation of triplicates. One-way ANOVA was used to calculate the statistical significance between vehicle control and treatment at each time point. (**** p < 0.0001).

Figure 5

Expression validation of selected genes (RNA-seq datasets) using qRT-PCR. Heatmap analysis showing the expression level of top 10 downregulated genes (a) in RNA-seq datasets. The expression for downregulated genes (b) was determined by quantitative RT-PCR in successive stages of transforming cells (0, 6, and 12 months). Expression levels were normalized to β-actin. Error bars indicate the standard deviation of triplicates. One-way ANOVA was used to calculate the statistical significance between vehicle control and treatment at each time point. * p < 0.05, *** p < 0.001, and **** p < 0.0001.

Figure 6

Sphere formation ability of transforming TRT-HU1 cells (0, 2, 4, 6, 8, 10, and 12 months). Control and transforming (2 × 10³ cells) were plated in 12-well ultra-low attachment plates, and allowed for 14 days, then the total number of spheres per well were counted. The horizontal bars represent the mean number of spheres, *** p < 0.001, and **** p < 0.0001.