Next Article in Journal
Characterization of Translationally Controlled Tumour Protein from the Sea Anemone Anemonia viridis and Transcriptome Wide Identification of Cnidarian Homologues
Previous Article in Journal
Acknowledgement to Reviewers of Genes in 2017
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide Methylation Patterns in Androgen-Independent Prostate Cancer Cells: A Comprehensive Analysis Combining MeDIP-Bisulfite, RNA, and microRNA Sequencing Data

Department of Laboratory Medicine, The First Affiliated Hospital of Wenzhou Medical University, ShangCai Village, Ouhai District of Wenzhou, Wenzhou 325000, China
*
Authors to whom correspondence should be addressed.
Submission received: 30 November 2017 / Revised: 21 December 2017 / Accepted: 30 December 2017 / Published: 11 January 2018
(This article belongs to the Section Technologies and Resources for Genetics)

Abstract

:
This study aimed to investigate the mechanisms underlying the development of the androgen-independent phenotype in prostate cancer. Methylation patterns were detected in androgen-independent and androgen-dependent lymph node carcinoma of the prostate (LNCaP) prostate carcinoma cells based on methylated DNA immunoprecipitation-bisulfite sequencing data and differentially methylated regions (DMRs) were identified. Differentially expressed genes (DEGs) and micro RNAs (miRNAs) with DMRs (named MDEGs and MDEmiRNAs) were identified by combining transcriptome and methylation data, and transcription factor (TF)-DEGs with DMRs in promoter (PMDEGs) and MDEmiRNA-MDEGs networks were constructed. Furthermore, a time-course analysis of gene transcription during androgen deprivation was performed based on microarray data and DMRs, MDEGs, and DEmiRNAs were validated. In total, 18,447 DMRs, 3369 MDEGs, 850 PMDEGs, and 1 MDEmiRNA (miR-429) were identified. A TF-target network (94 PMDEGs and 5 TFs) and a miRNA–target network (172 MDEGs and miR-429) were constructed. Based on the time-course analysis of genes in the networks, NEDD4L and PBX3 were targeted by SOX5, while GNAQ, ANLN, and KIF11 were targeted by miR-429. The expression levels of these genes and miR-429 were confirmed by quantitative real-time polymerase chain reaction. Additionally, 109 DMRs were confirmed using additional public datasets. The regulatory pathways SOX5-NEDD4L/PBX3, miR429-GNAQ/ANLNRHOA, and miR429-ANLNKIF11 may participate in the progression of the androgen-independent phenotype in prostate cancer.

1. Introduction

Prostate cancer (PC) is one of the leading causes of cancer death in males aged more than 60 years. There are estimated to have been 180,890 new cases of PC and 26,120 deaths caused by PC in the United States in 2016, accounting for 21% of total new cancer cases and 8% of total cancer deaths in males [1]. Currently, androgen ablation therapy (AAT) is a front-line therapy for PC. However, this therapeutic method would fail if the PC developed an androgen-independent phenotype resulting in castration-resistant PC (CRPC), which is a primary cause of refractory and lethal PC [2].
Previous studies have elucidated several molecular mechanisms underlying the development of the androgen-independent PC phenotype. Androgen receptor (AR) has been considered a driver of therapeutic resistance in PC via its aberrant activation [3,4] or the presence of AR gene mutations (e.g., AR F876L) or variants (e.g., ARv567es and AR-V7) [4,5]. Furthermore, PC stem cells, which are independent of androgens for survival and have the abilities of paracrine and autocrine androgen biosynthesis and overexpression of antiapoptotic proteins (e.g., B-cell lymphoma 2 (Bcl-2) and glucocorticoid receptor), can also promote resistance to AAT [4,6,7,8]. These genes or proteins that are involved in resistance to AAT are likely to participate in the development of the androgen-independent phenotype. In one of our previous studies, we identified eight novel fusion genes, 315 alternative splicing events, and 788 differentially expressed genes (DEGs) in an androgen-independent PC cell line, lymph node carcinoma of the prostate (LNCaP)-AI-F [9]. Additionally, microRNAs (miRNAs) have also been reported to be associated with the development of the androgen-independent phenotype. Several miRNAs, such as miR-1205, miR-373, miR-21, and miR-29a, are abnormally expressed in androgen-independent PC cells, compared to androgen-dependent PC cells [10,11,12,13]. Moreover, DNA methylation has been found to occur during PC progression and frequently occurs in CRPC [14]. Based on the Panomics gene array system (Panomics, Redwood, CA, USA), one previous study identified that the methylation patterns of several genes (e.g., CASP8, CD14, and cyclin dependent kinase inhibitor 2A) were significantly different between androgen-sensitive/dependent (LNCaP) and androgen-independent (DU145 and PC3) cells [15]. However, the global distribution of aberrant methylation patterns in the genome of androgen-independent PC cells remains unclear. There has been no comprehensive analysis combining RNA-sequencing, miRNA-sequencing, methylation-sequencing, and microarray data to investigate the global methylation patterns in the genome of androgen-independent PC cells.
Methylated DNA immunoprecipitation sequencing (MeDIP-seq) is a quantitative high-throughput DNA sequencing approach to identify methylated sites with low DNA concentrations based on methyl-DNA immunoprecipitation with an antibody against 5-methylcytosine (5mC) [16]. However, MeDIP-seq is unable to identify individual 5mC sites in captured reads or distinguish un-methylated reads captured by 5mC antibodies due to its non-specific binding, which increases the rate of false positives when detecting methylatedCGs (mCGs) [17]. Whole-genome bisulfite sequencing, the most powerful and complete strategy for the quantitative genome-wide detection of 5mC at a single-base resolution [17], has been reported to reduce the rate of false positives in MeDIP-seq [18]. Thus, a novel method, MeDIP-bisulfite sequencing (MB-seq), is generated by combining the MeDIP-seq protocol with bisulfite treatment. This method can detect the methylation status of each CpG and has a lower rate of false positives than comparable methods, showing more robust sensitivity and specificity for detecting methylatedCpGs (mCpGs) than MeDIP-seq and cytosine methylome sequencing (MethylC-seq) [18]. In MB-seq, the immunoprecipitation step enables the enrichment of methylated fragments of DNA, while the bisulfite treatment enables the distinction of methylated cytosines from un-methylated cytosines in the enriched fragments. Furthermore, in MB-seq, all the cytosines in the adapters are hydroxymethylated [18]. MB-seq has not yet been used to study the molecular mechanisms underlying the progression of the androgen-independent phenotype in PC cells.
In the present study, MB-seq was applied to investigate the global distribution of DNA methylation across the whole genomes of androgen-independent and androgen-dependent PC cells. Furthermore, RNA-sequencing, miRNA-sequencing, methylation-sequencing, and microarray data were comprehensively analyzed in combination with the MB-seq data. These results provide novel insights into the molecular mechanisms underlying the progression of the androgen-independent phenotype in PC cells.

2. Materials and Methods

2.1. Cell Lines

An androgen-dependent PC cell line, LNCaP, was purchased from the American Type Culture Collection (Rockville, MD, USA; passage 1). This was cultured in a routine medium and split once per week by trypsinization. The routine medium consisted of phenol red-positive F-12 medium (Gibco®, Life Technologies, Carlsbad, CA, USA) and 10% (v/v) fetal bovine serum (FBS) (Biowest, Nuaillé, France). LNCaP was chronically cultured in an androgen-deprived medium to establish an androgen-independent LNCaP line. First, LNCaP cells were cultured in an androgen-deprived medium composed of the phenol red-free dulbecco's modified eagle medium (DMEM)/F-12 medium (Gibco®) plus 10% (v/v) charcoal/dextran-treated FBS (Biowest) for 5 passages. As flutamide was able to competitively inhibit the binding of androgens to AR, LNCaP cells were then cultured in the above androgen-deprived medium plus 10−7 mol/L flutamide (Schering-Plough, Kenilworth, NJ, USA) [13] for 105 passages to obtain a complete androgen-independent cell line. The resulting cell line was named lymph node carcinoma of the cancer prostate—artificially induced—flutamide (LNCaP-AI-F) (passage 110), and its androgen-independence was confirmed by detecting cell morphological changes, cell proliferation rate, cell cycle, and the expression of BCL-2 and BAX in a previous study [13]. Meanwhile, the androgen-dependent LNCaP cells used as the control were cultured in the routine medium for 110 passages. All cells were cultured at 37 °C in an incubator with 5% CO2.

2.2. MeDIP-bisulfite Sequencing

