Hypoxia-Induced Suppression of Alternative Splicing of MBD2 Promotes Breast Cancer Metastasis via Activation of FZD1

The majority of cancer-related deaths can be attributed to the invasion–metastasis cascade of cancer cells migrated from the primary tumor to distant organs (1, 2). One hallmark of solid tumors is hypoxia and hypoxia-inducible factors (HIF), as transcription factors, were reported to reshape extracellular matrix and promote the epithelial-to-mesenchymal transition (EMT), local invasion, and the metastasis of breast cancer cells to the lungs by stimulating the expression of lysyl hydroxylases, lysyl oxidases, ANGPTL4, and L1CAM mRNA transcripts (3–5). In addition to regulating thousands of target genes at the transcriptional level, hypoxia can also affect posttranscriptional regulation, including alternative splicing of pre-mRNA (6). It is estimated that over 95% of human genes undergo alternative splicing, which can generate mRNAs encoding proteins with different, even opposite functions (7). However, whether and how hypoxia influences cancer metastasis by alternative splicing events is largely unknown.

Proteins of the methyl-CpG binding domain (MBD) family are main candidates for the readout of DNA methylation (8, 9). We previously reported that MBD2 represses reprogramming of human fibroblasts into induced pluripotent stem cells via downregulation of the glycolytic pathway (10). The MBD2 gene encodes at least three isoforms, including MBD2a, MBD2b, and MBD2c, among which, MBD2a and MBD2c have been the most studied splice variants (11). MBD2a and MBD2c play distinct roles in the self-renewal of human pluripotent stem cells and reprogramming of fibroblasts (12).

The long isoform MBD2a contains a GR-rich region and MBD domain near the N-terminus and a transcriptional repression domain (TRD) at the C-terminus, whereas the short isoform MBD2c shares the same MBD domain but differs at the C-terminus (9, 13). This structural difference leads to the interaction of MBD2a, but not MBD2c, with the nucleosome remodeling and deacetylase (NuRD) complex, which plays vital roles in transcriptional repression. As a reader of methylated DNA, MBD2a recognizes and binds to hypermethylated DNA promoters and recruits the NuRD complex to repress the expression of genes encoding tumor suppressors and thus plays an oncogenic role in several human cancers including breast, colorectal, prostate cancer, and glioblastoma (13–15). For instance, MBD2a stimulates tumorigenesis by binding to the aberrantly methylated promoter regions of the gene encoding the tumor suppressors p14/ARF and p16/Ink4A in colon cancer, and the GSTP1 gene in breast cancer (16, 17). Moreover, a previous study showed that MBD2a promotes breast cancer metastasis through activation of MMP2 (18). However, little is known regarding the functional role of MBD2c in tumorigenesis or metastasis.

In this study, we set out to determine the role of MBD2 splice variants in regulation of breast cancer metastasis. We discovered that MBD2a promotes, whereas MBD2c suppresses metastasis. Specifically, we found that HIF1 suppressed MBD2 alternative splicing by inducing miR-222, which inhibited the expression of SRSF2 under hypoxia. Moreover, our RNA-sequencing results showed that MBD2a and MBD2c inversely regulate the expression of β-catenin and Snail1, two EMT markers. Furthermore, chromatin immunoprecipitation sequencing (ChIP-seq) results revealed and ChIP-PCR confirmed that MBD2a promotes while MBD2c suppresses the expression of FZD1, a Wnt/β-catenin signaling pathway receptor that promotes EMT, by competitive binding to its CpG islands. Collectively, our results not only provide promising biomarkers to predict the invasiveness of breast cancer, but also establish a previously unappreciated link between MBD2 alternative splicing and breast cancer metastasis, which has important implications for understanding how breast cancer metastasis is triggered in response to hypoxia.

Short hairpin RNA (shRNA) targeting MBD2a, SRSF2, FZD1, HDAC1/2, MTA2, and RBBP7 were commercially purchased (Sigma-Aldrich). MBD2c shRNAs were inserted into PLKO vector. HIF1α shRNAs were inserted into pRPL-GFP-Puro vector as reported previously (19). shRNA targeting sequences and information of MBD2a, MBD2c, or SRSF2 overexpression plasmids were listed in Supplementary Table S1. Each lentiviral plasmid was cotransfected with plasmids encoding Δ8.9 and VSVG into HEK293T packaging cells using PEI (Invitrogen). Viral supernatant was collected 48 hours posttransfection, filtered (0.22-nm pore size), and added to breast cancer cells in the presence of 8 μg/mL polybrene (Sigma-Aldrich) and stable cells were selected with puromycin.

HEK293, HEK293T, MDA-MB-231, and MDA-MB-468 cells were cultured in DMEM supplemented with 10% FBS and 1% penicillin/streptomycin. MCF-7, T47D, BT-474, ZR-75–1, BT-549, and SKBR3 cells were cultured in RPMI1640 supplemented with 10% FBS, 1.5 g/L NaHCO₃, 2.5 g/L glucose, 0.11 g/L sodium pyruvate, and 1% penicillin/streptomycin. SUM149 cells were cultured in Ham F12 medium supplemented with 10% FBS, 1% penicillin/streptomycin, 5 μg/mL insulin, and 1 μg/mL hydrocortisone. All cell lines were authenticated using short tandem repeat (STR) profiling (by GENEWIZ Co. Ltd.). After authentication, large frozen stocks were made for future use. All cell lines were used within 15 passages (less than 2 months) after reviving from the frozen stocks. Cells were kept at 37°C in humidified 5% CO₂ in air. Trypan blue exclusion assay were used to assess cell growth and survival. Hypoxia condition was achieved by placing cells in a hypoxic working station (Whitley H35 Hypoxystation, Don Whitley Scientific/DWS), which contains 1% O₂, 5% CO₂, and 94% N₂. Detailed information for all reagents used is provided in Supplementary Table S2.

