Progesterone modulates the DSCAM-AS1/miR-130a/ESR1 axis to suppress cell invasion and migration in breast cancer

Whole transcriptome sequencing data from 30 breast tumors samples were re-analyzed. Ten tumors were derived from patients who were administered a single dose of 500 mg of hydroxyprogesterone within 15 days prior to surgery, with varying duration for individual patients, while 20 tumors were obtained from patients who were not exposed to hydroxyprogesterone [37]. Gene expression was quantified using Salmon [38]. Genes with expression > 5 reads in at least 20% of the cancer samples were retained. Design matrices were created based on progesterone treatment, and differential gene expression analyses were performed with progesterone-treated (n = 10) and control (n = 20) tumor samples, using DESeq2 [39]. Data were assessed in the R environment. ENSEMBL IDs were converted using bioconductor packages (org.Hs.eg.db), and gene names not matching the ENSEMBL IDs were obtained from LNCipedia.

T47-D, BT-474, MCF7 and MDA-MB-231 breast cancer cells were procured, confirmed, cultured, and maintained as explained previously [13, 40]. The human embryonic kidney 293FT cell line was purchased from Invitrogen (Cat No. R70007), cultured in DMEM with 10% FBS, and maintained at 37℃ with 5% CO₂.

Cells were serum-starved and treated with 10 nM progesterone (6 h), 100 nM mifepristone (2 h), or an equal volume of ethanol (vehicle control) as explained previously [13]. RNA isolation, DNaseI treatment, and cDNA synthesis for genes/lncRNAs and microRNAs were performed as explained previously [29, 40]. Further, the cDNAs were used to study gene/miRNA expression patterns by quantitative real-time PCR (qPCR) method using the KAPA SYBR real-time PCR master mix (Sigma, Cat No. KK4601) and QuantStudio 5 real-time PCR system (Applied Biosystems, Cat no. A34322). GAPDH or ACTB and U6 were used as internal controls to normalize the expression of genes and miRNAs, respectively. Differential gene expression changes were calculated as fold change values using the 2^−∆∆CT method. The sequences for qPCR primers were manually designed using SnapGene sequence viewer. Designed primers were tested and optimized using OligoCalc (Sigma), UCSC In Silico PCR, and NCBI blast. The primer sequences for the genes and miRNAs are listed in Additional file 2: Table S1.

Total RNA was isolated from progesterone treated and untreated T47-D and MDA-MB-231 cells. Good quality RNA samples (RNA integration number > 9) were used to prepare the sequencing library using TruSeq library prep kit v2 (Illumina) with ribosomal RNA depletion. Libraries were sequenced on HiSeq4000 with 100 bp pair-end chemistry. A minimum of 60 million paired-end reads were obtained for each RNA sample with good Phred scores (score > 30). Differential gene expression analysis was performed using the salmon-DeSeq2 pipeline. Briefly, all the raw reads were corrected using the trimmomatic version V0.32 [41], followed by alignment to the human reference pseudo-genome (GRCh38) using Salmon (version: 0.8.2) [38] and differential expression analysis using DeSeq2 [39]. Genes/lncRNAs with fold change > 2 and < 0.5 with p-value < 0.05 were considered to be significantly deregulated in response to progesterone.

ChIP-sequencing data with PR, ER, and p300 pulldown for progesterone-treated T47-D cell line were downloaded from the SRA database [18]. The raw data for these experiments were analyzed as described earlier [40]. Briefly, reads were aligned to the gencode (v30) human reference genome (GRCh38) using a BWA aligner (version 0.7.17). Peak calling was performed using the MACS tool (version 2.0) [42]. Aligned reads were used for differential protein binding in the genome using DiffBind (version 3.0) [43]. The 5 kb upstream and downstream regions for annotated genes/lncRNAs were analyzed for PR, ER, and p300 binding, and annotation of the peaks was performed using Uropa [44].

DIANA-LncBase v2 [45] database was used to predict the binding of miRNAs to DSCAM-AS1. Predicted miRNAs binding to DSCAM-AS1 were further determined using the "microRNA–ncRNA targets” module of MirWalk v2.0 [46], which includes prediction algorithms of miRanda [47], RNAHybrid [48], and Targetscan [49]. Further, miRTarBase [50] and MirWalk v2.0 were mined to extract miRNAs targeting 3’-UTR of ESR1.

Sense and antisense DNA oligonucleotides with T7 RNA promoter sequences were designed and synthesized by Sigma-Aldrich to prepare siRNAs targeting DSCAM-AS1, PGR, and ESR1. The complete method for synthesis of small RNA transcripts using T7 RNA polymerase has been described previously [51]. Briefly, sense and antisense strands of DNA oligonucleotides with T7 RNA promoter complementary sequences were annealed in a Thermocycler for synthesizing dsDNA. The dsRNAs were subjected to in vitro transcription reaction (37 ℃ for 2 h) using T7 RNA polymerase (Promega, Cat no. P2075) in 1 × T7 Transcription Buffer (Promega, Cat no. P118B). The single-stranded sense and antisense siRNAs were further annealed to prepare double-stranded siRNAs. The complete list of DNA oligonucleotides for siRNA synthesis is provided in Additional file 2: Table S1. Synthesized siRNAs were transfected in breast cancer cells using Lipofectamine 3000 kit (Invitrogen, Cat No. L3000015) in serum-free media. Progesterone treatment was given to transfected cells 48 h post-transfection and collected for downstream analysis.

For transient overexpression of miR-130a in breast cancer cells, the precursor miRNA sequence with 200 bp flanking gene sequence was amplified from T47-D genomic DNA and cloned in pJET1.2/blunt vector (Thermo Scientific, Cat no. K1232) followed by sub-cloning in pcDNA3.1(-) mammalian expression vector under CMV promoter. XbaI (NEB, Cat no. R0145) and HindIII (NEB, Cat no. R0104) recognition sequences in multiple cloning sites of pcDNA3.1(-) were used for cloning. For transient overexpression of DSCAM-AS1, the complete cDNA sequence was cloned in the pcDNA3.1( +) expression vector using XbaI and BamHI (NEB, Cat no. R0136). Cloned constructs were confirmed by restriction digestion and Sanger sequencing. Further, the overexpression plasmids were transfected into breast cancer cells using Lipofectamine 3000. Empty pcDNA3.1(-) vector was transfected as vector control. Cells were collected 48 h post-transfection and RNA was extracted. Overexpression was confirmed by qPCR analysis.

For stable overexpression of DSCAM-AS1 in breast cancer cells, a complete cDNA sequence was cloned in the pBABE-puro expression vector using BamHI and SalI (NEB, Cat no. R0138). Cloning was confirmed using restriction digestion and Sanger sequencing. The primer sequences used for cloning and Sanger sequencing are provided in Additional file 2: Table S1. The 293FT cells were used for transfection and retrovirus production. Transductions were performed for 16 h in T47-D cells, followed by a selection of positive clones using 1 µg/mL puromycin (HiMedia, Cat no.TC198-10MG). Puromycin-resistant clones were further confirmed for DSCAM-AS1 overexpression by qPCR analysis.

Transwell cell migration and invasion assays were performed as described previously [13]. Briefly, cell invasion assay was performed with Matrigel loaded onto the inserts in Boyden chambers; while, cell migration assay was performed without Matrigel. The number of cells that migrated or invaded through the membrane was counted and the total fraction of cells was plotted as percent cell migration or invasion, respectively.

Full-length DSCAM-AS1 cDNA sequence was amplified from T47-D cells and cloned in pJET1.2/blunt vector, followed by sub-cloning in pGL3-promoter vector (Promega, Luciferase expressing vector) downstream to firefly luciferase using XbaI. DSCAM-AS1 mutant construct was generated with mutated overlapping primers by site-directed mutagenesis using DpnI (NEB, Cat no. R0176) and PrimeSTAR GXL DNA polymerase (TaKaRa, Cat no. R050B). DSCAM-AS1 mutant construct contained mutations in miR-130a MRE in DSCAM-AS1 cDNA to prevent miR-130a binding. Wild-type and mutated DSCAM-AS1 constructs were confirmed by Sanger sequencing. A 620 bp fragment of 3’-UTR of ESR1 with miR-130a binding site was amplified from T47-D cDNA and cloned in pJET1.2/blunt vector. The cDNA was sub-cloned in the pGL3-promoter vector using XbaI. The primer sequences for cloning and generating mutant construct are provided in Additional file 2: Table S1.

