Deep sequencing and in silico analyses identify MYB-regulated gene networks and signaling pathways in pancreatic cancer

We have recently demonstrated that the transcription factor MYB can modulate several cancer-associated phenotypes in pancreatic cancer. In order to understand the molecular basis of these MYB-associated changes, we conducted deep-sequencing of transcriptome of MYB-overexpressing and -silenced pancreatic cancer cells, followed by in silico pathway analysis. We identified significant modulation of 774 genes upon MYB-silencing (p < 0.05) that were assigned to 25 gene networks by in silico analysis. Further analyses placed genes in our RNA sequencing-generated dataset to several canonical signalling pathways, such as cell-cycle control, DNA-damage and -repair responses, p53 and HIF1α. Importantly, we observed downregulation of the pancreatic adenocarcinoma signaling pathway in MYB-silenced pancreatic cancer cells exhibiting suppression of EGFR and NF-κB. Decreased expression of EGFR and RELA was validated by both qPCR and immunoblotting and they were both shown to be under direct transcriptional control of MYB. These observations were further confirmed in a converse approach wherein MYB was overexpressed ectopically in a MYB-null pancreatic cancer cell line. Our findings thus suggest that MYB potentially regulates growth and genomic stability of pancreatic cancer cells via targeting complex gene networks and signaling pathways. Further in-depth functional studies are warranted to fully understand MYB signaling in pancreatic cancer.

