Introduction

BRCA1 is one of the most important tumor suppressor genes. Carriers of BRCA1 germline mutations are at greatly increased risks for developing various tumors, including breast, ovarian, prostate, colon, pancreatic, and gastric cancers (1, 2).Notably, BRCA1 is most closely associated with breast cancer. However, why BRCA1-associated tumors exhibit significant tissue specificity is unknown, and the effects of BRCA1 defects in normal human mammary cells are largely obscure.

In addition to the tissue specificity of tumors caused by BRCA1 deficiency, breast cancers associated with BRCA1 deficiency exhibit a set of other distinctive features. Most BRCA1 mutation carriers develop triple-negative breast cancers (TNBCs), in which tumors show negative expression of estrogen receptor (ER), progesterone receptor (PR), and Her2 (human epidermal growth factor receptor 2), and these cancers are usually also classified as basal-like breast cancers. The currently recognized cells of origin of BRCA1-deficient breast cancers are the luminal progenitor cell (LPs) in the mammary gland (3-6). LPs show abnormal expansion in the normal mammary tissue of BRCA1 mutation carriers (3, 7), which might be caused by a differentiation blockade (8). Investigating the traits of LPs is important for tumor prevention and treatment in BRCA1 mutation carriers and a wider range of patients with TNBCs and basal-like breast cancer. However, the reason LPs are prone to transform into cancer cells when BRCA1 is deficient is unclear.

ELF3 is a member of the ETS family of transcription factors. Early reports suggested that ELF3 may regulate the development and differentiation of mammary tissue (9). Intriguingly, ELF3 has opposite functions in tumorigenesis in different tissues. In breast cancers, ELF3 shows high expression (10), plays an oncogenic role by promoting EMT (11) and forms a positive feedback loop with Her2 to maintain the transformation phenotype of Her2+ breast cancers (12). However, the mechanism of high ELF3 expression in certain types of breast cancer and its association with breast cancers beyond the Her2+ subtype is currently unknown.

In this study, utilizing single-cell sequencing of BRCA1 mutant breast cancers and time-resolved RNA-seq of BRCA1-deficient normal mammary cells, we revealed that LPs have a trait of higher replication stress compared to other cell populations in the normal mammary tissue, endowing them with the potential for transformation when BRCA1 is deficient, and ELF3 plays a key role in the process. Mechanistically,

ELF3 expression is significantly upregulated by replication stress and BRCA1 deficiency, dependent on the ATR-Chk1-E2F axis and the transcriptional regulation function of BRCA1, respectively. In addition, ELF3 can help suppress excessive genomic instability, promoting the tolerance of LPs to BRCA1 deficiency. Notably, BRCA1 deficiency causes cells to display LP transcriptional profile characteristics, leading to dedifferentiation and expansion of LPs, and ELF3 is a core transcription factor during this process. Our study reveals why BRCA1 deficiency is prone to result in tumorigenesis in LPs and elucidates the key role of replication stress and ELF3 during this process, suggesting that ELF3 could be a promising target for BRCA1-associated breast cancers.

Results

1. LPs have higher levels of replication stress compared with other mammary cell populations