For luciferase assay, 293FT cells (50,000 cells/well) were co-transfected with pGL3-DSCAM-AS1 (wild-type/mutant) along with pcDNA3.1(-)-miR-130a or pcDNA3.1 empty vector using Lipofectamine 3000 kit. Additionally, co-transfections were performed with pGL3-3’-UTR-ESR1, pcDNA3.1(-)-miR-130a, or pcDNA3.1 empty vector. pEGFP-N2 was transfected to measure transfection efficiency in all wells. Cells were lysed 48 h post-transfection, and luciferase activity was measured using a luminometer (Berthold Luminometer, Germany). Luminescence and fluorescence units were measured from each transfected well. The luciferase activity was calculated by normalization of luminescence units with fluorescence units from the same well and plotted as luciferase activity. Each experiment was performed in triplicates.

The total RNA and miRNA sequencing data for patients with breast cancer were downloaded from The Cancer Genome Atlas (TCGA). Data from 751 breast cancer samples sequenced for total RNA and miRNAs were considered for further analysis. The samples with normally distributed DSCAM-AS1 or ESR1 expression values were segregated into quartiles. The upper and lower quartile samples were compared. The miRNA levels were compared between patients with ESR1-high and -low expression (the upper and lower quartiles, respectively). A similar analysis was performed for miRNAs in patients with DSCAM-AS1-high and -low expression. The significance of differences between both the groups was calculated using the Wilcoxon–Mann–Whitney test.

The TCGA breast cancer samples with high and low miRNA expression were compared for survival outcomes. The KM plotter [52] and GEPIA [53] were used for Kaplan–Meier survival analysis within specified breast cancer groups. Overall and relapse-free survival of patients was calculated based on the levels of lncRNAs and miRNAs in the samples.

GraphPad Prism version 8 (GraphPad Software, La Jolla, CA) was used to calculate statistical significance between different experimental groups in qPCR, cell-based assays, and luciferase reporter assays. The student's unpaired t-test was used to investigate statistical significance. A p-value < 0.05 was considered to be statistically significant.

First, we analyzed 30 whole transcriptome datasets to identify differentially expressed lncRNAs upon progesterone treatment. Of the 30 tumor samples, 10 had received a single 500 mg dose of hydroxyprogesterone and 20 were controls [37, 54]. Sequencing of these samples generated 17.2–60.7 million reads per sample (median, 37.4 million), wherein > 94–96% reads aligned to the human genome. Differential gene expression analysis between the control and progesterone-treated patients aided in identifying 2,222 differentially expressed genes (FDR < 0.1; 764 up- and 1,458 down-regulated), containing 537 lncRNAs (287 up- and 250 down-regulated), while a majority of the deregulated genes were of protein-coding category (Fig. 1A, Additional file 2: Table S2).