With the advent of next generation sequencing (NGS), high resolution genomic and transcriptomic information can be retrieved through the whole genome sequencing (WGS) and RNA sequencing (RNA-seq) 12 . The transcriptome profiling is rapidly replacing the hybridization-based microarrays as it provides an unbiased, extensive and precise measurement of levels of transcripts and their isoforms 13 . Differentially expressed genes in two or more conditions can be identified through RNA-seq and the biological significance of the transcriptomic alterations can be analysed through a number of bioinformatics tools. Ingenuity Pathway Analysis (IPA) provides one such user-friendly interface that can translate the changes in gene expression to that of altered networks and pathways 14 .
In this study, we analyzed the differential expression of genes in MYB-silenced MiaPaCa cells, relative to the MYB-expressing parental cells. The genes modulated upon MYB-silencing were annotated through comparative analysis and the biological significance of the altered transcriptome was interpreted via IPA. EGFR and RELA were observed to be down-regulated upon MYB-silencing and were confirmed to be direct transcriptional targets of MYB. Moreover, the MYB-induced changes in gene expression were also verified by ectopic expression of MYB in BxPC3 cell line, further strengthening the role of MYB in pancreatic cancer. Analyses of the dataset also suggested other novel functions of MYB in pancreatic cancer that warrant in-depth investigation to comprehend their functional relevance and are subject of ongoing research. The expression of a few select genes from the resulting MiaPaCa (shMYB Vs. NTScr) RNA-Seq data was also validated using qPCR. In a converse approach, we also compared the (C) relative expression of these genes in ectopic MYB-expressing (BxPC3-MYB) vs. control BxPC3 (BxPC3-Neo) cell lines. (D,E) The expression at protein level of these genes was also analyzed by immunoblotting. PUMA and GLUT1 represent the proteins encoded by BBC3 and SLC2A1 genes, respectively.
Scientific RepoRts | 6:28446 | DOI: 10.1038/srep28446 c-MYB, ADM, ALDH1, BBC3 (PUMA), HDAC5, KLF4, LDHA, MDM2, SHH, SLC2A1 (GLUT1); relevant to PC pathobiology were examined by quantitative reverse transcription PCR (qPCR) using gene specific primer set. The change in the expression profile of these genes was consistent with the RNA-seq data (Fig. 1B). Furthermore, to ascertain the role of MYB in modulating the expression of these genes in PC, in a converse approach wherein we overexpressed MYB in low MYB expressing BxPC3 cells, the expression profile observed was vice versa to MYB-knockdown in MiaPaCa PC cells (Fig. 1C). The expression of these genes at protein level was also examined in high (MiaPaCa-NT-Scr and BxPC3-MYB) and low (MiaPaCa-shMYB and BxPC3-Neo) MYB-expressing PC cells by immunoblotting. The changes at the protein level were in accordance with the changes observed at the transcript level for the analyzed genes (Fig. 1D,E).
Identification of MYB-regulated gene networks in pancreatic cancer. Next, we performed the Ingenuity Pathways Analysis (IPA) to interpret our dataset in the context of MYB-regulated gene networks in PC cells. IPA transforms the list of genes into a set of networks built on the records maintained in their database which are largely based on published reports 15,16 . Each network is generated algorithmically based on their connectivity and assigned a score. IPA analysis revealed 25 different gene networks upon MYB-silencing in PC cells (Supplementary Table 2). Amongst them, the top three modulated networks are presented in Fig. 2. Network 1 was predicted with a score of 54 and comprised of 34 focus molecules ( Fig. 2A). In this network, thirty genes were found to be down-regulated, including several splicing factors such as SRSF1, SRSF5, SRSF7 and SRSF10, while four genes (KIFC2, CHD3, CUL7 and KCNK3) were upregulated. As revealed by the IPA software, the top function associated with this network was "RNA Post-Transcriptional Modification, Molecular Transport and RNA Trafficking". Network 2 (Score -49, focus molecules -32; Fig. 2B), consisted of three upregulated (ARFGEF3, HES4, PLD3) and twenty-nine downregulated genes including EGFR, GINS complex, MCM complexes and primases, suggesting an overall downregulation of this network. The top function of the genes in this network is "Cellular Assembly and Organization". Network 3 (Score -41, focus molecules -41; Fig. 2C) consisted of eleven up-regulated and eighteen down-regulated genes. NF-κ B complex, RPA and KIAA1524 occupied the central nodes of this network and were all down-regulated. As determined by IPA software, these molecules are mainly involved in "Cell cycle, DNA Replication, recombination and Repair".
MYB-modulated signaling pathways in pancreatic cancer. In order to identify the signaling pathways catalogued in the IPA library that are of significance to our dataset, "canonical pathway" analysis of the IPA was employed. Ten enriched canonical signaling pathways upon MYB-silencing are represented in Fig. 3. An upregulation of pathways associated with cell cycle checkpoint regulation including "G2/M DNA Damage Regulatory Pathway" (z score: 0.302), "G1/S Cell Cycle Regulatory Pathway" (z score: 0.333), along with the upregulation of "p53 signaling" (z score: 0.632) was observed. These pathways are involved in induction of cell cycle arrest and apoptosis, and up-regulation of these pathways in MYB-silenced cells is supportive of the de-repression of cellular regulatory mechanisms in the absence of MYB. Surprisingly, "BRCA1-mediated DNA Damage Signaling" and "ATM signaling" (z score − 2.53 and − 1.897, respectively) were observed to be down-regulated. These signaling pathways are mostly reported to be inversely correlated with the expression of oncogenes, and further investigations are required to interpret the significance of these signaling pathways in PC, particularly in the context of MYB signaling. Further, as expected, the inhibition of "Pancreatic Adenocarcinoma Signaling" with a z-score of −0.333 was also revealed by the IPA analysis. Supplementary Table 3 lists the enriched pathways, genes involved, ratio of the genes present in our dataset to that present in IPA database, z-score and p-values. However, while changes in expression of genes involved in "Chromosomal Replication Pathway, Homologous Recombination Repair Pathway, Mismatch Repair Pathway" and "HIF-1α Signaling" was observed, IPA did not predict any activity pattern for these networks and their activation state could not be inferred. Together, this data suggests the power of IPA analysis for the quick identification of perturbed signaling pathways, and also validates the role of MYB in the regulation of several canonical pathways that are of significance to PC pathogenesis. Having observed down-regulation of pancreatic adenocarcinoma signaling pathway upon MYB-silencing, we next focused on the genes of this pathway that were present in our dataset. We observed that in our dataset, a number of factors representing pancreatic adenocarcinoma signaling pathways were down-regulated. These included EGFR, NF-κ B, VEGF, PI3K, MDM2 and CDK2 (Fig. 4).