Genomic DNA was extracted from LNCaP and LNCaP-AI-F cells using the DNeasy® blood and tissue kit (Qiagen, Hilden, Germany). After the verification of DNA integrity and DNA quantification, 1 μg DNA was sonicated into fragments of 100–500 base pairs (bp). Afterwards, the end-repair of DNA fragments, addition of adenine, and ligation to C-hydroxymethylated Illumina multiplexing adapters (adapter 1: 5′-pho-GATCGGAAGAGCACACGTCT-3′; adapter 2: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-3′) [19] were performed using the TruSeq DNA sample prep kit (Illumina, San Diego, CA, USA). The DNA was then subjected to a standard MeDIP assay [20,21], and MeDIP-enriched DNA was recovered from the resultant Dynabead®-MeDIP-DNA mixture using phenol-chloroform and ethanol. Subsequently, MeDIP-DNA was bisulfite-treated twice using the EpiTect® bisulfite kit (Qiagen) and further amplified for 12 cycles with Illumina multiplexing polymerase chain reaction (PCR) primers 1.0 and 2.0 (final concentration: 0.5 mM; Kapa Biosystems, Wilmington, MA, USA) in an 8 × 25 µL reaction system. After purification of the PCR products and gel electrophoresis on a 2% (w/v) low melting agarose gel, libraries with insert sizes of 270–370 bp were selected and quantified using the Quant-iT™ PicoGreen® double stranded DNA (dsDNA) reagent and kits (Invitrogen, Carlsbad, CA, USA). Afterwards, the DNA was sequenced for 200 cycles on a next-generation sequencing instrument (Illumina Hiseq 2000 with a TruSeq SBS Kit v3-HS) [18]. The raw MB-seq data were deposited in the Sequence Read Archive (SRA) database (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/sra) of the National Center for Biotechnology Information (NCBI) with accession number of SRP067471.

2.3. Raw Data Preprocessing and Reads Mapping