LPs in the mammary gland are the currently recognized cells of origin of BRCA1-deficient breast cancers. However, the reason for LPs are prone to transform with BRCA1 deficiency has not been elucidated. To address this question, we attempt to investigate from two aspects that may influence the process of tumorigenesis. Firstly, we considered that LPs, which are stimulated by upstream progesterone paracrine signals through RANKL (5), may have higher levels of replication stress, which can trigger tumorigenesis by increasing genomic instability. We utilized both published data and our published single-cell sequencing data (8) to address this question. We analyzed the data of Lim et al. (3) by scoring multiple DNA replication pathway gene expression to infer the replication stress in different mammary cell populations as previously studied (13). The result showed that LPs have the highest level of replication stress among all subgroups of normal human mammary cells (Fig. 1A). Our single-cell sequencing data were obtained from the tumor cells and normal mammary cells of four BRCA1 mutation carriers (case #1-4) and normal mammary cells of three matched non-carriers (case #5-7). According to our previous analysis, case #3 BRCA1 mutant carrier developed basal-like breast cancer, which was derived from LPs. Importantly, in this case of BRCA1 mutant carrier, single-cell sequencing data revealed that tumor cells had the highest levels of replication stress, and LPs had the second-highest levels (Fig. 1B). Moreover, we applied Monocle to reconstruct the transformation trajectories from normal mammary cells to tumor cells in the BRCA1 mutation carrier and found that the levels of replication stress showed an increasing trend during the process of LP transformation into tumor cells (Fig. 1C and Supplemental Fig. S1A). These data demonstrate that replication stress is an important trait of LPs and might be a key driver of BRCA1-deficient tumorigenesis, which is consistent with the previous report that BRCA1 mutant mammary epithelial cells are defective in replication stress suppression (14).

Figure 1.

Next, we aimed to find other key regulators during BRCA1-deficient tumorigenesis in LPs. In addition to inducing replication stress, BRCA1 deficiency also profoundly affects the transcriptional profile of normal mammary cells by disrupting the transcriptional function of BRCA1, which may further boost the process of tumorigenesis. To investigate the transcriptional impact of BRCA1 deficiency in normal human mammary epithelial cells, we took advantage of the Tet-on system to build an inducible knockdown model in MCF10A cells. By adjusting the treatment condition of DOX (doxycycline), different time points of BRCA1 deficiency can be controlled. Treatment with DOX causes a constant and efficient knockdown of BRCA1, verified at mRNA and protein levels (Supplemental Fig. S1B). We showed that DOX treatment induces significant growth arrest (Supplemental Fig. S1C) and an increase in genomic instability (Supplemental Fig. S1D), consistent with the phenotypes of BRCA1-deficient cells (15), confirming that our system functioned as expected. We then performed time-resolved RNA-seq with different durations of DOX treatment (2, 5 and 10 days; Supplemental Fig. S1E). Notably, 10 days of DOX treatment resulted in growth arrest, thus cells collected for sequencing have gained growth advantage. Differential gene expression analysis revealed that BRCA1 was efficiently knocked down during the whole experiment (Fig. 1D). Importantly, the principal component analysis (PCA) revealed that samples showed a time-course distribution in the dimension of PC1, which has the dominant power of explanation (Fig. 1E and Supplemental Fig. S1F). These data demonstrate that the variations between samples can be largely attributed to time. To further investigate the modes of differentially expressed genes (DEGs) that change over time, we carried out the DEG cluster analysis using Mfuzz (16). Analysis of the Day 10 DEGs showed that genes could be divided into 5 clusters with different expression change patterns (Fig. 1F). We mainly focused on Cluster 1, for these sharply upregulated genes may play more important roles in cells gaining growth advantage with BRCA1 deficiency at Day 10, mimicking the initiation of BRCA1 deficiency associated tumorigenesis. To find out the most potent regulator upregulated both by replication stress and BRCA1 deficiency, we performed RNA-seq in HU-treated MCF10A cells and overlapped the upregulated genes in HU RNA-seq and Cluster1 genes in BRCA1 deficiency RNA-seq (Fig. 1G). The results showed that ELF3, a transcription factor, is one of the most dominant genes in this analysis. ELF3 has been reported to be involved in functions including EMT in human breast cancer and promotes the transformation of normal human mammary cells (11, 17). Therefore, we speculated that ELF3 could be an important driving factor during the tumorigenesis of BRCA1-deficient and basal-like breast cancer.

2. ELF3 is upregulated in BRCA1-associated breast cancer and is related to a worse prognosis

We next explored the expression of ELF3 in breast cancers in the TCGA and METABRIC databases to further investigate the association between ELF3 and BRCA1 and validate our previous analysis. In both clinical databases, ELF3 expression is significantly higher in BRCA1-associated breast tumors than in non-BRCA1-associated breast tumors (Fig. 2A). Moreover, the expression levels of BRCA1 and ELF3 showed a significant negative correlation, consistent with our findings (Fig. 2B). Since most BRCA1 mutation carriers develop basal-like breast cancers, we next investigated the association between ELF3 and basal-like breast cancers. We exploited the TCGA and METABRIC databases to compare ELF3 expression levels among different breast cancer subtypes. As we expected, ELF3 showed the highest expression levels in basal-like breast cancer in both databases (Fig. 2C). We took advantage of novel subgroups identified by the METABRIC (18) and further explored ELF3 expression in ten IntClusts of the METABRIC database. Among all subgroups, IntClust 5 was enriched with Her2+ cancers, and IntClust 10 was enriched with basal-like cancers. These two subgroups showed the highest expression levels of ELF3, with no significant difference between them (Fig. 2D). The higher levels of ELF3 expression in Her2 subtype breast cancer are consistent with previous studies showing that ELF3 may be a downstream transcription target of Her2 and that there might be a positive feedback loop between them (12, 19, 20). On the other hand, IntClust 4 shows a favorable outcome and has a ‘CNA-devoid’ feature, which has the lowest level of ELF3 expression. These data are consistent with our findings that ELF3 is upregulated by BRCA1 loss, indicate that ELF3 may play an important role in basal-like breast tumorigenesis, and moreover, the expression of ELF3 may be related to genomic stability status. Furthermore, among human breast cancer tissue samples, ELF3 expression levels were significantly higher in TNBCs than in other subtypes (Fig. 2E and F). Additionally, TNBCs have higher ELF3 expression levels than the Her2 subtype (Fig. 2E and F), revealing a tighter connection between ELF3 and BRCA1-associated breast cancers.

Figure 2.

ELF3 shows significant copy number amplification in breast cancer (21). Thus, we wondered that whether the high expression levels of ELF3 in BRCA1-associated breast cancers are due to copy number variation. Therefore, we analyzed the proportion of ELF3 copy number amplification in different subtypes (Supplement Fig. S2) and found that this was not the case. In TCGA basal-like breast cancer, the ELF3 amplification proportion was not higher than that in other subtypes but was only 3.5%, which was lower than the total level (7.3%). These data demonstrate that the upregulation of ELF3 in BRCA1-associated breast cancer is not caused by copy number amplification.

Furthermore, we investigated whether ELF3 expression could be a biomarker for breast cancer prognosis. Consistent with our hypothesis, higher ELF3 expression is correlated with worse prognosis in TNBC and basal-like breast cancer (Fig. 2G). In the Her2 subtype, higher ELF3 expression is also related to shorter overall survival, consistent with previous reports that as a Her2 target gene, ELF3 may help promote tumor invasion (20). Interestingly, ELF3 expression shows less or no prognostic prediction value in luminal A and luminal B subtypes, suggesting that ELF3 plays a more powerful role in TNBC and basal-like breast cancer. These data imply that in BRCA1-associated breast cancer, ELF3 tends to be upregulated and plays a “driver” role, promoting tumor progression and leading to a worse prognosis.

3. ELF3 upregulation is induced by replication stress via the ATR-Chk1-E2F axis and by BRCA1 deficiency via GATA3 transcription

We next focused on the mechanisms of ELF3 upregulation in the presence of replication stress and BRCA1 deficiency. Consistent with the HU RNA-seq analysis, our single-cell RNA-seq data showed that in both BRCA1 WT and mutant carriers, the levels of replication stress and ELF3 expression showed a significant positive correlation (Fig. 3A). ELF3 was also upregulated by HU-induced replication stress in MCF10A cells (Fig 3B). Since ATR is the key kinase that responds to replication stress, we treated cells with ATRi and found that it could efficiently block ELF3 upregulation induced by HU, while ATMi had no such effect (Fig. 3B). Moreover, downstream of ATR, ELF3 expression upregulation depended on Chk1 activation but not Chk2 activation (Fig. 3C). It has been reported that the transcriptional activation induced by replication stress-ATR-Chk1 pathway largely depends on E2F transcription factor activation (22). Under replication stress, ATR-Chk1-mediated phosphorylation of E2F6 blocks its inhibitory function on E2F and activates E2F transcription. Thus, we overexpressed E2F6 and found that E2F6 blocks the upregulation of ELF3 under HU treatment conditions (Fig. 3D), validating that ELF3 upregulation relies on E2F6 inactivation and E2F de-suppression. Moreover, utilizing ENCODE ChIP-seq data, we found that the ELF3 promoter is indeed bound by E2F1 (Fig. 3E). These data demonstrate that ELF3 is upregulated by the ATR-Chk1-E2F axis under replication stress conditions.

Figure 3.

Then we investigate the mechanism of ELF3 upregulation induced by BRCA1 deficiency. ELF3 expression gradually increased during 10 days of DOX treatment (Fig. 3F). ELF3 upregulation was further confirmed by siRNA-mediated BRCA1 knockdown in MCF10A cells to strictly exclude the unknown influence of the Tet-on system (Fig. 3G). BRCA1 has transcription regulation functions by interacting with other transcription factors to bind to gene promoters and regulate gene expression (23). Therefore, we explored transcription factors associated with BRCA1, including GATA3 (24) and Myc (25), for whether the binding motifs of these transcription factors appear in the ELF3 promoter sequence utilizing the JASPAR database (Supplement Fig. S3A). We found that E2F1 was present in the prediction list with a high relative score, consistent with our results. The prediction results show that the motif of GATA3 appeared in the ELF3 promoter multiple times with an overall high relative score. Indeed, knockdown of GATA3 in MCF10A cells also strongly induced ELF3 expression (Fig. 3H). Similarly, we found multiple GATA3 binding peaks at the promoter of ELF3 in two independent ChIP-seq datasets from the ENCODE database, consistent with two of the GATA3 binding sequences predicted in the JASPAR database (Fig. 3I). These data suggest that ELF3 is a GATA3 downstream target gene. Importantly, when both BRCA1 and GATA3 were knocked down, there was no significant superimposed effect on ELF3 upregulation (Fig. 3J and Supplement Fig. S3B), implying that BRCA1 and GATA3 act in the same pathway for ELF3 induction.

Notably, we found that and BRCA1 knockdown and HU treatment have a superimposed effect on upregulating ELF3 expression levels (Fig. 3K and Supplement Fig. S3C). Furthermore, GATA3 knockdown and HU treatment also had a significant superimposed effect in upregulating ELF3 expression (Fig. 3L and Supplement Fig. S3D). These data validate that the BRCA1-GATA3 axis and replication stress regulates ELF3 expression independently.

4. ELF3 helps suppress excessive genomic instability

Next, we investigated the impact of ELF3 upregulation under replication stress conditions and BRCA1 deficiency. Based on the association between ELF3 and the worse prognosis in BRCA1-associated breast cancer patients (Fig. 2G), we hypothesized that upregulation of ELF3 could help cells deal with replication stress and suppress excessive genomic instability, facilitating cancer initiation and cancer evolution. As expected, ELF3 knockdown greatly increased the sensitivity of MCF10A cells to HU and cisplatin (Fig. 4A and B, Supplement Fig. S4A). This finding suggests that ELF3 can suppress genomic instability and help normal cells deal with replication stress. Furthermore, ELF3 knockdown in HCC1937 and SUM149PT cells, both of which are BRCA1 mutant human breast cancer cell lines, also caused a significant decline in cell proliferation (Fig. 4C and D, Supplement Fig. S4B and C), suggesting that ELF3 is essential for BRCA1 mutant breast cancer cell survival. Moreover, HCC1937 and SUM149PT cells became more sensitive to HU and cisplatin when ELF3 was deficient (Fig. 4E and F, Supplemental Fig. S4D and E). Since PARP inhibitors can induce replication fork stalling by PARP1 trapping and lead to replication stress in BRCA1-deficient cells, we also tested olaparib sensitivity in HCC1937 cells and obtained similar results (Supplement Fig. S4F). These results demonstrate that in BRCA1 mutant breast cancer cells, ELF3 can help suppress excessive genomic instability.

Figure 4.

In addition, utilizing the CellMiner database, we found that among nearly sixty diverse human cancer cell lines, lower ELF3 expression is correlated with higher drug sensitivity to cisplatin and etoposide (Fig. 4G). Consistent results were obtained using the GDSC database (Fig. 4H, Supplement Fig. 4G). Moreover, in nude mice, SUM149PT tumors with ELF3 deficiency were significantly smaller than those in the control group (Fig. 4I). These data indicate that ELF3 can help suppress excessive genome instability in the presence of replication stress, helping normal cells tolerate replication stress and subsequent DNA damage, allowing the cells to transform into malignant tumors. During cancer progression, BRCA1 mutant cancer cells gradually become reliant on high ELF3 expression to preserve a certain degree of genomic stability, making ELF3 a promising therapeutic target and biomarker in BRCA1-associated breast cancers.

5. ELF3 helps maintain the stability of DNA replication

We then investigated the mechanism by which ELF3 preserves genomic stability. ELF3 deficiency leads to more endogenous DNA damage in MCF10A cells (Supplement Fig. S5). In HCC1937 cells, knockdown of ELF3 resulted in a significant increase in spontaneous γH2AX foci (Fig. 5A), indicating that ELF3 deficiency leads to more endogenous DNA damage in BRCA1-deficient breast cancer cells. Moreover, with HU treatment and release for repair, ELF3 deficiency resulted in a significant increase in the number of γH2AX and 53BP1 foci per cell and the proportion of cells with damage foci (Fig. 5B-D). Under replication stress conditions, γH2AX and 53BP1 foci represent DNA double-strand breaks formed by collapsed replication forks (26). These data suggest that ELF3 deficiency results in a defective ability to repair stalled replication forks; thus, more DNA damage is generated. The resulting upregulation of ELF3 expression during replication stress helps cells better manage stalled replication forks and prevents excessive DNA damage.

Figure 5.

To further explore how ELF3 contributes to replication stress management, RNA-seq on ELF3-knockdown MCF10A cells was performed. Gene set enrichment analysis (GSEA) showed that multiple pathways associated with DNA replication were significantly downregulated in ELF3-deficient cells (Fig. 5E), which was further verified by RT–qPCR (Fig. 5F). These data indicate that DNA replication pathways are affected by ELF3 deficiency.

DNA fiber assay was performed to more closely investigate the effects of ELF3 deficiency on DNA replication. ELF3 deficiency significantly slowed down the DNA replication progression (Fig. 5G), with significantly reduced symmetry of replication forks (Fig. 5H), suggesting that more replication forks were stalled or unstable. Moreover, with HU treatment, ELF3 knockdown undermined the stability of stalled replication forks (Fig. 5I). Notably, double-knockdown of ELF3 and BRCA1 led to a further decrease in replication fork stability (Fig. 5I). These data reveal that ELF3 expression upregulation when BRCA1 is deficient can help cells deal with replication stress, while simultaneous ELF3 deficiency leads to intolerable levels of genomic stability.

6. ELF3 is a core transcription factor inducing LP genes expression under BRCA1 deficiency

Finally, we explored the possible reasons for the tissue specificity of BRCA1-associated tumors. It has been reported that BRCA1 deficiency disturbs the development and differentiation of mammary tissue and causes abnormal expansion of the LP population (3, 8, 27). We hypothesized that BRCA1 deficiency might affect the normal differentiation process of LPs by disturbing the transcriptional profile, thus leading to abnormal expansion and transformation of LPs. We analyzed the time-course RNA-seq data using a published signature gene set from LPs (3) and found that over time with BRCA1 deficiency, normal mammary cells gradually show a higher signature score of LP gene expression (Fig. 6A). Similarly, using the sets of genes up-and downregulated in the LP population compared to in other mammary cell populations obtained from the data of Pellacani et al. (28), we found that genes upregulated in LPs also tended to be upregulated in BRCA1-deficient conditions, whereas genes downregulated in LPs tended to be downregulated (Fig. 6B). These results suggest that BRCA1 may broadly affect the transcription profile of mammary cells, and therefore, when BRCA1 is deficient, LPs may tend to maintain a dedifferentiated state, leading to the abnormal expansion and tumorigenic properties of these cells in BRCA1 mutant carriers.

Figure 6.

Multiple ETS family transcription factors were shown to be involved in human mammary development and differentiation in previous reports. However, the function of ELF3 in LPs has not yet been elucidated. We used GSEA to explore the transcriptional profiles of BRCA1-deficient mammary cells, ELF3-overexpressing mammary cells and human LPs from previous reports (28). Surprisingly, these three transcriptional profiles have significant correlations (Fig. 6C), indicating that ELF3 is not only a core transcription node under BRCA1 deficiency conditions but also, more importantly, that ELF3 plays significant roles in the transcriptional regulation of human LPs. In addition, in our single-cell sequencing data, the expression of ELF3 in LPs is the highest among all mammary cell populations (Fig. 6D and Supplement Fig. S6A). Same results were obtained analyzing the published data from Pellacani et al. (28) (Supplement Fig. S6B). To investigate ELF3 functions in regulating LP genes transcription, we used SEA (simple enrichment analysis) of MEME to perform motif enrichment of ESE and ELF subfamilies transcription factors in promoters of LP genes, as these transcription factors are more closely related to ELF3 in the ETS family, and several, such as ELF5 and EHF, have been reported to regulate LPs (28-30). ELF3 was at the top of the list (Fig. 6E), demonstrating that ELF3 has higher affinity for the promoters of LP genes. Moreover, our single-cell data revealed that ELF3 expression levels were the highest in LPs compared to other members of the ESE and ELF subfamilies (Fig. 6F), including ELF5, which has been reported to play vital roles in LP differentiation. In addition, we analyzed the ATAC-seq data in Brca1-deficient mouse LPs from Bach et al. (31). The ATAC-seq compared the chromatin accessibility of LPs in Brca1/p53-deficient mice with WT mice, and motif enrichment analysis of transcription factors was performed in regions with increased accessibility of Brca1-deficient mouse LPs. ELF3 was among the top genes (Supplement Fig. S6C) and is significantly enriched among all ETS transcription factors (Fig. 6G). These data further confirm our hypothesis that ELF3 plays an important role in the transcriptional regulation of LPs.

Discussion

In this study, we revealed that LPs have the trait of higher replication stress, an important factor for their potential to transformation. The difference of LPs from other cell populations in normal mammary tissue has received considerable attention in the field. Research in mouse mammary tissue has revealed that the overall luminal epithelial cell population is enriched in mitotic, cell cycle and DNA repair-related pathways under estrogen-and progestogen-stimulated conditions relative to the basal and stromal cell populations (32). Subsequently, Casey et al. (33) found that LPs of mouse mammary tissue are enriched in DNA repair and cell division pathways. A recent preprint from the same group (34) continued the analysis of mouse tissue sequencing data and found that in LPs, DNA repair and replication pathways are significantly upregulated, indicating that LPs are more capable of DNA damage repair and thus more tolerant to BRCA1 deficiency. Furthermore, Nolan E et al. (5) found that RANK+ LPs are highly proliferative and prone to DNA damage. These studies did not reveal the reason for LP tendency to transform. In our study, we focus on the characteristic of LPs compared to other cell populations and using sequencing data from normal human mammary tissues (3), we reveal that the DNA replication pathways and the cellular response to replication stress pathways are significantly upregulated in LPs compared to those in other cell populations. We argue that high levels of replication stress may be the major source of DNA damage in LPs, resulting in DDR gene upregulation in these cells, and the high proliferation feature of LP could be the main reason for the high replication stress.

Replication stress is a hallmark of cancer, fueling cancer evolution by increasing genomic instability (35). LPs have inherently higher replication stress levels, and when BRCA1 is deficient, replication stress is further elevated during the transformation from normal LPs to tumor cells, suggesting that replication stress is an important driver of BRCA1-deficient tumorigenesis and the suppression of replication stress could be a new approach for preventive intervention in BRCA1 mutation carriers with the advantages of being less damaging than prophylactic resection. Moreover, we found that both replication stress and BRCA1 deficiency itself can induce high expression of ELF3, which also plays an important role in the tendency of LPs to transform with BRCA1 deficiency. First, ELF3 can help suppress excessive genomic instability, which is crucial for cell tolerance and survival with BRCA1 deficiency. We speculate that high ELF3 expression can help LPs maintain moderate levels of genomic instability that allow cells to proliferate and is sufficient to fuel cancer evolution, including subsequent events such as p53 mutation and BRCA1 loss of heterozygosity (LOH). Second, ELF3 upregulation can boost LP gene transcription, leading to LP dedifferentiation and expansion. Finally, as an EMT promoter in mammary tissue, ELF3 can accelerate the transformation of LPs to malignant cells. Therefore, high levels of replication stress and ELF3 expression are important factors resulting in the high tendency for LPs to transform (Fig. 7). A key feature of BRCA1-associated tumors is the tissue specificity, the mechanism of which has long been a research focus. The transcriptional regulation functions of BRCA1 are critical for the development and differentiation of normal mammary tissue. LPs with abnormal expansion have been observed in BRCA1 mutation carriers. Several important differentiation pathways, such as WNT, Notch and Hedgehog, were downregulated in these abnormal LPs (8), suggesting that BRCA1 deficiency may lead to a differentiation blockade in LPs and, therefore, abnormal expansion. Our RNA-seq data revealed that the transcriptional profile of BRCA1-deficient MCF10A cells has a striking similarity to that of LPs in normal human mammary tissue, and ELF3 might be a core mediating transcription factor. This finding suggests that BRCA1 has a broad regulatory effect on gene expression in normal mammary cells and has specific functions in LPs; thus, its deficiency causes LPs to maintain a dedifferentiated state, leading to abnormal expansion of this cell population, and ELF3 may be a core transcription factor inducing these changes. In addition, our study reveals another layer of the mechanism underlying the tissue specificity of BRCA1-associated tumors—the capability of BRCA1 to regulate ELF3 expression. ELF3 is an epithelial-specific transcription factor, but the functions of ELF3 in different epithelial tissues vary. For example, ELF3 promotes tumor development in breast, prostate, and lung cancers (10, 36, 37), whereas it suppresses tumorigenesis in bladder and oral squamous cancers (38, 39). In addition, ELF3 itself plays an important role in the development and differentiation of normal mammary tissue; thus, dysregulation of ELF3 may disturb these processes. Therefore, we suggest that replication stress caused by BRCA1 deficiency leads to upregulation of ELF3 expression in LPs, thus promoting the tissue-specific oncogenic functions of ELF3 and breast tumorigenesis.

Figure 7.

This study also has some limitations. First, we used the MCF10A cell line, which is widely used for investigations into normal human mammary epithelial cells with the advantage of being easy to manipulate and modify, and compared to patient-derived cells, unknown genetic backgrounds in addition to BRCA1 mutations are avoided. However, the transcriptional profile and epigenetic characteristics of this cell line are distinct from those of normal human mammary cell populations and LPs (28). Although our data show that MCF10A cells can reflect the actual functions of BRCA1 in terms of maintaining genomic stability (Fig. 1B and C) and our RNA-seq data demonstrate enrichment in the known transcriptional target genes of BRCA1, the effect of BRCA1 defects in LPs needs to be further verified in primary human mammary cell populations. Second, the higher replication stress levels in LPs may lead to altered expression or modifications of molecules other than ELF3, affecting the development and progression of BRCA1-associated breast cancer, which also needs to be further explored. In addition, there are likely other DEGs in our BRCA1-deficient RNA-seq that play important roles in the progression of BRCA1-associated breast cancer, and these genes also have potential research value.

Taken together, we revealed that LPs have a trait of higher replication stress, endowing them with the potential for transformation with BRCA1 deficiency. Moreover, the levels of replication stress gradually increase during tumorigenesis, indicating that replication stress is a driving factor of BRCA1-associated malignant transformation. Mechanistically, ELF3 is significantly upregulated upon replication stress and BRCA1 deficiency, which can help suppress excessive genomic instability by stabilizing the DNA replication process, empowering LPs to tolerate and to transform with BRCA1 deficiency. Notably, BRCA1 deficiency causes a tendency for cells to display characteristics of the LP transcriptional profile, and ELF3 is a core transcription factor mediating this process. Our study reveals why BRCA1 deficiency is prone to result in tumorigenesis in LPs and elucidates the key role of replication stress and ELF3 during this process, and suggests promising targets for BRCA1-associated breast cancers.

Materials and Methods

Cell culture

MCF10A cells were cultured in DMEM/F12 supplemented with 5% horse serum, 1% penicillin/streptomycin, 20 ng/mL EGF, 10 μg/mL insulin, 0.5 μg /mL hydrocortisone and 100 ng/mL cholera toxin. For MCF10A-shBRCA1-Tet-on cells, except for the replacement of horse serum with 8% Tet-on system approved certified fetal bovine serum, all other components remain the same. HEK293T and SUM149PT cells were cultured in DMEM supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin. HCC1937 cells were cultured in RPMI1640 supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin.

Reagents and plasmids

Reagents (antibodies and chemicals) used in this study were listed in Table S1. For the construction of Tet-on-shBRCA1-hygro plasmid, shBRCA1 sequence was subcloned into the Tet-on-pLKO-hygro vector (SinaSun). For ELF3 knockdown by lentivirus in SUM149PT, shELF3 and shCtrl sequence was subcloned into the pLKO.1-puro vector. The shRNA sequences used in this study are listed in Table S2. For ELF3 overexpression by lentivirus in MCF10A, human ELF3 was subcloned into the pCDH-CMV-MCS-EF1-puro vector, which was a kind gift from Yin lab. Myc-tagged human E2F6 was subcloned into the pDEST-N-Myc vector. PCR primers used in this study were listed in Table S3.

Transfection of siRNA and infection of shRNA

All the siRNAs were from GenePharma. For siRNA transfection, cells were seeded in 6-well plate or 35 mm² dish at the day before transfection. 100 μL OPTI-MEM, 40 pmol siRNA and 3 μL Lipofectamine RNAiMAX reagent were mixed and added to the cells. After 48-72 h, cells were collected for subsequent experiments. The siRNA sequences used are listed in Table S3. For shRNA infection, lentivirus was generated by co-transfection of the lentiviral vector (Tet-pLKO.1-shBRCA1, pLKO.1-shCtrl and pLKO.1-shELF3) with envelope plasmids pMD2.G and psPAX2 in HEK293T cells. Lentivirus was collected after 48 h. Cells were infected with lentivirus with polybrene. After 24 h, cells were re-infected and harvested at 72 h for subsequent experiments.

Quantitative PCR (RT-qPCR)

Total RNA was extracted using Trizol reagent (Invitrogen) according to the manufacturer’s guidelines. cDNA was synthesized by reverse transcription PCR using Hifair III 1st Strand cDNA Synthesis SuperMix kit (Yeasen). qPCR was performed using SYBR Green Mix (Yeasen) on ABI StepOnePlus qPCR instrument (ThermoFisher). Each group was performed in triplicate. The relative expression mRNA levels were calculated by relative quantitative method (ΔΔCt). The primers used are listed in Table S5.

Immunofluorescence

Cells were seeded on coverslips and treated as indicated. Cells were than fixed in cold 4% paraformaldehyde for 15 min and permeabilized in 0.25% Triton X-100 solution for 5 min at room temperature. After blocking with 2% BSA for 10 min, the primary antibodies were added onto the coverslips, and then incubated for 1 h at room temperature, followed by secondary antibody incubation for 30 min. Finally, the coverslips were stained with DAPI for 5 min, and images were acquired using the ZEN software.

Immunohistochemistry

All tissue microarray chips were purchased from Shanghai Outdo Biotech Co.Ltd. IHC was performed as described previously, except that the following antibodies were used. The anti-ELF3 rabbit antibody were purchased from SIGMA. After incubation with primary antibody, the HRP-conjugated rabbit secondary antibody (ZSGB-Bio) was applied and IHC slides were observed with Olympus BX51 microscope and Olympus DP73 CCD photographic system.

Colony formation assay

MCF10A-shBRCA1-Tet-on cells were seeded in 6 mm² plates. 1 μg/mL DOX was added for the DOX group during the whole process of cell growth to knockdown BRCA1 expression. After 14 days, the cells were washed with phosphate-buffered saline (PBS) and stained with 0.1% Coomassie brilliant blue in 10% ethanol for 30 minutes at room temperature. The stained dishes were washed with water, and the colonies were counted.

DNA fiber

The DNA fiber assay was performed as described previously with slight modifications (40). Briefly, cells transfected with siRNA were labeled by 25 μM IdU (Yuanye Bio-Technology) for 20 min and subsequently 200 μM CldU (Abcam) for 20 min. Then cells were treated with 4 mM HU for 4 h for replication fork stability assay or collected for replication fork speed and symmetry assay. Cells were resuspended in lysis buffer (200 mM Tris-HCl, pH 7.5, 50 mM EDTA and 0.5% SDS) to extract DNA. Cell lysates were then dripped onto a glass slide and DNA fibers were stretched along the slide after 3 minutes incubation. The slides were air-dried and fixed in 3:1 methanol/acetic acid solution for 10 min, followed by 2.5 M HCl treatment for 1 h, then blocked in 3% BSA PBST buffer for 20 min. DNA fibers were incubated with anti-BrdU antibody (Abcam) binding CldU and anti-BrdU antibody (BD) binding IdU at 37°C for 1.5 h. After washing with PBST, the slides were incubated with Alexa Fluor 488- and Alexa Fluor 594-conjugated secondary antibody (Invitrogen) for 45 minutes. Finally, the slides were mounted in mounting reagent (Solarbio) and DNA fibers images were acquired using ZEN software.

Comet assay

The comet assay was performed as described previously with slight modifications (41). Briefly, cells were collected by trypsin digestion and resuspended in PBS. 0.5% normal melting agarose (NMA) was added to glass slides and quickly covered with coverslips. The slides were placed at 4 °C for 10 min. The coverslips were removed carefully. The cell suspensions were mixed gently with equal volume of 1% low melting agarose (LMA) and quickly added onto the solidified NMA, quickly covered with coverslips. The slides were placed at 4 °C for 10 min. After removing the coverslips, the slides were incubated overnight in the fresh cell lysis buffer at 4 °C overnight. Then the slides were rinsed with ddH₂O and electrophoresis at 4°C for 30 minutes at 20 V in neutral electrophoresis buffer. Then the slides were stained with 5 μg/ml of propidium iodide for 10 min. The comet was observed using a fluorescence microscope (Zeiss AxioCam 503 color). The pictures were analyzed using CASP software (http://casplab.com/).

Cell survival (CCK-8) assay

Cell survival assay was performed using CCK-8 kits (meilunbio, MA0218-3) according to the manufacturer’s protocols. Cells were seeded in 96-well plates. After 12 h, medium containing indicated concentration of drugs was changed and cells were incubated for 3-5 days. DMEM containing 0.1% DMSO was used as the control group. CCK-8 solution was diluted in 10% concentration and added to wells, and the cells were incubated for 45 min-1.5 h. DMEM containing 10% CCK-8 was used as the blank group. The absorbance at 450 nm was detected by a microplate reader. Cell survival was calculated according to the formula: (ODdrug – ODblank) /(ODcontrol – ODblank).

RNA sequencing

Total RNA of cells was extracted with Trizol reagent. RNA sequencing was performed by Novogene (China). The resulting data was mapped to human reference genome (hg38) and gene counts were calculated for each sample using STAR (v2.6.1) (42) with the parameter “--quantMode GeneCounts”. To identify differentially expressed genes (DEGs), genes with low expression levels (average counts per million reads < 0.05) were filtered out. DESeq2 was employed to detect DEGs with |log2 fold change| > 2 and FDR < 0.001. The raw RNA sequencing data have been deposited into GEO under the accession number GSExxxxxx. PCA analysis was performed using the “prcomp” function of stats (v3.6.1) R package.

Analysis of TCGA breast cancer and METABRIC data

TCGA (https://www.cancer.gov/tcga) breast cancer and METABRIC (18) data were downloaded using the cgdsr R package. For TCGA data, RSEM (batch normalized from Illumina HiSeq_RNASeqV2) was used as the expression levels of the genes. For METABRIC data, mRNA expression log intensity levels (Illumina Human v3 microarray) were used as the expression levels of the genes. BRCA1 associated cancers were identified as cancers with BRCA1 mutations, heterozygous loss or homozygous deletion. Pearson correlation coefficient between ELF3 expression and BRCA1 expression was calculated. Subtypes and IntClust types of breast cancers were identified according to the clinical data of TCGA and METABRIC databases.

Chip-Seq analysis of ENCODE databases

The Chip-Seq data were downloaded from the ENCODE database (https://www.encodeproject.org/) (43, 44) with the following identifiers: ENCFF858GLM, ENCFF692OYJ, ENCFF384CPN, ENCFF342GNN, ENCFF352QVM, and ENCFF437NQS. Images were generated using IGV (v2.7.2).

Drug sensitivity analysis of Cell Miner and GDSC data

Cell line drug sensitivity data and mRNA expression data were downloaded from the CellMiner (45, 46) and GDSC (47) databases. Pearson correlation coefficient between ELF3 expression and drug sensitivity was calculated for the indicated drugs.

Replication stress score and luminal progenitor score analysis

For replication stress score calculation, human mammary subgroups gene expression matrixes were obtained from the datasets of Lim et al. (3) (GSE16997). Replication stress associated pathways were selected referring to Dreyer et al. (13). and gene lists were downloaded from the GSEA website (www.gsea-msigdb.org). For each pathway, replication stress scores were calculated according to the formula below, where e_g is the expression of the gene in different mammary subgroups, and n_g is the number of genes in this pathway.

For luminal progenitor score calculation, the LP signature gene set with average log-fold change was obtained from data of Lim et al. (3), and log CPM matrix of MCF10A-shBRCA1-Tet-on DOX 10 days vs Ctrl RNA-seq data was used. For each sample, LP scores were calculated according to the formula below, which refers to Lim et al. (3), where x_g is the log-fold-change for the gene from LP gene set, and y_g is the log cpm for the same gene in BRCA1 knockdown RNA-seq.

GSEA analysis and barcode plot

GSEA analysis was performed using GSEA software (v4.1.0) (48, 49) and clusterProfiler (v3.14.3) (50). For the customized gene set enrichment analysis of the transcriptional profiles of BRCA1 KD, ELF3 OE and LP, BRCA1 KD data were from MCF10A-shBRCA1-Tet-on DOX 10 days vs Ctrl RNA-Seq data and ELF3 OE data were from MCF10A OE ELF3 vs Ctrl RNA-seq data. LP gene set was obtained from 2016 cell reports, and OE ELF3 gene sets were generated using MCF10A OE ELF3 vs Ctrl RNA-seq data, where the log-fold-change cutoff was 1.3 and P value cutoff was 0.01. For the barcode plot, LP gene sets were obtained from data of Pellacani et al. (28), and log-fold-change of all genes from MCF10A-shBRCA1-Tet-on DOX 10 days vs Ctrl RNA-seq data was used for the enrichment analysis. P value was calculated using the “roast” function from limma (v3.42.2) (51).

Motif enrichment analysis

Motif matrixes of transcription factors were downloaded from the JASPAR database (52). LP gene list was obtained according to the data of Pellacani et al. (28) and LP genes are identified to significantly upregulated in LPs compared with all three other subgroups (BC, LC and SC). Promoter sequences of LP genes were extracted from the genome from 3000 bp upstream to 200 bp downstream of the TSS sites. Motif enrichment analysis was performed using SEA (Simple Enrichment Analysis) of the MEME Suite (53).

Single-cell sequencing data analysis

Single-cell RNA-seq (scRNA-seq) expression profiles of breast cancer tissue and normal mammary tissue from four BRCA1 mutation carriers, as well as normal mammary tissue from three non-carriers were generated using the 10x Genomics Chromium platform and NovaSeq 6000 platform. The detail information of sample processing, scRNA-seq data processing and quality control were described in our previous study (8). Briefly, scRNA-seq data analyses was completed utilizing the 10× Genomics Cell Ranger software suite. Principal component analysis (PCA) (twenty-five principal components) and k-means algorithm (k was set as 9.0) were used to cluster cells. We used a series of cell specific markers to identify the cell type of each cluster, as described in our previous study1. Among these, luminal progenitors cluster was identified based on the expression of KRT8/18+, CD24+, GABRP+, KIT+, and ALDH1A3+. Triple negative tumor cell population was identified based on tumor sample specific cluster, ER/PR/ HER2 expression coincident with the clinical records and PAM50 subtypes analysis. Replication stress score of each cell was calculated based on the sum expression level of genes involving in replication stress-related pathways referring to Dreyer et al. (13). In this study, we selected luminal lineage cells and tumor cells from a BRCA1-mutated triple-negative breast cancer patient to construct tumorigenesis trajectories with the R package Monocle 2 with standard settings.

Animal experiments

5×10⁶ SUM149PT tumor cells were subcutaneously injected into the flanks of female BALB/c nude mice. Tumor growth was measured using a caliper every three days, and the tumor volume was calculated using the formula: tumor volume = length ×width²/2.

Statistics Analysis

Statistical analysis was performed using GraphPad Prism (GraphPad Software). The tests used included the paired and unpaired two-tailed Student’s t test, the Mann-Whitney U test, one-way and two-way ANOVA, and Pearson’s correlation statistical test. A P value of less than 0.05 was considered statistically significant. Details can be found in figure legends.

Study approval

The normal human breast cancer tissue and breast cancer samples used for single-cell RNA-seq in this study were collected from Peking University Cancer Hospital. This study was carried out in accordance with the ethical principles of the Declaration of Helsinki and was approved by the Research and Ethics Committee of Peking University Cancer Hospital. Written informed consent was obtained from all participants prior to participation. Animal experiments in this study were performed in accordance with the Guidelines of Peking University Animal Care and Use Committee.

Data availability

All the data supporting the findings of this study are available within the article and supplementary information file and from the corresponding author upon reasonable request. The RNA-sequencing raw data reported in this paper have been deposited in the Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (PRJCA010883) that are publicly accessible at https://ngdc.cncb.ac.cn.

Acknowledgements

We thank Dr. Jianming Zeng (University of Macau), Xiaojie Sun and all the members of the bioinformatics team, biotrainee, for generously sharing their experience and codes. We thank Ence Yang, Changyu Tao and Zelin Liu for the kind sharing of computing resources and guiding on the RNA-seq data analytical methodology. We thank Lim et al., Pellacani et al. and Bach et al. for their outstanding contribution on BRCA1 associated breast tumorigenesis research and mammary cell populations sequencing data. We thank the ENCODE Consortium and the ENCODE production laboratories for generating the ChIP-seq datasets.

Authors’ Contributions

J.Zhou., X.A.Z. and J.W. designed the research. J.Zhou. performed most of the experiments and data analysis. L.H. performed the single cell sequencing data analysis. X.A.Z. performed the mouse experiments. Y.M. contributed to the RT-qPCR experiments. J.Zhan. performed the immunohistochemistry analysis. Z.X. contributed to experiment ideas and data analysis. M.Z. and Q.S. assisted the experiments. Y.X. provided the single cell sequencing data. J.Zhou. wrote the manuscript. X.A.Z. and J.W. reviewed and edited the manuscript.

Funding Statement

National Key R&D Program of China (2017YFA0503900 and 2016YFC1302100) National Natural Science Foundation of China (81672981 and 81872282) Beijing Municipal Natural Science Foundation (7182082) Clinical Medicine Plus X-Young Scholars Project of Peking University (PKU2018LCXQ01) Special Projects for Strengthening Basic Research of Peking University (BMU2019JC006) The National Science Fund for Distinguished Young Scholars (82125031)