Characterization of EGFR and RELA as direct transcriptional targets of MYB. EGFR and
NF-κ B-p65 have been reported to be of significance to PC pathogenesis and our RNA-seq data reported a 2.214 and 4.638 fold decrease, respectively, in transcript levels of these genes. Therefore, we subsequently explored the mechanism through which the inhibition of both EGFR and NF-κ B (RELA) was attained in MYB-silenced cells. We first measured the mRNA expression of EGFR and RELA in MiaPaCa-shMYB and BxPC3-MYB cells, relative to their respective controls. The qRT-PCR data showed a 7.14 and 4.02 fold decrease in transcript levels of EGFR and RELA, respectively, upon MYB silencing in MiaPaCa cells (Fig. 5A, top panel), which is consistent with their decreased expression as obtained in RNA-seq. Conversely, the mRNA expression of these genes was enhanced by forced expression of MYB in BxPC3 cells (Fig. 5A, bottom panel) with EGFR and RELA reporting a 6.54 and 3.74 fold increase in transcript levels, respectively. These results were also reflected at protein expression as MYB-silencing reduced expression of both EGFR and RELA in MiaPaCa-shMYB cells, while an opposite effect was seen in MYB overexpressing BXPC3 cells (Fig. 5B left panel). The densitometric analysis of the bands revealed a 9.65 and 1.65 fold decrease in protein expression of EGFR and RELA, respectively, upon Ten PC-related canonical pathways identified have been illustrated. Pathways identified are represented on the x-axis. The left y-axis corresponds to the -log of the P-value (Fisher's exact test) and the right y-axis represents the ratio (orange points) of the number of genes in a given pathway that meet cutoff criteria, divided by the total number of genes that map to that pathway.
Scientific RepoRts | 6:28446 | DOI: 10.1038/srep28446 MYB-silencing in MiaPaCa cell while its overexpression resulted in 6.90 and 1.49 fold increase in expression of these proteins (Fig. 5B, right panel). Moreover, alteration of MYB-expression in the PC cells significantly altered NF-κ B transcriptional activity and was observed to directly correlate with the expression of MYB (Fig. 5C).
To examine if EGFR and RELA are direct transcriptional targets of MYB, we performed in silico analysis of ~1 kilobase DNA sequence 5′ -upstream of the coding sequence (RefSeq ID NM_005228 and NM_021975, respectively) using online tools (ALGGEN-PROMO, TFBind). Two putative MYB binding sites were predicted in RELA promoter region (at − 386 and − 510 bases with respect to the first codon) and three potential binding sites for MYB were predicted in the promoter region of EGFR (at − 499, − 644 and − 713 with respect to the first codon) (Fig. 5D). The direct binding of MYB to EGFR and RELA promoters was confirmed by ChIP assay. Pull-down of the EGFR and RELA promoter sequence was significantly decreased in MYB-silenced MiaPaCa-shMYB cells while opposite results were seen in MYB-overexpressing BxPC3 cells (Fig. 5E). Taken together, our data demonstrate that MYB enhances the expression of EGFR and RELA in PC by directly binding to their promoter region.

Discussion
We recently identified MYB as an important regulator of pancreatic cancer pathogenesis which modulates tumor growth, aggressiveness as well as metastasis 6 . The analysis of the RNA-seq data of MYB-silenced MiaPaCa PC cells identified known MYB-target genes, e.g., ADM, BRCA1, CBFB, KLF4, LDHA, PCNA etc. 6,[17][18][19] , along with several novel direct/indirect targets. These modulated genes belong to various protein families serving as cytokines, cellular transporters, enzymes, cell cycle regulators, DNA damage and repair molecules, transcription regulators etc. (Supplementary Tables 1 and 2). Being a transcription factor, silencing of MYB was expected to down-regulate a number of genes [17][18][19] , however, a simultaneous up-regulation of genes was also observed indicating a potential direct/indirect transcription repressor activity of MYB.
IPA analysis is a powerful tool to study interactions between molecular factors. IPA-network analysis illustrates the biological relationships between genes in a dataset at the molecular level. In this study, the Network 1 ( Fig. 2A) suggested down-regulation of several proteins involved in alternative splicing of genes. The SRSF proteins depicted in this network have been demonstrated to have oncogenic potential and are frequently over expressed in many types of human tumors 20,21 . Similarly, MALAT1 (metastasis associated lung adenocarcinoma transcript 1); a long non-coding RNA frequently overexpressed in PC 22 was also decreased upon MYB-silencing. MALAT1 has been demonstrated to act as transcriptional regulator of genes involved in cell growth, migration and invasion 23 . Apart from splicing factors, down-regulation of DDX39A (DEAD box RNA helicase 39A), a RNA helicase, that plays a role in RNA splicing/export was also observed. DDX39A up-regulation in lung squamous cell cancer has been demonstrated to promote tumor growth 24 . Indeed this network holds significant importance as the spliceosome is now being widely acknowledged as a novel target of oncogenic stress in cancers 25,26 . Alternative splicing and the differential expression of splicing factors has been reported in pancreatic cancer as well as other cancers to produce aberrant variants that facilitate tumor associated characteristics [27][28][29] . However, the factors governing aberrant splicing in cancers are largely unknown. There is evidence to suggest role of MYB in regulation of alternative splicing in normal hematopoietic cells 30 , but its role in cancer remains to be explored. Overall, the down-regulation of several oncogenic factors that are represented in the network 1, upon MYB silencing, supports the oncogenic potential of MYB.
The second network identified genes around the EGFR, GINS complex, MCM complex and primase nodes. EGFR, also known as ErbB1, is encoded by the c-ERBB-1 a proto-oncogene and is a transmembrane tyrosine kinase growth factor receptor overexpressed in pancreatic cancer 31,32 . Down-regulation of GINS complex subunits (GINS1-4) is observed in our dataset (Fig. 2(B) and Supplementary Table 1). The PSF1/GINS1 is reported to be over-expressed in breast tumor cells as well as in highly proliferating cells where it enhances cell proliferation via increased DNA replication and anchorage independent growth of breast cancer cells 33,34 . The other node in this network includes mini chromosome maintenance complex (MCM4, MCM5, MCM10, and MCMBP). The MCM complex is over-expressed in human cancers and also in pre-cancerous cells undergoing malignant transformations. Its deregulation results in chromosomal defects and contributes to tumorigenesis 35 . Another axis involves the enzymes or proteins involved in DNA replication such as POLA1, encoding the catalytic subunit of DNA polymerase, which plays an essential role in DNA replication. Prim1 and Prim2 encoding smaller and larger subunit of DNA primase respectively form heterodimer that functions as a DNA-directed RNA polymerase to synthesize small RNA primers that are used to create okazaki fragments on the lagging strand of the DNA.
The third-top network identified upon silencing of MYB centered majorly around the down-regulation of NF-κ B complex node. Constitutive activation of NF-κ B transcriptional factors is observed in pancreatic cancers and regulates tumor development and progression. RIPK4 expression is also down-regulated in this group. It activates the NF-κ B signaling 36 and is also an important regulator of Wnt/beta-catenin signaling 37 , thus facilitating tumor growth and proliferation. Silencing of RIPK4 in cervical cancer cells inhibits cell migration and invasion through down-regulation of vimentin, MMP2 and fibronectin 37 . Interestingly, MIB2 (mindbomb homolog 2, also known as skeletrophin) is one of the genes upregulated in this network. This gene encodes a RING (Really Interesting New Gene) finger-dependent ubiquitin ligase, and acts as a negative regulator of invasion in many malignant tumors through down-regulation of the Met oncogene 38 . The second central node in this network was identified around KIAA1524/CIP2A gene and the encoded protein stabilizes the Met oncogene, promoting epithelial mesenchymal transition (EMT) through the increased expression of vimentin and snail proteins 39 . Moreover, CIP2A also associates with H-Ras activating MEK/ERK pathways to promote EMT 39 . The third central node in this network involves genes that are known to regulate DNA replication and recombination, i.e. TIPIN, Rfc complex and ATM/ATR complex.
In addition to the gene networks, genes modulated upon MYB-silencing were assigned to several canonical signaling pathways indicating that MYB regulates different pathways to exert its oncogenic potential in PC. Some of these enriched canonical signaling pathways corresponded to the DNA replication and repair pathways. For example, MYB can directly/indirectly modulate the expression of CDC6, CDC7 and MCM group of proteins that are required for the initiation of replication. Further, activation of cell cycle regulatory pathways responsible for either G1/S arrest or G2/M arrest was also observed. Interestingly, up-regulation of the p53 signaling pathway was identified upon silencing of MYB. Mutations in the p53 gene, considered to inactivate the transcription factor, have been observed in nearly all pancreatic cancers and are also present in the MiaPaCa PC cells. However, in the MYB-silenced dataset, an up-regulation of PUMA suggests an additional degree of p53-independent regulation of the tumor.
In MYB-silenced cells, the pancreatic adenocarcinoma signaling was also negatively regulated due to the inhibition of EGFR and NF-κ B. Furthermore, we identified both RELA and EGFR to be the direct transcriptional targets of MYB. It is possible that their reported overexpression in PC is through the MYB-mediated increase in mRNA levels. Overexpression of RELA has been reported in pancreatic adenocarcinomas and correlated with the activation of NF-κ B pathway 40 . Interestingly, while NF-κ B-P65 can be activated in response to various cytokines, growth factors, activating mutations; EGFR-dependent activation of NF-κ B signaling has also been reported in pancreatic cancer 41 . Moreover, NF-κ B activation in response to EGFR-targeted therapy has been identified as an adaptive resistance mechanism, thereby limiting EGFR therapy 42 .
In summary, we report for the first time various molecular mechanisms regulated by MYB in pancreatic cancer on a systemic level. We performed RNA-Seq analysis in MiaPaCa cell line model to generate a general model of MYB-dependent pathways in pancreatic cancer. We then used IPA analysis to predict the gene networks and canonical pathways affected by MYB signaling. The effects on many critical genes were not only validated in MiaPaCa model, but also verified in a reciprocal model by forced-expression of MYB in non-expressing cell line, further strengthening the role of MYB in pancreatic cancer pathogenesis. Among the many promising factors from several key signaling pathways, we identified a novel mechanistic involvement of EGFR and NF-κ B in the MYB signaling which needs to be further explored, particularly in the context of therapeutic targeting of MYB signaling in PC patients.
Total RNA isolation. Total RNA was extracted using TRIzol reagent (Invitrogen) from pancreatic cancer cell lines. The integrity and quality of the RNA was checked by denaturing RNA gels and quantitated on the Nanodrop 1000 (Thermo Scientific). cDNA synthesis and quantitative reverse transcription-PCR (qRT-PCR). Two μ g of total RNA isolated from pancreatic cancer cells was used for cDNA synthesis using the High Capacity Complementary DNA Reverse Transcription Kit (Thermo Scientific) following manufacturer's instructions. cDNA constructed was used as a template with SYBR Green Master Mix on an iCycler system (Bio-Rad, Hercules, CA) along with specific primer pair sets to perform qRT-PCR in 96-well plates. Supplementary Table 3 lists the sequences of the primers used, their respective Tm, and GC content. The thermal conditions employed for real-time PCR assays were as follows: cycle 1: 95 °C for 10 min, cycle 2 (× 40): 95 °C for 10 sec and annealing temperature respective to each primer pair for 45 sec.
RNA-sequencing and bioinformatics analysis. RNA-sequencing and bioinformatics analysis were performed at the Genomics Core Facility at the Heflin Center for Genetics, University of Alabama at Birmingham, as previously described 43 . In brief, polyA + mRNA isolated from MiaPaCa-NT-Scr and MiaPaCa-shMYB cells was converted to cDNA using random primers. Thereafter, cDNA sequencing libraries were generated using Scientific RepoRts | 6:28446 | DOI: 10.1038/srep28446 TruSeq library generation kit (Illumina, San Diego, CA) and quantified by qPCR using Roche LightCycler 480 with the Kapa Biosystems kit for library quantification (Kapa Biosystems, Wilmington, MA). Sequencing was performed on the Illumina HiSeq 2500 platform employing latest versions of sequencing reagents and flow cells generating up to 300 gigabytes of sequence information per flow cell. TopHat was used to align the raw RNA-Seq fastq reads to the human hg19 genome using the short read aligner Bowtie 44 . Cufflinks (version1.3.0) was used to align reads from TopHat and assemble transcripts, estimate their abundances and test for differential expression and regulation. The resulting RNA-Seq data was submitted to NCBI's Gene Expression Omnibus (GEO) database and is accessible through GEO Series accession number GSE61290 (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc= GSE61290).
Immunoblotting analysis. Immunoblotting was performed by the standard procedures as described earlier 6 . Cell lysates resolved on 10% polyacrylamide gels were transferred to PVDF membranes and subjected to immunodetection procedure using specific antibodies of interest and visualized using the SuperSignal West Femto Maximum sensitivity substrate kit (Thermo Scientific) on Bio-Rad image analyzer (Bio-Rad). ImageJ (imagej.nih.gov) software was used for densitometric analysis. The values were normalized to the expression of β -actin and the experiments were performed in triplicates.
Ingenuity pathway analysis. Networks and canonical pathways were generated through the Ingenuity Pathway Analysis (IPA; Ingenuity Systems, Qiagen, http://www.ingenuity.com/products/ipa) on the list of differentially expressed genes with the cutoff p-value ≤ 0.05 and fold change ± 1.5. Gene symbols were used as the identifier and the IPA based Ingenuity Knowledge Base was used as a reference to perform core analysis. Network algorithm was used to generate the network of genes based on their connectivity and were assigned a score based on the number of focus genes. Canonical pathway analyses identified pathways from the IPA library that were of significance to the dataset. NF-κB transcriptional activity assay. In order to study the transcriptional activity of NF-κ B, pancreatic cancer cell lines grown in 6-well plates were transfected with 1 μ g of NF-κ B-firefly luciferase based promoter reporter plasmid (pGL4.32 [luc2P/NF-κ B-RE/Hygro]) and 0.5 μ g of control reporter plasmid containing Renilla reniformis luciferase gene downstream of the TK promoter (pRL-TK). After 48 h, the transfected cells were harvested in passive lysis buffer and luciferase activity was measured using the Dual Luciferase Assay System (Promega, Madison, WI).
Chromatin immunoprecipitation (ChIP) assay. ChIP assay was performed using a ChIP-IT enzymatic kit as previously described 42 . Briefly, DNA-protein cross-linking was done with paraformaldehyde (37%) followed by enzymatic DNA shearing. Sheared DNA was then subjected to immunoprecipitation using anti-MYB or normal rabbit IgG (as control). Subsequently, cross-linking reversed, proteins digested with proteinase K and DNA isolated. PCR was performed using specific primers sets flanking MYB-binding promoter regions (Supplementary Table 3) and amplification products resolved on a 1.5% agarose gel and visualized using ethidium bromide staining. Input DNA (without immunoprecipitation) and normal IgG-precipitated DNA were used as positive and negative controls, respectively.