An EMT-driven alternative splicing program occurs in human breast cancer and modulates cellular phenotype

Abstract

Epithelial-mesenchymal transition (EMT), a mechanism important for embryonic development, plays a critical role during malignant transformation. While much is known about transcriptional regulation of EMT, alternative splicing of several genes has also been correlated with EMT progression, but the extent of splicing changes and their contributions to the morphological conversion accompanying EMT have not been investigated comprehensively. Using an established cell culture model and RNA-Seq analyses, we determined an alternative splicing signature for EMT. Genes encoding key drivers of EMT-dependent changes in cell phenotype, such as actin cytoskeleton remodeling, regulation of cell-cell junction formation, and regulation of cell migration, were enriched among EMT-associated alternatively splicing events. Our analysis suggested that most EMT-associated alternative splicing events are regulated by one or more members of the RBFOX, MBNL, CELF, hnRNP, or ESRP classes of splicing factors. The EMT alternative splicing signature was confirmed in human breast cancer cell lines, which could be classified into basal and luminal subtypes based exclusively on their EMT-associated splicing pattern. Expression of EMT-associated alternative mRNA transcripts was also observed in primary breast cancer samples, indicating that EMT-dependent splicing changes occur commonly in human tumors. The functional significance of EMT-associated alternative splicing was tested by expression of the epithelial-specific splicing factor ESRP1 or by depletion of RBFOX2 in mesenchymal cells, both of which elicited significant changes in cell morphology and motility towards an epithelial phenotype, suggesting that splicing regulation alone can drive critical aspects of EMT-associated phenotypic changes. The molecular description obtained here may aid in the development of new diagnostic and prognostic markers for analysis of breast cancer progression.

Figure 1. Alternative mRNA isoform expression in EMT.