The raw data were preprocessed using GA pipeline software version 1.3 (Illumina) with the default settings. Sequencing reads were filtered using FASTX-Toolkit (http://hannonlab.cshl.edu/). Bases with a quality assessment score < 10 at the 5′-terminal or continuous “N” at the terminals and reads in which more than 20% of bases had a quality assessment score < 20 were removed. Only clean reads were subjected to the subsequent analyses.
The Bismark software [22], a flexible aligner and methylation caller for bisulfite-seq data, was applied to map the clean reads to the human genome sequence (hg19) that was downloaded from the Genome Bioinformatics Site at the University of California Santa Cruz (UCSC, http://genome.ucsc.edu). The methylation of C in CpG, CHG, and CHH sequences was detected using the Bismark software [22]. The software parameters were set as seed length = 30 and the number of mismatches within seed ≤ 3.

2.4. Differentially Methylated Regions Screening

Based on the clean reads, differential methylation patterns were detected between LNCaP-AI-F and LNCaP samples within a window of 50 bp, and the adjacent windows were merged (length of merged window ≤ 2000 bp). A t-test was used to assess the significance of differential methylation, and the raw p-value was adjusted into a false discovery rate (FDR) using the Benjamini–Hochberg method [23]. Only the regions with FDR < 0.05 and fold change > 1.25 were identified as differentially methylated regions (DMRs).

2.5. Annotation and Distribution of Differentially Methylated Regions

The annotation information about the transcription start site, transcription termination site, exon, intron, and repeats were downloaded from the UCSC database. Information about untranslated regions (UTRs) was obtained from UTRdb (http://utrdb.ba.itb.cnr.it/), and information about enhancers was generated from the VISTA Enhancer Browser (http://enhancer.lbl.gov/). Based on the above information, transcription start sites, transcription termination sites, exons, introns, repeats, UTRs, and enhancers that overlapped with the DMRs were identified.

2.6. Analysis of Transcription Factor-Binding Motifs in Differentially Methylated Regions -Promoter-Overlapping Regions

Methylation within CpG islands can impede transcription factor (TF)-DNA binding, affecting the transcription of corresponding genes. The DMRs that overlapped with promoters were obtained using the Perl program (https://www.perl.org/), and motif enrichment analysis was performed using the online tool MEME Suite [24]. Afterwards, the enriched motifs were aligned with the TF-binding motifs of vertebrates in JASPAR CORE 2014 (http://jaspar.genereg.net) [25].

2.7. Comprehensive Analysis of Differentially Methylated Regions and Transcriptome

The DEGs and differentially expressed miRNAs (DEmiRNAs) between LNCaP-AI-F and LNCaP cells were reported in our previous studies [9,13]. The chromosomal location information of the DEGs, DEmiRNAs, and DMRs, which was obtained from the Ensembl genome database project (http://asia.ensembl.org/index.html), was input into Linux shell programming [26], and the DEGs and DEmiRNAs that overlapped with DMRs were identified and defined as MDEGs and MDEmiRNAs, respectively. As methylation of the TF-binding region in a gene could inhibit the binding of TF to that region, we specifically defined the DEGs with DMRs in their promoters as PMDEGs.

2.8. Pathway and Functional Enrichment Analyses of Differentially Expressed Genes with DMRs (MDEGs)

To study the potential functions of MDEGs, the online tool Targetmine (http://targetmine.nibio.go.jp/) was used to perform pathway and functional enrichment analyses based on the KEGG [27], Reactome [28], WikiPathways [29], NCI-PID [30], and Gene Ontology (GO) [31] databases. The p-value was adjusted by the Holm–Bonferroni method, and the adjusted p-value < 0.05 was set as the cut-off criterion.

2.9. Construction of Transcription Factor-Target Network

Promoters of PMDEGs were aligned with TF-binding motifs in the DMR-promoter-overlapping regions identified above using the MotifDb R package [32] (criterion: matching value in pulsewidth modulation algorithm > 85%). Then the corresponding TF-target pairs and Protein-Protein interactions (PPIs) of targets obtained from the STRING database (http://string-db.org/) [33] were used to construct a TF-target network, which was visualized using Cytoscape software [34].

2.10. Construction of micro RNA–Target Network

Based on the information of miRNA–target gene pairs in the online database miRDB [35], the target genes of MDEmiRNAs were predicted (criterion: target score > 50), and the MDEGs among these target genes were selected to construct a miRNA–target network, which was visualized by Cytoscape [34].

2.11. Time-Course Analysis of Gene Transcription during Androgen Deprivation

As the regulation of gene expression is a dynamic process, it is important to identify and characterize changes in gene expression over time. Arrays collected over a time course allow for the study of the dynamic behavior of gene expression and clustering analysis has been commonly applied to time-course microarray data [36]. GSE8702 [37], a gene transcription profile during the androgen deprivation of LNCaP, was downloaded from the GEO database (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/). The raw data for this dataset were generated via the Affymetrix Human Genome U133 Plus 2.0 Array. GSE8702 contained six LNCaP samples cultured in normal medium (control group) and ten LNCaP samples cultured in androgen-deprived medium (deprivation group) for 3 weeks, 1 month, 5 months, 11 months, and 12 months. In our study, samples in the deprivation group were utilized to study the time-course changes in gene transcription during androgen deprivation. After data pre-processing using the Affy package [38], the fuzzy c-means clustering method [39] in Mfuzz software [39,40] was utilized to cluster genes based on their messenger RNA (mRNA) levels at the five time-points. The parameters in the clustering analysis were set as minimum standard deviation = 0.5 and score = 0.6, with a higher score representing a closer clustering. Genes in the Mfuzz-clusters were further clustered, and their transcription levels were shown using pheatmap in R [41].

2.12. Comprehensive Differentially Expressed Genes Analysis of Sequencing Data and Microarray Data

To identify the genes that were consistently up-regulated or down-regulated during androgen deprivation, the genes in the above Mfuzz-clusters were compared with the DEGs identified in a previous study [9], PMDEGs in the TF-target network, or MDEGs in the miRNA–target network. Venn diagrams were drawn using the online software Venny 2.0.2 [42].

2.13. Data Validation of the Identified Differentially Methylated Regions and Differentially Expressed micro RNAs

The methylation profiling dataset GSE41701 [43] deposited in GEO was used for the data validation of the above-identified DMRs. This dataset contained seven benign prostate tissue samples, seven primary localized prostate cancer tissue samples, and 6 metastatic CRPC tissue samples. Only the primary localized prostate cancer tissue samples and metastatic CRPC tissue samples were used for analysis. Firstly, DMRs between the primary and CRPC samples in GSE41701 were identified by a previously reported method [43]. Briefly, bisulfite reads were aligned to the bisulfite-converted hg19 reference genome using Bismark software. All samples had bisulfite conversion rates > 99.5%. Differentially methylated CpGs were identified using the Fisher exact test with correction for multiple testing by the Benjamini–Hochberg method. DMRs were defined as regions containing at least five differentially methylated CpGs (FDR = 0.05) with a fold change > 1.25. Afterwards, the DMRs in GSE41701 that overlapped with those in the present study were identified.
Furthermore, the non-coding RNA profiling GSE55829 was utilized to validate the expression of the above identified DEmiRNAs, which included four xenograft prostate cancer tissue samples at the androgen-dependent stage and four at the CRPC stage. DEmiRNAs between the androgen-dependent stage and the CRPC stage were identified using a previously reported method [13]. Subsequently, the DEmiRNAs in GSE55829 that overlapped with those in the previous study [13] were identified.

2.14. Quantitative Reverse Transcription Polymerase Chain Reaction Assay

Total RNA for mRNA expression analysis was extracted and purified from LNCaP-AI-F and LNCaP cells using TRIzol® reagent (Invitrogen) according to the manufacturer’s protocol. Afterwards, total RNA was reverse transcribed using Moloney murine leukemia virus reverse transcriptase (Promega, Madison, WI, USA). After complementary DNA (cDNA) synthesis, mRNA expression levels were measured using SYBR® Green quantitative polymerase chain reaction (qPCR) Supermix (Bio-Rad, Hercules, CA, USA). Furthermore, total RNA for miRNA expression detection was isolated using the mirVana™ miRNA isolation kit (Ambion, Austin, TX, USA) and purified by the miRNeasy® mini kit (Qiagen). miRNA expression levels were determined using the mirVana™ qPCR detection kit (Ambion). Primers were designed using Primer5 software (PRIMER-E, Plymouth, UK). Each PCR was repeated at least three times.
Expression levels of both mRNA and miRNA were analyzed by a CFX96™ real-time detection system (Bio-Rad). Relative gene expression was calculated using the 2−ΔΔCt method [44]. The expression levels of mRNAs were normalized against glyceraldehyde-3-phosphate dehydrogenase and the expression level of the miRNA was normalized against U6 small nuclear RNA.

2.15. Statistical Method

The data were presented as the mean ± standard deviation. Statistical analysis was performed using SPSS 19.0 (IBM, Armonk, NY, USA). Differences in the mRNA expression levels of genes and miRNA between LNCaP-AI-F and LNCaP cells were analyzed by two-tailed Student’s t-tests and were considered significant if p < 0.05.

3. Results

3.1. Overview of the MeDIP-bisulfite Sequencing Data

The MB-seq analysis of the LNCaP-AI-F and LNCaP cell lines generated approximately 48.0 and 50.6 million reads, respectively. After data filtering, approximately 80% clean reads were uniquely mapped to the human genome. The ratio of total methylated cytosine in LNCaP-AI-F cells was 8.45%, which was a little higher than that in LNCaP cells (8.20%). In addition, the methylation levels in CpG were much higher than those in CHG and CHH in both LNCaP-AI-F and LNCaP cells.

3.2. Annotation and Distribution of Differentially Methylated Regions

A total of 18,447 DMRs were identified between LNCaP-AI-F and LNCaP cells, including 10,135 hypermethylated and 8312 hypomethylated DMRs in LNCaP-AI-F cells. DMRs were present in 55.6% of genes and 24.54% of intergenic regions (Figure 1A). Additionally, 24.11% of the introns and 23.32% of the mRNAs contained DMRs (Figure 1B).

3.3. Transcription Factor-Binding Motifs in Differentially Methylated Regions -Promoter-Overlapping Regions

A total of 2816 DMRs were identified to overlap with the promoters of 2587 genes, and these genes were mainly associated with the functions of cellular processes (GO:0009987) and binding (GO:0005488) (Figure 1C). Additionally, these DMR-promoter-overlapping regions were enriched in ten TF-binding motifs, which were similar to the binding motifs of 20 TFs (Table 1).

3.4. Comprehensive Analysis of Differentially Methylated Regions and Transcriptome Data

Based on the obtained DMRs in LNCaP-AI-F cells, and DEGs (4653 up-regulated DEGs and 4074 down-regulated) identified in our previous study [9], a total of 3369 MDEGs were identified, including 1614 up-regulated MDEGs (involving 3840 DMRs) and 1775 down-regulated MDEGs (involving 4925 DMRs). Among the 3840 DMRs, 1989 were hypermethylated and 1850 were hypomethylated in LNCaP-AI-F cells compared with LNCaP cells. Meanwhile, among the 4925 DMRs, 2883 were hypermethylated and 2042 were hypomethylated. Furthermore, among the 2587 genes whose promoters overlapped with DMRs, 381 genes involving 385 DMRs were up-regulated in LNCaP-AI-F cells, and they were defined as up-regulated PMDEGs. In contrast, 469 genes involving 519 DMRs were down-regulated and were defined as down-regulated PMDEGs. Thus, a total of 850 PMDEGs were identified.
Meanwhile, based on the 86 up-regulated DEmiRNAs and 91 down-regulated DEmiRNAs found in our previous study [13], only one MDEmiRNA was identified (hsa-miR-429). Specifically, hsa-miR-429 was down-regulated in LNCaP-AI-F cells, and it was located within a hypomethylated DMR.

3.5. Pathway and Functional Enrichment Analyses of Differentially Expressed Genes with DMRs

The up-regulated MDEGs were significantly enriched in eight pathways and 119 GO terms. The pathways and the top five terms were mainly related to cell cycle, mitotic, nuclear part, and nucleic acid binding (Table 2). Meanwhile, down-regulated MDEGs were significantly enriched in three pathways and 35 GO terms. These pathways and the top five terms were mainly associated with phagocytosis, transport, cytoplasm, transferase activity, and nucleic acid binding (Table 2).

3.6. Analysis of the Transcription Factor-Target and micro RNA–Target Networks

To uncover the potential regulatory relationships among the identified PMDEGs, TFs were predicted from the PMDEGs. The promoters of the 850 PMDEGs were aligned with the ten TF-binding motifs that were identified above (Table 1). Five TFs targeting 94 PMDEGs and 57 protein-protein interactions of PMDEGs were identified, and they were used to construct a TF-target network (Figure 2A). Notably, in the TF-target network, SOX5 targeted 26 up-regulated PMDEGs and 36 down-regulated PMDEGs.
Furthermore, 768 potential target genes of hsa-miR-429 were predicted and 172 of them were MDEGs, including 98 up-regulated ones and 74 down-regulated ones. In the miRNA–target network of 172 MDEGs and hsa-miR-429 (Figure 2B), RHOA interacted with 15 MDEGs. ZEB1, shown in Table 1, was predicted to be a TF targeting hsa-miR-429, and is known to inhibit the translation of hsa-miR-429 [45].

3.7. Time-Course Analysis of Gene Transcription during Androgen Deprivation

Based on the mRNA levels at five time-points, a total of six gene-clusters involving 1367 genes were identified due to their differential expression patterns (Figure 3A). A total of 171 genes in cluster 1, 112 genes in cluster 2, and 313 genes in cluster 6 were up-regulated, whereas 348 genes in cluster 3 and 245 genes in cluster 5 were down-regulated in LNCaP cells during androgen deprivation. Additionally, 178 genes in cluster 4 were first up-regulated and then down-regulated during androgen deprivation. According to the hierarchical clustering analysis, genes in Mfuzz-clusters 1, 2, 3, 4, 5, and 6 were clearly distinguished via pheatmap-clusters a, b, c, d, e, and f (Figure 3B). These results suggested the accuracy of both clustering analyses and the identified genes were defined as time-course genes.

3.8. Comprehensive Differentially Expressed Genes Analysis of Sequencing Data and Microarray Data

Genes in clusters 1, 2, and 6 were first compared with the 4653 up-regulated DEGs found in the previous study [9], and respectively 91, 35, and 217 common genes were obtained (Figure 4A). Similarly, genes in clusters 3 and 5 were compared with the 4074 down-regulated DEGs found in the previous study [9], and 74 and 101 overlapped genes were identified, respectively (Figure 4B). These genes were consistently dysregulated in LNCaP cells during androgen deprivation and differentially expressed between LNCaP-AI-F and LNCaP cells.
Genes in clusters 1, 2, and 6 were compared with the up-regulated PMDEGs in the TF-target network (Figure 2A) and only one common gene DDX41 was found between cluster 1 and up-regulated PMDEGs (Figure 4C). Genes in clusters 3 and 5 were compared with the down-regulated PMDEGs in the TF-target network (Figure 2A) and only NEDD4L and PBX3 in cluster 5 were among those PMDEGs (Figure 4D). Moreover, both NEDD4L and PBX3 were targeted by SOX5, namely, SOX5-NEDD4L/PBX3.
Genes in clusters 1, 2, and 6 were compared with the up-regulated MDEGs in the miRNA–target network (Figure 2B). PKIA in cluster 2, as well as KIF11, WHSC1, ANLN, LBR, and NCAPG2 in cluster 6 were identified as MDEGs in that network (Figure 4E). Genes in clusters 3 and 5 were compared with the down-regulated MDEGs in the miRNA–target network (Figure 2B), and RELN, GNAQ, and KBTBD11 in cluster 3, as well as LAMC1, MBOAT2, ZNF532, and DPY19L3 in cluster 5 were among those MDEGs (Figure 4F). Additionally, both GNAQ and ANLN interacted with RHOA (namely, miR429-GNAQ/ANLNRHOA) and ANLN also interacted with KIF11 (namely, miR429-ANLN—KIF11).

3.9. Validated DMRs and Differentially Expressed micro RNAs in Public Datasets

A total of 109 DMRs identified in this study overlapped with those in the dataset GSE41701. The TF targets containing these confirmed DMRs (e.g., FBXO44, VPS72, and RIMS1) are also shown in Figure 2A.
Additionally, a set of 33 up-regulated differentially expressed (DE)miRNAs found in the previous study [13] overlapped with those in GSE55829. None of the down-regulated DEmiRNAs overlapped between those datasets. Strikingly, hsa-miR-429, which was down-regulated in the previous study [13], was up-regulated in GSE55829

3.10. Validation of Differentially Expressed Genes with DMRs and Differentially Expressed micro RNAs by Quantitative Real-Time Polymerase Chain Reaction

The mRNA expression levels of NEDD4L, PBX3, and GNAQ, as well as hsa-miR-429 were significantly lower in LNCaP-AI-F cells than in LNCaP cells (p < 0.05). Furthermore, the mRNA expression levels of ANLN, RHOA, and KIF11 were significantly higher in LNCaP-AI-F cells than in LNCaP cells (p < 0.05, Figure 5). The results of qRT-PCR were consistent with those obtained via bioinformatics analysis.

4. Discussion

In this study, we utilized a novel method, MB-Seq, to detect DNA methylation, and then combined it with RNA-sequencing, miRNA-sequencing, and microarray data to investigate the mechanisms of androgen-independence in PC cells. We identified three potential regulatory pathways, namely SOX5-NEDD4L/PBX3, miR429-GNAQ/ANLNRHOA, and miR429-ANLNKIF11, which involve the time-course down-regulated PMDEGs of NEDD4L and PBX3, the time-course down-regulated MDEG of GNAQ, the time-course up-regulated MDEGs of ANLN, KIF11, and RHOA, the TF of SOX5, and the down-regulated miRNA of hsa-miR-429.
For the SOX5-NEDD4L/PBX3 pathway, NEDD4L and PBX3 were consistently down-regulated in LNCaP cells during androgen deprivation and in androgen-independent LNCaP-AI-F cells compared with those in androgen-dependent LNCaP cells. In addition, they were predicted to be targeted by SOX5. NEDD4L encodes a ubiquitin ligase, which has been reported to be up-regulated after androgen treatment in LNCaP cells [46]. Decreased expression of NEDD4L has been reported to correlate with poor prognosis in gastric cancer patient [47]. Moreover, NEDD4L is suggested to be down-regulated in colorectal and lung cancers [48,49]. In the present study, the down-regulation of NEDD4L during androgen deprivation and in LNCaP-AI-F cells was consistent with the report by Qi et al. [46]. SOX5 and DNA methylation may be additional regulators of NEDD4L expression during the progression of the androgen-independent phenotype. PBX3 can regulate genes involved in steroidogenesis and the differentiation of urogenital organs, and steroidogenesis is related to the development of castration-resistant PC. Ramberg et al. reported that PBX3 is up-regulated in PC tissue and post-transcriptionally regulated by androgen in PC cells in an AR-independent manner [50]. Additionally, PBX3 is overexpressed in other cancers, such as gastric and colorectal cancers [51,52]. SOX5 encodes sex-determining region Y-box 5, and its overexpression is associated with the progression and distant metastasis of PC, nasopharyngeal carcinoma, and hepatocellular carcinoma [53,54,55]. SOX9, a homologue of SOX5, regulates the expression of AR in PC cells and enhances PC invasion [56,57]. In the present study, the expression levels of NEDD4L, PBX3, and SOX5 were confirmed by qRT-PCR. These findings suggested that NEDD4L, PBX3, and SOX5 may play crucial roles in the development of the androgen-independent phenotype in PC cells. This pathway may be used as a potential therapeutic target in androgen-independent PC.
For the miR429-GNAQ/ANLNRHOA pathway, GNAQ was down-regulated, whereas ANLN and RHOA were up-regulated in androgen-independent LNCaP-AI-F cells, compared with those in androgen-dependent LNCaP cells, which was confirmed by qRT-PCR. GNAQ mutations have been found in thyroid cancer [58] and GNAQ/GNA11 mutations are able to initiate human uveal melanoma [59]. Although there are few reports on the function of GNAQ in androgen-independent PC, our results suggest the involvement of GNAQ in the transition from androgen-dependent LNCaP to androgen-independent LNCaP-AI-F. ANLN encodes anillin, an actin binding protein, which has been overexpressed in non–small cell lung and breast cancers, serving as a potential target candidate [60,61]. Importantly, a previous microarray analysis found that ANLN was up-regulated in hormone-refractory PC progression [62], which was in agreement with our results. Interestingly, in the present study, both GNAQ and ANLN were predicted to interact with RHOA. RHOA encodes a RhoA-GTPase, the activation of which is associated with metastasis in various malignancies. For instance, Kamai et al. [63] found that overexpression of RHOA is associated with progression of testicular cancer. The inactivation of RhoA-GTPase mediates the inhibitory effects of FTY720 on the invasion of androgen-independent PC cells [64], indicating that the activation of RHOA may be related to the invasion of androgen-independent PC cells. This result supports our finding that the up-regulation of RHOA in androgen-independent LNCaP-AI-F cells coincided with the metastatic property of androgen-independent PC. miR-429 belongs to the miR-200 family, and the miR-200 family of miRNAs are found to be down-regulated in PC cells [65], which is consistent with the results of our study. Interestingly, miR-429 was up-regulated in ovarian and colorectal cancers [66,67], suggesting the different roles of miR-429 in different cancer types. Collectively, GNAQ, ANLN, and RHOA might participate in the androgen-independent progression of PC cells and their expression might be post-transcriptionally regulated by miR-429. They may serve as biomarkers in androgen-independent PC.
For the miR429-ANLNKIF11 pathway, ANLN and KIF11 were consistently up-regulated in LNCaP cells during androgen deprivation, as well as in androgen-independent LNCaP-AI-F cells. The up-regulation of the two genes in LNCaP-AI-F cells was validated by qRT-PCR. KIF11, also known as KSP or Eg5, encodes the human kinesin Eg5. A high expression level of KIF11 has been reported in many types of human cancers, such as oral and gastric cancers [68,69]. KIF11 is required for bipolar spindle formation, which is essential for cell division during mitosis. A previous study found that androgen-independent PC3 cells expressed more Eg5 protein compared with LNCaP cells [70]. These findings suggested that KIF11 is involved in the development of androgen independence in PC cells and its expression might be regulated by miR-429. Thus, we speculated that the miR429-ANLNKIF11 pathway may serve as novel prognostic biomarker and therapeutic target for androgen-independent PC.
In the present study, miR-429 was the only identified MDEmiRNA. Although a previous study reported the abnormal expression of miR-429 in LNCaP-AI cells compared to LNCaP cells [13], the methylation status of miR-429 in androgen-independent PC cells has not previously been studied to the best of our knowledge. Strikingly, miR-429 was confirmed by qRT-PCR to be down-regulated in LNCaP-AI-F cells compared with LNCaP cells. However, it was found to be up-regulated in xenograft prostate cancer tissue samples at the CRPC stage compared with the androgen-dependent stage, based on the data validation using GSE55829. We speculate that the differences in sample types (cells were used in the previous study [13] whereas tissues were used in the dataset GSE55829) may contribute to the conflicting results. Further experimental validation in CRPC tissues will be conducted to resolve the conflicting observation.
Despite the aforementioned results, there are several limitations to this study. The regulatory relationships of regulators (e.g., miR-429 and SOX5) and their targets need to be confirmed by experiments. Furthermore, the driving effect of DNA methylation on gene expression and the associations between DNA methylation and the development of androgen-independent PC need to be further experimentally investigated. These questions will be addressed in our future studies.

5. Conclusions

In conclusion, the regulatory pathways, SOX5-NEDD4L/PBX3, miR429-GNAQ/ANLNRHOA, and miR429-ANLNKIF11 may participate in the transition from androgen-dependent LNCaP cells to androgen-independent LNCaP-AI-F cells. Our results may provide novel information for the study of the molecular mechanisms underlying the progression of the androgen-independent phenotype. Those genes and miR-429 may serve as therapeutic targets for delaying the development of androgen-independent PC.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Project number: 81401737), the Research Project of Zhejiang Provincial Education Department (Project number: Y201430425), and the Wenzhou Municipal Science and Technology Bureau of China (Project number: Y20140024).

Author Contributions

Y.W., T.Q., W.H., and B.C. contributed to the sample acquisition, analysis and interpretation of data and drafting of the manuscript. M.D. and G.X. contributed to the conception and design of the study.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Siegel, R.L.; Miller, K.D.; Jemal, A. Cancer statistics, 2016. CA Cancer J. Clin. 2016, 66, 7–30. [Google Scholar] [CrossRef] [PubMed]
  2. De La Taille, A.; Vacherot, F.; Salomon, L.; Druel, C.; De Medina, S.G.D.; Abbou, C.; Buttyan, R.; Chopin, D. Hormone-refractory prostate cancer: A multi-step and multi-event process. Prostate Cancer Prostatic Dis. 2001, 4, 204–212. [Google Scholar] [CrossRef] [PubMed]
  3. Kahn, B.; Collazo, J.; Kyprianou, N. Androgen receptor as a driver of therapeutic resistance in advanced prostate cancer. Int. J. Biol. Sci. 2014, 10, 588–595. [Google Scholar] [CrossRef] [PubMed]
  4. Karantanos, T.; Evans, C.P.; Tombal, B.; Thompson, T.C.; Montironi, R.; Isaacs, W.B. Understanding the mechanisms of androgen deprivation resistance in prostate cancer at the molecular level. Eur. Urol. 2015, 67, 470–479. [Google Scholar] [CrossRef] [PubMed]
  5. Nelson, W.G.; Yegnasubramanian, S. Resistance emerges to second-generation antiandrogens in prostate cancer. Cancer Discov. 2013, 3, 971–974. [Google Scholar] [CrossRef] [PubMed]
  6. Amaral, T.M.S.; Macedo, D.; Fernandes, I.; Costa, L. Castration-resistant prostate cancer: Mechanisms, targets, and treatment. Prostate Cancer 2012, 2012, 327253. [Google Scholar] [CrossRef] [PubMed]
  7. Sonpavde, G.; Matveev, V.; Burke, J.; Caton, J.R.; Fleming, M.T.; Hutson, T.E.; Galsky, M.D.; Berry, W.R.; Karlov, P.; Holmlund, J.T.; et al. Randomized phase II trial of docetaxel plus prednisone in combination with placebo or AT-101, an oral small molecule Bcl-2 family antagonist, as first-line therapy for metastatic castration-resistant prostate cancer. Ann. Oncol. 2012, 23, 1803–1808. [Google Scholar] [CrossRef] [PubMed]
  8. Arora, V.K.; Schenkein, E.; Murali, R.; Subudhi, S.K.; Wongvipat, J.; Balbas, M.D.; Shah, N.; Cai, L.; Efstathiou, E.; Logothetis, C.; et al. Glucocorticoid receptor confers resistance to antiandrogens by bypassing androgen receptor blockade. Cell 2013, 155, 1309–1322. [Google Scholar] [CrossRef] [PubMed]
  9. Wang, Y.; Liu, Q.; Xu, G.; Mao, F.; Qin, T.; Teng, H.; Cai, W.; Yu, P.; Cai, T.; Zhao, M.; et al. Comparative RNA-Seq analysis reveals potential mechanisms mediating the conversion to androgen independence in an LNCaP progression cell model. Cancer Lett. 2014, 342, 130–138. [Google Scholar] [CrossRef] [PubMed]
  10. Durojaiye, V.; Ilboudo, A.; Levine, F.; Osborne, J.; Park, J.Y.; Ogunwobi, O.O. miR-1205/FRYL as a novel regulatory mechanism in androgen-independent prostate cancer. Cancer Res. 2015, 75 (Suppl. 15). [Google Scholar] [CrossRef]
  11. Yang, K.; Handorean, A.M.; Iczkowski, K.A. MicroRNAs 373 and 520c are downregulated in prostate cancer, suppress CD44 translation and enhance invasion of prostate cancer cells in vitro. Int. J. Clin. Exp. Pathol. 2009, 2, 361–369. [Google Scholar] [PubMed]
  12. Ribas, J.; Ni, X.; Haffner, M.; Wentzel, E.A.; Salmasi, A.H.; Chowdhury, W.H.; Kudrolli, T.A.; Yegnasubramanian, S.; Luo, J.; Rodriguez, R.; et al. miR-21: An androgen receptor–regulated microRNA that promotes hormone-dependent and hormone-independent prostate cancer growth. Cancer Res. 2009, 69, 7165–7169. [Google Scholar] [CrossRef] [PubMed]
  13. Xu, G.; Wu, J.; Zhou, L.; Chen, B.; Sun, Z.; Zhao, F.; Tao, Z. Characterization of the small RNA transcriptomes of androgen dependent and independent prostate cancer cell line by deep sequencing. PLoS ONE 2010, 5, e15519. [Google Scholar] [CrossRef] [PubMed]
  14. Friedlander, T.W.; Roy, R.; Tomlins, S.A.; Ngo, V.T.; Kobayashi, Y.; Azameera, A.; Rubin, M.A.; Pienta, K.J.; Chinnaiyan, A.; Ittmann, M.M.; et al. Common structural and epigenetic changes in the genome of castration-resistant prostate cancer. Cancer Res. 2012, 72, 616–625. [Google Scholar] [CrossRef] [PubMed]
  15. Mishra, D.K.; Chen, Z.; Wu, Y.; Sarkissyan, M.; Koeffler, H.P.; Vadgama, J.V. Global methylation pattern of genes in androgen-sensitive and androgen-independent prostate cancer cells. Mol. Cancer Ther. 2010, 9, 33–45. [Google Scholar] [CrossRef] [PubMed]
  16. Taiwo, O.; Wilson, G.A.; Morris, T.; Seisenberger, S.; Reik, W.; Pearce, D.; Beck, S.; Butcher, L.M. Methylome analysis using MeDIP-seq with low DNA concentrations. Nat. Protoc. 2012, 7, 617–636. [Google Scholar] [CrossRef] [PubMed]
  17. Harris, R.A.; Wang, T.; Coarfa, C.; Nagarajan, R.P.; Hong, C.; Downey, S.L.; Johnson, B.E.; Fouse, S.D.; Delaney, A.; Zhao, Y.; et al. Comparison of sequencing-based methods to profile DNA methylation and identification of monoallelic epigenetic modifications. Nat. Biotechnol. 2010, 28, 1097–1105. [Google Scholar] [CrossRef] [PubMed]
  18. Cai, W.; Mao, F.; Teng, H.; Cai, T.; Zhao, F.; Wu, J.; Sun, Z.S. MBRidge: An accurate and cost-effective method for profiling DNA methylome at single-base resolution. J. Mol. Cell Biol. 2015, 7, 299–313. [Google Scholar] [CrossRef] [PubMed]
  19. Zheng, Z.; Advani, A.; Melefors, O.; Glavas, S.; Nordström, H.; Ye, W.; Engstrand, L.; Andersson, A.F. Titration-free 454 sequencing using Y adapters. Nat. Protoc. 2011, 6, 1367–1376. [Google Scholar] [CrossRef] [PubMed]
  20. Ruike, Y.; Imanaka, Y.; Sato, F.; Shimizu, K.; Tsujimoto, G. Genome-wide analysis of aberrant methylation in human breast cancer cells using methyl-DNA immunoprecipitation combined with high-throughput sequencing. BMC Genom. 2010, 11, 137. [Google Scholar] [CrossRef] [PubMed]
  21. Weber, M.; Davies, J.J.; Wittig, D.; Oakeley, E.J.; Haase, M.; Lam, W.L.; Schuebeler, D. Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. Nat. Genet. 2005, 37, 853–862. [Google Scholar] [CrossRef] [PubMed]
  22. Krueger, F.; Andrews, S.R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 2011, 27, 1571–1572. [Google Scholar] [CrossRef] [PubMed]
  23. Benjamini, Y.; Hochberg, Y. On the adaptive control of the false discovery rate in multiple testing with independent statistics. J. Educ. Behav. Stat. 2000, 25, 60–83. [Google Scholar] [CrossRef]
  24. Bailey, T.L.; Johnson, J.; Grant, C.E.; Noble, W.S. The MEME Suite. Nucleic Acids Res. 2015, 43, W39–W49. [Google Scholar] [CrossRef] [PubMed]
  25. Mathelier, A.; Zhao, X.; Zhang, A.W.; Parcy, F.; Worsley-Hunt, R.; Arenillas, D.J.; Buchman, S.; Chen, C.; Chou, A.; Ienasescu, H.; et al. JASPAR 2014: An extensively expanded and updated open-access database of transcription factor binding profiles. Nucleic Acids Res. 2014, 42, D142–D147. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Tansley, D.S. Linux and Unix Shell Programming; Addison-Wesley Professional: Boston, MA, USA, 2000. [Google Scholar]
  27. Kanehisa, M.; Goto, S.; Sato, Y.; Furumichi, M.; Tanabe, M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012, 40, D109–D114. [Google Scholar] [CrossRef] [PubMed]
  28. Jupe, S.; Akkerman, J.W.; Soranzo, N.; Ouwehand, W.H. Reactome-a curated knowledgebase of biological pathways: Megakaryocytes and platelets. J. Thromb. Haemost. 2012, 10, 2399–2402. [Google Scholar] [CrossRef] [PubMed]
  29. Kelder, T.; van Iersel, M.P.; Hanspers, K.; Kutmon, M.; Conklin, B.R.; Evelo, C.T.; Pico, A.R. WikiPathways: Building research communities on biological pathways. Nucleic Acids Res. 2012, 40, D1301–D1307. [Google Scholar] [CrossRef] [PubMed]
  30. Schaefer, C.F.; Anthony, K.; Krupa, S.; Buchoff, J.; Day, M.; Hannay, T.; Buetow, K.H. PID: The pathway interaction database. Nucleic Acids Res. 2009, 37, D674–D679. [Google Scholar] [CrossRef] [PubMed]
  31. Gene Ontology Consortium. Gene ontology annotations and resources. Nucleic Acids Res. 2013, 41, D530–D535. [Google Scholar]
  32. Shannon, P.; Richards, M. MotifDb: An Annotated Collection of Protein-DNA Binding Sequence Motifs. R package version 1.2. 2014. Available online: http://bioconductor.org/packages/release/bioc/html/MotifDb.html (accessed on 30 November 2017).
  33. Franceschini, A.; Szklarczyk, D.; Frankild, S.; Kuhn, M.; Simonovic, M.; Roth, A.; Lin, J.; Minguez, P.; Bork, P.; von Mering, C.; et al. STRING v9.1: Protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013, 41, D808–D815. [Google Scholar] [CrossRef] [PubMed]
  34. Smoot, M.E.; Ono, K.; Ruscheinski, J.; Wang, P.-L.; Ideker, T. Cytoscape 2.8: New features for data integration and network visualization. Bioinformatics 2011, 27, 431–432. [Google Scholar] [CrossRef] [PubMed]
  35. Wong, N.; Wang, X. miRDB: An online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 2015, 43, D146–D152. [Google Scholar] [CrossRef] [PubMed]
  36. Storey, J.D.; Xiao, W.; Leek, J.T.; Tompkins, R.G.; Davis, R.W. Significance analysis of time course microarray experiments. Proc. Natl. Acad. Sci. USA 2005, 102, 12837–12842. [Google Scholar] [CrossRef] [PubMed]
  37. D’Antonio, J.M.; Ma, C.; Monzon, F.A.; Pflug, B.R. Longitudinal analysis of androgen deprivation of prostate cancer cells identifies pathways to androgen independence. Prostate 2008, 68, 698–714. [Google Scholar] [CrossRef] [PubMed]
  38. Gautier, L.; Cope, L.; Bolstad, B.M.; Irizarry, R.A. affy—Analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 2004, 20, 307–315. [Google Scholar] [CrossRef] [PubMed]
  39. Futschik, M.E.; Carlisle, B. Noise-robust soft clustering of gene expression time-course data. J. Bioinform. Comput. Biol. 2005, 3, 965–988. [Google Scholar] [CrossRef] [PubMed]
  40. Kumar, L.; Futschik, M.E. Mfuzz: A software package for soft clustering of microarray data. Bioinformation 2007, 2, 5–7. [Google Scholar] [CrossRef] [PubMed]
  41. Kolde, R. Pheatmap: Pretty Heatmaps. R package version 0.7.7. Available online: https://cran.r-project.org/src/contrib/Archive/pheatmap/ (accessed on 30 November 2017).
  42. Oliveros, J.C. VENNY—An Interactive Tool for Comparing Lists with Venn Diagrams. 2007. Available online: http://bioinfogp.cnb.csic.es/tools/venny/index.html (accessed on 30 November 2017).
  43. Lin, P.C.; Giannopoulou, E.G.; Park, K.; Mosquera, J.M.; Sboner, A.; Tewari, A.K.; Garraway, L.A.; Beltran, H.; Rubin, M.A.; Elemento, O.; et al. Epigenomic alterations in localized and advanced prostate cancer. Neoplasia 2013, 15, 373–383. [Google Scholar] [CrossRef] [PubMed]
  44. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  45. Bracken, C.P.; Gregory, P.A.; Kolesnikoff, N.; Bert, A.G.; Wang, J.; Shannon, M.F.; Goodall, G.J. A double-negative feedback loop between ZEB1-SIP1 and the microRNA-200 family regulates epithelial-mesenchymal transition. Cancer Res. 2008, 68, 7846–7854. [Google Scholar] [CrossRef] [PubMed]
  46. Qi, H.; Grenier, J.; Fournier, A.; Labrie, C. Androgens differentially regulate the expression of NEDD4L transcripts in LNCaP human prostate cancer cells. Mol. Cell. Endocrinol. 2003, 210, 51–62. [Google Scholar] [CrossRef] [PubMed]
  47. Gao, C.; Pang, L.; Ren, C.; Ma, T. Decreased expression of Nedd4L correlates with poor prognosis in gastric cancer patient. Med. Oncol. 2012, 29, 1733–1738. [Google Scholar] [CrossRef] [PubMed]
  48. Tanksley, J.P.; Chen, X.; Coffey, R.J. NEDD4L is downregulated in colorectal cancer and inhibits canonical WNT signaling. PLoS ONE 2013, 8, e81514. [Google Scholar] [CrossRef] [PubMed]
  49. Qu, M.-H.; Han, C.; Srivastava, A.K.; Cui, T.; Zou, N.; Gao, Z.-Q.; Wang, Q.E. MiR-93 promotes TGF-β-induced epithelial-to-mesenchymal transition through downregulation of NEDD4L in lung cancer cells. Tumor Biol. 2016, 37, 5645–5651. [Google Scholar] [CrossRef] [PubMed]
  50. Ramberg, H.; Alshbib, A.; Berge, V.; Svindland, A.; Tasken, K.A. Regulation of PBX3 expression by androgen and Let-7d in prostate cancer. Mol. Cancer 2011, 10, 50. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Li, Y.; Sun, Z.; Zhu, Z.; Zhang, J.; Sun, X.; Xu, H. PBX3 is overexpressed in gastric cancer and regulates cell proliferation. Tumor Biol. 2014, 35, 4363–4368. [Google Scholar] [CrossRef] [PubMed]
  52. Han, H.-B.; Gu, J.; Ji, D.-B.; Li, Z.-W.; Zhang, Y.; Zhao, W.; Wang, L.M.; Zhang, Z.Q. PBX3 promotes migration and invasion of colorectal cancer cells via activation of MAPK/ERK signaling pathway. World J. Gastroenterol. 2014, 20, 18260. [Google Scholar] [CrossRef] [PubMed]
  53. Ma, S.; Chan, Y.P.; Woolcock, B.; Hu, L.; Wong, K.Y.; Ling, M.T.; Bainbridge, T.; Webber, D.; Chan, T.H.M.; Guan, X.Y.; et al. DNA fingerprinting tags novel altered chromosomal regions and identifies the involvement of SOX5 in the progression of prostate cancer. Int. J. Cancer 2009, 124, 2323–2332. [Google Scholar] [CrossRef] [PubMed]
  54. Huang, D.Y.; Lin, Y.T.; Jan, P.S.; Hwang, Y.C.; Liang, S.T.; Peng, Y.; Huang, C.Y.; Wu, H.C.; Lin, C.T. Transcription factor SOX-5 enhances nasopharyngeal carcinoma progression by down-regulating SPARC gene expression. J. Pathol. 2008, 214, 445–455. [Google Scholar] [CrossRef] [PubMed]
  55. Wang, D.; Han, S.; Wang, X.; Peng, R.; Li, X. SOX5 promotes epithelial-mesenchymal transition and cell invasion via regulation of Twist1 in hepatocellular carcinoma. Med. Oncol. 2015, 32, 461. [Google Scholar] [CrossRef] [PubMed]
  56. Wang, H.; McKnight, N.C.; Zhang, T.; Lu, M.L.; Balk, S.P.; Yuan, X. SOX9 is expressed in normal prostate basal cells and regulates androgen receptor expression in prostate cancer cells. Cancer Res. 2007, 67, 528–536. [Google Scholar] [CrossRef] [PubMed]
  57. Wang, H.; Leav, I.; Ibaragi, S.; Wegner, M.; Hu, G.F.; Lu, M.L.; Balk, S.P.; Yuan, X. SOX9 is expressed in human fetal prostate epithelium and enhances prostate cancer invasion. Cancer Res. 2008, 68, 1625–1630. [Google Scholar] [CrossRef] [PubMed]
  58. Murugan, A.K.; Dong, J.; Xie, J.; Xing, M. Uncommon GNAQ, MMP8, AKT3, EGFR, and PIK3R1 mutations in thyroid cancers. Endocr. Pathol. 2011, 22, 97–102. [Google Scholar] [CrossRef] [PubMed]
  59. Feng, X.; Degese, M.S.; Iglesias-Bartolome, R.; Vaque, J.P.; Molinolo, A.A.; Rodrigues, M.; Zaidi, M.R.; Ksander, B.R.; Merlino, G.; Sodhi, A.; et al. Hippo-independent activation of YAP by the GNAQ uveal melanoma oncogene through a trio-regulated rho GTPase signaling circuitry. Cancer Cell 2014, 25, 831–845. [Google Scholar] [CrossRef] [PubMed]
  60. Suzuki, C.; Daigo, Y.; Ishikawa, N.; Kato, T.; Hayama, S.; Ito, T.; Tsuchiya, E.; Nakamura, Y. ANLN plays a critical role in human lung carcinogenesis through the activation of RHOA and by involvement in the phosphoinositide 3-kinase/AKT pathway. Cancer Res. 2005, 65, 11314–11325. [Google Scholar] [CrossRef] [PubMed]
  61. Zhou, W.; Wang, Z.; Shen, N.; Pi, W.; Jiang, W.; Huang, J.; Hu, Y.; Li, X.; Sun, L. Knockdown of ANLN by lentivirus inhibits cell growth and migration in human breast cancer. Mol. Cell. Biochem. 2015, 398, 11–19. [Google Scholar] [CrossRef] [PubMed]
  62. Tamura, K.; Furihata, M.; Tsunoda, T.; Ashida, S.; Takata, R.; Obara, W.; Yoshioka, H.; Daigo, Y.; Nasu, Y.; Kumon, H.; et al. Molecular features of hormone-refractory prostate cancer cells by genome-wide gene expression profiles. Cancer Res. 2007, 67, 5117–5125. [Google Scholar] [CrossRef] [PubMed]
  63. Kamai, T.; Yamanishi, T.; Shirataki, H.; Takagi, K.; Asami, H.; Ito, Y.; Yoshida, K.I. Overexpression of RhoA, Rac1, and Cdc42 GTPases is associated with progression in testicular cancer. Clin. Cancer Res. 2004, 10, 4799–4805. [Google Scholar] [CrossRef] [PubMed]
  64. Zhou, C.; Ling, M.T.; Kin-Wah Lee, T.; Man, K.; Wang, X.; Wong, Y.C. FTY720, a fungus metabolite, inhibits invasion ability of androgen-independent prostate cancer cells through inactivation of RhoA-GTPase. Cancer Lett. 2006, 233, 36–47. [Google Scholar] [CrossRef] [PubMed]
  65. Mongroo, P.S.; Rustgi, A.K. The role of the miR-200 family in epithelial-mesenchymal transition. Cancer Biol. Ther. 2010, 10, 219–222. [Google Scholar] [CrossRef] [PubMed]
  66. Chen, J.; Wang, L.; Matyunina, L.V.; Hill, C.G.; McDonald, J.F. Overexpression of miR-429 induces mesenchymal-to-epithelial transition (MET) in metastatic ovarian cancer cells. Gynecol. Oncol. 2011, 121, 200–205. [Google Scholar] [CrossRef] [PubMed]
  67. Li, J.; Du, L.; Yang, Y.; Wang, C.; Liu, H.; Wang, L.; Zhang, X.; Li, W.; Zheng, G.; Dong, Z. MiR-429 is an independent prognostic factor in colorectal cancer and exerts its anti-apoptotic function by targeting SOX2. Cancer Lett. 2013, 329, 84–90. [Google Scholar] [CrossRef] [PubMed]
  68. Daigo, K.; Takano, A.; Thang, P.M.; Yoshitake, Y.; Shinohara, M.; Tohnai, I.; Murakami, Y.; Maegawa, J.; Daigo, Y. Characterization of KIF11 as a novel prognostic biomarker and therapeutic target for oral cancer. Int. J. Oncol. 2018, 52, 155–165. [Google Scholar] [CrossRef] [PubMed]
  69. Imai, T.; Oue, N.; Nishioka, M.; Mukai, S.; Oshima, T.; Sakamoto, N.; Sentani, K.; Matsusaki, K.; Yoshida, K.; Yasui, W. Overexpression of KIF11 in gastric cancer with intestinal mucin phenotype. Pathobiology 2017, 84, 16–24. [Google Scholar] [CrossRef] [PubMed]
  70. Hayashi, N.; Koller, E.; Fazli, L.; Gleave, M.E. The role of the miR-200 family in epithelial-mesenchymal transition. Prostate 2008, 68, 1283–1295. [Google Scholar] [CrossRef] [PubMed]
Figure 1. DMRs between lymph node carcinoma of the cancer prostate—artificially induced—flutamide (LNCaP-AI-F) and LNCaP cells. (A) Distribution of DMRs by gene structure; (B) Distribution of DMRs by functional elements; (C) Functions of genes whose promoters overlapped with DMRs. DMRs: Differential methylated regions; CGI: CpG island; UTR: untranslated region; CDS: coding sequence; Downstream: 2 kilobase downstream region of genes; TSS: transcription start site; TTS: transcription termination site; Upstream: 2-kilobase upstream region of genes. In human cells, exons include the CDS and 3′UTR or 5′UTR.
Figure 1. DMRs between lymph node carcinoma of the cancer prostate—artificially induced—flutamide (LNCaP-AI-F) and LNCaP cells. (A) Distribution of DMRs by gene structure; (B) Distribution of DMRs by functional elements; (C) Functions of genes whose promoters overlapped with DMRs. DMRs: Differential methylated regions; CGI: CpG island; UTR: untranslated region; CDS: coding sequence; Downstream: 2 kilobase downstream region of genes; TSS: transcription start site; TTS: transcription termination site; Upstream: 2-kilobase upstream region of genes. In human cells, exons include the CDS and 3′UTR or 5′UTR.
Genes 09 00032 g001aGenes 09 00032 g001b
Figure 2. TF target and miRNA–target regulatory networks. (A) TF-target network. Red diamond: up-regulated TF; green diamond: down-regulated TF; grey diamond: TF that is not differentially expressed between LNCaP-AI-F and LNCaP cells; red circle: up-regulated PMDEG; green circle: down-regulated PMDEG; red line with arrow: a TF targets a PMDEG and the TF-binding motif overlaps with a hypermethylated DMR; green line: a TF targets a PMDEG and the TF-binding motif overlaps with a hypomethylated DMR; dotted line: the confirmed regulatory relationships by data validation using the GSE41701 dataset; (B) bmiRNA–target network. Yellow rectangle: down-regulated MDEmiRNA; red circle: up-regulated MDEG; green circle: down-regulated MDEG; yellow line with T-shape: a MDEmiRNA targets a MDEG; grey line: protein-protein interaction between two MDEGs. TF: transcription factor; miRNA: microRNA; PMDEG: differentially expressed gene whose promoter overlapped with DMR(s); MDEmiRNA: differentially expressed microRNA that overlapped with DMR(s); MDEG: differentially expressed gene that overlapped with DMR(s); DMR: differentially methylated region.
Figure 2. TF target and miRNA–target regulatory networks. (A) TF-target network. Red diamond: up-regulated TF; green diamond: down-regulated TF; grey diamond: TF that is not differentially expressed between LNCaP-AI-F and LNCaP cells; red circle: up-regulated PMDEG; green circle: down-regulated PMDEG; red line with arrow: a TF targets a PMDEG and the TF-binding motif overlaps with a hypermethylated DMR; green line: a TF targets a PMDEG and the TF-binding motif overlaps with a hypomethylated DMR; dotted line: the confirmed regulatory relationships by data validation using the GSE41701 dataset; (B) bmiRNA–target network. Yellow rectangle: down-regulated MDEmiRNA; red circle: up-regulated MDEG; green circle: down-regulated MDEG; yellow line with T-shape: a MDEmiRNA targets a MDEG; grey line: protein-protein interaction between two MDEGs. TF: transcription factor; miRNA: microRNA; PMDEG: differentially expressed gene whose promoter overlapped with DMR(s); MDEmiRNA: differentially expressed microRNA that overlapped with DMR(s); MDEG: differentially expressed gene that overlapped with DMR(s); DMR: differentially methylated region.
Genes 09 00032 g002aGenes 09 00032 g002b
Figure 3. Clustering analyses of genes in LNCaP cells during androgen deprivation. (A) Clustering analysis based on Mfuzz software. Red: most close to the clustering center; blue: moderately close to the clustering center; green: least close to the clustering center; (B) Hierarchical clustering analysis based on pheatmap software and Mfuzz clusters. Blue: low messenger RNA (mRNA) level; red: high mRNA level. Scale bar represents the duration of androgen deprivation and pheatmap-clusters a, b, c, d, e, and f correspond to Mfuzz clusters 1, 2, 3, 4, 5, and 6, respectively.
Figure 3. Clustering analyses of genes in LNCaP cells during androgen deprivation. (A) Clustering analysis based on Mfuzz software. Red: most close to the clustering center; blue: moderately close to the clustering center; green: least close to the clustering center; (B) Hierarchical clustering analysis based on pheatmap software and Mfuzz clusters. Blue: low messenger RNA (mRNA) level; red: high mRNA level. Scale bar represents the duration of androgen deprivation and pheatmap-clusters a, b, c, d, e, and f correspond to Mfuzz clusters 1, 2, 3, 4, 5, and 6, respectively.
Genes 09 00032 g003
Figure 4. Comprehensive analysis of sequencing and microarray data. (A) Genes in clusters 1, 2, and 6 were compared with the 4653 up-regulated DEGs; (B) Genes in clusters 3 and 5 were compared with the 4074 down-regulated DEGs; (C) Genes in clusters 1, 2, and 6 were compared with the up-regulated PMDEGs in the TF-target network; (D) Genes in clusters 3 and 5 were compared with the down-regulated PMDEGs in the TF-target network; (E) Genes in clusters 1, 2, and 6 were compared with the up-regulated MDEGs in the miRNA–target network; (F) Genes in clusters 3 and 5 were compared with the down-regulated MDEGs in the miRNA–target network. DEGs: differentially expressed genes; TF: transcription factor; miRNA: microRNA; PMDEGs: differentially expressed genes whose promoters overlapped with DMRs; MDEGs: differentially expressed genes that overlapped with DMRs; DMRs: differentially methylated regions.
Figure 4. Comprehensive analysis of sequencing and microarray data. (A) Genes in clusters 1, 2, and 6 were compared with the 4653 up-regulated DEGs; (B) Genes in clusters 3 and 5 were compared with the 4074 down-regulated DEGs; (C) Genes in clusters 1, 2, and 6 were compared with the up-regulated PMDEGs in the TF-target network; (D) Genes in clusters 3 and 5 were compared with the down-regulated PMDEGs in the TF-target network; (E) Genes in clusters 1, 2, and 6 were compared with the up-regulated MDEGs in the miRNA–target network; (F) Genes in clusters 3 and 5 were compared with the down-regulated MDEGs in the miRNA–target network. DEGs: differentially expressed genes; TF: transcription factor; miRNA: microRNA; PMDEGs: differentially expressed genes whose promoters overlapped with DMRs; MDEGs: differentially expressed genes that overlapped with DMRs; DMRs: differentially methylated regions.
Genes 09 00032 g004
Figure 5. Relative mRNA expression levels of MDEGs and miR-429 determined by quantitative real-time PCR. The expression levels of genes were normalized against glyceraldehyde-3-phosphate dehydrogenase, and the expression level of miR-429 was normalized against U6 small nuclear RNA. * p < 0.05. MDEGs: differentially expressed genes that overlapped with differentially methylated regions.
Figure 5. Relative mRNA expression levels of MDEGs and miR-429 determined by quantitative real-time PCR. The expression levels of genes were normalized against glyceraldehyde-3-phosphate dehydrogenase, and the expression level of miR-429 was normalized against U6 small nuclear RNA. * p < 0.05. MDEGs: differentially expressed genes that overlapped with differentially methylated regions.
Genes 09 00032 g005
Table 1. Transcription factor (TF)-binding motifs in differentially methylated regions (DMR)-promoter-overlapping regions.
Table 1. Transcription factor (TF)-binding motifs in differentially methylated regions (DMR)-promoter-overlapping regions.
MotifSiteWidthE-ValueKnown or Similar Motifs
167303.00 × 10−69IRF1 (MA0050.2)
FOXP1 (MA0481.1)
2170952.30 × 10−18STAT1 (MA0137.3)
STAT4 (MA0518.1)
STAT2::STAT1 (MA0517.1)
397291.40 × 10−15ZNF263 (MA0528.1)
EGR1 (MA0162.2)
4120751.60 × 10−9ZEB1 (MA0103.2)
T (MA0009.1)
SMAD2::SMAD3::SMAD4 (MA0513.1)
5161551.30 × 10−5GATA4 (MA0482.1)
6100062.90 × 10−5NHLH1 (MA0048.1)
MYOG (MA0500.1)
PAX5 (MA0014.2)
710481.20 × 10−3EGR1 (MA0162.2)
813291.70 × 10−3SOX5 (MA0087.1)
94784.20 × 10−3RUNX2 (MA0511.1)
1051057.50 × 10−3ARID3A (MA0151.1)
Letters in bold represent TFs.
Table 2. Pathways (all) and functions (top 5) enriched by MDEGs.
Table 2. Pathways (all) and functions (top 5) enriched by MDEGs.
MDEGsCategoryTerm IDTerm DescriptionAdjusted p-ValueGene Count
UpPathwayR-HSA-69278Cell Cycle, Mitotic [R-HSA-69278]4.29 × 10−685
PathwayR-HSA-1640170Cell Cycle [R-HSA-1640170]1.91 × 10−595
PathwayR-HSA-453274Mitotic G2-G2/M phases [R-HSA-453274]3.59 × 10−328
PathwayR-HSA-1266738Developmental Biology [R-HSA-1266738]6.80 × 10−3102
PathwayR-HSA-69275G2/M Transition [R-HSA-69275]8.10 × 10−327
Pathwayvegfr1_2_pathwaySignaling events mediated by VEGFR1 and VEGFR2 [vegfr1_2_pathway]1.13 × 10−220
PathwayR-HSA-68877Mitotic Prometaphase [R-HSA-68877]1.27 × 10−226
Pathwayhsa04520Adherens junction [hsa04520]1.43 × 10−220
GO BPGO:0000278mitotic cell cycle [GO:0000278]4.75 × 10−13153
GO BPGO:0007049cell cycle [GO:0007049]1.61 × 10−9199
GO BPGO:0022402cell cycle process [GO:0022402]2.19 × 10−9170
GO BPGO:1903047mitotic cell cycle process [GO:1903047]6.39 × 10−9119
GO BPGO:0090304nucleic acid metabolic process [GO:0090304]2.52 × 10−7454
GO CCGO:0005622intracellular [GO:0005622]1.99 × 10−181144
GO CCGO:0044424intracellular part [GO:0044424]1.31 × 10−171124
GO CCGO:0031981nuclear lumen [GO:0031981]1.42 × 10−17448
GO CCGO:0005654nucleoplasm [GO:0005654]2.48 × 10−17392
GO CCGO:0044428nuclear part [GO:0044428]4.85 × 10−17472
GO MFGO:0005515protein binding [GO:0005515]1.75 × 10−121042
GO MFGO:0005488binding [GO:0005488]1.61 × 10−101141
GO MFGO:0003723RNA binding [GO:0003723]3.90 × 10−7198
GO MFGO:0003676nucleic acid binding [GO:0003676]8.01 × 10−7310
GO MFGO:0044822poly(A) RNA binding [GO:0044822]3.15 × 10−6172
DownPathwayhsa04144Endocytosis [hsa04144]4.00 × 10−452
Pathwayhdac_classi_pathwaySignaling events mediated by HDAC Class I [hdac_classi_pathway]4.83 × 10−320
Pathwayhsa04666Fc gamma R-mediated phagocytosis [hsa04666]9.38 × 10−324
GO BPGO:0016192vesicle-mediated transport [GO:0016192]7.26 × 10−5154
GO BPGO:0051179localization [GO:0051179]9.44 × 10−5511
GO BPGO:0006810transport [GO:0006810]5.10 × 10−4419
GO BPGO:0051234establishment of localization [GO:0051234]3.12 × 10−3423
GO BPGO:0048518positive regulation of biological process [GO:0048518]3.96 × 10−3490
GO CCGO:0044424intracellular part [GO:0044424]3.47 × 10−101156
GO CCGO:0005622intracellular [GO:0005622]2.75 × 10−81167
GO CCGO:0005737cytoplasm [GO:0005737]4.37 × 10−6901
GO CCGO:0030054cell junction [GO:0030054]3.72 × 10−5104
GO CCGO:0012505endomembrane system [GO:0012505]1.57 × 10−3327
GO MFGO:0016773phosphotransferase activity, alcohol group as acceptor [GO:0016773]1.58 × 10−390
GO MFGO:0016772transferase activity, transferring phosphorus-containing groups [GO:0016772]2.30 × 10−3108
GO MFGO:0016301kinase activity [GO:0016301]2.52 × 10−396
GO MFGO:0000166nucleotide binding [GO:0000166]1.29 × 10−278
GO MFGO:1901265nucleoside phosphate binding [GO:1901265]1.40 × 10−278
MDEGs: differentially expressed genes with differentially methylated region; GO: gene ontology; BP: biological process; CC: cellular component; MF: molecular function; ID: identifier.

Share and Cite

MDPI and ACS Style

Wang, Y.; Qin, T.; Hu, W.; Chen, B.; Dai, M.; Xu, G. Genome-Wide Methylation Patterns in Androgen-Independent Prostate Cancer Cells: A Comprehensive Analysis Combining MeDIP-Bisulfite, RNA, and microRNA Sequencing Data. Genes 2018, 9, 32. https://0-doi-org.brum.beds.ac.uk/10.3390/genes9010032

AMA Style

Wang Y, Qin T, Hu W, Chen B, Dai M, Xu G. Genome-Wide Methylation Patterns in Androgen-Independent Prostate Cancer Cells: A Comprehensive Analysis Combining MeDIP-Bisulfite, RNA, and microRNA Sequencing Data. Genes. 2018; 9(1):32. https://0-doi-org.brum.beds.ac.uk/10.3390/genes9010032

Chicago/Turabian Style

Wang, Yumin, Tingting Qin, Wangqiang Hu, Binghua Chen, Meijie Dai, and Gang Xu. 2018. "Genome-Wide Methylation Patterns in Androgen-Independent Prostate Cancer Cells: A Comprehensive Analysis Combining MeDIP-Bisulfite, RNA, and microRNA Sequencing Data" Genes 9, no. 1: 32. https://0-doi-org.brum.beds.ac.uk/10.3390/genes9010032

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop