RET Ligands Mediate Endocrine Sensitivity via a Bi-stable Feedback Loop with ERα

The molecular mechanisms of endocrine resistance in breast cancer remain poorly understood. Here we used PRO-seq to map the location of hundreds of genes and thousands of distal enhancers whose transcriptional activities differ between endocrine sensitive and resistant MCF-7 cells. Our genome-wide screen discovered increased transcription of the glial-cell line derived neurotrophic factor (GDNF), a RET tyrosine kinase receptor ligand, which we validate as both necessary and sufficient for resistance in MCF-7 cells. GDNF caused endocrine resistance by switching the active state of a bi-stable feedback loop in the MCF-7 regulatory network from ERα signaling to GDNF-RET signaling. To cause this switch, GDNF downregulated ERα transcription and activated the transcription factor EGR1, which, in turn, induced GDNF. Remarkably, both MCF-7 cells and ER+ primary tumors appear poised for endocrine resistance via the RET signaling pathway, but lack robust RET ligand expression and only develop resistance upon expression of GDNF or other RET ligands.


Introduction 20
Estrogen receptor alpha (ERα) is the major driver of ~75% of all breast cancers. ERα is 21 a transcription factor whose genomic actions are induced upon binding its cognate ligand, 17β-22 estradiol (E2). E2-liganded ERα activates and represses thousands of ERα target genes and 23 non-coding RNAs (Carroll et al., 2006;Hah et al., 2011Hah et al., , 2013. Genes whose transcription is 24 directly activated by ERα promote a mitogenic response in breast cancer cells, resulting in entry To directly test the hypothesis that the genomic actions of ERα are substantially altered 178 in the TamR lines, we treated B7 TamS and G11 TamR MCF-7 cells for 40 minutes with either E2 or 179 Tamoxifen, and monitored transcriptional changes using PRO-seq. As expected, RNA . Surprisingly, E2 also activated 183 transcription in G11 TamR lines (Figure 2A bottom), strongly suggesting that E2 signaling 184 continued to function in TamR lines despite the almost complete lack of GREB1 and PGR. 185 Likewise, direct E2 target genes defined in a previous GRO-seq study (Hah et al., 2011) were 186 largely up-or down-regulated as expected in both B7 TamS and G11 TamR MCF-7 cells ( Figure 2B).
Notably, however, we observed a much more muted effect of E2 on both enhancer and gene 188 transcription in G11 TamR compared with B7 TamS (Figures 2A and 2B), explaining the enrichment 189 in E2 target genes and ERE motifs in differences between TamS and TamR lines, as described 190 above. The reduced effect of E2 on transcription may reflect that the lack of GREB1 or PGR in 191 these lines reduces the effect that ERα has on transcriptional activation. Additionally, however, 192 we also observed a 2.44-fold reduction in the abundance of ERα mRNA ( Figure 2C). Thus,it 193 appears that, while E2 signaling becomes less responsive in G11 TamR MCF-7 cells, the E2 194 signaling pathway remains largely functional and able to affect gene transcription in a stimulus-195 dependent manner. 196 One current model of tamoxifen resistance posits that tamoxifen can function as an ERα 197 agonist in resistant breast cancer cells (Osborne et al., 2003). If this hypothesis is correct, then 198 tamoxifen should promote the activation of ERα target genes in the G11 TamR cells. However, our 199 results showed that tamoxifen had no effect on either enhancer or gene transcription in either 200 B7 TamS or G11 TamR lines (Figures 2A and 2B). Looking genome-wide, the tamoxifen treated 201 B7 TamS and G11 TamR MCF-7 cells are very highly correlated with untreated controls (Spearman's 202 rank correlation ρ > 0.99; Figure S2). The lack of transcriptional differences in either line is 203 consistent with ERα signaling having already been largely shut down under these conditions by 204 three-days of growth in charcoal-stripped FBS, which depletes endogenous E2 from the media. 205 Importantly, these results demonstrate that tamoxifen does not appear to function as an agonist 206 in G11 TamR cells contrary to one current model for endocrine resistance (Osborne et al., 2003). 207 Given that our findings also suggested that E2 signaling remains functional, but muted in 208 the TamR line, we next tested whether ERα was required for the growth of our tamoxifen 209 resistant cells. We found that the viability of both G11 TamR and H9 TamR MCF-7 cells was largely 210 unaffected by treatment with the ER degrader, fulvesterant ( Figure 2D). Therefore, endocrine 211 resistance in G11 TamR and H9 TamR MCF-7 cells appears to occur independently of ERα signaling, suggesting that these TamR lines are likely using an alternative pathway for cell survival and 213 proliferation when grown in the presence of tamoxifen. RET ligands, GDNF, NRTN, ARTN, and PSPN (Sariola and Saarma, 2003). Remarkably, one of 220 these ligands, glial-cell derived neurotrophic factor (GDNF), was among the most highly up-221 regulated genes in both G11 TamR and H9 TamR MCF-7 lines ( Figure 3A). We confirmed the 222 transcriptional differences in GDNF between B7 TamS and G11 TamR MCF-7 cells using qPCR and 223 found that GDNF mRNA levels were increased by ~25 fold in the resistant line ( Figure 3B). 224 Thus both GDNF transcription and mRNA abundance correlate with endocrine resistance in 225 MCF-7 cells, suggesting that GDNF may contribute to the endocrine resistance phenotype. 226 We directly tested this hypothesis by manipulating GNDF levels in our MCF-7 model. We 227 first examined the effects of 10 ng/mL of recombinant GDNF protein on the growth of B7 TamS 228 cells in the presence of antiestrogens. Remarkably, GDNF completely rescued B7 TamS MCF-7 229 cells when challenged with both tamoxifen ( Figure 3C) and fulvestrant ( Figure S3A). Moreover, 230 GDNF treatment without tamoxifen increased the proliferation rate of B7 TamS MCF-7 cells by 231 ~20% (Figure 3C), suggesting that the growth pathways activated by GDNF can work 232 independently of ERα. Next we tested whether GDNF was necessary to confer endocrine 233 resistance in our model system by using short hairpin RNAs (shRNA) to knockdown GDNF in 234 G11 TamR MCF-7 cells. Results show that GDNF depletion (GDNF-KD) reduced GDNF mRNA 235 levels by 57.38% ( Figure 3D) and that these cells were significantly more sensitive to tamoxifen 236 treatment than G11 cells transfected with a scrambled control ( Figure 3E). Moreover, endocrine 237 resistance could be restored to GDNF-KD G11 cells by the addition of 5 ng/ mL recombinant GDNF protein (Figure 3E), demonstrating that growth inhibition does not reflect an off-target 239 effect of the GDNF shRNA. Taken together, these data demonstrate that GDNF plays a central 240 and causal role in establishing endocrine resistance in G11 TamR MCF-7 cells. 241 Having shown that GDNF expression promotes endocrine resistance in our MCF-7 cell 242 model, we next asked whether GDNF mRNA abundance predicts poor relapse free survival 243 (RFS) using publicly available microarray data (Györffy et al., 2010). Indeed, high GDNF 244 expression significantly predicted poor RFS with a hazard ratio of 2.2 (p = 0.028) in one cohort 245 of 88 breast cancer patients ( Figure 3F). GDNF remained significantly correlated with RFS after 246 controlling for expression of ESR1 (ERα), MKI67, and HER2 (ERBB2) using a multivariate 247 analysis (HR = 2.27; p = 0.027). Across 10 sufficiently powered cohorts of patients, GDNF had 248 hazards ratios greater than 1 (i.e., high expression predicts poor RFS outcomes) in seven of 249 these cohorts (mean = 1.758; p = 0.03, two-sided Wilcoxon rank sum test). Moreover, the three 250 studies with significant or borderline significant p-values all had hazards ratios greater than 1 251 (1.62, 1.75, and 2.2; Supplementary Table 2). Taken together, these results suggest a trend in 252 which high transcription of GDNF predicts poor RFS in breast cancer patients, possibly 253 suggesting that GDNF plays a role in endocrine resistance in the clinic. 13 direct transcriptional target of ERα in vivo as well. GFRA1 mRNA encodes the GDNF co-265 receptor, GFRα1, and, together with RET, activates RET-ligand signaling. Further analysis of 266 the mRNA-Seq data set found that GFRA1 is also strongly correlated with ESR1 mRNA in 267 breast cancers (Spearman's ρ = 0.67, p < 2.2e-16; Figure S4A), suggesting that it is also a 268 direct target of E2 signaling. In our MCF-7 endocrine resistance model, GFRA1 transcription is 269 5-higher in TamS MCF-7 cells compared to TamR lines and RET transcription is not significantly 270 different (Figures 4B and 4C), demonstrating that neither factor is overexpressed in TamR 271 MCF-7 cells. These observations suggest that additional mechanisms beyond a high RET or 272 GFRA1 expression level cause endocrine resistance in cell models and in vivo. 273 Our finding that recombinant GDNF was sufficient for endocrine resistance in B7 TamS 274 MCF-7 cells demonstrates that GDNF is a key limiting factor, whose absence prevents TamS 275 cells from taking on a resistant phenotype. To extend this hypothesis to primary breast cancers, 276 we asked whether GDNF expression is low in general, such that it might limit RET pathway 277 activation in most ER+ breast cancers. Indeed, GDNF expression was detectible in only 565 of 278 1,177 primary breast cancers (48%) analyzed by TCGA ( Figure S4B). In principal, RET 279 signaling may be activated by any of the four RET ligands (GDNF, NRTN, ARTN, and PSPN). 280 However, only low levels of NRTN, ARTN, or their co-receptors were detected in primary breast 281 tumors (Figures 4D and 4E; Figure S4B). Thus, we conclude that RET ligand expression is 282 low compared with cell surface receptors, especially RET and GFRα1, which are activated in 283 part by ERα. This contrast between RET receptors and ligands supports a model in which the 284 RET signaling pathway is 'poised' for endocrine resistance by expression of the receptors and 285 that limiting levels of GDNF expression, or possibly of other RET ligands, ensures endocrine 286 sensitivity in most tumors. 287 Next we asked whether high RET ligand expression in a subset of ER+ tumors may 288 explain some cases of endocrine resistance. A careful examination of the GDNF expression 289 distribution in TCGA breast cancers revealed a long tail, indicating high GDNF expression in a 290 handful of cases in the TCGA dataset ( Figure 4E). Our hypothesis that GDNF expression limits 291 RET-dependent endocrine resistance implies that these GDNF-high samples should be prone to 292 endocrine resistance. We devised a simple non-parametric computational approach, which we 293 call the 'outlier score', to quantify the degree to which GDNF is highly expressed based on the 294 symmetry of the empirical probability density function (see methods; Figure 4E, blue line). 295 Based on this score, we conservatively estimate that, of 925 ER+ breast cancer patients in the 296 TCGA dataset, 122 have high expression of at least one of the RET ligands (13%), 57 of which 297 had high levels of GDNF ( Figure 4F). If our proposed model that RET ligands are the limiting 298 factor for endocrine resistance is accurate, cases with this long tail are those that are more likely 299 to be resistant to endocrine therapies. To test this hypothesis, we analyzed expression 300 microarray data collected prospectively by biopsies of patients that either respond, or do not 301 respond, to the aromatase inhibitor letrozole (Miller et al., 2012). A score comprised of the sum 302 of the outlier scores from all four RET ligands is significantly higher in patients that do not 303 respond to letrozole treatment (p= 0.016, one-sided Wilcoxon rank sum test; Figure 4G). By 304 contrast, RET shows no significant difference between patients that respond or do not respond 305 to letrozole. These results suggest that RET ligand expression, but not RET itself, explain the 306 differences in response to letrozole in this cohort of patients. 307 To further explore whether RET ligands contribute to endocrine resistance in primary 308 breast cancers, we asked whether high expression of RET ligands predict RFS. We have  We set out to identify the transcriptional targets activated by GDNF-induced RET 322 signaling. To identify both direct and indirect target genes that respond to GDNF-RET, and to 323 distinguish between them, we collected kinetic PRO-seq data following 0, 1, and 24 hours of 324 GDNF treatment in B7 TamS , C11 TamS , G11 TamR , and H9 TamR MCF-7 cells. We sequenced PRO-325 seq libraries to a high read depth (Table S1) and verified that biological replicates (B7 TamS and 326 C11 TamS ; G11 TamR and H9 TamR ) have highly correlated transcriptional patterns across the time 327 course (Spearman's rank correlation ρ > 0.95; Figures S5A and S5B). 328 We first compared transcriptional changes induced by GDNF between TamS and TamR 329 MCF-7 cells. Because GDNF is both sufficient for resistance in B7 TamS and necessary for 330 resistance in G11 TamR (Figures 3C and 3E), we hypothesized that its effects on gene 331 transcription are also likely to be highly similar in TamS and TamR MCF-7 cells. Consistent with 332 this expectation, transcriptional changes induced by GDNF were highly correlated between 333 TamS and TamR cell lines (Pearson's R > 0.73, p < 2.2e-16; Figures S5C and S5D). As 334 expected, transcriptional responses were lower in magnitude in TamR MCF-7 cells following 335 both 1 and 24 hours of GDNF treatment (Figures S5C and S5D), likely reflecting a dampened 336 GDNF response in TamR lines due to higher basal levels of GDNF. Given these observations, 337 we focused our downstream analyses on TamS MCF-7 cells. 338 We found that GDNF treatment changed the transcription of 4,921 genes, covering 339 rank sum test). Likewise, down-regulated genes were observed to have slightly but significantly 383 higher pausing indices (p < 2.2e-16; Wilcoxon rank sum test). These results suggest that GDNF 384 treatment activates and represses genes in part by changing the rate at which Pol II is 385 transitions from a paused state into productive elongation. 386

