Transcriptomic analysis reveals that combinations of vitamins A, D2, and D3 have synergistic effects in HCT-116 colon cancer cells by altering the expression of genes involved in multiple canonical pathways including apoptosis, regulation of the epithelial mesenchymal transition and immunity
April 12, 2021
Gail Mahady
Pinal Kanabar
Nina Los
Temitope Lawal
Shitalben Patel
Mark Maienschein-Cline
Zarema Arbieva
Research and medical notice. The vitamin A and D combinations described here have been studied only in cells and animals, not in people. This research is shared for education and is not medical advice or a recommendation to take vitamins A or D, which can be harmful at high doses. Talk with your doctor before starting a supplement or changing any treatment.
Funding and acknowledgements: We would like to express our sincere appreciation for the generous educational sponsorship from the First Analysis Institute of Integrative Studies (GBM); a research grant from the Regenstein Foundation (GBM); a Postdoctoral fellowship award from the Schlumberger Foundation (TOL/GBM); and a Raman Post-Doctoral Fellowship by the University Grants Commission, Govt. of India to (NAR/GBM). PK and MMC were supported in part by University of Illinois at Chicago (UIC) Center for Clinical and Translational Science (CCTS) through Grant UL1TR002003. The contents are solely the responsibility of the authors and do not necessarily represent the official views of the funding agencies. We would like to express our sincere thanks and appreciation to F. Oliver Nicklin at the First Analysis Institute of Integrative Studies for the theoretical discussions, as well as his assistance with hypothesis generation and the direction of the work.
How to cite: Kanabar P., Los N., Lawal T., Patel S., Maienschein-Cline M., Arbieva Z., Mahady G. Transcriptomic analysis reveals that combinations of vitamins A, D2, and D3 have synergistic effects in HCT-116 colon cancer cells by altering the expression of genes involved in multiple canonical pathways including apoptosis, regulation of the epithelial mesenchymal transition and immunity. Functional Foods in Health and Disease 2021; 11(4): 154-178. https://doi.org/10.31989/ffhd.v11i4.784
Abstract
Integrated systems biology approaches suggest that combinations of nutrients may be more effective against cancer due to the large number of signaling pathways associated with cancer initiation and promotion. In a previous work, we have reported that combinations of vitamins A (as all trans-retinoic acid, ATRA), D2, and D3 act synergistically to induce apoptosis in colon and gastric cancer cells.
In this work, we use whole-genome transcriptomic profiling to detect gene expression changes using RNA-seq to more comprehensively investigate the biological pathways affected by the combination of vitamin D2, D3 and ATRA. HCT-116 colon cancer cells were harvested, RNA was isolated and RNA-seq libraries were prepared using a Universal Plus mRNASeq kit. Sequencing was carried out on NovaSeq 6000. General quality-control metrics were obtained using FastQC and raw reads were aligned to human reference genome hg38 using STAR and BWA MEM. ENSEMBL genes were quantified using FeatureCounts, and differential expression statistics were computed using EdgeR.
Specific gene expression was validated using qPCR. Transcriptomic analysis showed that of 26,313 genes analyzed, the expression of 8,402 genes was significantly altered (4030 up-regulated and 4373 down-regulated, FDR<0.05) in the treated cells, of which, 3,621 genes were differentially expressed (fold change <-1 or >+1 and an FDR <0.05). IngenuityÒ Pathway analysis revealed the involvement of 97 canonical pathways, with the top pathways including: mechanisms of cancer, apoptosis, myc-mediated apoptosis, regulation of the epithelial mesenchymal transition, and immunity.
INTRODUCTION
Colorectal carcinoma (CRC) is a leading cause of cancer mortality worldwide and is responsible for ~900,000 deaths annually [1-3]. Symptoms usually occur later in disease progression; therefore, early detection is important for overall survival. Advancements in the understanding of CRC pathology have significantly increased the range of treatment options for patients and increased survival rates, especially in patients without metastatic disease [1, 2]. Colon cancer is also one of the best genetically characterized cancers, and specific somatic mutations in oncogenes and tumor suppressor genes have been identified [2].
Epidemiological and experimental studies suggest that populations living in the South and Southwest have reduced cancer incidence and mortality as compared with their Northern counterparts, due in part to more extensive exposure to solar ultraviolet (UV) B radiation and vitamin D levels [1-5]. Over the years, this hypothesis has since gained considerable support due to the results of geographical, ecological, observational, mechanistic, and clinical studies [2-10]. The risk of breast, colorectal, gastric, prostate, and other cancers is increased in patients with a vitamin D deficiency, and for gastric cancer, lower vitamin D levels are associated with a poor prognosis [9-15]. Colorectal cancer patients with an elevated serum vitamin D level and intake had a better prognosis, and showed improved survival and reduced recurrence [7,13- 19]. Vitamin A and all trans retinoic acid (ATRA) are reported to inhibit the growth of breast, gastric and colon cancers, and leukemia [20-30].
Over the past 10 years, new hypotheses have integrated systems biology approaches and suggest that combinations of nutrients may be more effective against cancer due to the large number of signaling pathways associated with cancer initiation and promotion [31-33]. In particular, whole-genome transcriptomic profiling to detect gene expression changes using RNA-seq has been used to develop a more comprehensive view of the biological pathways involved in the effects of natural compounds and nutrients in vitro. To this end, we have been investigating the combinations of vitamins A and D on colon, breast and gastric cancers. In our previous work, we reported that vitamins D2 and D3 inhibited the growth of cultured colon and gastric cell lines, but, when combined with all-trans-retinoic acid (ATRA), the IC50 was significantly reduced, indicating synergistic effects [34]. In the present study, we utilized the human tumor cell line HCT116 as we have previously reported that this colon cancer cell line was most susceptible to the combination of vitamins A and D, with an IC50 of 1 µM in vitro. We have used this human colon cell line as it can be used for future xenograph studies in rodent models and to investigate CRC-induced cachexia, and liver metastases. In vitro, HCT-116 cells treated with the combination of vitamins A and D exhibited increased activities of caspase 3/7 and 8 at 4-8 hours after treatment, and significant increases in the Bax/Bcl-2 ratio, indicating induction of apoptosis. In addition, the combination reduced the time to cytotoxicity from 72 hours to ~24 hrs [34].
In the current study, we analyzed the transcriptome of HCT-116 colon cancer cells treated with the combination of vitamin A and D (all-trans-retinoic acid, ATRA+D2+D3; 1:1:1) at the IC50 concentration (1 µM) using RNA-seq to measure differentially expressed genes (DEGs) and possible novel canonical pathways to determine molecular mechanisms of action. Combined bioinformatics and Ingenuity® Pathway Analyses (IPA; Qiagen, CA) were used to determine the effects of this combination on differential gene expression (DEG) in HCT-116 cells, as well as canonical pathways and causal networks involved in these mechanisms.
MATERIALS AND METHODS
Cell culture and maintenance: The HCT-116 (ATCC® CCL-247™) colon cancer cells were obtained from American Type Culture Collection and cultured and maintained as we have previously described [34]. Briefly, HCT 116 cells were cultured in McCoy’s 5a Medium (Gibco, Life Technology, Grand Island, NY, USA) supplemented with 10% FBS and 1% Penicillin/Streptomycin. All cells were incubated at 37°C in a humidified atmosphere of 5% CO2 and 95% air. At 80% confluence, cells were harvested by adding 0.25% trypsin/EDTA and counted by means of Trypan blue and hemocytometer. The cells were then re-suspended at an appropriate concentration and plated in 96 or 6 well plates for cellular assays. Cell were treated with ATRA+D2+D3 at the IC50 concentration (1µM) for 4 hours, then harvested and RNA isolated. The cells were treated for 4 hours as our previous study suggested that Bcl-2 protein gene expression and caspase activation were already increased by 4 to 8 hours after treatment [34].
RNA Isolation and Quality Control: HCT-116 cells were seeded at a density of 1.0 – 1.2 x 10⁶ cells in 1 ml per well of a 6-well plate and incubated overnight for cell attachment. Total RNA and proteins were extracted from cells after treatment for 4 hours at the IC50 concentration using Trizol (ThermoFisher Scientific, Waltham, MA, USA). The RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Tech., Palo Alto, CA). Each of the six RNA samples was used to generate parallel RNA-Seq and qPCR data. Vehicle-treated control cells were included as negative controls. RNA samples were quantified using NanoDrop™ One Spectrophotometer (Thermo Scientific) and analyzed for integrity using Agilent 4200 TapeStation and RNA Screen Tape (PN: 5067-5576). Levels of remaining DNA were checked on a subset of samples using a Qubit fluorometer (Invitrogen). DNA amounts did not exceed 10% of the total amount of RNA.
RNAseq Library Preparation: Libraries were prepared with the Universal Plus mRNASeq kit (Tecan, PN: 0520B-A01) according to manufacturer's instructions. In brief, 250 ng of total RNA was used to make each library. Library construction steps included poly(A) RNA selection, RNA fragmentation and double-stranded cDNA generation using a mixture of random and oligo(dT) priming, followed by end repair to generate blunt ends, ligation of UDI adaptors, strand selection, and PCR amplification to produce the final library. The number of amplification cycles was determined by qPCR using a small pre-amplification library aliquot for a subset of samples, in this case it was determined to be 15 PCR cycles. Final amplified libraries were purified using Agencourt RNAClean XP Bead (Beckman Coulter, PN: A63987).
RNAseq Library Validation and Quantification: Library fragment size distribution was confirmed by electrophoresis the 2200 TapeStation system using D1000 ScreenTape (Agilent, PN: 5067-5582, 5067-5589), and was determined to range between 264 and 294 nt. Individual libraries were quantified on a Qubit 2.0 Fluorometer with the Qubit dsDNA HS Assay Kit (Life Technologies, PN: P7581). Libraries were normalized and pooled to a final 10 nM concentration.
Sequencing: Concentration of the final library pool was confirmed by qPCR using KAPA Library Quantification Kit (Roche, KAPA Biosystems, PN: KK4873). The library pool was sequenced on a MiniSeq (Illumina) in order to check and adjust accordingly the proportions of individual libraries. Sequencing was carried out on NovaSeq 6000, SP flowcell, 2x50 nt reads.
Bioinformatics Statistic and Database Annotation: Processing of the raw data was performed by the University of Illinois at Chicago Core for Research Informatics (UICCRI). General quality-control metrics for next-generation sequencing data was obtained using FastQC. Raw reads were aligned to human reference genome hg38 in a splice-aware manner using STAR and BWA MEM [35]. The ENSEMBL database (www.ensembl.org) was used to analyze the molecular function and biological processes of differentially expressed genes (DEGs) detected from RNA-seq. ENSEMBL genes were quantified using FeatureCounts [36-37]. Differential expression statistics were computed using edgeR [38-39] on raw expression counts with the exactTest function. P-values were adjusted for multiple testing using the false discovery rate (FDR; q value) correction of Benjamini and Hochberg [40]. DEGs were categorized into functional clusters based on Gene Ontology (GO) function enrichment analysis using edgeR.
Ingenuity® Pathway Analysis (IPA): All processed array files were inspected for the quality metrics such as average signal present, signal intensity of species-specific housekeeping genes, relative signal intensities of labeling controls, absolute signal intensities of hybridization controls, and across-array signal distribution plots [41]. To identify the molecular signaling pathways associated with the observed changes in transcription levels, DEGs were further analyzed according to the predicted protein function using the Ensembl database. Expression and statistical data were uploaded into the Ingenuity® Pathway Analysis software (Qiagen, USA). Up-regulation or down-regulation of genes, filtered by fold change of <-1 and > 1 and a false discovery rate (q value) < 0.05, were used to identify canonical pathways and mapped to IPA’s knowledge bases, and their relevant biological functions, networks and pathways. The IPA knowledge base contains > 5M observations extracted from the biomedical literature and integrated third party-databases [41].
Quantitative polymerase chain reaction to validate candidate genes: QPCR was performed to validate the transcriptome levels of fragments per kilobase of transcript per million mapped reads (FPKM) generated by RNA-Seq. Briefly, total RNA was reverse transcribed and amplified using Power SYBR Green RNA- to-CT 1- step kit (Applied Biosystems, Foster City, CA, USA) as described by manufacturer with the Step One Plus Real Time PCR System (Applied Biosystems, Foster City, CA, USA). Briefly, each reaction was performed in triplicate in a 10 μl volume containing Power SYBR Green RT-PCR Mix (2X), 200 nM of each primer, RT Enzyme Mix (125X) and 100 ng RNA. The cycling conditions were as follows: 48°C for 30 min, 95°C for 10 min, followed by 50 cycles of 95°C for 15 sec and 60°C for 1 min. PCR reaction specificity was confirmed by melt curve analysis at 95°C for 15 sec, 60°C for 15 sec, and 95°C for 15 sec. The qRT- PCR assays were conducted in triplicate and ΔCt values (Ct for GADPH or β-actin - Ct for the test gene) were calculated for each RNA sample. The Student t-test was used to analyze whether there was a significant difference between the mean ΔCt for the control vs. the treated HCT-116 groups, with a threshold significance level p < 0.05. The fold change in gene expression was calculated as 2ΔΔCt (ΔΔCt = ΔCt of the treatment group - ΔCt of control group).
The specific primers used were (forward (F) and reverse (R) primers, respectively): GADPH: (F) 5′- ACCACAGTCCATGCCATCAC-3′ and (R) 5′-TCCACCACCCTGTTGCTGTA-3′; and Caspase 3: (F) 5’- AGGAGGGACGAACACGTCT-3’ and (R) 5’-CAAAGAAGGTTGCCCCAATCT-3’. Caspase 9: (F) 5’- AGGAGGGACGAACACGTCT-3’ and (R) 5’-CAAAGAAGGTTGCCCCAATCT-3’. Caspase 7: (F) 5’- GGACCGAGTGCCCACTTATC-3’ and (R) 5’-TCGCTTTGTCGAAGTTCTTGTT-3’. Caspase-10: (F) 5′- GCCTGTTCCTGTGATGTGGAG-3’ and (R) 5′-TGCCCACAGACATTCATACAGTTTC-3′. Caspase 8: (F) 5’- CTGGGAAGGATCGACGATTA-3’ and (R) 5’-CATGTCCTGCATTTTGATGG-3’. Bcl-2: (F) 5’-CGCATCAGGAAGGCTAGAGT-3’ and (R) 5’-AGCTTCCAGACATTCGGAGA-3’; Bax: (F) 5’-TGCCAGCAAACTGGTGCTCA-3’ and (R) 5’- GCACTCCCGCCACAAAGATG-3’; β-actin: (F) 5’TGACGTGGACATCCGCAAAG-3’ and (R) 5’-CTGGAAGGTGGACAGCGAGG-3’; HDAC3: (F) 5’-CAAGCCATACCAGGCCTCCCAGC-3’ and (R) 5’-GAGATGCGCCTGTGTAACGCGAG-3’; HDAC1: (F) 5’-CTACTACGACGGGGATGTTGG-3’ and (R) 5’-GAGTCATGCGGATTCGGTGAG-3’; SIRT1: (F) 5’- TGCTGGCCTAATAGAGTGGCA-3’ and (R) 5’-CTCAGCGCCATGGAAAATGT-3’; SIRT3: (F) 5’- GTCGGGCATCCCTGCCTCAAAGC-3’ and (R) 5’-GGAACCCTGTCTGCCATCACGTCAG-3’; p53: (F) 5’- AAGTCTGTGACTTGCACGTACTCC-3’ and (R) 5’-GTCATGTGCTGTGACTGCTTGRTAG-3’; BIM: (F) 5′- ATGTCTGACTCTGACTCTCG-3′ and (R) 5′-CCTTGTGGCTCTGTCTGTAG-3′; NOXA: (F) 5′-TCCTGAGCAGAAGAGTTTGG-3′ and (R) 5′-GGAGATGCCTGGGAAGAAGG-3′; PUMA: (F) 5′-ATGGCGGACGACCTCAAC-3′ and (R) 5′- AGTCCCATGAAGAGATTGTACATGAC-3′; SMAD7: (F) 5’-GGCTGTGTTGCTGTGAATC-3’ and (R) 5’- GGTATCTGGGAGTAAGGAGGAG-3’; NOTCH: (F) 5′CACCCATGAC-CACTACCCAGTT3′ and (R) 5′CCTCGGACCAATCA-GAGATGTT3′; JAG1: (F)5’-TGCCAAGTGCCAGGAAG-3’ and (R) 5’-GCCCCATCTGGTATCACACT-3’; CRMP1 (F) 5′- ATGCCCTGAGCAGACCTGAAGAGC-3′ and (R) 5′-AGTAATGGGTGCCATCGGTCCCCAG-3′; ADAM23: (F) 5′- TATGAGCAGCTGTCCACTCG-3′, and (R) 5′-CCCCAGCCTGTGCCCCCAAG-3’; AXIN2: (F) 5’- TTATGCTTTGCACTACGTCCCTCCA-3’and (R) 5’-CGCAACATGGTCAACCCTCAGAC-3’; IL12B: (F) 5’- GGCCAGTACACCTGTCACAA-3’ and (R) 5’-CAGAGAGTGTAGCAGCTCCG’3’ using published primers and NIH Primer BLAST [42-50]. PCR products were visualized by addition of ethidium bromide on 1% agarose gel. GADPH and β-actin were used to control and calibrate cDNA synthesis.
Availability of data: The RNA-seq raw datasets supporting the conclusions of this work are deposited and available at the National Center for Biotechnology Information Gene Expression Omnibus repository, number GSE160109, https://www.ncbi.nlm.nih.gov/geo.
Statistical analysis: All experiments were performed in triplicate and data are presented as mean plus standard deviation (S.D.). Statistical significance was determined using the Student’s t test using GraphPad version 8.4 (San Diego, CA). P-values were corrected for multiple testing using the false discovery rate (FDR; q value) correction of Benjamini and Hochberg [40]. FDR values of < 0.05 were considered significant.
RESULTS
RNA seq analysis: To develop a more comprehensive view of the biological pathways involved in the antiproliferative effects of the combination treatment ATRA+D2+D3 in HCT-116 colon cancer cells, we used whole-genome transcriptomic profiling using RNA-seq with validation by qPCR. This combination treatment inhibited HCT-116 colon cancer cell proliferation and induced apoptosis, as well as reduced kill times [34]. Cell samples from HCT-116 colon cancer cells treated with vehicle solvent (Control, DMSO 0.01%) or ATRA+D2+D3 (IC50 1 µM) in triplicate were used to isolate RNA and generate the RNA-seq data. Isolated RNA quality was measured using an Agilent 2100 Bioanalyzer and RNA samples had an RNA Integrity Number score ≥ 9.6-10 (Figure 1).
The RNA-seq libraries were prepared using a Universal Plus mRNASeq kit and library construction included poly(A) RNA selection, RNA fragmentation and double-stranded cDNA generation using a mixture of random and oligo(dT) priming, followed by end repair to generate blunt ends, ligation of UDI adaptors, strand selection, and PCR amplification to produce the final library. Amplified libraries were purified, library fragment size distribution was confirmed by electrophoresis and was determined to range between 264 and 294 nt. Libraries were normalized and pooled to a final 10 nM concentration.
Figure 1. RNA QC of total RNA. A1, B1, C1 = control vehicle solvent cells; D1, E1, F1 = treated HCT-116 colon cancer cells. RNA samples were quantified using NanoDrop™ One Spectrophotometer (Thermo Scientific) and analyzed for integrity using Agilent 4200 TapeStation and RNA Screen Tape. Levels of remaining DNA did not exceed 10% as determined Qubit fluorometer (Invitrogen).
Transcriptomic profiles of control versus treated HCT-116 colon cancer cells: RNA-seq data were used to prepare a principal component analysis (PCA) to determine the differences observed between the control (DMSO 0.01%) and treated (ATRA+D2+D3) HCT-116 cells. As can be seen in Figure 2A, the differences between the control and treated cells were distinct from each other with the treatment group distantly clustered from each other in the same PCA plot. These data confirm that there was a different gene expression pattern for a wide range of genes between the control and treatment groups. Transcriptomic analysis of 26,313 genes, showed that the expression of 8,402 genes was significantly altered, 4,030 were up-regulated, and 4,372 were down-regulated in treated versus control cells (FDR < 0.05; Figure 2B volcano plot). The two-dimensional volcano plot shows the pattern of 8,402 genes scattered outward away from the center of origin (Figure 2B). The stratified volcano plot visually emphasizes differentially expressed genes using colored points with external information [51], in this case, the up-regulated DEGs are represented in red and the down-regulated DEGs in blue. Of these, 3,621 DEGs had a fold change (FC) of <-1 and >1, (2,028 were up-regulated and 1,593 were down-regulated; FDR <0.05). Statistical analysis was performed using a modified t-test with the p-value corrected by the Benjamini-Hochberg algorithm [40].
Figure 2A and B. Overview of the differentially expressed genes from RNA-seq data from control HCT-116 colon cancer cells (0.01% DMSO) as compared with HCT-116 cells treated with 1.0 µg/ml of ATRA+D2+D3. (A) Principal component analysis (PCA) of the six RNA-seq data from HCT-116 colon cancer cells. Control (red dots; DMSO 0.01%) samples; treated samples (blue dots) treated with ATRA+D2+D3. (B) A stratified volcano plot showing the relationship between the magnitude of gene expression change (Log2 fold change, x axis) and the statistical significance of this change (-log10 of the false discovery rate (FDR), Q-value; y axis) and the average abundance of the 8,402 significantly expressed genes (FDR < 0.05) in treatment versus control HCT-116 cells. Each point represents the average value of one transcript and three replicates. The expression difference is considered significant at a q value of ≤ 0.05 (FDR). Points are colored according to their average expression in all data sets. Blue = down-regulated transcripts and red = up-regulated transcripts. FDR = false-discovery rate. Significantly up and down differentially expressed genes are highlighted in red and blue, respectively.
Top differentially expressed genes (DEGs) in treated HCT-116 cells: In addition to the volcano plot depicting the 8,402 DEGs, a heatmap was generated for the top 500 DEGs in HCT-116 colon cancer cells treated with ATRA+D2+D3, as compared with control HCT-116 cells (0.01% DMSO only; Figure 3). Using a false discovery rate (q < 0.01) and Log2FC of <-2 or > 2; 500 genes were differentially expressed (Figure 3 heatmap of the top 500 DEGs and hierarchical clustering of DEGs). Furthermore, 50 of the top DEGs are presented in Tables 2 and 3. Different gene expression patterns were observed between control and treated HCT-116 cells in the heatmap (Figure 3), indicating the reliability of the DEGs.
Figure 3. Heatmap visualization of the top 500 differentially expressed genes with gene ontology across data sets. Each column represents the value for one replicate of the experiment, three for control HCT-116 cells, and three for HCT-116 cells treated with 1 µg/ml (IC50) of ATRA+D2+D3 for 4 hours. The colors represent Z-scored Log2 CPM. To interpret the biological function of DEGs, an enrichment analysis of gene ontology to display their functional distribution is also depicted.
Table 1. Some of the top differentially expressed genes down-regulated after treatment of HCT-116 colon cancer cells with ATRA+D2+D3.
Table 2. Top differential expressed genes up-regulated after treatment of HCT-116 colon cancer cells with ATRA+D2+D3.
Table 2 (continued).
Canonical pathways affected in the treated HCT-116 colon cancer cells as compared with controls: We analyzed DEG profiles of ATRA+D2+D3 treated HCT-116 cells versus control cells using the Ingenuity Pathway Analysis (IPA) software (Qiagen, USA) to determine the relationship between highly significant genes and the most significant canonical pathways and biological networks involved. A total of 3,621 DEGs showed an FC of < -1 or >1 with an FDR < 0.05. The analysis was based on the ratio of the number of DEGs in the dataset to the total number of reference genes in corresponding pathways in the IPA knowledge bases. The IPA analysis used the Fisher’s exact test (p < 0.05) to determine the significant canonical pathways correlated with the DEGs from the RNA-seq data. Enrichment of DEGs in a particular pathway was determined on the basis of the Ensembl database. IPA core analysis showed that a total of 97 canonical pathways were significantly (FDR < 0.05) correlated with the changes in gene expression in treated HCT-116 colon cancer cells versus control cells. Figure 4 is a graph of the top 26 interconnected canonical pathways affected after treatment of HCT-116 cells with the ATRA+D2+D3 combination, the number of genes in each pathway, and p-values. Table 3 describes the top five most significant canonical pathways, along with their -log (P-value), and the ratio of affected genes over the total number of genes in a particular canonical pathway. The most significant canonical pathways impacted after treatment of HCT-116 cells with ATRA+D2+D3 were molecular mechanisms of cancer, apoptosis signaling, myc-mediated apoptosis signaling, death receptor signaling, and regulation of the epithelial-mesenchymal transition (EMT) and immunity (Table 3), indicating that this combination treatment affected multiple signaling pathways associated with cancer. Table 4 presents the most significant diseases, biological functions, as well as cellular and molecular functions impacted by treatment of the HCT-116 cells. Tables 4 and 5 also includes p-values and the number of molecules in the pathway. Not surprisingly, cancer was the most affected disease state, with genes related to cell death, proliferation, survival and movement significantly impacted.
Table 3. List of the top five most significant canonical pathways impacted in HCT-116 colon cancer cells treated with ATRA+D2+D3, their p-value, percentage overlap and the number of genes affected in the pathway over the total number of genes in the pathway. Data was analyzed using RNA-seq DEGs in Ingenuity® Pathways Analysis (IPA). Of the DEGs, 3,621 had an FC cutoff of < -1 and > 1 and a false discovery rate (FDR, q value) < 0.05. IPA analysis used the Fisher’s exact test (p < 0.05) to determine the significant canonical pathways correlated with the DEGs from the RNA-seq data. Enrichment of DEGs in a particular pathway was determined on the basis of the Ensembl database.
Table 4. Analysis of the top, most significant diseases and disorders, and molecular and cellular functions impacted after treatment of HCT-116 colon cancer cells with ATRA+D2+D3 using the Ingenuity® Pathway Analysis (IPA). Up-regulation or down-regulation of genes was filtered by fold change of > 1 and a false discovery rate (FDR) of <0.05, were used to identify and construct functional networks. These data were mapped to IPA knowledge bases, and their relevant biological functions, networks and pathways. IPA analysis used the Fisher’s exact test (p < 0.05) to determine the significant canonical pathways correlated with the DEGs from the RNA-seq data. Enrichment of DEGs in a particular pathway was determined on the basis of the Ensembl database.
Differentially expressed genes in IPA canonical pathways: To further investigate the molecular mechanisms involved in the effects of ATRA+D2+D3 in HCT-116 colon cancer cells, we analyzed DEGs involved in specific canonical pathways associated with mechanism of cancer, apoptosis, myc-mediated apoptosis, death-receptor signaling, and the EMT. In the mechanisms of cancer pathway, a total of 75 of the 391 genes in this pathway were differentially expressed in HCT-116 cells after ATRA+D2+D3 treatment (Table 3; Figure 5).
Of the 75 DEGs in this pathway, 31 were down-regulated and 44 genes associated with apoptosis were up-regulated. DEGs in the apoptosis and myc-mediated apoptosis signaling pathways were similarly affected. These two pathways represented 149 genes of which 44 were differentially expressed, 27 genes were up-regulated and 17 were down-regulated (Table 3). In the regulation of the epithelial-mesenchymal transition (EMT) canonical pathway, 42 of the 192 genes in the pathway were differentially expressed (Table 3; Figure 6).
Figure 4. Graphical analysis of the top 26 interconnecting canonical pathways impacted by treatment of HCT-116 colon cancer cells with ATRA+D2+D3 (IC50 1µM). Data was analyzed using RNA-seq DEGs in Ingenuity® Pathways Analysis (IPA). Data shows the number of genes in the pathway and the percentage of genes either upregulated or down-regulated in that pathway and the associated p-value representing 3,621 analyzed DEGs with an FC cutoff of < -1 or > 1 and a false discovery rate (FDR, q value) < 0.05 (2,028 were up-regulated and 1,593 were down-regulated) in control (vehicle solvent) versus treated HCT-116 colon cancer cells. IPA analysis used the Fisher’s exact test (p < 0.05) to determine the significant canonical pathways correlated with the DEGs from the RNA-seq data. Enrichment of DEGs in a particular pathway was determined on the basis of the Ensembl database. Green represents down-regulated genes and red represents up-regulated genes.
Figure 5. Mechanisms in cancer was the top canonical pathway impacted in HCT-116 colon cancer after treatment with ATRA+D2+D3. Genes highlighted in red/pink show significant up-regulation and genes highlighted in green/light green show significant down-regulation. Seventy-five of the 391 genes in this pathway were differentially expressed after treatment, however only those genes with a false discovery rate FDR < 0.01 are included in this figure.
Figure 6. Treatment of HCT-116 cells with ATRA+D2+D3 impacted the epithelial-mesenchymal transition (EMT) canonical pathway. Genes highlighted in red/pink show significant up-regulation and genes highlighted in green/light green show significant down-regulation. Forty-two of the 192 genes in this pathway were differentially expressed after treatment, however only those genes with a false discovery rate FDR < 0.01 are included in this figure
Validation of DEGs in the apoptosis canonical pathway by qPCR: Because the majority of alterations of genes were observed in the canonical pathway of mechanisms of cancer and apoptosis in ATRA+D2+D3 treated HCT-116 cells, qPCR was used to validate the changes in DEG expression from RNA-seq. Differentially expressed genes (associated with apoptosis) in treated versus control HCT-116 cells are presented in Figures 7A and B. The down-regulated and up-regulated DEGs were ranked based on fold change of < -1 or > 1 and an FDR < 0.01. A fold change lower < 1 indicated that gene expression was down-regulated by the treatment and a fold change > 1 suggested an increase in gene expression. Genes that showed a significant up-regulation included those of the Bcl-2 protein family, p53 signaling, NOTCH1, SMAD7 and multiple caspases (Casp). As can be seen in Figures 7A and B, gene expression was similar for the expression of specific DEGs from RNA-seq, with only Bcl-2 and Casp3 showing a difference between RNA-seq and qPCR. As can be seen, genes in both the extrinsic and intrinsic apoptotic pathways are up-regulated including Bim, Bax, p53, Puma, Noxa and caspases 3, 7, 8, 9 and 10. In as much as we observed a reduction in the expression of HDACs in HCT-116 cells treated with ATRA+D2+D3, we further compared the data from RNA-seq with qPCR results to determine the effect of treatment on the expression of specific HDACs and sirutins (SIRTs). Results from RNA-seq data showed that HDAC1 and 3, as well as SIRT1 and 3 were significantly down-regulated, supporting our previous observations, and these data corresponded with the qPCR data (Figure 7B). Interestingly, Adam23, CRMP1 and IL-12 were also significantly up-regulated in RNA-seq and validated by qPCR (Figures 7A).
In Figure 8, IPA analysis also showed an increase in gene expression in the canonical pathway for Communication between Innate and Adaptive Immune Cells in HCT-116 colon cancer after treatment with ATRA+D2+D3. Genes highlighted in red/pink show significant up-regulation and included Toll-like receptors, IL-15 and IL-12β (FDR < 0.01).
Figures 7A and B: A. RNA-seq and qPCR comparison analysis of relative gene expression of genes in the canonical pathways of in HCT-116 after treatment with the IC50 of the combination of ATRA+D2+D3 or vehicle control. A. Upregulated expression of genes in apoptosis pathways including Bcl-2 proteins, caspases and CRMP1, ADAM23 and IL12. B. Down-regulated gene expression of HDACs, SIRTs, SMAD7, and NOTCH1 in HCT-116 was determined after treatment with the IC50 of the combination of D2+D3+ATRA as compared with β-actin. RNA was isolated from the cells 4 hr after treatment using the Trizol method. Total RNA was reverse transcribed and amplified using Power SYBR Green RNA-to-CT one-step kit for gene quantification. For qPCR, the quantitation of gene expression was performed using β-actin gene as an endogenous control and relative to the control using the ΔΔCT calculation.
Figure 8. Communication between Innate and Adaptive Immune Cells was one of the canonical pathways impacted in HCT-116 colon cancer after treatment with ATRA+D2+D3. Genes highlighted in red/pink show significant up-regulation and genes highlighted in green/light green show significant down-regulation. Analysis of the top 1400 DEGs showed that Toll-like receptors, IL-15 and IL-12β were significantly (FDR < 0.01) upregulated, and the downstream effects that impact the adaptive, humoral and cellular immune responses.
DISCUSSION
The antiproliferative activities of vitamins A and D have been reported in a range of breast, gastric, leukemia, pancreatic, and colon cancer cells [20-30, 34, 53-70]. For example, Muto et al., reported that calcitriol (1,25-D3) induced apoptosis in colorectal adenomas by up-regulating pro-apoptotic proteins including Bax, down-regulating the anti-apoptotic protein Bcl-2, and inducing p53-independent apoptosis [67]. Hu et al., reported that all-trans-retinoic acid (ATRA) inhibited the proliferation of leukemia cells and induced cell cycle arrest through the induction of p53 signaling [70]. Furthermore, combinations of vitamin D3 with 5-fluorouracil or metformin had synergistic effects and increased cytotoxicity and apoptosis in breast and colon cancer cells [61-62]. In our previous work, we reported that vitamin D2 and D3 inhibited the growth of colon and gastric cancer cells with IC50s ranging from 19-56 μM, however when combined with ATRA (ATRA+D2+D3) the active concentration was reduced to IC50 1µM, indicating synergistic effects at very low concentrations in both gastric and colon cancer cells [34]. To put this concentration into a clinical perspective, the median serum concentrations of vitamin A in adults was reported at 1.7 to 2.2 μM/L, and biologically 1 μg of vitamin D is equivalent to 40 IU [71, 72]. In HCT-116 cells, the combination significantly inhibited the proliferation, induced apoptosis, and increased the activities of caspase 3/7 and 8, as well as the Bax/Bcl-2 ratio as compared with the vitamins alone [34]. The data further suggested that the mechanisms associated with this combination were pleiotropic, involving the regulation of multiple signaling pathways involved in cell survival/apoptosis (Bcl-2, Bax, caspase 3/7), and autophagy (mTOR) in colon and gastric cancer cells [34].
In this study, we have established a global transcriptomic profile to further elucidate the molecular and cellular mechanisms by which the combination of ATRA+D2+D3 inhibited the proliferation of HCT-116 colon cancer cells and to identify novel genes and biological pathways involved. We have used the clinically relevant HCT-116 human colon cancer cell line, as it can be used in rodent xenograph studies, and to investigate CRC-induced cachexia, as well as liver metastases. In the generated results, Transcriptomic analysis of 26,313 genes, showed that 3,470 had a significant fold change (FC) of < -1 and > 1 (2,406 were upregulated and 1,064 were down-regulated, FDR < 0.05), in treated versus control HCT-116 cells. Approximately 500 DEGs showed a FC of < -2, > 2, or -2, with an FDR of < 0.01. Ingenuity Pathway Analysis revealed that treatment impacted multiple canonical pathways including molecular mechanisms in cancer, apoptosis and myc-mediated apoptosis signaling supporting previous works [34, 67, 73]. IPA analysis further showed the involvement of the regulation of the EMT, and interestingly, communication between innate and adaptive immune cells.
In terms of canonical pathways, apoptosis signaling pathways were significantly up-regulated with the expression of 75 of the possible 391 genes altered. IPA data analysis showed that ATRA+D2+D3 treatment of HCT-116 cells increased the expression of Bax, Bim, p53, Noxa, Puma, and caspases 3/7, 8, 9 and 10 (Figure 5), all involved in both intrinsic and extrinsic apoptosis [68-69, 73]. The expression of Bcl-2 and Bax (members of the B-cell CLL/lymphoma 2 family of proteins), was significantly altered, supporting our previous work [34]. However, transcriptomic analysis further showed that BIM (a pro-apoptotic Bcl-2-like 11 protein, Figure 5) was also significantly up-regulated in ATRA+D2+D3 treated HCT-116 cells. BIM is known to activate other pro-apoptotic proteins, namely Bax/Bak, thereby increasing the permeability of the mitochondrial outer membrane leading to the release of cytochrome c and activating a cascade of aspartate specific cysteine proteases (caspases) to induce apoptosis [68-69, 73]. In addition, the data also showed significant up-regulation of multiple caspases, including caspase 3, 7, 8, 9 and 10. Previous studies have reported that ATRA alone, or in combination with vitamins D2 and D3 induced apoptosis by activation of cleaved caspases 3 and 7 [34, 70]. Our data showed that ATRA+D2+D3 treatment of HCT-116 cells up-regulated the expression of caspase 3/7 supporting this previous study, but further showed significant up-regulation of the expression of caspases 8, 9, and 10 that are involved in extrinsic apoptotic signaling [74-76]. Activation of caspase 8 directly initiates extrinsic apoptosis by cleaving and activating the executioner caspases 3, 6, and 7, but is also involved in the intrinsic apoptotic pathway, showing convergence between these two pathways [75-76]. In terms of caspase 10, while its role is not well understood, some evidence suggests that it may also be involved in intrinsic apoptosis triggered by cytotoxic drugs [76]. Further IP analysis showed the involvement that the Myc-mediated apoptosis canonical pathway, particularly up-regulation of the expression of the p53 upregulated modulator of apoptosis (PUMA), and Noxa, also pro-apoptotic proteins, and also members of the Bcl-2 family [77-80]. Furthermore, our data showed the down-regulation of SMAD7 gene expression. SMAD7 expression is associated with increased inflammation and is associated with growth and survival of colorectal cancer (CRC) cells, and SMAD7 is over-expressed by tumor cells in both sporadic CRC and colitis-associated CRC [81-82]. Thus, ATRA+D2+D3 treatment significantly alters the expression of important genes in multiple apoptotic pathways in HCT-116 colon cancer cells in the favor of apoptosis.
Transcriptomic analysis further showed that ATRA+D2+D3 treatment of HCT-116 cells down-regulated the expression of histone deacetylases (HDACs). HDACs are involved in the epigenetic regulation of gene expression via histone modification and chromatin remodeling, and are categorized into class I (HDAC1, 2, 3 and 8), class II (HDAC4, 5, 6, 7, 9 and 10), class III (Sirutins 1-7) and class IV (HDAC11) [83-85]. Compounds that inhibit HDACs are emerging as novel drugs for the treatment of cancer [83]. Our data demonstrate that ATRA+D2+D3 treatment of HCT-116 cells inhibited gene expression of HDACs (confirmed by qPCR) including HDAC 1 and 3, as well as SIRT1 and 3, suggesting that these compounds may reduce cancer cell growth and increase apoptosis by impacting the epigenome and acting as HDAC antagonists in HCT-116 cells.
Beyond apoptosis, ATRA+D2+D3 treatment of HCT-116 cells significantly impacted genes in the regulation of the epithelial-mesenchymal transition (EMT) canonical pathway (Table 1; Figure 6), In the EMT, epithelial cancer cells undergo significant morphological changes characterized by a transition from an epithelial phenotype to a mesenchymal phenotype allowing cancer cells to migrate, invade, and be more resistant to apoptosis [86-88]. Our data show that ATRA+D2+D3 treatment of HCT-116 cells induced the differential expression of 42 of the 192 genes in the EMT pathway (Table 1), including AXIN2 (-7.1 FC), NOTCH1 (-2.5 FC), JAGGED (-2.8 FC), all important proteins involved in the EMT. High-levels of expression of NOTCH1 and its ligand JAGGED-1 are associated with a poorer prognosis in breast, bladder, colon, leukemia, and prostate cancers, and high levels of JAGGED-1 expression in a subset of clinically localized tumors were significantly associated with recurrence, and metastatic disease [86-87]. NOTCH signaling is a key regulator of the EMT, and its activation in endothelial cells results in the phenotypical and functional changes that are consistent with mesenchymal transformation [86-88]. Our data indicate that NOTCH1 and JAGGED are down-regulated by the ATRA+D2+D3 treatment of HCT-116 cells, suggesting that the combination may reduce NOTCH signaling and the EMT, thereby being a potential preventative treatment for metastatic disease. ATRA+D2+D3 treatment of HCT-116 cells also up-regulated the expression of genes in the EMT, including the collapsin response mediator protein-1 (CRMP1, +10 FC) that acts as an invasion and metastasis suppressor [89]. Over-expression of CRMP1 inhibited the EMT in vivo, indicating that up-regulation of CRMP1 may be a preventative strategy for reducing cancer metastasis [89]. Furthermore, ADAMs (A Disintegrin And Metalloproteases), particularly ADAM23 was also significantly up-regulated (+7.5 FC) in HCT-116 cells after treatment. Reduced expression of ADAM23 was reported to increase cell migratory capability, increase cancer cell proliferation, and was associated with a poor prognosis and reduced survival rates [90-92]. Thus, our data suggest that ATRA+D2+D3 treatments of HCT-116 cells alter the expression of multiple genes that play a role in the regulation of the EMT by down-regulating the expression of genes that promote metastasis and up-regulating genes involved in the suppression of metastasis.
Finally, ATRA+D2+D3 treatments appeared to significantly alter the expression of genes in the “Communication between Innate and Adaptive Immune Cells” canonical pathway, particularly a significant increase in the gene expression of interleukin -12 (IL-12; Figures 7A and 8). IL-12 is one of the members of the interleukins, a group of cytokines that mediate communication between cells, regulate cell growth, differentiation, and stimulate immune responses [93]. The IL-12 family of proteins provide a bridge between innate and adaptive immune systems by priming CD4+ T cells to differentiate into cytokine-producing T-helper subsets and memory T cells [93]. IL-12 cytokines further modulate the cellular pathways needed for proper functioning of the immune system. Thus, these data suggest that ATRA+D2+D3 may further impact cancer cell growth by improving immune functions.
CONCLUSIONS
Our previous study, and that of others now support the hypothesis that combinations of vitamins A and D are more effective than either vitamin alone in reducing the proliferation of epithelial cell cancers. At the cellular level, vitamins D2 and D3 bind to the vitamin D receptor (VDR) and vitamin A (retinoids) binds to the retinoic acid receptor (RXR), then these two receptors form a heterodimer, the VDR-RXR complex, that increases or reduces the transcription of target genes. This suggests that the synergistic effects of vitamins A and D combinations may be mediated through increased formation of the VDR-RXR complex and enhanced gene transcription. In this work, transcriptomic analysis using RNA sequencing provided an efficient high-throughput method to robustly characterize differential gene expression and provide molecular insight into the mechanisms by which compounds impact multiple signaling pathways in cancer cells. RNA-seq used in combination with qPCR showed that the combination of ATRA+D2+D3 significantly altered the transcriptome of HCT-116 colon cancer cells by impacting multiple canonical pathways associated with apoptosis, EMT, and immunity. The genes for members of the Bcl-2 family of proteins, multiple caspases, NOTCH1, SMAD7, CRMP1, and ADAM23 were differentially expressed. Interestingly, our data further indicate that this combination may be useful for the prevention of metastatic disease by reducing the EMT, as well as by increasing immune function, all of which have important clinical significance. Our future work will investigate this combination in a xenograph mouse model and determine the cellular and molecular mechanisms in vivo using RNA-seq, qPCR and Nanostring technologies in excised tumors to confirm these in vitro results.
List of abbreviations: ADAM: A Disintegrin And Metalloproteases, ATRA: all trans-retinoic acid, AXIN2: axis inhibition protein 2, Bcl-2: B-cell lymphoma 2, Bax: bcl-2-like protein 4, BIM: Bcl-2-like protein 11, CRC: colorectal cancer, CRMP1: Collapsin Response Mediator Protein 1, D2:ergocalciferol, D3: cholecalciferol, DGE: differential gene expression, EMT: epithelial mesenchymal transition, FC: Fold change, FDR: false discovery rate; GADPH: Glyceraldehyde 3-phosphate dehydrogenase, HDAC: histone deacetylase, IL-12: interleukin-12; IPA: Ingenuity Pathway Analysis, Jag1: Jagged-1, NOTCH: notch protein, Noxa: Phorbol-12-myristate-13-acetate-induced protein, PCA: principal component analysis, PUMA: p53 upregulated modulator of apoptosis, qPCR: quantitative polymerase chain reaction, SIRT: sirutin, SMAD7: SMAD Family Member 7
Competing interests: The authors declare that there are no conflicts of interest.
Authors’ contributions: GBM, TOL, SP, designed and conducted the research, NSL and ZA performed the RNA purification and RNA-seq; PNK and MMC performed the Bioinformatics and statistical analyses. GBM, TOL, SP grew, treated and harvested the HCT-116 cells, isolated RNA, analyzed data, performed IPA analyses, wrote and edited the manuscript.
For a full list of references, please download the pdf.