TRPS1 is associated with breast cancer incidence and promotes breast cancer cell number accumulation

A recent genome-wide association study (GWAS) identified 32 novel single nucleotide polymorphisms (SNPs) associated with breast cancer susceptibility [28]. When we queried the NHGRI EBI GWAS Catalog to find published associations with genetic variants within the TRPS1 genomic locus, we found one of the lead SNPs from this study [29]. Furthermore, in the authors’ subtype-specific analysis of the results, the association with this variant was strongest for luminal breast cancers relative to other subtypes [28].

We used LocusZoom to plot the data within this locus (Fig 1A). A plot of these data indicates that two sets of SNPs have low linkage disequilibrium with one another, indicating that they are inherited independently and each confer risk. One set of variants is within an intronic region of the TRPS1 gene, and one is about 400 kilobases upstream from the transcription start site (TSS) of the TRPS1 gene. There are often many genes in close proximity to the lead SNP in a GWAS, making it difficult to predict which gene mediates the effect on the associated phenotype. However, in this case the nearest gene is almost a megabase upstream from TRPS1, suggesting that TRPS1 contributes to the breast cancer susceptibility associated with one or both of these sets of genetic variants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. TRPS1 is associated with breast cancer incidence and promotes breast cancer cell fitness.

A: LocusZoom plot of the TRPS1 genomic locus depicting the location and significance of SNPs associated with breast cancer susceptibility. TRPS1 is the closest gene to two sets of genetic variants in low linkage disequilibrium with one another. Data from [28], generated with summary statistics downloaded from the NHGRI-EBI GWAS Catalog, using LocusZoom [29, 30]. B: Scatter plot of TRPS1 and ESR1 knockout scores for each gene tested. Scores are normalized such that knockout of a gene with a score of 0 has no effect on cell number, and knockout of a gene with a score of -1 has an effect equal to that of knocking out one of a set of universally essential genes. Luminal breast cancer cell lines are colored in red. The data are from the Cancer Dependency Map project [31]. C) Violin and box and whisker plots of TRPS1 knockout scores from (B) for luminal breast cancer cell lines versus all other cell lines. Wilcoxon rank sum test p-value of 3.2*10^-5.

https://doi.org/10.1371/journal.pgen.1011159.g001

Based on this result, we hypothesized that perturbation of TRPS1 in luminal breast cancer cell lines would affect cell fitness. Using data from the Cancer Dependency Map project, we found that sensitivity to TRPS1 knockout correlated with sensitivity to knockout of ESR1, the gene encoding ER (Fig 1B) [31]. Furthermore, while TRPS1 knockout led to an increase in cell number for most cancer cell lines, luminal breast cancer cell lines were significantly enriched for TRPS1 dependency (Fig 1C).

Taken together, these data indicate that TRPS1 influences breast cancer incidence and is required for maximal breast cancer cell fitness. Next we sought to determine how TRPS1 regulates its target genes to mediate these breast cancer patient and cellular outcomes.

Endogenously degron-tagged TRPS1 is rapidly degraded in T47D cells