ESR1 and GDNF-EGR1 form a bi-stable feedback loop 388
We set out to define the transcriptional regulatory network associated with GDNF-389 dependent endocrine resistance. The dynamics of gene transcription can rigorously separate 390 direct and indirect target genes following a stimulus (Danko et al., 2013). Genes up-regulated 391 during the first 1 hour following GDNF treatment are largely assumed to be direct targets 392 because not enough time has elapsed for transcription, translation, and successive rounds of 393 transcriptional activation. We therefore defined all direct targets of E2 and GDNF signaling as those genes responding by 40 min or 1 hr of treatment, respectively. Secondary targets, defined 395 as transcriptional changes following 24 hrs of GDNF treatment, were assigned to TFs whose 396 transcription changed following 1 or 24 hrs using ChIP-seq data in MCF-7 cells (Euskirchen et 397 al., 2007). 398 The resulting transcriptional regulatory network inferred from the data shows extensive 399 crosstalk between E2 and GDNF signaling programs ( Figure 6A). We predict that GDNF/RET 400 and E2/ERα form a bi-stable feedback loop in which GDNF immediately (1 hr) inactivates the 401 transcription of ERα and activates transcription of EGR1, which, in turn, activates GDNF 402 transcription at 24 hrs ( Figures 6A and 6B). Thus, GDNF is an indirect target of GDNF/ RET 403 signaling that reinforces its own activity through a positive feedback loop dependent on EGR1. 404 In turn, EGR1 transcription is directly down-regulated following 40 min of E2. Thus, GDNF-RET 405 and ERα form a bi-stable feedback loop dependent on EGR1, in which either ERα or 406 GDNF/RET signaling can remain at a high level. 407 408

