METTL3 regulates breast cancer-associated alternative splicing switches

Abstract

Alternative splicing (AS) enables differential inclusion of exons from a given transcript, thereby contributing to the transcriptome and proteome diversity. Aberrant AS patterns play major roles in the development of different pathologies, including breast cancer. N⁶-methyladenosine (m⁶A), the most abundant internal modification of eukaryotic mRNA, influences tumor progression and metastasis of breast cancer, and it has been recently linked to AS regulation. Here, we identify a specific AS signature associated with breast tumorigenesis in vitro. We characterize for the first time the role of METTL3 in modulating breast cancer-associated AS programs, expanding the role of the m⁶A-methyltransferase in tumorigenesis. Specifically, we find that both m⁶A deposition in splice site boundaries and in splicing and transcription factor transcripts, such as MYC, direct AS switches of specific breast cancer-associated transcripts. Finally, we show that five of the AS events validated in vitro are associated with a poor overall survival rate for patients with breast cancer, suggesting the use of these AS events as a novel potential prognostic biomarker.

Fig. 1. Identification of AS events in non-tumorigenic and breast cancer cell lines.

A Venn diagrams showing the number of common alternative splicing events and genes in the non-tumorigenic mammary epithelial cell line MCF10-A and the breast cancer cell lines MCF7 and MDA-MB-231. Violin plots of changes of the significant percent splicing inclusion (∆PSI) in the breast cancer cell lines (B) MCF7 and (C) MDA-MB-231 related to the normal mammary epithelial cell line MCF10-A. D Dot plots representing the Gene Ontology enrichment (GO) analysis of the common spliced genes in MCF7 and MDA-MB-231. BP biological process, KEGG KEGG pathways. The size and the color of the dots are proportional to the number of genes enriched in each GO term and the significance of the enrichment (1.10⁻⁷⁵< P-value < 1.10⁻¹), respectively. E RT-PCR showing the different splicing events between the non-tumorigenic MCF10-A and breast cancer MCF7 and MDA-MB-231 cell lines. The number of the skipped exons are depicted for each transcript. The PSI was calculated in percentage for each gene. Non-specific bands are indicated with an asterisk.

Fig. 2. METTL3 promotes cell growth in breast cancer cell lines.

A RT-qPCR analysis of METTL3 mRNA level (upper) in the non-tumorigenic mammary epithelial cell line MCF10-A and the breast cancer cell lines MCF7 and MDA-MB-231. METTL3 is normalized to β-ACTIN. Western blot of METTL3 (lower) on whole cell extracts (WCE) from MCF10-A, MCF7 and MDA-MB-231 cell lines. β-ACTIN (ACTB) is used as the loading control. B LC-MS/MS quantification of m⁶A, m¹A and m⁶A_m in mRNA of MCF10-A, MCF7 and MDA-MB-231 cell lines. Methylated adenosines are normalized to the total of unmodified adenosines. RT-qPCR analysis and western blot of METTL3 mRNA and protein levels, respectively in METTL3 knockdown (sh1 and sh2) and scramble (scr) control in MCF10-A (C), MCF7 (D) and MDA-MB-231 (E) cell lines. METTL3 is normalized to β-ACTIN. Cell proliferation rate of scramble (scr) and METTL3 knockdowns (sh1 and sh2) in MCF10-A (F), MCF7 (G) and MDA-MB-231 (H) assessed over 4 days. Colony formation assay in MCF10-A (I), MCF7 (J) and MDA-MB-231 (K) cells in scramble (scr) and upon depletion of METTL3 (sh1 and sh2) at 7 days after seeding. Quantification of the relative number of colonies was calculated with scr set as 1. Percentage of apoptotic cells in control and METTL3 depleted cells in MCF10-A (L), MCF7 (M), and MDA-MB-231 (N) cell lines. + indicates Annexin V positive cells; and – indicates Annexin V negative cells. Volcano plots representing the Log₂ fold change of differentially expressed genes upon METTL3 knockdown (sh1 and sh2) in MCF10-A (O), MCF7 (P) and MDA-MB-231 (Q) cell lines in comparison to control cells (scr). The significant up- and down-regulated genes are shown in red and blue, respectively. P-value < 0.05. Data are mean ± SEM; n = 3; ****p < 0.0001; ***p < 0.001; **p < 0.01; *p < 0.05.

Fig. 3. METTL3 modulates AS in breast cancer cell lines.

A Venn diagrams showing the common alternative splicing events between knockdowns of METTL3 in MCF10-A (blue), MCF7 (yellow) and MDA-MB-231 (red). The corresponding total number of genes is indicated in brackets. P-value < 0.05, statistically significant. B, C Violin plots of the significant percent splicing inclusion (∆PSI) in knockdowns of METTL3 in MCF7 and MDA-MB-231 cells related to control cells. Dot plots of the GO enrichment analysis of the differentially expressed (DEG) and spliced (DSG) genes in (D) MCF10-A, (E) MCF7 and (F) MDA-MB-231 upon depletion of METTL3. The size and the color of the dots are proportional to the number of genes enriched in each GO term and to the significance of the enrichment (1.10⁻⁹ < P-value < 1.10⁻¹), respectively.