(A) Schematics of the in vitro EMT induction experiment. Immortalized human mammary epithelial cells (HMLE) expressing Twist fused to Estrogen Receptor (ER) were induced to undergo EMT by addition of tamoxifen into the culture media. mRNA was collected before EMT induction (epithelial sample) and after EMT induction (mesenchymal sample). cDNA pools from both samples were deep sequenced (RNA-Seq) and analyzed (See Materials and Methods). (B) Western blot analysis of N-cadherin, E-cadherin, fibronectin and vimentin expression with antibodies as indicated in cell lysates that were obtained before (1- untreated) and after (2- tamoxifen-treated) induction of EMT in HMLE/Twist-ER cells. α- tubulin was used as a loading control. (C) Gene ontology enrichment analysis bar graph of changes in alternative splicing events with |ΔΨ|> = 10% between samples. Gene ontology ‘biological process’, GO_BP_FAT, annotation is indicated in red on the y axis. KEGG Pathway (http://www.genome.jp/kegg/) annotation is indicated in blue on y axis. Benjamini FDR (−log10) is indicated on the x axis. Vertical dotted line marks Benjamini FDR = 0.05. (D) Column 1 shows different kinds of splicing events that have been analyzed. Columns 2–5 show the number of events of each type: (2) all known events based on AceviewAceView annotation; (3) events with both isoforms supported by RNA-Seq reads; (4) events detected at a False Discovery Rate (FDR) of 5% with ΔΨ > = 10% between samples; (5) events detected at an FDR of 5% with ΔΨ> = 30% between samples.

Figure 2. Motif analysis reveals splicing factors that are involved in the regulation of EMT–specific splicing.

(A) Pentamer motifs significantly enriched (FDR<0.1) in the 4 flanking 250-nt intronic regions of EMT-regulated skipped exons. Statistics of motifs resembling known binding sites of splicing factors are annotated as described in the key. Motifs that are not recognized as known binding sites are grouped into the “Other” group. * = at least one known motif of that splicing factor has an FDR<0.05. (B) Scatter plot of expression levels of RNA binding proteins and mRNA splicing regulators in epithelial and mesenchymal cells. Some splicing factors whose motifs were enriched in (A) are highlighted. Asterisks mark splicing factors which are also regulated by alternative splicing of their mRNA transcripts. Genes encoding components of cleavage/polyadenylation machinery are also highlighted. (C) A venn diagram showing potential regulation of EMT-associated skipped exon events by ESRP1, PTP and RBFOX2 splicing factors based on the microarray analysis of ESRP1 depleted MDA-MB-231 cells and CLIP-Seq analysis of FOX and PTB , (See Text S1). The universe of the Venn diagram consists of all EMT-regulated SE events by FDR of 5% and |ΔΨ|> = 10%. P (FOX2) = 8.58e-05; p (PTB) = 0.0013; p (ESRP1) = 9.27e-16.

Figure 3. EMT–associated alternative splicing events are confirmed in breast cancer cell lines.

(A) Alternative exon inclusion in four mRNA transcripts, as indicated, in eight breast cancer cell lines determined by a qRT-PCR analysis and depicted as a fold change relative to exon inclusion in T47D luminal cell line. (B) Alternative exon inclusion in five mRNA transcripts, as indicated, in eight breast cancer cell lines depicted as a fold change relative to exon inclusion in BT549 basal B cell line. (C) Distribution of all epithelial inclusion events combined. Each event is depicted as a fold change relative to inclusion in T47D. (D) Distribution of all mesenchymal inclusion events combined. Each event is depicted as a fold change relative to inclusion in BT549 cells. For (C) and (D), *** = p<0.001.

Figure 4. EMT–associated alternative mRNA isoforms classify breast cancer subtypes.

(A) Unsupervised hierarchical clustering of NCI-60 breast cancer cell lines , as indicated, using the 307 EMT-associated SE events (|ΔΨ|>0.1, FDR<0.05). Luminal cell lines are marked in blue. Basal B cell lines are marked in red. Shades of yellow indicate positive pearson correlation. Shades of blue indicate negative pearson correlation. The color scale is shown on the right. The dotted line separates luminal and basal cell line clusters. (B) Clustering correlation matrix of NCI-60 breast cancer cell lines , as indicated, using the 24 coherent events.

Figure 5. Alternative mRNA isoforms are expressed in FNA samples from breast cancer patients.

(A) An example of a fine needle aspiration (FNA) spread from a benign and an invasive human breast tumor. (B) A table describing gene names, gene functions, change in inclusion levels during EMT (ΔΨ) and proposed functions of six SE events used in the FNA qRT-PCR analysis in (C). (C) Heat plot of pairwise Pearson correlation coefficients obtained from correlation analysis of exon inclusion ratios in 40 IDC samples normalized to an average fibroadenoma sample exon inclusion for six alternative splicing events. Yellow indicates a correlation of 1, black indicates a correlation of 0, light blue indicates a correlation of −1. Shades of yellow and blue mark correlation in-between.

Figure 6. Expression of ESRP1 confers epithelial migration properties to mesenchymal cells.

(A) Western blot analysis of cell lysates from HMLE/pBP, HMLE/pBP-Twist and HMLE/pBP-Twist/ESRP1 cells probed with antibodies as indicated. α- tubulin was used as a loading control. (B) Still images from a live cell-tracking experiment of cells, as indicated, labeled with a cellular dye CMFDA and plated in the monolayer mixed 1∶20 with unlabelled cells. Cells were tracked for 12 hours. Cell tracks were generated using semi-automated cell tracking and represent single cell tracks over 12 hours with 10 minutes intervals. Centroids of fluorescent cells are indicated by grey circles. Videos are provided as Videos S1, S2, S3. Windrose plots of the range of motion of individual cells of each cell type are shown next to still images from a live-cell imaging experiment, as indicated. Windrose plots were generated by placing starting points of all cell tracks obtained in the cell tracking experiment into the same spot. (C) The box plot depicts speed distribution of individual cells inferred from live-cell imaging of cells in (A) and analyzed by the Imaris software. Edges of the boxes indicate 25^th and 75^th percentile and the whiskers 5^th and 95^th percentile. The line in the box indicates the median of the distribution. n = 138 cells for HMLE/pBP; n = 125 cells for HMLE/pBP-Twist; n = 113 cells for HMLE/pBP-Twist/ESRP1-EGFP. *** = p<0.001.

Figure 7. Expression of ESRP1 changes actin organization and localization of junctional markers in mesenchymal cells towards epithelial morphology.

(A) Immunofluorescence of cells, as indicated, using anti-ZO-1 antibody and Alexa350-phalloidin. Scale bar, 20 µm. Insets were 5× magnified. (B) Immunofluorescence of cells, as indicated, using anti-p120catenin antibody and Alexa405-phalloidin. Red arrows mark peripheral actin. Yellow arrows mark stress fibers. Blue arrows mark p120catenin at cell junctions. Scale bar, 5 µm. (C) Bar graph depicting movement of Texas Red–dextran across confluent monolayers of HMLE cells, as indicated, at 2 hrs and 4 hrs after addition of dextran compared to control cells expressing pBP (*, P<0.05; n = 6). Error bars represent SD.