GDNF-RET signaling down-regulates the E2/ ERα regulatory program 409
To validate the transcriptional regulatory network inferred to underlie endocrine 410 resistance, we first focused on validating ESR1, which encodes the ERα protein, as a direct and 411 immediate GDNF target gene. PRO-seq data found that ESR1 undergoes a two-fold decrease 412 in transcription following 1 hr of GDNF treatment and that this transcriptional change is stable 413 through 24 hrs (Figures 6C and 6D). These changes in ESR1 transcription lead to a 2-fold 414 decrease in ESR1 mRNA abundance following 4 and 24 hours of GDNF treatment ( Figures 6E). 415 Although changes following 1 hour of GNDF treatment are unlikely to reflect indirect effects, it is 416 nevertheless plausible that a transcription factor encoded by a short gene, such as FOSL1, is 417 transcriptionally activated, translated, and responsible for inactivating ESR1. 418 To determine whether changes in ESR1 transcription are a primary or secondary effect 419 of GDNF signaling, we set out to estimate the time at which ESR1 is down-regulated following the addition of GDNF to the cell culture media. To estimate the time at which the ESR1 421 promoter decreases transcriptional activity, we identified the end of the retreating wave of RNA 422 polymerase II 104,000 bp from the transcription start site at 60 min following GDNF treatment 423 To explore the dynamics with which changes in ESR1 transcription lead to differences in 432 ERα protein level, we used Western blotting to track the abundance of ERα and 433 phosphorylated-ERα protein following the addition of GDNF in B7 TamS MCF-7 cells. We found a 434 noticeable effect on ERα protein level as early as two hours after the addition of GDNF and that 435