To rapidly deplete TRPS1 and isolate primary TRPS1-regulated genes, we employed the dTAG system for targeted protein degradation [32, 33]. We inserted an inducible degron tag into the endogenous TRPS1 locus in the luminal breast cancer cell line T47D. We generated three independent clones that express the tagged TRPS1 protein that can be degraded by the addition of the small molecules dTAG-13 and dTAG^V-1 at 50nM each (dTAG) (Fig 2A). Of note, we depleted around 50% of TRPS1 from whole cell lysates in 10 minutes of treatment with dTAG, as determined by quantitative western blot, with less than 10% detected as soon as 20 minutes and as late as 48 hours after treatment (Fig 2B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Endogenously degron-tagged TRPS1 is rapidly degraded in T47D cells.

A: Quantitative Western blot with a serial dilution of the parental T47D cells followed by three independent dTAG-TRPS1 clones treated with DMSO or dTAG-13 and dTAG^V-1 at 50nM each (dTAG) for 2 hours. Membranes were probed with anti-TRPS1 or anti-ACTB antibodies. B: Quantitative Western blot with a serial dilution of dTAG-TRPS1 Clone 28 followed by a time course of treatment with dTAG. Membranes were probed as in (A). C: Heatmap of TRPS1 ChIP-seq peaks, in rows ranked by intensity, in cells treated with DMSO or dTAG for 30 minutes. D: MA plot of TRPS1 ChIP-seq peaks, with fold change values representing binding intensity in the dTAG condition relative to the DMSO condition. All points are colored blue to indicate they are significantly decreased at an FDR of 0.1.

https://doi.org/10.1371/journal.pgen.1011159.g002

To ensure this treatment depleted TRPS1 from chromatin, we performed chromatin immunoprecipitation with sequencing (ChIP-seq) using an anti-HA antibody to recognize the 2xHA tag within the degron tag. Our ChIP-seq libraries were of high quality, as measured using the quality control metrics of fraction of reads in peaks (S1 Fig) and the peak calling-independent strand cross-correlation (S2 Fig), generated using the ChIPQC R package [34]. We called peaks using MACS2 [35], using all DMSO samples together and all dTAG samples as the control. We observed a genome-wide decrease in TRPS1 binding intensity, with over 80 percent of TRPS1 depleted from chromatin after 30 minutes of dTAG treatment (Fig 2C and 2D).

TRPS1 directly represses regulatory element activity

With this system in hand, we set out to test the effects of TRPS1 depletion on regulatory element (RE) (i.e., enhancer or promoter) activity. To capture the dynamics of chromatin accessibility after TRPS1 depletion, we conducted a time course analysis using the Assay for Transposase-Accessible Chromatin with sequencing (ATAC-seq). The time points included dTAG treatments at 30 minutes, 1 hour, 2 hours, 4 hours, and 24 hours, while a DMSO treatment served as the vehicle control, assigned as the zero minute time point. Our ATAC-seq libraries were of high quality, as determined by plotting the distribution of fragment sizes (S3 Fig) and the enrichment of signal around TSS’s (S4 Fig), generated using the ATACseqQC R package [36]. We generated a consensus peak set using MACS2 [35] for all samples together. At the earliest time point after degradation, our best estimate of the primary effects of TRPS1 depletion, we identify 472 peaks that increased in intensity and 36 that decreased in intensity, at a false discovery rate (FDR) of 0.1. (Fig 3A and S1 and S2 Tables). We hypothesized that the increased peaks result from loss of a direct TRPS1 reduction of chromatin accessibility, with the decreased peaks an indirect effect of the redistribution of limiting cofactors.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. TRPS1 directly represses regulatory element activity.

A: MA plot of ATAC-seq peaks, with fold change values representing accessibility in the 30 minute dTAG-13 and dTAG^V-1 at 50nM each (dTAG) treatment condition relative to the DMSO condition. B: De novo motif identified in increased peaks from (A) (below), matched to the TRPS1 motif (above). C: Bar charts of prevalence of a representative GATA motif in increased, unchanged, and decreased peaks from (A). Chi-square test p-value < 2.2*10^-16. D: Bar charts of prevalence of a representative ER half-site in increased, unchanged, and decreased peaks from (A). Chi-square test p-value 4.5*10^-15. E: Heat map with hierarchical clustering of chromatin accessibility in ATAC-seq peaks that are significantly changed over the time course at an FDR of 0.1. F: Heat map of prevalence of the representative GATA motif in increased, unchanged, and decreased peaks, as in (C), for each time point relative to the DMSO condition. G: Density plot for composite PRO-seq signal across all TRPS1 ChIP-seq peaks, separated by strand and treatment condition.

https://doi.org/10.1371/journal.pgen.1011159.g003

To test this hypothesis, we performed de novo motif identification in the increased and decreased peaks. We identified a GATA motif in the increased peaks but not the decreased peaks (Fig 3B). To explicitly calculate the motif prevalence in each class of peaks, we found individual motif occurrences genome-wide and intersected the peaks with these motif instances. We found a significant enrichment of a representative GATA motif in the increased peaks over the unchanged and decreased peaks (Fig 3C).

To identify additional TFs that might be contributing to changes in chromatin accessibility, we analyzed motif enrichment in both increased peaks relative to unchanged peaks as well as decreased peaks relative to unchanged peaks (S3 and S4 Tables). To our surprise, we identified nuclear receptor motifs in both classes of peaks. To follow up this analysis, we again calculated motif prevalence in each class of peaks for a representative ER half-site. For this motif, we found a significant enrichment in both increased and decreased peaks relative to unchanged peaks (Fig 3D).

We predicted that the enrichment of the GATA motif within increased peaks would wane over time. To test this prediction, we turned to our time course data. We identified ATAC-seq peaks that were significantly changed over the time course at an FDR of 0.1, using a likelihood ratio test within DESeq2 [37] to identify peaks for which including a variable for the time point increased the predictive power of the model over one without this information. We performed hierarchical clustering of these peaks and found that the replicates for each time point clustered together and that the majority of dynamic peaks changed gradually over the time course, with additional clusters displaying different kinetics (Fig 3E). We called increased, unchanged, and decreased peaks as in Fig 3A for each time point relative to the control condition. We then calculated GATA motif prevalence in each set of peaks as in Fig 3C. As we expected, the GATA motif prevalence in the increased peaks decreased over time, but the majority of increased peaks at 24 hours still contained a GATA motif (Fig 3F). This is consistent with the primary effects of TRPS1 depletion driving a large proportion of the changes in chromatin accessibility even as late as 24 hours after TRPS1 depletion.

As an orthogonal readout of RE activity, we measured bidirectional transcription around TRPS1 binding sites using precision genomic run-on with sequencing (PRO-seq) [38]. Our libraries were of high quality using several quality control metrics (S5 Fig) [39]. Using a window centered on each summit of TRPS1 ChIP-seq intensity, we observed an increase in bidirectional transcription 30 minutes after TRPS1 depletion (Fig 3G).

Since composite profiles have limitations, we chose to look more closely at each TRPS1 peak and determine if increasing bidirectional transcription is a reproducible trend. We counted PRO reads in the window around the TRPS1 peak summits shown in Fig 3G and used DESeq2 to identify differentially transcribed regions. We normalized our counts based on the size factors generated from the DESeq2 analysis of the reads in genes from Fig 4. As read counts in small putative regulatory elements tend to be fewer than in large genes, it is not surprising that only 9 of these regions have individually statistically significantly higher bidirectional transcription and none with significantly lower transcription, at a false discovery rate of 0.1 (S6 Fig). However, we did find a TRPS1 cistrome-wide increase in PRO signal, with over 62% of these regions increasing in bidirectional transcription upon TRPS1 depletion. We performed ANOVA on a linear model predicting the logarithm of the normalized PRO reads for each region based on the DMSO or 30 minute dTAG treatment condition across the four replicates, and the p-value for the F-test was < 2.2*10^-16.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. TRPS1 directly represses transcription of target genes.

A: Kinetic traces of the two major clusters of activated genes over the time course. B: Kinetic traces of the two major clusters of repressed genes over the time course. C: Over-representation analysis of the activated genes from (A). D: Over-representation analysis of the repressed genes from (B). E: Cumulative distribution function plot of the distance from the TSS of each gene to the nearest TRPS1 ChIP-seq peak overlapping an increased ATAC-seq peak, by gene class. Kolmogorov–Smirnov test between activated and unchanged genes: p-value = 0.011. F: Cumulative distribution function plot of the distance from the TSS of each gene to the nearest TRPS1 ChIP-seq peak overlapping an decreased ATAC-seq peak, by gene class. KS test between repressed and unchanged genes: p-value > 0.1.

https://doi.org/10.1371/journal.pgen.1011159.g004

The biological interpretation is consistent with our other results that indicate a role of TRPS1 in chromatin compaction. We would not necessarily expect an increase in bidirectional transcription at each TRPS1 peak upon TRPS1 degradation, but the increase in chromatin accessibility could facilitate downstream factors that promote transcription initiation. As chromatin accessibility and bidirectional transcription can affect one another, we next isolated ATAC-seq peaks without bidirectional transcription, as identified using dREG [40]. As we saw in Fig 3A, chromatin accessibility predominantly increased upon TRPS1 depletion (S7 Fig). Along with our accessibility data and motif analysis, these data indicate that TRPS1 directly represses RE activity primarily via its effect on chromatin accessibility.

TRPS1 directly represses transcription of target genes

Downstream from the changes in RE activity, we measured changes in nascent transcription within genes with PRO-seq over the same time course of TRPS1 depletion as in the ATAC-seq experiment. As with our ATAC-seq time course analysis, we identified genes that were significantly changed over the time course at an FDR of 0.1, using a likelihood ratio test within DESeq2 [37] to identify genes for which including a variable for the time point increased the predictive power of the model over one without this information. We identified 1,425 dynamic genes over the time course and performed hierarchical clustering to classify the genes based on their expression kinetics (Fig 4A and 4B and S5, S6, S7 and S8 Tables). Over-representation analysis (ORA) of the activated genes identified several enriched Hallmark gene sets [41], most prominently cholesterol homeostasis genes (Fig 4C). ORA on the repressed genes revealed that the two estrogen response gene sets were the most significantly enriched of the Hallmark gene sets.

We hypothesized that TRPS1 regulates these gene sets by distinct mechanisms. Specifically, we predicted that TRPS1 directly represses the dTAG-activated genes by repressing the activity of REs proximal to these genes. To test this prediction, we measured the distance from the transcription start site (TSS) of each gene to the nearest TRPS1 ChIP-seq peak overlapping an increased ATAC-seq peak. By constructing a cumulative distribution function, we found that the activated genes are significantly closer to these activated REs (Fig 4E).

In contrast to the genes activated by TRPS1 depletion, we predicted that the effects on the genes repressed by TRPS1 depletion are indirect and distal to TRPS1 binding. To test whether TRPS1 directly activates a subset of REs to activate transcription of proximal genes, we measured the distance from the TSS of each gene to the nearest TRPS1 ChIP-seq peak overlapping a decreased ATAC-seq peak. There are few examples of this class of RE, so the distances were much farther, and there was no significant enrichment of repressed genes proximal to these peaks (Fig 4F). These data suggest that, while TRPS1 positively and negatively regulates hundreds of primary response genes, TRPS1 only represses transcription of its direct target genes.

TRPS1 redistributes ER binding to modulate ER target gene transcription

Based on the correlation between cancer cell line sensitivity to TRPS1 knockout and ESR1 knockout (Fig 1B), the enrichment of an ER binding motif in both increased and decreased ATAC-seq peaks (Fig 3D), and the over-representation of estrogen response gene sets in the genes repressed by TRPS1 depletion (Fig 4D), we focused on ER target genes to explore a possible mechanism by which TRPS1 indirectly activates transcription. We first defined direct ER target genes using our previously-generated PRO-seq data from parental T47D cells that were hormone starved for three days and then acutely stimulated with estrogen or a DMSO vehicle control for 90 minutes (S8 Fig) [39]. We exclusively focused on estrogen-activated genes because ER directly activates these genes [8, 42, 43]. 65 ER target genes were activated, and 58 were repressed by acute TRPS1 depletion (Fig 5A). To test the robustness of this change in ER target gene transcription, we additionally performed PRO-seq in each of the three independent TRPS1-dTAG clones generated from the parental T47D cells. Indeed, the genes we identified as TRPS1-regulated ER target genes in the one clone used for the time course experiment tend to be regulated in the same direction upon TRPS1 degradation across these three clones, suggesting that this effect is robust across the cell lines in which we can acutely deplete TRPS1 (S9 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. TRPS1 redistributes ER binding to modulate ER target gene transcription.

A: Kinetic traces of activated and repressed ER target genes over the time course. B: Violin and box and whisker plots for ER binding intensity fold change upon TRPS1 depletion at ER ChIP-seq peaks within 100kb of the TSS of each gene, grouped by gene class defined by change in expression upon TRPS1 depletion. (In (B) and (C), *** represents a significant one-sample t-test p-value < 10^-3. N.S. represents a non-significant p-value > 0.1.) C: Violin and box and whisker plots for ER binding intensity fold change upon TRPS1 depletion at ER ChIP-seq peaks, grouped by summit-to-summit distance to the nearest TRPS1 ChIP-seq peak. D: Model of TRPS1-mediated ER redistribution and modulation of ER target gene transcription. Transparent boxes indicate reduced binding intensity or attenuated transcription. Above, at baseline, TRPS1 directly decreases ER binding intensity proximal to TRPS1, attenuating ER activation of proximal ER target genes. Distal to TRPS1, ER binding intensity is not directly affected by TRPS1, and ER fully activates proximal ER target genes. Below, after TRPS1 depletion, ER binding proximal to TRPS1 increases in intensity, augmenting ER activation of proximal ER target genes. Distal to TRPS1, ER binding intensity is indirectly decreased, as limiting ER molecules are redistributed to TRPS1-proximal regulatory elements, attenuating ER activation of proximal ER target genes. E: ChIP, ATAC, and PRO density around an example increased ER binding site near an activated ER target gene. At this TRPS1-proximal ER binding site, upon dTAG treatment, TRPS1 binding intensity decreases, ER binding intensity increases, chromatin accessibility increases, and gene expression increases. F: ChIP, ATAC, and PRO density around an example decreased ER binding site near an repressed ER target gene. At this TRPS1-distal ER binding site, upon dTAG treatment, ER binding intensity decreases, chromatin accessibility decreases, and gene expression decreases. In (F) and (G), dTAG refers to dTAG-13 and dTAG^V-1 at 50nM each.

https://doi.org/10.1371/journal.pgen.1011159.g005

We hypothesized that these changes in ER target gene transcription are mediated by changes in the genomic distribution of ER binding. We performed ER ChIP-seq to test the prediction that ER binding intensity proximal to dynamic ER target genes would change in concordance with the change in gene transcription. As before, our ChIP-seq libraries were of high quality (S1 and S2 Figs). We called peaks using MACS2 [35], using all ER samples together and all IgG samples as the control. Consistent with our hypothesis, we found that ER ChIP-seq peaks within a 100 kilobase (kb) window around the TSS of activated genes tended to increase in intensity, and ER ChIP-seq peaks within a 100kb window around the TSS of repressed genes tended to decrease in intensity (Fig 5B).

We further hypothesized that only the increased ER binding sites represent a direct effect of TRPS1 activity. Consistent with this hypothesis, we found that ER binding proximal to TRPS1 tends to increase in intensity, and ER binding distal to TRPS1 tends to decrease in intensity (Fig 5C). Together, these data suggest a model in which TRPS1 depletion redistributes ER binding from TRPS1-distal sites to TRPS1-proximal sites and modulates ER target gene transcription proximal to the dynamic ER binding sites (Fig 5D). To illustrate this phenomenon, we provide an example of a TRPS1-proximal, increased ER peak near an activated ER target gene in (Fig 5E), and a TRPS1-distal, decreased ER peak near a repressed ER target gene in (Fig 5F).

TRPS1 activity is associated with breast cancer patient outcomes

We next sought to connect the primary TRPS1-responsive genes with downstream cellular and patient-related outcomes. We defined a new steady state of transcription with the 24 hour time point after TRPS1 depletion. We ranked genes based on their shrunken fold change in PRO signal (Fig 6A). Using this ranking, we performed gene set enrichment analysis with the Hallmark gene sets and found multiple cell-cycle-related gene sets to be negatively enriched, including E2F Targets (Fig 6B) [44, 45]. Consistent with this, we observed a significant decrease in cell number doubling rate of T47D dTAG-TRPS1 cells upon TRPS1 depletion (Fig 6C). Importantly, the isogenic parental T47D cells do not display a cell number defect with dTAG treatment, so we attribute this effect to TRPS1 depletion and not a non-specific effect of the compounds (S10 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. TRPS1 activity is associated with breast cancer patient outcomes.

A: MA plot of PRO signal in each gene, with shrunken log fold change values representing transcription in the 30 minute dTAG-13 and dTAG^V-1 at 50nM each (dTAG) treatment condition relative to the DMSO condition. B: Mountain plot of the Hallmark E2F Targets gene set, using genes ranked by shrunken fold change from (A). A negative enrichment score indicates an enrichment of the gene set among repressed genes. Adjusted p-value 2.5*10^-19. C: Cell number over time of dTAG-TRPS1 cells treated with dTAG or DMSO. Analysis of variance for the coefficient corresponding to the difference in doubling rates between the conditions in a linear model of the logarithm of cell number versus time: p-value 1.1*10^-5. D: Differential enrichment score (dES) calculation for an example patient with the highest TRPS1 activity. Above, genes ranked by scaled expression in this patient relative to all other patients in the METABRIC cohort [46, 47]. Below, gene set enrichment analysis of TRPS1-repressed and TRPS1-activated genes defined by response after 30 minutes of TRPS1 depletion. E: Patients from the METABRIC cohort, ranked by dES as calculated in (D), with classifications of the tumors on the right. F: Kaplan-Meier curves for patients in the METABRIC cohort, stratified by TRPS1 activity as in (E). Logrank p-value 4.99*10^-4.

https://doi.org/10.1371/journal.pgen.1011159.g006

Finally, we calculated a TRPS1 activity score by adapting methods developed by [48, 49]. We used our PRO-seq data to determine a primary TRPS1 regulon based on the differentially expressed genes 30 minutes after TRPS1 depletion. We classified breast cancer patients from the METABRIC cohort as having high TRPS1 activity if both a) TRPS1-repressed genes are negatively enriched and b) TRPS1-activated genes are positively enriched, relative to all other patients in the cohort (example patient in (Fig 6D)) [46, 47]. Similarly, we classified patients as having low TRPS1 activity if both a) TRPS1-repressed genes are positively enriched and b) TRPS1-activated genes are negatively enriched. We classified the remaining patients as having intermediate TRPS1 activity. We ranked patients based on their TRPS1 activity and found no association with other clinical covariates (Fig 6E). When we stratified patients by TRPS1 activity, we found TRPS1 activity to be significantly associated with shorter survival time (Logrank p-value 4.99*10^-4) (Fig 6F). When we first separated tumors by ER-positivity or by intrinsic subtype, we found that this association was specific for ER-positive tumors and Luminal A tumors (S11 Fig). We also performed this analysis using genes differentially expressed after 24 hours of TRPS1 depletion (S12 Fig). However, by this time cell cycle genes dominate, and we speculate that the association with breast cancer patient survival is due to these genes. We believe using the primary response genes offers a unique measure of TRPS1 activity not achievable using previous surrogates.