Total cellular protein was isolated from cells using RIPA buffer (50 mmol/L Tris-HCl, pH 8.0, 150 mmol/L NaCl, 5 mmol/L EDTA, 0.1% SDS, and 1% NP-40) supplemented with cocktail. Protein concentration was measured using the Bradford assay kit. Equal amount of proteins was loaded and separated by SDS-PAGE. Actin served as loading control. Detailed information for all antibodies used is provided in Supplementary Table S2.

Total RNA was isolated using TRIzol reagent (Thermo Fisher Scientific) followed by DNA-free (Invitrogen) treatment and reverse transcription with the HiScript II 1st Strand cDNA Synthesis Kit (Vazyme). qPCR was performed using SYBR Green master mix (Vazyme) on a Bio-Rad iCycler. Primer sequences used were shown in Supplementary Table S1. All samples were normalized to 18S rRNA.

miRNA agomirs and antagomirs were purchased from GenePharma Company. 100 nmol/L miRNA agomirs or 100 nmol/L antagomirs were transfected into MDA-MB-231 cells in a 60-mm dish with Lipofectamine 2000 (Thermo Fisher Scientific) according to the manufacturer's instruction.

Full-length 3′UTR of SRSF2, or the SRSF2 3′UTR harboring wild-type or potential miR-222 binding site mutants were inserted into the pSI-CHECK-2 dual-luciferase reporter vector. HEK293 cells were cotransfected with the reporter plasmids along with control agomir or miR-222 agomir using Lipofectamine 2000 (Invitrogen) in a 48-well plate. Luciferase activity was measured after transfection for 48 hours using the Dual-Luciferase Reporter Assay System (Promega). Renilla luciferase was normalized to firefly luciferase activity.

RNA immunoprecipitation (RIP) assay was performed following the procedure as we reported previously (20). Primary antibody against GFP (ProteinTech) was used in the RIP assay. The primers sequences used in these experiments are listed in Supplementary Table S1. In brief, cells were cross-linked in 0.1% formaldehyde prior to lysis. Cell lysates were sonicated, immunoprecipitated, and the eluates were reverse-crosslinked. Relative occupancy values were calculated by determining the IP efficiency and normalized to the level observed by immunoprecipitation using nonspecific IgG.

Total RNA was extracted and checked for a RIN number to inspect RNA integrity by an Agilent Bioanalyzer 2100. Detailed RNA-sequencing and analysis were determined on the basis of our previous protocol (21). A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEB Next Ultra RNA Library Prep Kit for Illumina (NEB). Reads were aligned to the human genome hg19 with tophat2 v2.1.1. Enrichment pathway analysis of genes was compiled from Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases by DAVID Bioinformatics Resources (22). Original data are available in the NCBI Gene Expression Omnibus (GEO; accession number GSE151952).

For the cell migration assay, 1 × 10⁵ MDA-MB-231 cells with different expression of MBD2a or MBD2c were cultured in the top chamber of the 24-well transwell plates (Corning Inc.) in serum-free DMEM, 600 μL of DMEM containing 20% FBS was added to the bottom chamber. Cells were then cultured, fixed, and stained. The migrated cells were counted under a microscope in a blinded manner.

For monitoring of cell migration in real time, the xCELLigence Real Time Cell Analyzer Dual Plate (RTCA DP) instrument was used according to the manufacturer's recommendations (Roche Applied Science). Briefly, 4 × 10⁴ cells in serum free DMEM were seeded per well of a 16-well CIM plate, and the bottom chamber was loaded with DMEM, supplemented with 20% FBS. Quantitative measurements were monitored and plotted using the RTCA software 1.2.1 of the RTCA xCELLigence system.

Genomic DNA from cells was extracted and then underwent bisulfite modification utilizing the Bisulfite Conversion Kit (Active Motif). Modified DNA was amplified by PCR. Sequences of the primers designed to detect the methylation status of CpG sites were listed as in Supplementary Table S1. For bisulfite DNA sequencing, PCR products were gel-purified and cloned into the TA-cloning Vector (Takara) according to the manufacturer's protocol. Plasmid DNA was purified and then sequenced. Results were analyzed by BIQ-Analyzer software (23).

The ChIP assay was performed with an EZ-ChIP kit (Millipore) following the manufacturer's instructions. Briefly, cells were fixed with 1% formaldehyde and quenched with 0.125 mol/L glycine. Cells were sonicated by Scientz JY99-IIDN. DNA was immunoprecipitated by either control IgG or Flag (Sigma-Aldrich) antibody, followed by qPCR analysis. The oligos used for this analysis are listed in Supplementary Table S1.