Further, to better understand the underlying mechanisms of action of progesterone, we performed whole transcriptome sequencing of T47-D (PR + /ER + /Her2-) and MDA-MB-231 (PR-/ER-/Her2-) breast cancer cells in response to progesterone treatment. A minimum of 60 million pair-end reads were obtained for each sample with > 90% of reads with a Phred score > 30, suggesting good quality of the data. Of 382 and 206 differentially expressed genes in T47-D and MDA-MB-231 cells, respectively, 18 lncRNAs were significantly deregulated in response to progesterone (-1 < log2FC > 1; p-value < 0.05) (Fig. 1A; Additional file 1: Figure S1; Additional file 2: Tables S3, S4). MDA-MB-231, a PR-negative cell line, also showed active transcriptional response to progesterone treatment, likely due to the PR-independent mode of action of progesterone mediated by glucocorticoid receptor (GR) [40]. Interestingly, expression of a few lncRNAs was consistently deregulated in the progesterone-treated breast tumor and cell line transcriptome data, viz.., DSCAM-AS1, PCED1B-AS1, RP11-21L23.2, RP11-363E7.4, and AC012358.8 (Fig. 1A). Of these, expression of DSCAM-AS1 was considerably downregulated in progesterone-treated breast cancer patients transcriptome data. Moreover, the normalized DSCAM-AS1 expression in progesterone untreated samples range from 3–860, compared to 3–450 in progesterone treated samples (Additional file 1: Figure S2). This suggests the variable DSCAM-AS1 expression across progesterone treated and untreated primary breast tumor samples. Further, consistent with our previous study, SGK1 was found to be significantly upregulated [13, 29, 40], in addition to deregulated expression of some lncRNAs, in progesterone-treated breast cancer samples (Fig. 1A). Taken together, the transcriptome analyses of tumor and cell lines identified novel progesterone-responsive lncRNAs in breast cancer.

An orthologous validation of the differentially expressed lncRNAs by real-time PCR identified a long non-coding RNA Down syndrome cell adhesion molecule antisense, DSCAM-AS1, as downregulated in ER/PR-positive T47-D and BT-474 cells upon progesterone treatment compared to ER/PR-negative MDA-MB-231 cells (Fig. 1B–D). In contrast to T47D and BT474 cells, we detected no significant change in the expression of DSCAM-AS1 in MCF7 cells (Additional file 1: Figure S3), consistent with its distinct transcriptional landscape, as described earlier, in response to progesterone treatment [18, 55]. The downregulation of DSCAM-AS1 could be effectively blocked by mifepristone, an antagonist of progesterone receptor (PR) and glucocorticoid receptor (GR) (Fig. 2A, B). However, the siRNA-mediated knockdown of PR, but not GR, rescued the down-regulation of DSCAM-AS1 in response to progesterone treatment, suggesting that PR mediates the downregulation of DSCAM-AS1 in response to progesterone in PR-positive cells (Fig. 2C). Next, we tested whether DSCAM-AS1 affects the inhibition of migration and invasion ability of breast cancer cells in response to progesterone [13]. Interestingly, DSCAM-AS1 knockdown could mimic the effect of progesterone by inhibiting breast cancer cell migration and invasion comparable to the extent obtained following treatment of the PR-positive T47-D and BT474 cells with progesterone (Fig. 2D–G).

Estrogen receptor (ER) has previously been shown to regulate DSCAM-AS1 expression via binding near the promoter region [27]. Consistent with the literature, we observed a significantly higher expression of DSCAM-AS1 transcript in TCGA breast cancer patient samples and ER/PR-positive T47-D and BT-474 cells than in MDA-MB-231 cells (Fig. 3A, B, Additional file 1: Figure S4). We hypothesized that ER/PR could modulate the DSCAM-AS1 expression in response to progesterone by binding to its upstream regulatory or distal regions. We analyzed chromatin immunoprecipitation ChIP-sequencing data following progesterone treatment in PR-positive T47-D cells, as described earlier [18, 40]. We identified enrichment of PR, ER, and p300 binding peak upon progesterone treatment at the “region 3” regulatory sequence of DSCAM-AS1 (Additional file 1: Figure S5). This suggests that progesterone alters the binding occupancy of PR and ER near DSCAM-AS1. Surprisingly, siRNA-mediated knockdown of DSCAM-AS1 in turn led to a significant decrease in the expression of ESR1 transcript, comparable to progesterone treatment, suggesting a possible feedback mechanism by which DSCAM-AS1 regulates the expression of ESR1 in T47-D and BT-474 cells (Fig. 3C, D). In contrast, overexpression of DSCAM-AS1 in T47-D cells, but not MDA-MB-231 cells, led to overexpression of ESR1 (Fig. 3E–H), suggesting that progesterone reduces expression of DSCAM-AS1 that further suppresses expression of ESR1 to inhibit cell migration and invasion in PR-positive breast cancer cells.