GWAS and DepMap data visualization

Summary statistics from [28] were downloaded from the NHGRI-EBI GWAS Catalog [29]. SNPs in the TRPS1 locus were plotted using LocusZoom [30]. Knockout scores and luminal breast cancer identifiers were downloaded from the Cancer Dependency Map project [31] and plotted using the statistical programming language R [62].

Cell culture

T47D cells (RRID:CVCL_0553) (ATCC) were cultured in RPMI 1640 medium (Gibco) supplemented with 10% fetal bovine serum (Gemini) and 10μg/ml insulin from bovine pancreas (Sigma, made as a 1000x solution in 1% aqueous glacial acetic acid).

Plasmid generation for gene editing

DNA for transfection was prepared as previously described [63, 64]. A CRISPR sgRNA (TTATCTTTGCAGATATGGTC) targeting the 5’ end of the TRPS1 coding sequences was designed using Benchling. The sgRNA was cloned into hSpCas9 plasmid PX458 (Addgene #48138) as previously described [65], using the following primers:

1. 5’ -CACCGTTATCTTTGCAGATATGGTC-3’ and
   5’ -AAACGACCATATCTGCAAAGATAAC-3’.

A plasmid harboring a synthetic HygR-P2A-2xHA-FKBP_F36V insert was generated with Cold Fusion (System Biosciences), starting with the HygR-P2A-AID cassete in pMGS58 (Addgene #135311) [64] and the Puro-P2A-2xHA-FKBP_F36V casette in (Addgene #91793) [32]. The linear donor was generated by PCR using primers (IDT) that contain 50-nucleotide homology tails and gel-purified. The primers contained 5’ phosphorothioate modifications to increase PCR product stability in the cell [66]. The primers used for making PCR donor fragments were:

1. 5’ -G*T*AACTTTCAGATAACACTGTATCTGCCTTTT
   CCCTTTATCTTTGCAGATATGAAAAAGCCTGAACT
   CACCG-3’
   and
   5’ -T*T*CACTTGCAACGTTTCTCAGAGGGGGGTTCT
   TTTTCCGGACACCTGAACCTGAACCTCCAGATCCA
   CCAGATCTTTCCAGTTTTAGAAGCTCCACATCG-3’
   with asterisks representing the phosphorothioate modifications.

dTAG-TRPS1 clone generation

Clones were generated as previously described [63, 64], with modifications. An initial round of cloning was performed using puromycin selection, but upon genomic DNA sequencing this clone did not appear to have a dTAG insertion event within TRPS1. Nevertheless, this clone was used for a second round of cloning using hygromycin selection. 3*10⁶ cells were plated in 10cm plates. The next day, cells were cotransfected with 15μg of CRISPR/Cas9-sgRNA plasmid and 1.85μg of linear donor PCR product using Lipofectamine 3000 (Thermo Fisher Scientific) in Optimem (Gibco). One day after transfection, the media was replaced. Starting four days after transfection, cells were selected for two weeks with 200μg/mL of Hygromycin B (Invitrogen) with 20% conditioned media, replaced twice per week. Colonies were then grown in 20% conditioned media, replaced twice per week, until they were large enough to be picked and passaged to a 24-well plate. Clones were expanded and frozen at 8 passages after transfection. Integration was tested with Western blotting, PCR, and Sanger sequencing. In each of the three clones, two to three of the four genomic copies of TRPS1 are knocked out, and only tagged TRPS1 protein is expressed. Details about each of the determined alleles are available at https://guertinlab.github.io/TRPS1_ER_analysis/Vignette.html#allele-sequencing.

Western blotting

8*10⁵ cells per sample were plated in each well of a 6-well plate. Cells were treated with DMSO or 50nM dTAG-13 and 50nM dTAG^V-1 in DMSO at various time points and collected simultaneously. At the time of harvest, cells were scraped and lysed in RIPA buffer (1% Nonidet P-40, 1% sodium deoxycholate, 0.1% sodium dodecyl sulfate, 2mM EDTA, 150mM NaCl, 10mM sodium phosphate, 50mM NaF, 50mM Tris pH 7.5), with 100μM benzamidine, 5μg/mL aprotinin, 5μg/mL leupeptin, 1μg/mL pepstatin, 1mM phenylmethylsulfonyl fluoride, and 2mM sodium orthovanadate added fresh. Lysates were sonicated in a Biorupter UCD-200 (Diagenode) on high for 30 seconds on and 30 seconds off for 5 cycles, and clarified by centrifugation at 14,000rpm for 15 min in 4°C. Protein concentration was measured by BCA assay and diluted to the same concentration. 10× Laemmli buffer was added to a final concentration of 1x, and 2-mercaptoethanol was added to a final concentration of 1%. Samples were boiled at 95°C for 10 minutes, and 30μg of each was loaded into a 10% polyacrylamide gel. Samples were separated by gel electrophoresis and transferred to nitrocellulose membranes. Membranes were incubated in blocking buffer (3% bovine serum albumin, 1X Tris buffered saline) for 1 hour at room temperature with rocking. Primary antibodies (anti-TRPS1, Cell Signaling #17936S, and anti-ACTB, Cell Signaling #3700S) were diluted 1:1,000 in primary buffer (3% bovine serum albumin, 0.1% sodium azide, 0.1% Tween-20, 1X Tris buffered saline) at 4°C with rocking overnight. Fluorescent secondary antibodies were diluted 1:10,000 in secondary buffer (5% bovine serum albumin, 0.1% sodium azide, 0.1% Tween-20, 1X Tris buffered saline) and incubated for 1 hour at room temperature with rocking, and fluorescence was measured (Odyssey, Licor).

ChIP-seq library preparation

2.4*10⁷ cells per sample were plated across 3 15cm dishes 2 days before harvest. Cells were treated with DMSO or 50nM dTAG-13 and 50nM dTAG^V-1 in DMSO for 30 minutes and collected simultaneously. At the time of harvest, cells were fixed with 1% formaldehyde (Sigma) for 10 minutes at 37°C and quenched with 125mM Glycine (Fisher) for 10 minutes at 37°C. Plates were moved to ice, and cells were washed and scraped into ice cold PBS containing Complete EDTA-free Protease Inhibitor Cocktail (Roche). Cells were pelleted in aliquots of 3.6*10⁷ cells, snap frozen in liquid nitrogen, and stored at -80°C. Pellets were thawed, and cells were lysed in 1mL Cell Lysis Buffer (85mM KCl, 0.5%NP40, 5mM PIPES pH 8.0), with protease inhibitor cocktail added fresh, for 10 minutes with rotation at 4°C. Nuclei were pelleted at 3300g at 4°C for 5 minutes and resuspended in 500μL ChIP lysis buffer (0.5% SDS, 10mM EDTA, 50mM Tris-HCl pH 8.1), with protease inhibitor cocktail added fresh, for 10 minutes with rotation at 4°C. Lysates were moved to 15ml polystyrene conical tubes (Falcon) and sonicated in a Biorupter UCD-200 (Diagenode) on high for 30 seconds on and 30 seconds off for 4 sets of 5 cycles. Before each set, ice in the water bath was replaced, and samples were gently vortexed to mix. Sonicated lysates were then move to 1.5ml tubes and clarified by centrifugation at 14,000rpm for 15 min in 4°C. 500μL of the supernatant was diluted into 6.5mL Dilution Buffer (0.01% SDS, 1.1% Triton X-100, 1.2mM EDTA, 167mM NaCl, 16.6mM Tris-HCl pH 8.0), with protease inhibitor cocktail added fresh (1*10⁶ cells in 200μL). 1ml (5*10⁶ cells) was aliquoted into each of 3 tubes with antibody (1.25 μg anti-HA, Cell Signaling #3724S, 2.5μg anti-ER, Millipore #06–935, or 2.5μg IgG control, Cell Signaling #2729S), and incubated with end-over-end rotation at 4°C overnight.

50μL Protein A/G Magnetic Beads (Pierce) per sample were washed with bead washing buffer (PBS with 0.1% BSA and 2mM EDTA) and then incubated with samples for 2 hours with rotation at 4°C. The samples were washed once each with low salt immune complex buffer (0.1% SDS, 1% Triton x-100, 2mM EDTA, 150mM NaCl, 20mM Tris HCl pH 8.0), high salt immune complex buffer (0.1% SDS, 1% Triton x-100, 2mM EDTA, 500mM NaCl, 20mM Tris Hcl pH8.0), LiCl immune complex buffer (0.25M LiCl, 1% NP-40, 1% deoxycholate, 1mM EDTA, 10mM Tris-HCl pH8.0), and 1xTE (10mM Tris-HCl, 1mM EDTA pH8.0). Immune complexes were eluted in elution solution, (1% SDS, 0.1M sodium bicarbonate) in a thermomixer for 30 min at 65°C at 1,200rpm. Crosslinks were reversed and proteins were digested with the addition of 200mM NaCl and 2ul Proteinase K in a thermocycler at 65°C for 16 hours. DNA was purified with a Qiaquick PCR cleanup (Qiagen), and libraries were prepared with a NEBNext Ultra II Library Prep Kit (New England Biolabs).

ChIP-seq analysis

Adapters were removed using cutadapt [67]. Reads were aligned to the hg38 genome assembly with bowtie2 [68]. Duplicate reads were removed, and the remaining reads were sorted into BAM files and converted to bed format for counting with samtools [69]. Reads were also converted to bigWig format with deeptools [70]. Peaks were called with MACS2 [35]. Reads were counted in peaks using bedtools, and differentially bound peaks were identified with DESeq2 [37, 71]. Heatmaps were generated with deeptools. Peak proximity to and overlap with other features were calculated with bedtools.

ATAC-seq library preparation

ATAC-seq libraries were prepared as previously described [72], with modifications. 4 replicates were performed from cells treated and collected at different times in the same day. 4*10⁵ cells per sample were plated in each well of a 6-well plate 2 days before harvest. Cells were treated with DMSO or 50nM dTAG-13 and 50nM dTAG^V-1 in DMSO at various time points and collected simultaneously. At the time of harvest, cells were moved to ice and scraped in 1mL ice cold PBS, and 100μL (∼5 x 10⁴ cells) were transferred to 1.5 mL tubes. Cells were centrifuged at 500 x g for 5 minutes at 4°C, and the pellets were resuspended in 50 μL cold lysis buffer (10mM Tris-HCl, 10mM NaCl, 3mM MgCl₂, 0.1% NP-40, 0.1% Tween-20, 0.01% Digitonin, adjusted to pH 7.4) and incubated on ice for 3 minutes. Samples were washed with 1 mL cold wash buffer (10mM Tris-HCl, 10mM NaCl, 3mM MgCl₂, 0.1% Tween-20). Cells were centrifuged at 500 x g for 10 minutes at 4°C, and pellets were resuspended in the transposition reaction mix (25 μL 2X TD buffer (Illumina), 2.5 μL TDE1 Tn5 transposase (Illumina), 16.5 μL PBS, 0.5 μL 1% Digitonin, 0.5 μL 10% Tween-20, 5 μL nuclease-free water) and incubated in a thermomixer at 37°C and 100rpm for 30 minutes. DNA was extracted with the DNA Clean and Concentrator-5 Kit (Zymo Research). Sequencing adapters were attached to the transposed DNA fragments using NEBNext Ultra II Q5 PCR mix (New England Biolabs), and libraries were amplified with 8 cycles of PCR. PEG-mediated size fractionation [73] was performed on the libraries by mixing SPRIselect beads (Beckman) with each sample at a 0.5:1 ratio, then placing the reaction vessels on a magnetic stand. The right side selected sample was transferred to a new reaction vessel, and more beads were added for a final ratio of 1.8:1. The final size-selected sample was eluted into nuclease-free water. This size selection protocol was repeated to further remove large fragments.

ATAC-seq analysis

Adapters were removed using cutadapt [67]. Reads aligning to the mitochondrial genome with bowtie2 [68] were removed. The remaining reads were aligned to the hg38 genome assembly with bowtie2. Duplicate reads were removed, and the remaining reads were sorted into BAM files with samtools [69]. Reads were converted to bed format with seqOutBias and bigWig format with deeptools [70, 74]. Accessibility peaks were called with MACS2 [35]. Reads were counted in peaks using bedtools, and differentially accessible peaks were identified with DESeq2 [37, 71]. de novo motif identification was performed on dynamic peaks with MEME, and TOMTOM was used to match motifs to the HOMER, Jaspar, and Uniprobe TF binding motif databases [75–78]. AME was used to identify motifs enriched in increased or decreased peaks relative to unchanged peaks [79]. FIMO and bedtools were used to assess motif enrichment around peak summits [80]. Dynamic peaks were clustered into response groups using DEGreport [81].

PRO-seq library preparation

Cell permeabilization was performed as previously described [82], with modifications. 4 replicates were performed from cells treated and collected at different times in the same day. For the time course experiment, 8*10⁶ dTAG-TRPS1 Clone 28 cells per sample were plated in 15cm dishes 2 days before harvest. Cells were treated with DMSO or 50nM dTAG-13 and 50nM dTAG^V-1 in DMSO at various time points and collected simultaneously. For the three clone experiment, 4*10⁶ cells per sample were plated in 10cm dishes 1 day before harvest. Cells were treated with DMSO or 100nM dTAG-13 in DMSO for 90 minutes and collected simultaneously.

At the time of harvest, cells were scraped in 10mL ice cold PBS and washed in 5mL buffer W (10mM Tris-HCl pH 7.5, 10mM KCl, 150mM sucrose, 5mM MgCl₂, 0.5mM CaCl₂, 0.5mM DTT, 0.004U/mL SUPERaseIN RNase inhibitor (Invitrogen), Complete protease inhibitors (Roche)). Cells were permeabilized by incubating with buffer P (10 mM Tris-HCl pH 7.5, KCl 10 mM, 250 mM sucrose, 5 mM MgCl₂, 1 mM EGTA, 0.05% Tween-20, 0.1% NP40, 0.5 mM DTT, 0.004 units/mL SUPERaseIN RNase inhibitor (Invitrogen), Complete protease inhibitors (Roche)) for 3 minutes on ice. Cells were washed with 10 mL buffer W before being transferred into 1.5mL tubes using wide bore pipette tips. Finally, cells were resuspended in 50μL buffer F (50mM Tris-HCl pH 8, 5mM MgCl₂, 0.1mM EDTA, 50% Glycerol, 0.5 mM DTT). Cells were snap frozen in liquid nitrogen and stored at -80°C.

PRO-seq libraries were prepared as previously described [83], with modifications. RNA extraction after the run-on reaction was performed with 500μL Trizol LS (Thermo Fisher) followed by 130μL chloroform (Sigma). The equivalent of 1μL of 50μM for each adapter was used. A random eight base unique molecular identifier (UMI) was included at the 5’ end of the adapter ligated to the 3’ end of the nascent RNA. 37°C incubations were performed with rotation with 1.5mL tubes placed in 50mL conical tubes in a hybridization oven. For the reverse transcription reaction, RP1 was used at 100μM and dNTP mix was used at 10mM each. Libraries were amplified by PCR for a total of 8 cycles in 100μL reactions with Phusion polymerase (New England Biolabs). No PAGE purification was performed to ensure that our libraries were not biased against short nascent RNA insertions.

PRO-seq analysis

Adapters were removed using cutadapt [67]. Libraries were deduplicated using fqdedup and the 3’ UMIs [84]. UMIs were removed, and reads were reverse complemented with theseqtk. Reads aligning to the rDNA genome with bowtie2 [68] were removed. The remaining reads were aligned, sorted, and convert to bed and bigWig files with bowtie2, samtools, seqOutBias, and deeptools, respectively [69, 70, 74]. Composite profiles around TRPS1 peaks were generated with deeptools. Reads were counted in genes using bedtools, and differentially expressed genes were identified with DESeq2 [37, 71]. Dynamic genes were clustered into response groups using DEGreport [81]. Over-representation analysis was performed with enrichr [85], and gene set enrichment analysis was performed with fgsea [86], both using the Hallmark gene sets [41].

Genome browser visualization

Genome browser [87] images were taken from the following session:https://genome.ucsc.edu/s/tgscott/dTAG_TRPS1_ChIP_PRO_ATAC.

Cell number enumeration

1.25*10⁴ cells per sample were plated in a 24-well plate. The next day (day 0), cells were treated with DMSO or 50nM dTAG-13 and 50nM dTAG^V-1 in DMSO. Media was replaced, maintaining the treatment condition, every 2 days. Cells were enumerated using a hemocytometer. 2 technical replicates were used on day 0, 2 for each treatment on day 3, and 3 for each treatment on days 7, 10, 14. The technical replicates were merged, and the experiment was performed in 4 biological replicates from different cell passages. The data were imported into R [62] for visualization and statistical analysis. A linear model was fit for the log-transformed cell number and the time. A second linear model was fit that included an interaction term between the time and the treatment condition, representing the effect of treatment on the doubling rate. Analysis of variance was performed on the two models to test for the significance of the interaction term.

TRPS1 activity score and patient outcome stratification

Primary TRPS1-regulated genes were defined based on the 30 minute time point using DESeq2 [37]. This TRPS1 regulon was then used in RTN [48, 49] to define a TRPS1 activity score for each patient within the METABRIC cohort [46, 47].