Dysregulation of mRNA expression by hsa-miR-186 overexpression in arsenic-induced skin carcinogenesis

Abstract

Dysregulated miRNA expression contributes to development of arsenic-induced cutaneous squamous cell carcinoma (cSCC). hsa-miR-186 (miR-186) is overexpressed in arsenical cSCC tissues as well as in preclinical cell line model of arsenical cSCC. Simultaneous miR-186 overexpression and chronic inorganic trivalent arsenite (iAs; 100 nM) exposure transformed human HaCaT cell line preferentially over miR-186 overexpression or iAs exposure alone. Both iAs and miR-186 regulate the expression of wide range of mRNA targets. However, how their interaction impacts the transcriptome-wide mRNA expression landscape ushering in cancer is unknown. We performed longitudinal RNA-seq analysis in passage-matched HaCaT cell clones (±miR-186 overexpression) with simultaneous chronic iAs exposure (0/100 nM) at 12 and 29 weeks. We determined the impact of each factor and their interaction towards differential gene expression and pathway dysregulation employing two different statistical approaches (t-statistic and 2-factor ANOVA). We show that a core set of pathways are dysregulated deterministically irrespective of the statistical approach chosen, possibly representing necessary changes for transformation. The data suggest that each clonal line could take a unique route to dysregulate this core set of pathways necessary for transformation, highlighting the possible role of stochasticity in cancer development. Evidence is presented to sift the strengths and weaknesses of each statistical methodology in providing biological understanding of events that play crucial roles in carcinogenesis in large datasets with multiple contributing variables.

Keywords: Arsenic; Carcinogenesis; Skin cancer; mRNA; miRNA.

Fig. 1.

Experimental setup showing four experimental groups (±miR-186; ±iAs) at each time point (12 and 29 weeks) as well as the pairwise t-statistic comparisons performed. Details about each pairwise comparison is provided in Table S1.

Fig. 2.

Differential gene expression analysis by pairwise t-statistic comparison. A. Venn diagram showing the number of differentially expressed mRNAs (p < 0.05) in each group and overlaps at 12 weeks. B. Venn diagram showing the number of differentially expressed mRNAs (p < 0.05) in each group and overlaps at 29 weeks. C. Bar graph representing the number of induced and suppressed mRNAs (p < 0.05) for each pairwise comparison at 12 and 29 weeks. Numbers are also provided in Table S2. Venn diagrams were generated using Bioinformatics & Evolutionary Genomics venn diagram tool (https://bioinformatics.psb.ugent.be/webtools/Venn/). Detailed description about DEG1-DEG10 is provided in Table S1.

Fig. 3.

Differential gene expression analysis by two-factor ANOVA analysis. A. Venn diagram showing the number of differentially expressed mRNAs (p < 0.05) explained by each factor and overlaps at 12 weeks. B. Venn diagram showing the number of differentially expressed mRNAs (p < 0.05) explained by each factor and overlaps at 29 weeks. C. Bar graph representing the number of induced and suppressed mRNAs (p < 0.05) explained by each factor at 12 and 29 weeks. Numbers are also provided in Table S3. Venn diagrams were generated using Venn Diagram Plotter program (https://pnnl-comp-mass-spec.github.io/Venn-Diagram-Plotter/).

Fig. 4.

Predicted dysregulated pathways populated by differentially expressed mRNAs (p < 0.01; fold change > ± 10 %) identified by pairwise t-statistic comparison. A. Bar graph representing the number of dysregulated pathways at 12 and 29 weeks. Dysregulated pathways are defined as -log(p-value) > 1.3. Dysregulated pathways without a z-score or with z-score between − 1 and + 1 are considered “Undetermined”. Activated pathways defined as -log(p-value) > 1.3; z-score > 1; inhibited pathways defined as -log(p-value) > 1.3; z-score < − 1 B. Heat Map of predicted activated/inhibited pathways based on differentially expressed mRNAs at 12 weeks stratified by specific pairwise comparisons (DEG1, DEG2 and DEG5). C. Heat Map of predicted activated/inhibited pathways based on differentially expressed mRNAs at 29 weeks stratified by specific pairwise comparisons (DEG6, DEG7 and DEG10). The color code bar on the top refers to the z-score values. Absence of a bar (represented by white) signifies that the pathway was not predicted to be activated or inhibited at that time point. The pathways are presented in the same order as in Table S9. Pathway numbers do not signify the same pathways at both time points.

Fig. 5.

Predicted dysregulated pathways populated by differentially expressed mRNAs (p < 0.01; fold change > ± 10 %) identified by two-factor ANOVA analysis. A. Bar graph representing the number of dysregulated pathways at 12 and 29 weeks. Dysregulated pathways are defined as -log(p-value) > 1.3. Dysregulated pathways without a z-score or with z-score between − 1 and + 1 are considered “Undetermined”. Activated pathways defines as -log(p-value) > 1.3; z-score > 1; inhibited pathways defined as -log(p-value) > 1.3; z-score < − 1 B. Heat Map of predicted activated/inhibited pathways based on differentially expressed mRNAs at 12 weeks stratified by factors (iAs, miR-186, interaction). C. Heat Map of predicted activated/inhibited pathways based on differentially expressed mRNAs at 29 weeks stratified by stratified by factors (iAs, miR-186, interaction). The color code bar on the top refers to the z-score values. Absence of a bar (represented by white) signifies that the pathway was not predicted to be activated or inhibited at that time point. The pathways are presented in the same order as in Table S10.