changes reached their lowest level at 4 hrs (Figures 6G). 436
After 24 hrs of GDNF treatment, we found that the down-regulation of ERα protein likely 437 results in the transcriptional down-regulation of E2 target genes. For example, PGR, GREB1, 438 ELOVL2, and NOS1AP are unaffected at 1 hr, but transcriptionally down-regulated between two 439 and four fold following 24 hrs of GDNF treatment (Figures 6H). We conformed by qPCR that 440 the GDNF-induced decrease in PGR mRNA occurs at 24 hrs but not at 4 hrs ( Figures S6A). 441 Several lines of genome-wide evidence support the indirect effects of GDNF on ERα target 442 genes as well. First, we find that E2 target genes are more than three-fold enriched in the set of 443 genes responding to GDNF at 24 hrs, but not at 1 hr ( Figure 6I), and that transcriptional 444 changes at 24 hours of GDNF negatively correlate with 40 min of E2 (Pearson's R= -0.14; p = 445 4.2e-3). Second, the binding motif that was most enriched in TREs that change Pol II abundance following 24 hrs of GDNF treatment was the ERα binding site (p < 1e-9, Fisher's 447 exact test; Figure S6B). Taken together, these results demonstrate that GDNF-RET signaling 448 down-regulates the E2 regulatory program within 6 hours of treatment by immediate effects on 449 the transcriptional activity of ESR1 during the first 10 min of GDNF treatment. 450 451