LncRNAs are known to sponge miRNAs, and thus, reduce the availability of miRNAs for target gene suppression [56, 57]. We thus tested whether DSCAM-AS1 could sponge miRNAs targeting the 3’-UTR of ESR1. Using the DIANA-LncBase v2 database prediction module, we identified 167 miRNAs that could bind DSCAM-AS1 with miTG-score > 0.7 (Additional file 2: Table S5). Concomitantly, we identified 72 miRNAs predicted to target the 3’-UTR of ESR1 from the miRTarBase database (Additional file 2: Table S6), with 9 overlapping miRNAs, viz. miR-548x, miR-548aj, miR-335, miR-129, miR-4422, miR-3121, miR-193b, miR-130a, and miR-301a (Fig. 4A). A real-time PCR-based validation of these 9 miRNAs in response to progesterone or genetic knockdown of DSCAM-AS1 identified miR-130a as significantly upregulated in T47-D and BT-474 cells (Fig. 4B, Additional file 1: Figure S6). Interestingly, in BT474 cells, a greater number of miRNAs were downregulated in response to silencing DSCAM-AS1 than in response to progesterone treatment. This may be related to the de-repression of miRNAs upon silencing DSCAM-AS1, a miRNA sponge, as well as the activation and inhibition of various pathways in response to progesterone, such as the up-regulation of miR-129–2, which in turn regulates the expression of PR, as demonstrated before [29].

Next, to investigate the function, we ectopically expressed miR-130a in T47-D and BT-474 cells. The ectopic expression of miR-130a led a significant decrease in ESR1 transcript than in vector control (Fig. 4C, D), with a concomitant decrease in invasion and migration of T47-D and BT-474 cells. miR-130a overexpression could mimic progesterone treatment or DSCAM-AS1 knockdown in PR-positive T47-D and BT-474 cells (Fig. 4E–H). Furthermore, to test a direct interaction between miR-130a and DSCAM-AS1, DSCAM-AS1 cDNA was cloned downstream to the luciferase reporter gene. The findings revealed a decrease in normalized luciferase activity upon overexpression of wild-type miR-130a, but not with miR-130a construct with a mutated binding site. As a positive control, the 3’-UTR of ESR1 cloned downstream to the luciferase gene showed similar inhibition of luciferase activity (Fig. 4I). In contrast, overexpression of DSCAM-AS1 in T47-D, but not MDA-MB-231, cells led to a significant reduction in miR-130a levels than that in vector control (Fig. 4J, K, Additional file 1: Figure S7 A-C). Interestingly, miR-130a also showed a significant inverse correlation to DSCAM-AS1 and ESR1 expression in 752 TCGA breast cancer samples (Additional file 1: Figure S8 A-B). Taken together, these results validate the association between DSCAM-AS1 and miR-130a to maintain ESR1 levels in PR-positive breast cancer cells with a consistent inverse correlation of miR-130a with the expression of DSCAM-AS1 and ESR1 in the TCGA patient samples.

The prognostic value of DSCAM-AS1 and miR-130a expression in survival prediction was further tested in TCGA breast cancer datasets (n = 1062) generated by whole transcriptome sequencing to perform the Kaplan–Meier (KM) survival analysis. Patients in the datasets were divided into high- and low-expression classes by the median expression value of DSCAM-AS1 and miR-130a, and a log-rank test was performed for stratifying patients with different prognoses. The analysis showed a significantly better overall survival in patients with breast carcinoma with high miR-130a expression than those with low miR-130a expression (log-rank p = 0.02). Patients with high expression of miR-130a survived better (87 months) than those with low expression of miR-130a (69 months). Overall, we observed a survival benefit of 18 months in the miR-130a high expression cohort. Similar results were observed in patients with ER-positive subtype cancer (log-rank p = 0.05) (Fig. 5A, B). In contrast, KM analysis of patients with breast cancer did not show statistically significant change in overall survival in patients who exhibit high and low levels of DSCAM-AS1 (Fig. 5C, D). These findings imply that a high expression of miR-130a influence survival of patients with breast cancer.