Identification of mutations, gene expression changes and fusion transcripts by whole transcriptome RNAseq in docetaxel resistant prostate cancer cells

Docetaxel has been the standard first-line therapy in metastatic castration resistant prostate cancer. The survival benefit is, however, limited by either primary or acquired resistance. In this study, Du145 prostate cancer cells were converted to docetaxel-resistant cells Du145-R and Du145-RB by in vitro culturing. Next generation RNAseq was employed to analyze these cell lines. Forty-two genes were identified to have acquired mutations after the resistance development, of which thirty-four were found to have mutations in published sequencing studies using prostate cancer samples from patients. Fourteen novel and 2 previously known fusion genes were inferred from the RNA-seq data, and 13 of these were validated by RT-PCR and/or re-sequencing. Four in-frame fusion transcripts could be transcribed into fusion proteins in stably transfected HEK293 cells, including MYH9-EIF3D and LDLR-RPL31P11, which were specific identified or up-regulated in the docetaxel resistant DU145 cells. A panel of 615 gene transcripts was identified to have significantly changed expression profile in the docetaxel resistant cells. These transcriptional changes have potential for further study as predictive biomarkers and as targets of docetaxel treatment. Electronic supplementary material The online version of this article (doi:10.1186/s40064-016-3543-0) contains supplementary material, which is available to authorized users.


Background
Most metastatic prostate cancers respond to androgen deprivation therapy (ADT) but eventually develop castration resistance and become metastatic castration resistant prostate cancers (mCRPC) about 24-36 months after the treatment start (Harris et al. 2009;Attar et al. 2009;Watson et al. 2010). mCRPC is the major cause of cancer death in prostate cancer patients. Median survival time of patients with mCRPCs is 16-18 months from the start of progression (Amaral et al. 2012). Docetaxel chemotherapy can further prolong the median overall survival by 3-5 months (Galsky et al. 2012). However, docetaxel resistance is a critical problem because half of patients will not respond to docetaxel treatment (intrinsic resistance), while the other half, which responds initially, become resistant ultimately (acquired resistance) (Tannock et al. 2004). Failure of docetaxel treatment has been thought to be caused by either intrinsic or acquired resistance.
Docetaxel is a member of taxane family and widely been used to treat mCRPC patients. Docetaxel induces cancer cell death by binding β-tubulin, stabilizing microtubule assembly, suppressing dynamics of individual micro-tubules in G2-M phase tumor cells and preventing disassembly (Yvon et al. 1999;Eisenhauer and Vermorken 1998). Despite a decade of clinical use, the mechanism of resistance to docetaxel has not been fully investigated and there are no clinically reliable biomarkers to predict the drug resistance. Limited data suggests that the resistance may be caused by the following mechanisms: Open Access *Correspondence: chunde.li@ki.se 1 Department of Oncology-Pathology, Karolinska Institutet, Stockholm, Sweden Full list of author information is available at the end of the article (1) decreased drug concentration due to high expression of drug export pump proteins ABCB1, ABCB4, ABCC1 (Gottesman et al. 2002); (2) mutations in the drug targets (Berrieman et al. 2004); (3) inhibition of apoptotic pathways (Bhalla 2003); (4) altered expression profile of tubulins or microtubule-associated proteins (MAPs) (Seve and Dumontet 2008;Verrills et al. 2006). So far, only a few drugs have been developed with modest survival benefit in docetaxel resistant mCRPC.
This study applied next generation RNA sequencing (RNAseq) technology in combination with specific software (Ozsolak and Milos 2011) to determine gene expression changes, mutations and fusions in docetaxel sensitive cell lines versus docetaxel resistant cell lines. The comparison between these cell lines identified a panel of genes potentially involved in the development of docetaxel resistance. The clinical importance was further addressed by comparing with published RNA sequencing results in prostate cancer samples from patients.

Mutations acquired in docetaxel resistant cell lines
We generated docetaxel resistant variants of Du145 prostate cancer cells as described in M&M. We used triplicates of each cell line (Du145, Du145-R and Du145-RB) for whole transcriptome RNA-sequencing and found 4864 mutations totally (Additional file 1). We compared TaxR (docetaxel resistant) and TaxS (docetaxel sensitive) cell lines to find mutations acquired after docetaxel treatment. Only mutations, which were absent in TaxS (Du145) but present in all Du145-R triplicates and Du145-RB triplicates, were chosen as "stably acquired mutations". Forty-two such mutations were identified (Table 1) and 4 randomly selected mutations were validated by PCR followed by SANGER sequencing (Fig. 1).
By matching with previously published whole transcriptome analyses, we could identify that 34 of these genes had mutations in prostate cancer samples from patients (Table 1) (Robinson et al. 2015). For many genes, e.g. ABCB2, there are published data to support their importance in the development of drug resistance (Aberuyi et al. 2014;Rahgozar et al. 2014).