GDNF-EGR1 feedback loop results in GDNF activation 452
We next focused on validating the activation loop between GDNF and the transcription 453 factor EGR1. Whereas GDNF transcription increased by 16-fold after 24 hrs of GDNF treatment 454 ( Figure 6B), no changes were found in any of the earlier time points we examined, strongly 455 suggesting that GDNF transcription is indirectly activated by GDNF-induced RET signaling. 456 Regarding how GDNF induces its own expression, we predict that GDNF treatment promotes 457 the upregulation of the transcriptional activator, EGR1, which, in turn, binds to the GDNF 458 promoter thus activating GDNF expression. In support of this hypothesis, we found that EGR1 459 transcription was upregulated more that 30-fold following 1 hr of recombinant GDNF treatment 460 ( Figure 7A). These changes in EGR1 transcription led to an 83-fold increase in EGR1 mRNA 461 abundance following 4 hrs of GDNF treatment ( Figure 7B). In further support of our hypothesis, 462 we identified a TRE in the first intron of GDNF that is bound by EGR1 in MCF-7 cells ( Figure  463 6B). Our model also predicts that EGR1 is directly activated by GDNF signaling, which is likely 464 mediated by an SRF binding site in the EGR1 promoter ( Figure 7A) the abundance of GDNF mRNA following a time course of tamoxifen treatment. As predicted, 477 tamoxifen significantly increased both EGR1 and GDNF mRNA levels following 24 hours of 478 tamoxifen treatment of B7 TamS MCF-7 cells, but not at 40 min or 4 hrs (Figures 7C and 7D). 479 Moreover, several lines of evidence suggest that this mutual suppressive relationship between 480 ERα and EGR1 also holds in primary breast cancers. First, we note a highly significant negative 481 correlation between EGR1 and ESR1 mRNA abundance among ER+ breast cancers analyzed 482 using TCGA RNA-seq data (Pearson's R= -0.21; p = 2.7e-10; Figure 7E). Second, we found 483 that EGR1 transcription increases substantially in primary tumor biopsies following treatment 484 with the aromatase inhibitor letrozole (p = 1.775e-06, Wilcoxon rank sum test; Figure 7F

Discussion 488
Here we have used genomic tools to reconstruct a gene regulatory network that we 489 demonstrate is responsible for endocrine resistance in an MCF-7 breast cancer model. Our 490 approach is uniquely able to distinguish primary from secondary target genes by using PRO-seq 491 to measure nascent transcription over short (≤1 hr) and long (24 hrs presence of ERα protein in breast cancer cells, and would not be resistant to ERα degradation by fulvestrant. In addition, other lines of evidence rule out these models as well, including our 514 direct genome-wide experimental observations demonstrating that E2 remains an agonist and 515 tamoxifen an antagonist in our MCF-7 model (Figures 2A and 2B), as well as a complete lack 516 of genetic changes in ERα protein-coding sequence that are unique to either G11 TamR or H9 TamR 517 cell lines (data not shown). 518 We are also the first to propose that RET-mediated endocrine resistance occurs when 519 Based on these findings, we hypothesize that increased expression of any one of the 531 four RET ligands, GDNF, ARTN, NRTN, or PSPN confers endocrine resistance on cells 532 expressing the RET receptor. In support of this model, we report here that our scoring system 533 based on RET ligand overexpression in tumors significantly separates breast cancer patients 534 that respond to letrozole from those who do not ( Figure 4G). Moreover, we found that RET 535 ligands are predictive of relapse free survival in other cohorts of patents, even after accounting 536 for the expression of other prognostic markers such as ER, PR, and HER2 ( Figure 3F). Several 537 findings also strongly support the involvement of GDNF in endocrine resistance in our MCF-7 538 model, most notably the observations that GDNF rescues B7 TamS lines and that GDNF knockdown in G11 cells restores sensitivity to tamoxifen ( Figure 3E). These observations are 540 also supported by existing studies showing that another RET ligand, ARTN, contributes to 541 tamoxifen resistance in MCF-7 cells (Kang et al., 2010), extending and supporting the findings 542 reported here. However, there is one RET ligand that is notably an outlier. PSPN does not 543 appear to have any predictive value in patients, and thus may not play the same role in 544 resistance as the other three RET ligands. This may reflect the extremely low expression of its 545 co-receptor, GFRA4, in primary breast cancers (Figure S4B), preventing PSPN from having 546 much effect on breast cancer cells. Taken together, these findings suggest that RET ligand 547 expression, especially GDNF, ARTN, and NRTN, explain endocrine resistance in many cases. 548 One finding that our current study cannot yet explain is that our proposed bi-stable 549 regulatory network between ERα and GDNF/ EGR1 leads to the activation of GDNF in TamS as 550 well as in TamR MCF-7 cells. Under our model, tamoxifen treatment in either TamS or TamR 551 lines leads to the transcriptional activation of GDNF within 24 hours, a prediction of our model 552 which we were able to validate by qPCR ( Figure 7D). Thus, it remains unclear why endogenous 553 transcription at the GDNF locus is not sufficient to confer endocrine resistance in B7 TamS cells. 554 One potential explanation is that higher basal GDNF expression in TamR MCF-7 cells grown in 555 estrogen containing media (Figures 3A and 3B) gives TamR lines a "head start" when 556 switching growth signaling programs from ERα to GDNF/RET. Testing this model will require It is also unclear how RET ligand expression is activated in primary tumors. The 560 abundance of GDNF mRNA appears to be extremely low in primary breast tumors analyzed by 561 TCGA (Figures 4D, 4E and, S4B), which were largely collected before therapeutic intervention 562 prospective clinical studies targeting larger cohorts of patients starting endocrine therapies will 578 be required to fully validate our proposed mechanism of endocrine resistance.

Author Contributions 580
The project was conceived by CGD, SAC, and SH. All cell culture and molecular experiments 581 were done by SH, HZ, LJA, CM, and BAM. PRO-seq experiments were conducted by EJR and 582 SH. Genomic data was analyzed by CGD, TC, and SH. Data collection, experiments, and 583 analysis was supervised jointly by CGD and SAC. The paper was written by SH, CGD, and SAC 584 with input from all authors. 585 586 Acknowledgements 587 We thank X. Yao, L. Lan, as well as all members of the Danko and Coonrod labs for valuable 588 insights and discussions. We also thank J. Lewis for input on the manuscript draft. Work in this