Androgen-dependent alternative mRNA isoform expression in prostate cancer cells

Cell culture

Cell culture was as described previously^25,36. All cells were grown at 37°C in 5% CO₂. LNCaP cells (CRL-1740, ATCC) were maintained in RPMI-1640 with L-Glutamine (PAA Laboratories, R15-802) supplemented with 10% Fetal Bovine Serum (FBS) (PAA Laboratories, A15-101). For androgen treatment of cells, medium was supplemented with 10% dextran charcoal stripped FBS (PAA Laboratories, A15-119) to produce a steroid-deplete medium. Following culture for 72 hours, 10 nM synthetic androgen analogue methyltrienolone (R1881) (Perkin-Elmer, NLP005005MG) was either added (Androgen +) or absent (Steroid deplete) for the times indicated.

RNA-Seq analysis

RNA-seq transcript expression analysis of previously generated data²⁵ was performed according to the Tuxedo protocol³⁷. All reads were first mapped to human transcriptome/genome (build hg19) with TopHat³⁸/Bowtie³⁹, followed by per-sample transcript assembly with Cufflinks⁴⁰. The mapped data was processed with Cuffmerge, Cuffdiff and Cuffcompare, followed by extraction of significantly differentially expressed genes/isoforms; expression changes between cells grown with androgen and cells grown without androgens were assessed. Reference files for the human genome (UCSC build hg19) were downloaded from the Cufflinks pages: (UCSC-hg19 package from June 2012 was used.). The software versions used for the analysis were: TopHat v1.4.1, SAM tools Version: 0.1.18 (r982:295), bowtie version 0.12.8 (64-bit) and cufflinks v1.3.0 (linked against Boost version 104000). The Tuxedo protocol³⁷ was carried out as follows: For steps 1–5, no parameters (except for paths to input/output files) were altered. In step 15, additional switches -s, -R, and -C were used when running cuffcompare. Steps 16–18 (extraction of significant results) were performed on the command line.

RNA extraction, RT–PCR and real-time PCR

Cells were harvested and total RNA extracted using TRIzol (Invitrogen, 15596-026) according to manufacturer's instructions. RNA was treated with DNase 1 (Ambion, AM2222) and cDNA was generated by reverse transcription of 500ng of total RNA using the Superscript VILO cDNA synthesis kit (Invitrogen, 11754-050). Alternative events were analysed by either reverse transcriptase PCR or real-time PCR. Exon profiles were monitored and quantified using the Qiaxcel capillary electrophoresis system (Qiagen) and percentage inclusion was calculated as described previously¹⁰. Real time PCR was performed in triplicate on cDNA using SYBR® Green PCR Master Mix (Invitrogen, 4309155) and the QuantStudio 7 Flex Real-Time PCR System (Thermo Fisher Scientific). Samples were normalised using the average of three reference genes, GAPDH, β -tubulin and actin. Ct values for each sample were calculated using SDS 2.4 software (Applied Biosystems) and relative mRNA expression was calculated using the 2-ΔΔCt method. All primer sequences are listed in Supplementary Table 1. Raw Ct values are given in Dataset 1⁴¹.

Antibodies

The following commercial antibodies were used in the study: anti-RLN2 rabbit monoclonal (Abcam, ab183505 1:1000 dilution), anti-TACC2 rabbit polyclonal antibody (11407-1-AP, Proteintech 1:500 dilution), anti-NDUFV3 rabbit polyclonal antibody (13430-1-AP, Proteintech 1:500 dilution), anti-actin rabbit polyclonal (A2668, Sigma 1:2000 dilution), anti-α-Tubulin mouse monoclonal (Sigma, T5168 1:2000 dilution), normal rabbit IgG (711-035-152, Jackson labs 1:2000 dilution) and normal mouse IgG (715-036-150, Jackson labs 1:2000 dilution).

Gene ontology analysis