Fig. 4. METTL3 influences AS via m⁶A deposition.

A Dot plot representing the level of statistical significance of AS transcripts harboring m⁶A (red) or non-m⁶A modified AS transcripts (blue), in chromatin-bound transcripts dataset from HEK293T, MCF7, MDA-MB-231, and transcripts with premature termination codons (PTC) dataset from human glioblastoma [48]. The size of the dots is proportional to the frequency of the events. P-value < 0.01, statistically significant. Motif density of m⁶A peaks in the –80 to +80 nt region around the m⁶A peak in intronic or random regions (upper panels) and the corresponding HOMER motifs outputs (lower panels) in (B) MCF7 and in (C) MDA-MB-231. D RT-PCR of AS genes in METTL3 knockdown in breast cancer cell lines MCF7 and in MDA-MB-231 cells. The number of the skipped exons are depicted for each transcript. The PSI was calculated in percentage for each gene. Non-specific bands are indicated with an asterisk.

Fig. 5. m⁶A motifs in MYC 3´UTR promotes the translation of MYC mRNA.

A m⁶A peak distribution in MYC mRNA in MCF7 (left panel) and MDA-MB-231 (right panel) visualized in IGV. Input reads are represented in darker colors and the enriched RNA immunoprecipitated in yellow (MCF7) or red (MDA-MB-231). The amplified region by qPCR is depicted with a red line below MYC gene body. B RT-qPCR of m⁶A RNA immunoprecipitation (MeRIP) showing the enrichment of m⁶A in MYC relative to GAPDH in MCF7 (left) and MDA-MB-231 (right). C Relative level of SELECT products specific to m⁶A site in MYC 3´UTR, using total RNA from DMSO treated or STM2457 treated MDA-MB-231 cells. D RT-qPCR analysis of MYC after FLAG-METTL3 immunoprecipitation performed in control cells (+Dox) or in cells overexpressing Tet-off FLAG-METTL3 (-Dox) in MDA-MB-231 cells. E RT-qPCR analysis of MYC mRNA (upper panel) and western blot for MYC (lower panel) in MDA-MB-231 upon STM2457 treatment. βACTIN is used as loading control. F RT-qPCR analysis of MYC mRNA after treatment with actinomycin D at the time points 0, 10, 30 and 60 min in MDA-MB-231 control and treated with STM2457. G Relative Renilla luciferase activity of the psiCHECK2-MYC 3´UTR in MDA-MB-231 cells treated with DMSO (control) or with STM2457 for 48 h. Control cells were transfected with psiCHECK2 empty vector. Renilla luciferase activity was measured and normalized to Firefly luciferase. Data are mean ± SEM; n = 3 or 4; ****p < 0.0001; ***p < 0.001; *p < 0.05. In A, D, E, and G P-values were determined by two-tailed t-test; in C P-values were determined by one-tailed t-test. H Western blot showing the overexpression of SRSF11 in MDA-MB-231 in comparison to MCF10-A and MCF7. HDAC1 is used as loading control. I Western blot assessing the expression of SRSF11 in MDA-MB-231 upon STM2457 treatment. βACTIN is used as loading control. J Overlaps between AS events of knockdown of METTL3 in MCF7 and MYC-associated AS events (left panel); P-value < 0.0001. GO analysis of the common genes between AS events between knockdown of METTL3 in MCF7 and MYC-associated AS events (right panel); P-value < 0.05. K Overlaps between AS events of knockdown of METTL3 in MDA-MB-231 and MYC-associated AS events (left panel); P-value < 0.0001. GO analysis of the common genes between AS events in knockdown of METTL3 in MDA-MB-231 and MYC-associated AS events (right panel); P-value < 0.05.

Fig. 6. Identification of breast cancer prognosis-related AS events.

PSI values were analyzed in breast cancer patients (1094 samples) and normal samples (113 samples) for the AS events tested in (A) COMMD4_AS2, (B) GNAS, (C) MATR3, (D) RHOC, (E) COMMD4_AS1, (F) MARK3, (G) POLDIP3, (H) FASTK, (I) EXOC7, (J) BAX. Data were taken from the TCGA SpliceSeq database. Kaplan–Meier plots of overall survival (OS) for breast cancer patients classified according to the AS events expression (low or high) for (K) COMMD4_AS2, (L) EXOC7, (M) RHOC (N) BAX, (O) FASTK, (P) GNAS, (Q) MARK3, (R) MATR3, (S) COMMD4_AS1, (T) POLDIP3. U OS rate for the combination of COMMD4_AS1, COMMD4_AS2, EXOC7, RHOC and PODLIP3. p < 0.05, statistically significant.