The ChIP was performed with an EZ-ChIP kit (Millipore) following the manufacturer's instructions. For the MDA-MB-231 cells’ ChIP-seq experiments with spike-in control, Drosophila spike-in chromatin (Active Motif) and a Drosophila-specific spike-in antibody (Active Motif) were used for each reaction. The Drosophila melanogaster–specific antibody is used to immunoprecipitate Drosophila melanogaster chromatin from each reaction, which is then used as a reference for normalization during data analysis (21). ChIP-seq was performed at Novogene on the Illumina novaseq6000 platform. Sequenced reads were separately aligned to either the human genome hg19 or the Drosophila melanogaster genome (dm6) using Bowtie2 software (version 2.2.6; ref. 24). Peaks were classified on the basis of the location (UCSC annotation data) and shown in the following genome regions: intergenic, introns, downstream, upstream, and exons (25). The gene traces were visualized using the Integrative Genomics Viewer (IGV), while showing the enrichment of CHIP sequencing reads on FZD1 gene sequence (26). Original data are available in the NCBI GEO (accession number GSE151952).

Gene set enrichment analysis (GSEA; version 4.0.1; http://software.broadinstitute.org/gsea/index.jsp) was used to perform gene set enrichment analysis on various functional and/or characteristic gene signatures obtained from the Msig database (https://www.gsea-msigdb.org/gsea/index.jsp; refs. 27, 28). Statistical significance was assessed by comparing the ES to enrichment results generated from 1,000 random permutations of the gene set to obtain P values (nominal P value).

All animal studies were conducted with approval from the Animal Research Ethics Committee of the University of Science and Technology of China. For xenograft experiments, MDA-MB-231 cells with differential MBD2a or MBD2c expression were injected into the mammary fat pad (MFP) of 5-week-old female nude mice (SJA Laboratory Animal Company). Eight days after injection, tumor volumes were measured every 2 or 3 days with a caliper and calculated using the equation: volume = width × depth × length × 0.52.

All human tissue specimens used in this article were obtained from The First Affiliated Hospital of Medical University of Anhui (Hefei, China). For using of these clinical materials for research purposes, prior patients’ written informed consents and approval from the Ethics Committee of Anhui Medical University (Hefei, China) were obtained. The studies were conducted in accordance with ethical guidelines of the Declaration of Helsinki.

Following the procedure we reported previously (20), IHC and hematoxylin and eosin (H&E) staining were performed with formalin-fixed paraffin-embedded tissue. MBD2a (1:100; Abcam) antibody was used for IHC. All images were taken with a Zeiss AxioImager Z1, and quantification was performed with HistoQuest (TissueGnostics GmbH; www.tissuegnostics.com). Immunoreactivity was also graded according to the percentage of positive tumor cells (0, < 2%; 1, 2%–5%; 2, 6%–50%; 3, > 50%), and the intensity of staining (0, no staining; 1, weak; 2, moderate; and 3, strong). Cut-off levels for the sum of scores were 0 to 2 for low expression, 3 to 4 for moderate expression, and 5 to 6 for high expression.

Human in situ hybridization (ISH) was performed by the RNAscope 2.0 High Definition assay (29) according to the manufacturer's instructions (Advanced Cell Diagnostics, Inc.). In brief, cells in the clinical specimens are fixed and permeabilized using xylenes, ethanol, and protease to allow for probe access. Slides are boiled in pretreatment buffer for 15 minutes and rinsed in water. Next, two independent target probes are hybridized to the MBD2a/MBD2c RNA at 40°C for 2 hours, with this pair of probes creating a binding site of a preamplifier. The preamplifier is hybridized to the target probes and amplified. Cells are counterstained to visualize signal, followed by H&E staining and mounting in a xylene-based mounting media (29). Quantification was performed with HistoQuest.

Statistical analysis was determined using Fisher exact probability for Figs. 3A and 4C and Supplementary Fig. S4C. Spearman test for Supplementary Fig. S3A, and Student t test for other data. Data were presented as the mean (± SD) (Figs. 1–4) or mean (± SEM) (Fig. 5) of at least three independent experiments. * P < 0.05 or ^#P < 0.05 was considered significantly different.

Both MBD2a and MBD2c function to regulate self-renewal of human pluripotent stem cells, but only MBD2a has been shown to regulate cell proliferation in a tumor metastasis context. Interestingly, we found that expression of MBD2a was increased, whereas MBD2c was decreased, in highly invasive breast cancer cell lines at both the mRNA and protein levels (Fig. 1A–C). These results indicate that MBD2a and MBD2c may exert inverse functions in regulating tumor metastasis. Pursuing this idea of inverse functions for the two isoforms, we used the highly metastatic triple-negative breast cancer MDA-MB-231 and MDA-MB-468 cells (which have low expression of MBD2c but have relatively high expression of MBD2a). We overexpressed MBD2c or knocked down MBD2a using confirmed MBD2a-specific shRNA constructs. Immunoblotting and quantitative reverse transcriptase PCR (qPCR) confirmed the expected increase in the MBD2c level and the expected MBD2a decrease in these cells (Fig. 1D and E; Supplementary Fig. S1A). Transwell assays demonstrated that the cell migration capacity of MDA-MB-231 and MDA-MB-468 cells is significantly inhibited by MBD2c overexpression and by MBD2a knockdown (Fig. 1D and E; Supplementary Fig. S1B).

We also successfully overexpressed MBD2a and specifically knocked down MBD2c in MDA-MB-231 and MDA-MB-468 cells (Supplementary Fig. S1C–S1E). Transwell assays with these inversely manipulated MBD2 isoforms revealed significantly increased cell migration capacities for both cell lines (Supplementary Fig. S1D–S1F). Although depletion of MBD2a led to a decrease in cell growth, there were no differences in the proliferation rates of the other three cell lines (Supplementary Fig. S1G and S1H). We also monitored cell migration in real time using an xCELLigence Real Time Cell Analyzer Dual Plate (RTCA DP) instrument, which quantifies cell migration and yields the data as an index (30). Consistent with our transwell assays, the MBD2c overexpression cells and the MBD2a knockdown cells had significantly reduced basal migration rates as compared with their respective empty vector and NTC controls. Again, continuing the trend, the MBD2a overexpression cells and the MBD2c knockdown cells exhibited significantly increased migration rates (Fig. 1F and G). Thus, multiple lines of in vitro evidence support that MBD2a promotes and MBD2c represses the cell migration capacity of breast cancer cells.

To assess the clinical relevance of these findings, we examined the expression levels of MBD2a and MBD2c in human breast cancer samples, including the less invasive ductal carcinoma in situ (DCIS) tissues and the more invasive carcinoma tissues. Immunoblotting and qPCR results showed that MBD2a was expressed at significantly higher levels in the invasive breast cancer tissues as compared with DCIS tissues; consistently, the inverse trends were observed for MBD2c (Fig. 1H and I). Next, IHC was employed to analyze MBD2a expression in a retrospective cohort of 116 clinicopathologically characterized breast cancer samples (assessed for pathologic grades), including 9 cases of DCIS, 84 cases of invasive carcinoma, and 7 normal breast control tissues; we also examined lung tissues from 16 cases of breast cancer–related lung metastasis (Supplementary Fig. S1I). The IHC results revealed that MBD2a protein expression was generally not detectable in normal breast tissue, weak in DCIS, abundant in invasive carcinomas, with the highest detected levels in lung metastases (Fig. 1J). Quantitative analysis of the IHC results demonstrated that MBD2a expression of invasive samples was significantly increased compared with DCIS samples and normal breast control tissues, and MBD2a was significantly upregulated in lung metastasis compared with primary breast carcinoma (without or with metastasis), suggesting that the MBD2a protein levels increase during breast cancer progression (Fig. 1K).

Given the lack of a commercially available MBD2c antibody for the IHC detection, we successfully adopted an RNA ISH strategy that enabled specific detection of the two MBD2 isoforms (Supplementary Fig. S1J). Consistent with the protein data in our IHC results, tissues from the 22 invasive breast cancer cases we examined had obviously increased levels of MBD2a transcripts compared with the 23 cases of DCIS tissues. Moreover, the MBD2c transcript levels were obviously lower in the invasive tissues (Fig. 1L). Specific ISH signals of MBD2a and MBD2c mRNA in morphologically intact cells were assessed quantitatively using HistoQuest software, which revealed that the mean intensity of the signal for MBD2a transcripts was obviously higher in the invasive breast cancer samples compared with DCIS samples, and a significant inverse trend was detected for MBD2c transcripts (Fig. 1M). Taken together, these results support the clinical relevance of differential MBD2 isoform accumulation, specifically indicating that increased MBD2a, accompanied with decreased MBD2c mRNA levels are consistently associated with more aggressive, invasive, and metastatic capacity.

Next, we investigated the mechanism for the contextual changes of MBD2 isoforms in tumor metastasis. It is well known that intratumoral hypoxia is a hallmark of advanced cancers. We observed that exposure of MDA-MB-231 or MDA-MB-468 cells to hypoxic culture conditions induced MBD2a and repressed MBD2c mRNA (Fig. 2A). Similar results were obtained in other breast cancer cell lines, including MCF7, T47D, and SUM-149 cells, under hypoxic conditions (Supplementary Fig. S2A–S2C). Moreover, MBD2a protein levels were upregulated, whereas MBD2c was downregulated under hypoxia (Fig. 2A; Supplementary Fig. S2D). Notably, specific primers were used to amplify only MBD2 pre-mRNA, and qPCR results showed that MBD2 pre-mRNA remained the same under hypoxia (Fig. 2A), indicating that hypoxia triggers a shift between MBD2 splice isoforms.

A previous study demonstrated that splicing factor SRSF2 (also known as SC35 or SFRS2) was required for the alternative RNA splicing of MBD2 in human pluripotent stem cell (12). In MDA-MB-231 cells, results of RIP (RNA immunoprecipitation) assay showed that SRSF2 bound to MBD2a pre-mRNA at intron 2, which is consistent with that in IPS cells (Supplementary Fig. S2E; ref. 12). Our qPCR and immunoblotting assays further demonstrated that overexpression of SRSF2 led to markedly increased MBD2c and decreased MBD2a mRNA and protein (Fig. 2B). Opposite results were observed by specific shRNAs targeting SRSF2 in MDA-MB-231 cells (Supplementary Fig. S2F). More interestingly, hypoxia suppressed SRSF2 expression in different breast cancer cell lines at both the mRNA and protein levels (Fig. 2C; Supplementary Fig. S2A–S2D). Furthermore, we established MDA-MB-231 or MDA-MB-468 cells stably expressing HIF1α specific shRNAs (Supplementary Fig. S2G), and found that hypoxia-induced downregulation of SRSF2 and MBD2c, and upregulation of MBD2a, were reversed by shHIF1α at both the mRNA and protein levels (Fig. 2D). Furthermore, SRSF2 overexpression blocked hypoxia-induced switching from MBD2c to MBD2a in MDA-MB-231 and MDA-MB-468 cells (Fig. 2E), demonstrating that HIF1α–mediated SRSF2 inhibition serves as upstream regulator of MBD2 switching in cancer cells.

Given the role of miRNAs in RNA silencing and posttranscriptional regulation of gene expression, we next carried out a comprehensive bioinformatics analysis and found that miR-200c and miR-222 were candidate miRNAs that not only potentially bind the 3′UTR of SRSF2 mRNA, but also possess potential hypoxia response elements (HREs) in their promoters (Supplementary Fig. S2H). qPCR assay confirmed the increased expression of both miR-200c and miR-222 in hypoxic MDA-MB-231 and MDA-MB-468 cells (Fig. 2F). Further experiments using their specific agomirs or antagomirs revealed that miR-222, but not miR-200c, inhibited expression of SRSF2 and MBD2c, while enhancing expression of MBD2a at both the mRNA and protein levels in MDA-MB-231 cells (Fig. 2G; Supplementary Fig. S2I). Dual-luciferase reporter assay demonstrated that miR-222 bound to the 3′UTR of SRSF2 mRNA (Supplementary Fig. S2H; Fig. 2H). More importantly, hypoxia-induced suppression of SRSF2 and MBD2c, as well as MBD2a accumulation, were abolished by introducing miR-222 antagomir into MDA-MB-231 cells (Fig. 2I). Taken together, our data demonstrate that hypoxia inhibits production of the MBD2c splice isoform via miR-222-mediated SRSF2 suppression.

To pursue the mechanism of opposing action of MBD2a and MBD2c in breast cancer metastasis, RNA-sequencing was performed in MDA-MB-231 cells with MBD2a knockdown or MBD2c overexpression. RNA-seq data showed that 1,390 genes were similarly regulated by MBD2a knockdown and MBD2c overexpression, which account for 20.1% of MBD2a-knockdown and 30.6% of MBD2c-overexpression–regulated genes, respectively (Fig. 3A). Heatmap showed all regulated genes by shMBD2a or MBD2c overexpression (up- and downregulated genes are defined as having a log₂-fold change greater than 0.5 or less than -0.5, respectively; Supplementary Fig. S3A). Gene ontology term enrichment analysis (DAVID: https://david.ncifcrf.gov/) showed that genes similarly affected by MBD2a knockdown and MBD2c overexpression play prominent roles in multiple biological processes involved in tumor metastasis, including cellular component organization or biogenesis, locomotion, and biological adhesion (Fig. 3B; Supplementary Table S3). GSEA of RNA-seq data revealed that MBD2c overexpression or shMBD2a was positively associated with three kinds of gene signatures: HALLMARK_PROTEIN_SECRETION, HALLMARK_TGF_BETA_SIGNALING, and GO_BIOLOGICAL_ADHESION, but was negatively correlated with RICKMAN_METASTASIS_DN (Supplementary Fig. S3B). Therefore, these analyses confirmed the negative enrichment of metastasis-related gene signatures in MBD2c overexpressing or MBD2a knockdown cells compared with control cells.

Seventeen genes with altered expression according to RNA-seq data (Supplementary Table S4) involved in cell migration and metastasis were studied in MDA-MB-231 cells expressing shMBD2a or MBD2c (Fig. 3C and D; Supplementary Fig. S3C and S3D), which confirmed our RNA-seq data. Of note, MBD2a enhanced while MBD2c suppressed the expression of β-catenin and Snail1, two EMT markers, at both the mRNA and protein levels (Supplementary Fig. S3A; Fig. 3C–H). Similar results were observed in other breast cancer cell lines including MCF7, T47D, MDA-MB-468, and SUM-149 (Supplementary Fig. S3E–S3L). Further experiments showed that hypoxia-induced β-catenin and Snail1 accumulation was abolished by MBD2a knockdown or MBD2c overexpression (Fig. 3I and J). Transwell assays demonstrated that the hypoxia-induced migration of MDA-MB-231 cells was significantly suppressed by MBD2a knockdown or MBD2c overexpression (Fig. 3K and L), documenting that MBD2a and MBD2c play opposing roles on the expression of β-catenin and Snail1, and subsequent EMT and breast cancer metastasis under hypoxia.

More interestingly, MBD2a-stimulated and MBD2c-suppressed expression of β-catenin and Snail1 was markedly reversed by forced expression of MBD2c and MBD2a, respectively (Fig. 3M and N). On the other hand, overexpressed MBD2c and MBD2a potentiated the effect of shMBD2a and shMBD2c on β-catenin and Snail1 expression, respectively (Fig. 3O and P). These data prove that MBD2a and MBD2c have competitively antagonistic roles in regulating β-catenin and Snail1 expression, which provides a molecular correlate to our findings that MBD2a and MBD2c play opposite roles in breast cancer metastasis.

MBD2 has been shown to mainly mediate gene silencing of tumor suppressors by recruiting NuRD complex to hypermethylated gene promoters. Thus, next, ChIP-seq was performed in MDA-MB-231 cells with Flag antibody against Flag-tagged MBD2a without or with overexpression of GFP-MBD2c, and Flag-tagged MBD2c without or with overexpression of GFP-MBD2a, respectively. As a result, ChIP-seq data revealed that both MBD2a and MBD2c were enriched on DNA at MBD2a or MBD2c binding sites (Fig. 4A, lanes 2 and 4; Supplementary Fig. S4A, lanes 2 and 4), suggesting highly similar binding sequences of MBD2a and MBD2c. Most importantly, overexpression of MBD2c reduced the enrichment of MBD2a on DNA (Fig. 4A, lane 2 vs. lane 3; Supplementary Fig. S4A, lane 2 vs. lane 3), and vice versa (Fig. 4A, lane 4 vs. lane 5; Supplementary Fig. S4A, lane 4 vs. lane 5). Similar results were observed in aggregation plots analysis (depicting the average ChIP signal across a variety of genomic features) of the distribution of MBD2a and MBD2c tags around the MBD2a peaks or MBD2c peaks (Fig. 4B; Supplementary Fig. S4B), suggesting a mutual inhibition competition between MBD2a and MBD2c in their binding to DNA.

Further analysis of the ChIP-seq data showed that MBD2a binds to 6,445 genes and MBD2c binds to 4,405 genes, among which, 3,138 genes are overlapped (Supplementary Fig. S4C). On the other hand, RNA-seq data showed that 1,390 genes were simultaneously regulated by MBD2a and MBD2c (Fig. 3A). Combined analysis of the ChIP-seq and RNA-seq data showed only 174 genes were cobound and coregulated simultaneously by both MBD2a and MBD2c (Fig. 4C). DAVID analysis revealed that these 174 genes are mainly involved in locomotion, cellular component organization or biogenesis, and biological adhesion (Supplementary Fig. S4D and E; Supplementary Table S5). Among which, FZD1, an important upstream regulator of Wnt/β-catenin pathway and EMT, is on our top list of the interest (Supplementary Fig. S3A). Binding of G protein–coupled receptor FZD1 to WNTs activates the WNT/catenin pathway, which activates the gene transcription and mediates a variety of functions including cell growth, metastasis, and survival (31, 32). Paracrine activation of FZD1/WNT/β-catenin pathway supports the proliferation and growth of uterine leiomyomas (33). In addition, hypermethylation of FZD1 was associated with the inactivation of Wnt/catenin signal (34). Integrative Genomics Viewer (IGV) analysis revealed that the binding of MBD2a on the high CpG density regions of FZD1 gene were markedly reduced by MBD2c overexpression, and inverse trend was observed (Fig. 4D). Moreover, we selected three potential MBD2-binding sites upon CpG islands of FZD1 (Supplementary Fig. S4F). Our ChIP experiment showed that both MBD2a and MBD2c bound to these three sites, and the binding exhibited a mutual inhibitory competition pattern between MBD2a and MBD2c (Fig. 4E), confirming that MBD2a and MBD2c colocalized on chromatin and competitively bound to the same sequences.

Next, qPCR and immunoblotting assays showed that FZD1 expression was repressed by MBD2a depletion or MBD2c overexpression, and enhanced by MBD2c depletion, MBD2a overexpression, or hypoxia in MDA-MB-231 cells (Fig. 4F and G). FZD1 inhibition markedly suppressed the expression of β-catenin and Snail1 at both mRNA and protein levels in breast cancer cells (Fig. 4H; Supplementary Fig. S4G and S5H). More interestingly, suppression of FZD1 dramatically abolished the elevated expression of β-catenin and Snail1 by MBD2a overexpression, shMBD2c, or hypoxia exposure in breast cancer cells (Fig. 4I and J), suggesting that FZD1 is involved in MBD2a and MBD2c regulated EMT under hypoxic conditions.

MBD2a has long been considered as a transcriptional repressor by recognizing and binding to hypermethylated DNA promoters together with other complexes (35), but our results showed that knockdown of histone deacetylase complexes HDAC1/2, metastasis-associated gene 1 (MTA2) or retinoblastoma-binding protein 7 (RBBP7), the partners of MBD2a when represses gene transcription, exhibited no effect on the mRNA levels of MBD2a, MBD2c, FZD1, β-catenin, and Snail1 (Supplementary Fig. S4I), indicating that MBD2a- and MBD2c-regulated FZD1 transcription is independent of this mechanism.

Previous studies reported that MBD2a positively regulates gene transcription by acting as a DNA demethylase by removing repressive methyl residues (36, 37). Thus, next, we performed bisulfite sequencing PCR (BSP) and methylation-specific PCR (MSP) assays in MDA-MB-231 cells (Supplementary Fig. S4F). Results of BSP (Supplementary Fig. S4J–S4L) and MSP (Fig. 4K and L) assays showed that methylation of the FZD1 promoter was enhanced by shMBD2a or MBD2c overexpression, which resulted in a significant reduction of mRNA transcripts and protein levels of FZD1 (Fig. 4K and L). In the presence of 5-Aza-2′-deoxycytidine (decitabine; DAC), an inhibitor of DNA methyltransferase activity, the hypermethylated status of FZD1 promoter was repressed and, FZD1 mRNA levels increased (Fig. 4K and L). These results suggest that MBD2a promotes FZD1 transcription by reducing the methylation levels of the CpG islands of FZD1, whereas MBD2c functions oppositely. Taken together, MBD2a and MBD2c inversely regulate the expression of EMT markers by competitively binding to FZD1.

To further evaluate the effect of MBD2 splice variants on tumor metastasis, orthotopic mouse models were employed by mammary fat pad injection of MDA-MB-231 cells. Tumor growth tracking showed that MBD2a overexpression or MBD2c knockdown showed no significant effect on primary tumor growth (Supplementary Fig. S5A, S5C, S5F, and S5H), although MBD2c overexpression or MBD2a knockdown marginally delayed tumor growth (Supplementary Fig. S5B, S5D, S5E, and S5G). But importantly, both the number and size of lung metastases were obviously increased in mice bearing MBD2a-overexpressing compared with control cells, which were attenuated by MBD2c overexpression (Fig. 5A). The inverse trends were observed for MBD2c overexpression, which were abolished by MBD2a overexpression (Fig. 5B). On the other hand, MBD2c overexpression exacerbated the inhibition of breast cancer lung metastasis caused by MBD2a knockdown (Fig. 5C), while MBD2a overexpression potentiated the increased breast cancer lung metastasis caused by MBD2c knockdown (Fig. 5D). To obtain a more sensitive estimate of overall lung metastatic burden, genomic DNA was isolated from the contralateral lung and qPCR was performed using primers that only amplify human DNA. As a result, the metastatic burden in the lungs of mice confirmed the histologic analysis results (Fig. 5A–D, right).

Immunoblot assays using protein lysates from breast tumor tissues at the end of the experiment confirmed the expression of MBD2a and MBD2c in different groups (Fig. 5E–H). Consistent with the data from cultured cells, we also observed the inverse roles of MBD2a and MBD2c in regulating the expression of FZD1, β-catenin, and Snail1 in mouse tumor tissues (Fig. 5E–H). Collectively, these data provide evidence for the opposing roles of MBD2a and MBD2c in regulating spontaneous metastasis of primary breast tumors to the lungs in vivo.

Metastasis accounts for the great majority of breast cancer–related deaths worldwide, yet the biological underpinnings of metastatic diseases remain the least understood in contrast to the large body of evidence that have uncovered the detailed mechanisms leading to primary tumor formation. Metastasis is facilitated by epigenetic inheritance and alternative splicing of multiple genes involved in hypoxic microenvironment and methylation status (38). In this study, we demonstrated that suppression of MBD2 alternative splicing promoted breast cancer metastasis under hypoxia conditions. Hypoxia-induced miR-222 inhibited the expression of SRSF2, a splicing factor, and resulted in increased MBD2a and decreased MBD2c expression, which consequently promoted tumor metastasis through FZD1-regulated EMT both in vitro and in vivo. Our findings illustrate a novel link and the underlying mechanisms between MBD2 switching and breast cancer metastasis under hypoxia context.

Hypoxic microenvironment is a notable hallmark of solid tumor formation and progression and about 25%–40% of invasive breast cancers exhibit hypoxic regions (39). HIFs are comprehensively involved in tumor metastasis cascade including invasion, intravasation, extravasation, and metastatic niche formation (38, 40). In addition, hypoxia is involved in the alternative splicing of HIF and non-HIF target genes (6, 41). However, up to date, little is known about whether and how hypoxia-regulated alternative splicing affects breast cancer metastasis. Here, we show that HIF1 suppresses alternative splicing of MBD2 to promote breast cancer metastasis. We found that HIF1, as an upstream regulator, suppresses MBD2 alternative splicing at both the mRNA and protein levels without affecting the expression of MBD2 pre-mRNA (Fig. 2A). Our results further showed that SRSF2, critical for MBD2 splicing, was repressed by HIF1-induced miR-222 in tumor cells (Fig. 2B and C, F–I), and hypoxia-induced switching of MBD2c to MBD2a was dependent on HIF1–mediated SRSF2 inhibition (Fig. 2E). Furthermore, transwell assays showed that MBD2a promotes, whereas MBD2c suppresses breast cancer metastasis by regulating β-catenin and Snail1 expression in reverse, suggesting the difference between MBD2a and MBD2c in affecting breast cancer metastasis (Fig. 3).

Significantly, we observed that MBD2a and MBD2c competitively bound to the CpG islands of FZD1 and inversely regulate the expression of downstream EMT markers (Fig. 4A, B, D, and E). MBD2 has been considered as a transcriptional repressor; however, our results showed that MBD2a positively regulates FZD1, β-catenin, and Snail1 expression at transcriptional level and knockdown of the HDAC components did not influence their mRNA transcripts (Supplementary Fig. S4I). Previous reports documented that MBD2a, as a DNA demethylase, upregulates uPA and Foxp3 expression by removing repressive methyl residues and thereby giving rise to promoter-specific gene transcriptional activation (36, 37, 42). Moreover, it was speculated that MBD2 recruits “activators” to turn on gene expression. The complexes formed by MBD2 with other proteins such as cAMP-responsive factor (CEBPA), MBD2-interacting protein (MBDin), transforming-acid-coiled–coil (TACC3), focal adhesion kinase (FAK/PYK2) and nerve growth factor–inducible protein A (NGFI-A), in most cases, are mutually exclusive with HDAC-containing complexes thus relieving the repression potential of MBD2 even prior to eventual demethylation (43–47). In addition, it is not yet clear whether MBD2a directly promotes DNA demethylation, or performs this function by recruiting or regulating other DNA demethylation enzymes such as ten eleven translocation DNA demethylases (TETs) or TETs partner MLL observed in Treg cells (42, 48). Because of the structural difference at C-terminus of MBD2a and MBD2c, it is worthwhile to study whether MBD2c binds to HAT, TACC3, HTLV-1 TAX1, or other unknown MBD2a-binding transcriptional activators (43, 49). Here, we found that the methylation level of FZD1 was increased when depletion of MBD2a or overexpression of MBD2c by BSP or MSP assays (Fig. 4K and L; Supplementary Fig. S4J–S4L). Intriguingly, after addition of DAC, an inhibitor of DNA methyltransferase, into the culture of these cells, the increased methylation level of FZD1 was reversed, which resulted in an elevated FZD1 mRNA level (Fig. 4K and L). Thus, our results indicate that altered methylation status of FZD1 promoter when manipulating the expression of MBD2 splice variants could be responsible for the change of FZD1 mRNA expression. Our results significantly extended the distinct functions of MBD2 splice variants from self-renewal of human pluripotent stem cells and somatic reprogramming (12) to tumorigenesis and metastasis.

Compared with MBD2a, there are relatively few studies on the role of MBD2c in breast cancer metastasis or how MBD2c regulates gene expression. In this regard, we found that MBD2c suppresses breast cancer metastasis by decreasing FZD1-mediated β-catenin and Snail1 expression (Figs. 3 and 4). Strikingly, for the first time, our clinical breast cancer sample analysis through IHC and RNA ISH revealed that high-MBD2a but low-MBD2c expression correlate significantly with metastasis degree of human breast carcinoma, suggesting that MBD2a and MBD2c with differential expression may serve as potential biomarkers for the prediction of breast cancer invasiveness (Fig. 1J–M). Furthermore, we presented multiple evidence here to support the hypothesis that MBD2c and MBD2a competitively bind to the same gene area and affect the methylation status of target genes of MBD2a (Fig. 4; Supplementary Fig. S4). Thus, our results provide previously unsuspected insights on the roles of MBD2c in breast cancer metastasis and how MBD2c regulates gene expression. Nonetheless, this study did not provide the mechanistic details such as the sequential order of how these events occur in cancer cells. Crystal structural analysis and numerous other experiments will be needed to delineate the detailed and distinct mechanisms of how MBD2a and MBD2c regulate gene expression, which would be of great potential significance for clinical intervention and thus deserving further investigation.

Taken together, we have provided evidence here to demonstrate that suppression of MBD2 alternative splicing promotes breast cancer metastasis under hypoxia. It's well known that HIF1 not only plays important roles in pathogenesis of multiple diseases, but also has very important physiological significance. Though many clinical studies have shown that intratumoral HIFs is closely related to patient mortality and a large body of preclinical data suggests that HIF inhibitors will improve outcome when added to current therapy in a wide range of cancer types, currently available pan-HIF inhibitors including digoxin, acriflavine, echinomycin, and YC-1 all have dose-limiting side effects that will prevent their use as clinical agents (50). As a result, it seems unrealistic to target HIF1 alone for the treatment of cancers and metastasis. Our findings herein raise exciting possibilities of targeting MBD2 switching together with HIF inhibitors as a potential therapeutic strategy against breast cancer metastasis, particularly as mediated by hypoxia and HIF1.

No disclosures were reported.

Z. Liu: Conceptualization, data curation, software, validation, investigation, methodology, writing-original draft, writing-review and editing. L. Sun: Conceptualization, data curation, funding acquisition, validation, methodology, writing-original draft, writing-review and editing. Y. Cai: Resources, data curation. S. Shen: Software, investigation. T. Zhang: Data curation, investigation. N. Wang: Investigation. G. Wu: Investigation. W. Ma: Investigation. S.-T. Li: Data curation, investigation. C. Suo: Investigation, methodology. Y. Hao: Investigation. W.-D. Jia: Resources, investigation. G.L. Semenza: Formal analysis, writing-review and editing. P. Gao: Conceptualization, data curation, supervision, funding acquisition, methodology, writing-original draft, project administration, writing-review and editing. H. Zhang: Conceptualization, data curation, supervision, funding acquisition, validation, investigation, writing-original draft, project administration, writing-review and editing.

This work is supported in part by National Key R&D Program of China (2018YFA0107103, 2018YFA0800300, 2017YFA0205600), National Natural Science Foundation of China (91957203, 81525022, 81930083, 81530076, 81821001, 81702361, 81874060), the Chinese Academy of Sciences (XDB39020100), the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2017ZT07S054), Outstanding Scholar Program of Guangzhou Regenerative Medicine and Health Guangdong Laboratory (2018GZR110102001), and the Fundamental Research Funds for the Central Universities (YD2070002008, 2020ZYGXZR038).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.