Gene ontology (GO) analysis of RNA-Seq data was carried out as described previously⁴². Enrichment of GO terms (with b500 annotations) was calculated using the goseq R package (version 1.18.0). Genes were considered significant at a p-value threshold of 0.05 after adjustment using the Benjamini-Hochberg false discovery rate.

Bioinformatic analysis of patient transcriptome data

Available clinical and processed RNA-Seq data from The Cancer Genome Atlas (TCGA) prostate adenocarcinoma (PRAD) cohort, comprising 497 tumour samples from as many patients with different stages / Gleason grades and 52 matched samples taken from normal prostate tissue (were downloaded from the Broad Institute TCGA Genome Analysis Center (Firehose 16/01/28 run https://doi.org/10.7908/C11G0KM9⁴³). Transcriptome data from the TCGA PRAD cohort were analysed for alternative isoform expression, with transcript models relying on TCGA GAF2.1, corresponding to the University of California, Santa Cruz (UCSC) genome annotation from June 2011 (hg19 assembly). This annotation encompassed 42 of the 73 androgen-regulated alternative mRNA isoform pairs identified. These were studied using two types of analysis: 1) differential transcript expression between tumour and normal prostate tissue and 2) correlation between isoform expression in tumour samples and Gleason score or tumour stage.

Differential isoform and gene expression analysis was performed on estimated read counts using the limma software R package (version 3.7) following its RNA-Seq analysis workflow⁴⁴. This workflow was also used for differential isoform ratio analysis, relying on logit-transformed ratio (see below). An FDR-adjusted p-value of 0.05 for the moderated t-statistics was used as threshold for significance of differential expression. Individual isoform expression was estimated in TPM (transcripts per million mapped reads). The expression ratio, henceforth called PSI (percent spliced-in), of each annotated androgen-regulated isoform pair in each TCGA sample was calculated as the ratio between the expression of isoform 1 and the total expression of isoforms 1 and 2 combined, i.e. the sum of their expressions. For each isoform pair, ΔPSI is the difference of median PSI between the tumour and the normal groups of samples.

Two-tailed Spearman’s rank correlation tests were used to study the association between isoform expression and both Gleason score and tumour stage (these were used herein as numeric variables). An FDR-adjusted p-value of 0.05 was used as threshold for significance. Isoform expression differences between tumour and normal samples were considered equivalent to those detected in LNCaP cells under androgen stimulation when there was a statistically significant consistent change in the levels of the expected induced or repressed isoform (1 or 2), concomitant with no contradictory change in the PSI. Isoform “switches” were considered equivalent when there was a minimum (ΔPSI > 2.5%) and statistically significant consistent change in the PSI. Equivalent criteria were used to evaluate the equivalence between androgen-dependence and the associations with Gleason score and tumour stage.

Global identification of androgen-dependent mRNA isoform production in prostate cancer cells predicts a major role for alternative promoter utilisation

We analysed previously published RNAseq data from LNCaP cells²⁵ to globally profile how frequently androgens drive production of alternative mRNA isoforms in prostate cancer cells. This analysis identified a group of 73 androgen regulated alternative mRNA isoforms, which could be validated by visualisation on the UCSC Genome Browser⁴⁵ (Table 1). 64 AR regulated mRNA isoforms were novel to this study. Experimental validation in an independent RNA sample set using RT-PCR confirmed 17/17 of these alternative events at the mRNA level (Supplementary Figure 1). 73% of genes (53/73) with identified alternative androgen regulated mRNA isoforms also changed their overall expression levels in response to androgens (Table 2). Some of the androgen regulated alternative events are in genes are already implicated in in either prostate cancer or other cancer types (summarised in Table 3). However, Gene Ontology analysis of these 73 genes did not identify any significantly enriched biological processes.

The 73 identified mRNA isoforms were generated via androgen-regulated utilisation of 56 alternative promoters, 4 alternative 3′ ends and 13 alternative splicing events (Figure 1A). Of the 56 androgen regulated alternative promoters that were identified, 23 alternative promoters were induced by androgens (including LIG4, Figure 1B), 26 promoters were repressed by androgens, and for 7 genes there was a switch in usage from one promoter to another (Table 1). The alternative splicing events that were under androgen control included 12 alternative exons and one androgen-regulated intron retention (Table 1). 10 of these are novel to this study, including exclusion of an alternative exon in ZNF678 (Figure 1C). Of the alternative exons, six genes contained switches in previously unannotated protein-coding exons in response to androgen-exposure. We also identified four androgen regulated alternative mRNA 3' end isoform switches, including a switch in the 3’ end of the mRNA transcript for the MAT2A gene (Figure 1D).

Figure 1. Global identification of androgen-dependent mRNA isoform production in prostate cancer cells predicts a major role for alternative promoter utilisation.

(A) Analysis of RNAseq data from LNCaP cells grown with (A+) or without androgens (R1881) (steroid deplete, SD) for 24 hours identified 73 androgen regulated alternative mRNA isoforms. The 73 alternative events were generated via androgen-regulated utilisation of 56 alternative promoters, 4 alternative 3' ends and 13 alternative splicing events. (B) Androgens drive a promoter switch in the LIG4 gene, which produces an mRNA isoform with an alternative 5’UTR. Visualisation of our LNCaP cell RNA-seq reads for the LIG4 gene on the UCSC genome browser identified a switch from promoter 1 to alternative promoter 2 in cells grown in the presence of androgens. Promoter 2 is predicted to produce a different 5’UTR without influencing the protein sequence (left panel). Quantitative PCR using primers specific to each promoter indicate that in response to androgens there is repression of promoter 1 and induction of promoter 2 (right panel). (C) Androgens drive alternative splicing of the ZNF678 gene. Visualisation of our LNCaP cell RNA-seq reads for the ZNF678 gene on the UCSC genome browser identified a switch to inclusion of a cassette exon in the presence of androgens. Inclusion of the alternative cassette exon in the ZNF678 gene is predicted to induce a switch to an alternative non-coding mRNA isoform (left panel). Quantitative PCR using primers in flanking exons confirmed increased inclusion of the alternative exon in LNCaP cells exposed to androgens (right panel). (D) Androgens promote selection of an alternative 3’ end for the MAT2A gene. Visualisation of our LNCaP cell RNA-seq reads for the MAT2A gene on the UCSC genome browser indicates a switch to reduced usage of an alternative 3’ end in the presence of androgens (left panel). Quantitative PCR using primers specific to each isoform confirmed down-regulation of an alternative 3’ end (p<0.01). Alternative 3’ ends for the MAT2A gene are predicted to produce proteins with different amino acid sequences and to influence a known Pfam domain (right panel).

Androgen regulated events control the production of alternative protein isoforms, non-coding RNAs and alternative 5' UTRs

48/73 (66%) of the androgen regulated alternative events detected in response to androgen stimulation are predicted to change the amino acid sequence of the resulting protein (Table 1). Some of these are already known to have a well characterised role in prostate cancer progression, including an alternative promoter in the oncogene TPD52 that produces a protein isoform called PrLZ (Figure 2A)^46–49. Others are not so well characterised. Using western blotting we could detect a novel shorter protein isoform corresponding to androgen-driven selection of an alternative promoter in the TACC2 gene (Figure 2B); and exclusion of a cassette exon in the NDUFV3 gene, which we show also produces a novel shorter protein isoform (Figure 2C). We also detected a switch in the 3' end of the mRNA transcript for the MAT2A gene, which is predicted to produce a protein isoform with a shorter C-terminal domain (Figure 1D); and induction of an alternative 3' isoform of CNNM2, which is predicted to be missing a conserved CBS domain (Table 1 and Supplementary Figure 1).

Figure 2. Androgen regulated mRNA isoform switches control alternative protein isoforms and non-coding RNAs.

(A) Androgens induce an alternative promoter in the oncogene TPD52 that produces an isoform called PrLZ. Visualisation of our LNCaP cell RNA-seq reads for the TPD52 gene on the UCSC genome browser identified a switch from promoter 1 to alternative promoter 2 in cells grown in the presence of androgens. Promoter 2 is known to produce an alternative protein isoform of TPD52 known as PrLZ (left panel). Quantitative PCR using primers specific to each promoter indicate an induction of the PrLZ isoform in response to androgens (middle panel). PrLZ has an alternative N-terminal amino acid sequence which results in an alternative protein isoform and disrupts a known Pfam domain (right panel). (B) Androgens induce an alternative promoter in the TACC2 gene that produces a novel alternative protein isoform. Visualisation of our LNCaP cell RNA-seq reads for the TACC2 gene on the UCSC genome browser identified a switch from promoter 1 to alternative promoter 2 in cells grown in the presence of androgens. Promoter 2 is predicted to produce an alternative shorter protein isoform of TACC2 (isoform 2) (left panel). Quantitative PCR using primers specific to each promoter indicate a switch from isoform 1 to isoform 2 in response to androgens (middle panel). Detection of TACC2 protein in LNCaP by western blotting (cells were grown with or without androgens for 24 or 48 hours). Tubulin was used as a loading control. Exposure to androgens for 48 hours induces expression of the alternative TACC2 protein isoform (right panel). (C) Androgens drive alternative splicing of the NDUFV3 gene. Visualisation of our LNCaP cell RNA-seq reads for the NDUFV3 gene on the UCSC genome browser identified a switch to exclusion of a cassette exon in the presence of androgens (left panel). Quantitative PCR using primers in flanking exons confirmed less inclusion of the alternative exon in LNCaP cells exposed to androgens (middle panel). Exclusion of the alternative cassette exon is predicted to produce an alternative protein isoform. Detection of NDUFV3 protein in LNCaP cells using western blotting (right panel). (D) Androgens suppress an alternative promoter in the RLN2 gene, which produces a shorter non-coding mRNA isoform. Visualisation of our LNCaP cell RNA-seq reads for the RLN2 gene on the UCSC genome browser identified a switch from promoter 1 to alternative promoter 2 in cells grown in the presence of androgens. Promoter 2 is predicted to produce an untranslated non-coding mRNA isoform (left panel). Quantitative PCR using primers specific to each promoter indicated a significant switch in promoter usage in response to androgens (middle panel). Detection of RLN2 protein in LNCaP by western blotting (cells were grown with or without androgens for 48 hours). Actin was used as a loading control. As seen previously⁵⁵, androgens suppress RLN2 protein levels.

11 of the remaining identified androgen-regulated alternative events change the expression of mRNAs from coding to non-coding or untranslated (not predicted to produce a protein) (Table 1). These included promoter switches for the RLN1 and RLN2 genes which encode peptide hormones that may be important in prostate cancer^5,50–55. Androgens drive a promoter switch in both RLN1 and RLN2 to produce predicted non-coding or untranslated mRNA isoforms, reducing expression of protein-coding RLN1 and RLN2 mRNA isoforms. To test whether prostate cancer cells turn off gene expression by switching between utilisation of promoters that generate coding and noncoding mRNAs, we analysed RLN2 protein levels. Consistent with our hypothesis and a previous study⁵⁵, RLN2 protein production was negatively regulated by androgens in parallel to the switch to the non-coding mRNA isoform (Figure 2D).

14 of the identified androgen-dependent mRNA isoforms lead to/result in coding mRNAs with altered 5’ untranslated regions (5′ UTR) with no impact on the coding sequence. These include a promoter switch in the LIG4 gene (Figure 1B).

Differential expression of androgen-dependent mRNA isoforms in prostate adenocarcinoma versus normal tissue

To investigate potential links between androgen-dependent mRNA isoforms and tumourigenesis, we analysed the expression of 41 androgen-regulated mRNA isoform pairs in clinical prostate adenocarcinoma and normal prostate tissues. This analysis utilised transcriptomic data from 497 tumour samples and 52 normal samples in the PRAD TCGA cohort¹⁰⁴. The remaining isoform pairs identified within our dataset have not been previously annotated by UCSC, therefore it was not possible to include them in our comparison. A description of the cohort used is summarised in Table 4.

33 of the 42 mRNA isoform pairs exhibited significant differences in the expression of at least one of the isoforms, or in the isoform expression ratio between tumour and normal tissues (Table 5). 13 of those tumour-specific alterations mimicked the effect of androgen stimulation in LNCaP cells: the changes were in form of alternative promoters for TACC2, TPD52, NUP93, PIK3R1, RDH13, ZFAND6, CDIP1, YIF1B, LIMK2, and FDFT1; an alternative 3´ end in CNNM2; and alternative exons in NDUFV3 and SS18 (Figure 3, Table 5 & Supplementary Figure 2). Two of the alternative promoters (ZFAND6 and CDIP1) are predicted to introduce a change in the 5′UTR, whereas all the others are predicted to alter the resulting protein isoform. A number of mRNA isoforms that were androgen responsive in LNCaP cells showed tumour specific alterations opposite to the effect of androgen stimulation. These were LIG4, MAPRE2, OSBPL1A, SEPT5, NR4A1, and RCAN1 (all predicted to alter the resulting protein isoform except LIG4). For the remaining 14 mRNA isoform pairs, the data was inconclusive according to the consistency conditions listed in the methods section (Table 5).

Figure 3. Differential expression of androgen dependent mRNA isoforms in prostate cancer versus normal tissue within the PRAD TCGA cohort for TPD52, TACC2, NDUFV3 and CNNM2.

Violin-boxplots of expression in transcripts per million mapped reads (TPM) of Isoforms 1 (left panel) and 2 (central panel), and of their expression ratio in PSI (right panel) in normal and tumour samples. The mean log2 fold-change (logFC) in expression between tumour and normal samples and the associated FDR-adjusted p-value for the moderated t-statistic of differential expression are shown for both isoforms (left and central panels). The mean difference in PSI (deltaPSI) between tumour and normal samples and the associated FDR-adjusted p-value for the Mann-Whitney U test of differential splicing are shown (right panel).

Changes in androgen-dependent mRNA isoform expression during tumour progression

We next investigated whether the identified androgen-dependent mRNA isoforms are differentially expressed during prostate cancer progression by correlating the expression levels of each isoform with Gleason scores and prostate tumour grades within the PRAD TCGA cohort (Figure 4 & Figure 5, Table 6 & Table 7 and Supplementary Figure 3 & Supplementary Figure 4). For 6 of the alternative mRNA isoforms responsive to androgens (made from alternative promoters in LIG4, OSBPL1A, CLK3, TSC22D3 & ZNF32 and utilising an alternative exon in ZNF121), the expression changed significantly with Gleason score and showed specific alterations consistent with the effect of androgen stimulation. Conversely, 9 alternative isoforms (which were androgen responsive in LNCaP cells) showed tumour specific alterations opposite to the effect of androgen stimulation (including an alternative promoters in NUP93 and the alternative 3´end of MAT2A). 3 androgen regulated mRNA isoforms (OSBPL1A, CLK3 and TSC22D3) change significantly with both Gleason grade and tumour stage.

Figure 4. Differential alternative mRNA isoform expression in the TGCA PRAD cohort across different Gleason grades for OSBPL1A, CLK3, TSC22D and ZNF121.

Violin-boxplots of expression in transcripts per million mapped reads (TPM) of Isoforms 1 (left panel) and 2 (central panel), and of their expression ratio (right panel) by Gleason grade. Their respective Spearman’s correlation coefficient (Rho) with grade and associated FDR-adjusted p-value are shown.

Figure 5. Differential alternative mRNA isoform expression in the TGCA PRAD cohort across different tumour stages for OSBPL1A, CLK3 and TSC22D3.

Violin-boxplots of expression in transcripts per million mapped reads (TPM) of Isoforms 1 (left panel) and 2 (central panel), and of their expression ratio (right panel) by tumour stage. Their respective Spearman’s correlation coefficient (Rho) with stage and associated FDR-adjusted p-value are shown.