Fusion transcript detection and validation
ChimeraScan software was employed to find fusions from RNAseq data. Selecting gene-gene pairs supported by two or more unique alignment reads provided an initial list of 48, 75 and 66 fusion candidates in DU145, Du145-R and Du145-RB cell lines respectively (Additional file 2). We validated all fusion candidates that had a ChimeraScan score above 5 in at least 1 out of 3 cell lines. Of 16 fusion candidates selected (Table 2), 13 (81.25 %) were verified by Reverse Transcription PCR (RT-PCR) with primers covering the fusion break points (Additional file 3), and 5 of validated genes were further verified by Sanger sequencing in the Du145, Du145-R and Du145-RB cell lines (Fig. 2). Two gene fusions had been found by previous studies: UBE2L3-KRAS  expressed in all three cell lines and TAF15-AP2B1 (http://54.84.12.177/PanCanFusV2/Fusions!fusion) specific expressed in Du145 (Additional file 3). The other fourteen fusions were novel discovered. Figure 3a showed that all chromosomes were involved in gene fusion except chr 21 and chr Y. The two largest fusion groups were distributed in chr1 (27, 14.3 %) and chr6 (44, 23.3 %), and most of the fusions were intra-chromosomal (26 out of 27 in chr1; 39 out of 44 in chr6). Of the 16 chosen fusion candidates (3 of them could not be validated by PCR), 10 of them were commonly expressed in all 3 cell lines, one expressed only in TaxR cell line (MYH9-EIF3D) and 2 specifically in TaxS cell line (TAF15-AP2B1, VCL-ADK) ( Table 2; Fig. 3b). Among 10 commonly expressed fusions, two were upregulated in TaxR cell lines compared to TaxS cell line (LDLR-RPL31P11, SRGAP2P2-SRGAP2). Eight out of 16 were predicted to be in-frame suggesting their potential to produce functional fusion proteins (Table 2).
Four fusion candidates were validated by qPCR in Du145, Du145-R and Du145-RB cells, as well as at the protein level by western blot in plasmid transfected HEK293 cells (Fig. 3c, d), but their translation into protein could not be validated by western blot in Du145, Du145-R and Du145-RB cell lines.
Interestingly, when we validated the VCL-ADK fusion candidate by PCR, we found that there were two bands in the same PCR lane (Fig. 3e). Sanger sequencing results showed that both of the two bands were VCL-ADK fusions. The upper band was a fusion between VCL and ADK variant 1, 2 and 3, while the lower band was another fusion with ADK variant 4. Western blot showed that both fusions (VCL-ADK variant 1, 2, 3 and VCL-ADK variant 4) could be detected as protein in plasmid transfected HEK293 cell lines (Fig. 3d).

Identification of stably up-or down-regulated genes in the TaxR cell lines
Using gene expression of parental Du145 (TaxS) cells as a baseline, we identified 453 up-regulated and 473 downregulated genes in the Du145-R cells, and 483 up-and 365 down-regulated genes in the Du145-RB cells (Additional file 4). In addition, we found 216 genes with significantly different expression levels between DU145-RB and DU145-R. These 216 genes were presumably not related to the development of docetaxel resistance. By matching the three gene lists we further identified 615 (329 up-regulated and 286 down-regulated) genes that were shared by both DU145-R and DU145-RB as compared with DU145 (Additional file 4). These genes were thought to have stable expression changes after acquiring resistance to docetaxel. Of the 615 genes, the 40 most upand down-regulated in the TaxR cell lines were chosen for verification by RT-PCR and 37/40 (92.5 %) were confirmed (Additional file 5). Information about the most-differentially-expressed genes is shown in Additional file 6. The second most up-regulated gene was ABCB1, which encodes an ATPdependent drug efflux pump that mediates the development of resistance to anticancer drugs (Gottesman et al. 2002). The average fold changes (log) in TaxR cell lines were up to 8.9 and 10.2 in up-regulated genes and down-regulated genes, respectively. The largest functional group was transcription factors (Additional file 7). Twenty-one oncogenes and 16 translocated cancer genes were also among the enriched functional groups in the set of 615 stably differentially-expressed genes.
The 615 most significantly deregulated genes were put into the Panther Online tool (www.pantherdb.org), which yielded 528 functional hits distributed on 11 GOterms, where the two largest groups were Binding (GO: 0005488) and Catalytic Activity (GO: 0003824) (Fig. 4b).
Thomson Reuters was employed to analyze enriched networks of expression changing genes and showed that the NF-κb, EGR1 (Early Growth Response 1) and ETS (ETS family of transcription factors) were the three most enriched networks in the docetaxel resistant cells (Fig. 4c and Additional file 8). PLAU and PLAUR (Plasminogen Activator, Urokinase Receptor), a ligand-membrane receptor pair, are the only 'Convergence hubs' and MDR1 (ABCB1) was connected to all three pathways.
Next generation sequencing data of PC3 and LNCaP, two docetaxel-sensitive cell lines similar to Du145, were added into further analysis. Multivariate modeling with SIMCA resulted in a model, which separated all TaxS and TaxR cell lines into two classes and extracted those genes that contributed most to the model (Table 3).
When comparing the list of 615 stably up-or downregulated gene lists with the fusion gene list, we found 6 genes in common (Table 4), all of which were up-regulated. Table 5 summarizes the prostate cancer cell lines used in this study. LNCaP, PC3 and DU145 cell lines were originally ordered from the ATCC (American Type Culture Collection). Du145 was cultured in medium containing docetaxel (from low concentration to high concentration, increased gradually) for one year, until Du145 acquired docetaxel resistance (Du145-R). We also cultured Du145-R in normal medium without docetaxel for one month (Du145-RB) to see if it would revert to docetaxel sensitive again (Kharaziha et al. 2015). DU145-RB was frozen Mutation and amino acid changing were showed in column 3 and 4. Mutation Type was analyzed by program SIFT and gained by setting cutoff as 0.05. 'Damaging' means that the substitution is predicted to affect protein function. NA not analyzed. The last column shows the connection between cell lines and tumor samples. It is labeled "yes" if mutation can be found in both cell line sequence and tumor sample sequence

Analysis of differentially expressed genes
We analyzed RNAseq data according to a published TopHat and Cufflinks protocol (Trapnell et al. 2012). In summary, we used TopHat to align reads to the reference genome, Cufflinks to assemble and obtain expression values for all transcripts, Cuffdiff for testing differential expression of genes and transcripts and finally the CummeRbund R package for downstream analysis and visualization.

Fusion detection method
We used ChimeraScan, which aligns paired-end reads to a reference genome-transcriptome with Bowtie in an iterative process where read pairs that could not be aligned were trimmed into smaller fragments and realigned (Iyer et al. 2011). ChimeraScan uses a filter to avoid false-positive chimeras.

Statistical analysis
The online services Panther (http://www.pantherdb. org) and Thomson Reuters were applied for functional enrichment analysis (Mi et al. 2013  and 2 classes (TaxS and TaxR) were set in the model to obtain VIP scores by whichvariables (genes) are sorted based on importance (contribution to the model) of genes (Bylesjo et al. 2006).

PCR and qPCR validation
Total RNA was isolated from cell lines by TRIzol (Invitrogen, Catalog #15596018) according to the manufacturer's instructions. Cloned AMV First-Strand Synthesis Kit (Life Technologies, Catalog #12328) was used to transcribe mRNA to cDNA.
PCR primers for fusion validation were designed according to the sequence of fusion transcripts. Forward primer was located on the 5′ gene of the fusion gene and reverse primer on the 3′ gene of the fusion gene. Primers for validation of mutations covered mutation points. PCR was conducted using Platinum Taq DNA polymerase (Life Technologies, Catalog #10966018) and was followed by Sanger sequencing (conducted by Eurofins Genomics). To differentiate gene expression levels of selected genes, 20-32 amplifying cycles were used based on gene expression level.

Plasmid construction and western blot validation
PCR product was ligated into multiple cloning sites of pCMV-AC-GFP after digestion of restriction enzymes, Sgf I and Mlu I. pCMV-AC-GFP was purchased from ORIGENE (Catalog #PS100010), Sgf I enzyme from NEB (Catalog #R0630S), Mlu I enzyme from NEB (Catalog #R0198S), and T4 ligase from Promega (Catalog #M180A). We transfected HEK293 cells with constructed plasmid 48 h before collecting cells in lysis buffer. Western blot experiments were conducted using these cell lysates. Anti-TurboGFP antibody was purchased from Evrogen (Catalog #AB513).

Discussion
We have identified 42 genes with specific and stable mutations in TaxR cells. The functions of these genes may support their importance in the development of docetaxel resistance. Among these genes, SMAD4 is a co-activator and mediator of signal transduction by TGF-beta and acts as a tumor suppressor. Experiments have shown that SMAD4 inactivation promotes drug  (Zhang et al. 2014;Raz et al. 2014). ABCA2 is a member of ATP-binding cassette (ABC) transporters that transports many kinds of small molecules through membranes and is involved in drug resistance in leukemia cell lines (Dharmapuri et al. 2015).
Approximately 50 % of prostate cancer has primary resistance to docetaxel treatment. The other half is sensitive to docetaxel but eventually develops secondary (acquired) resistance (Marin-Aguilera et al. 2012). In this study, 34 out of the 42 mutations discovered in the resistant cell lines can be found in tumor samples from patients (Table 1), implicating that primary and acquired resistance may share the same molecular mechanism(s). In the case of primary resistance, most cancer cells carry the resistant genomic changes before the treatment, whereas for acquired resistance, just a few cancer cells carry these resistant genomic changes before treatment. By treatment selection or new mutational events, most cancer cells become carriers of resistant genomic changes. This hypothesis can be further tested in studies using tumor samples from patient cohorts with data of docetaxel treatment.
The four fusion transcripts (listed in Fig. 3c: MYH9-EIF3D, LDLR-RPL31P11, TAF15-AP2B1, VCL-ADK) could be detected by PCR and qPCR in the cell lines, but their translation into protein could not be validated by western blot in Du145, Du145-R and Du145-RB, probably due to the low expression of the fusion proteins. Fusion transcripts could be translated into protein in stably transfected HEK293 cells analyzed by western blot. Moreover, several genes involved in the fusion events have shown important functions in cancer development. TAF15, a member of the FET family, has been found rearranged with various transcription factors with cancer promoting functions in sarcomas as well as in rare hematopoietic and epithelial cancers (Kovar 2011). MYH9 is a member of the myosin superfamily and its function is related to migration, invasion and metastasis of cancer cells. EIF3D is associated with cell cycle regulation and motility of prostate cancer cells (Gao et al. 2015). MYH9 fusion proteins have been found in anaplastic large cell lymphoma and one example is the MYH9-ALK fusion protein that has tyrosine kinase activity in vivo (Lamant et al. 2003). The MYH9-USP6 detected by a previous study and MYH9-EIF3D found in the present study have the same fusion point in MYH9. MYH9, which is located in the 5′ part of the   When we compared the expression of DU145-RB and DU145-R, we found 216 genes that were differently expressed. We tested that DU145-RB was still docetaxel resistant, indicating these genes were not involved in maintaining docetaxel resistance of the two resistant cell lines. As expected, ABCB1 (MDR1) was confirmed as one of the top 10 differentially expressed genes that could separate TaxR from TaxS cells. Its functional importance was further supported by its connection with the NF-κb, EGR1 and ETS pathways (Fig. 4). ABCB1, which shows overexpression in some cancers, is involved in a common resistance mechanism. However, limited studies showed significant connection between ABCB1 and clinical outcomes, such as survival (Shaffer et al. 2012), indicating the importance of other molecular and biological changes. Researchers and pharmaceutical companies are trying to circumvent this strategy and find new potential genes or pathways to overcome resistance in cancer.
TGPSM2 and GRK5 are members of G-protein signaling pathway important in cancer progression. SKAP1 encodes a src kinase associated phosphoprotein 1 and is a member of the Ras signaling pathway and B cell receptor signaling pathway. LIMK1 is a serine/threonine kinase associated with the cytoskeletal structure in many cellular processes, and may have importance in the sensitivity of lung cancer and osteosarcoma cells to chemotherapy treatment (Chen et al. 2013;Zhang et al. 2011). The analysis further showed that PLAU and PLAUR (Plasminogen Activator, Urokinase Receptor), a pair of ligand and membrane receptor, constituted the only 'Convergence hub' by statistical analysis using the Thomson Reuters software. This novel finding may suggest that they may play a unique role in docetaxel resistance. It would be interesting to further study if they alone or, together with other important genomic findings in this study, can be further verified as important biomarkers to predict primary docetaxel resistance. Most importantly, they can even become attractive targets for the development of new drugs to overcome both primary and acquired docetaxel resistance.

Conclusion
The present study found both previous and novel mutations, genes with altered expression levels, and fusion proteins in docetaxel resistant prostate cancer cell lines, and provide some understanding of acquired docetaxel resistance at the gene transcription level. If some of these changes can be further verified with importance in primary resistance, they can be considered as predictive biomarkers for docetaxel treatment as well as targets for the development of new treatments to overcome the docetaxel resistance.