DOI: 10.5061/dryad.v41ns1s0s
title: ETV4 mediates dosage-dependent prostate tumor initiation and cooperates with p53 loss to generate prostate cancer

description: The mechanisms underlying ETS-driven prostate cancer initiation and
 progression remain poorly understood due to a lack of model systems that
 recapitulate this phenotype. We generated a genetically engineered mouse
 with prostate-specific expression of the ETS factor, ETV4, at lower and
 higher protein dosages through mutation of its degron. Lower-level
 expression of ETV4 caused mild luminal cell expansion without histologic
 abnormalities and higher-level expression of stabilized ETV4 caused
 prostatic intraepithelial neoplasia (mPIN) with 100% penetrance within 1
 week. Tumor progression was limited by p53-mediated senescence and Trp53
 deletion cooperated with stabilized ETV4. The neoplastic cells expressed
 differentiation markers such as Nkx3.1 recapitulating luminal gene
 expression features of untreated human prostate cancer. Single-cell and
 bulk RNA-sequencing showed stabilized ETV4 induced a novel luminal-derived
 expression cluster with signatures of the cell cycle, senescence, and
 epithelial to mesenchymal transition. These data suggest that ETS
 overexpression alone, at sufficient dosage, can initiate prostate
 neoplasia.
description: <strong>Mouse prostate digestion</strong>:
 Intraperitoneal injection of tamoxifen was administered in 8-week-old
 mice. 2 weeks after tamoxifen treatment, the mouse prostate was digested 1
 hour with Collagenase/Hyaluronidase (STEMCELL, #07912), and then 30
 minutes with TrypLE<sup>TM</sup> Express Enzyme (Thermo
 Fischer, # 12605028) at 37°C to isolate single prostate cells. The
 prostate cells were stained with PE/Cy7 conjugated anti-mouse CD326
 (EpCAM) antibody (BioLegend, 118216) and then, CD326 and EYFP double
 positive cells were sorted out by flow cytometry, which are luminal cells
 mainly from the anterior prostate and dorsal prostate. The mRNA or genomic
 DNA were extracted from these double-positive cells and then were used for
 ATAC-sequencing and RNA-sequencing analysis.
 <strong>ATAC-seq and primary data
 processing</strong>: ATAC-seq was performed as previously described.
 Primary data processing and peak calling were performed using ENCODE
 ATAC-seq pipeline (https://github.com/kundajelab/atac_dnase_pipelines).
 Briefly, paired-end reads were trimmed, filtered, and aligned against mm9
 using Bowtie2. PCR duplicates and reads mapped to mitochondrial chromosome
 or repeated regions were removed. Mapped reads were shifted +4/-5 to
 correct for the Tn5 transposase insertion. Peak calling was performed
 using MACS2, with p-value &lt; 0.01 as the cutoff. Reproducible peaks
 from two biological replicates were defined as peaks that overlapped by
 more than 50%. On average 25 million uniquely mapped pairs of reads were
 remained after filtering. The distribution of inserted fragment length
 shows a typical nucleosome banding pattern, and the TSS enrichment score
 (reads that are enriched around TSS against background) ranges between 28
 and 33, suggesting the libraries have high quality and were able to
 capture the majority of regions of interest.
 <strong>Differential peak accessibility</strong>:
 Reads aligned to peak regions were counted using R package
 GenomicAlignments_v1.12.2. Read count normalization and differential
 accessible peaks were called with DESeq2_v1.16.1 in R 3.4.1. Differential
 peaks were defined as peaks with adjusted p-value &lt; 0.01 and
 |log2(FC)| &gt; 2.  For visualization, coverage bigwig files were
 generated using bamCoverage command from deepTools2, normalizing using the
 size factor generated by DESeq2. The differential ATAC-seq peak density
 plot was generated with deepTools2, using regions that were significantly
 more or less accessible in ETV4<sup>AAA</sup> samples relative
 to EYFP samples. <strong>Motif
 analysis</strong>: Enriched motif was performed using MEME-ChIP
 5.0.0 with differentially accessible regions in
 ETV4<sup>AAA</sup> relative to EYFP. ATAC-seq footprinting was
 performed using TOBIAS. First, ACACCorrect was run to correct Tn5 bias,
 followed by ScoreBigwig to calculate footprint score, and finally
 BindDetect to generate differential footprint across regions.
 <strong>RNA-seq analysis</strong>: The extracted RNA
 was processed for RNA-sequencing by the Integrated Genomics Core Facility
 at MSKCC.  The libraries were sequenced on an Illumina HiSeq-2500 platform
 with 51 bp paired-end reads to obtain a minimum yield of 40 million reads
 per sample. The sequenced data were aligned using STAR v2.3 with GRCm38.p6
 as annotation. DESeq2_v1.16.1 was subsequently applied on read counts for
 normalization and the identification of differentially expressed genes
 between ETV4<sup>AAA</sup> and EYFP groups, with an adjusted
 p-value &lt; 0.05 as the threshold. Genes were ranked by
 sign(log2(FC)) * (-log(p-value)) as input for GSEA analysis using ‘Run
 GSEA Pre-ranked’ with 1000 permutations (48). The custom gene sets used in
 GSEA analysis are shown in Table S2.
 <strong>Unsupervised hierarchical
 clustering</strong>: To get an overall sample clustering as part of
 QC, hierarchical clustering was performed using pheatmap_v1.0.10 package
 in R on normalized ATAC-seq or RNA-seq data. It was done using all peaks
 or all genes, with Spearman or Pearson correlation as the distance metric.
 To have an overview of the differential gene expression from the RNA-seq
 data, unsupervised clustering was also performed on a matrix with all
 samples as columns and scaled normalized read counts of differentially
 expressed genes between ETV4<sup>AAA</sup> and EYFP as
 rows. <strong>Integrative analysis of ATAC-seq,
 RNA-seq, and ChIP-seq data</strong>: ERG ChIP-seq peaks were called
 using MACS 2.1, with an FDR cutoff of q &lt;
 10<sup>-3</sup> and the removal of peaks mapped to blacklist
 regions. Reproducible peaks between two biological replicates were
 identified as ETV4<sup>AAA</sup> ATAC-seq peaks. ERG ChIP-seq
 peaks and ETV4<sup>AAA</sup> ATAC-seq peaks were considered as
 overlap if peak summits were within 250bp. To determine whether the
 overlap was significant, enrichment analysis was done using
 regioneR_v1.8.1 in R, which counted the number of overlapped peaks between
 a set of randomly selected regions in the genome (excluding blacklist
 regions) and the ERG-ChIP seq peaks or ETV4<sup>AAA</sup>
 ATAC-seq peaks. A null distribution was formed using 1000 permutation
 tests to compute the p-value and z-score of the original
 evaluation. To assign ATAC-seq peaks to genes,
 ChIPseeker_v1.12.1 in R was used. Each peak was unambiguously assigned to
 one gene with a TSS or 3’ end closest to that peak. Differential gene
 expression between ETV4<sup>AAA</sup> and EYFP was evaluated
 using log2(FC) calculated by DESeq2. p-values were estimated with Wilcoxon
 rank t-test and Student t-test.
 <strong>scRNA-sequencing</strong>:
 <em>Tmprss2-CreER<sup>T2</sup>, EYFP;
 Tmprss2-CreER<sup>T2</sup>, ETV4<sup>WT</sup>;
 Tmprss2-CreER<sup>T2</sup>,
 ETV4<sup>AAA</sup></em>; and
 <em>Tmprss2-CreER<sup>T2</sup>,
 ETV4<sup>AAA</sup>; Trp53<sup>L/L</sup></em>
 mice were euthanized 2 weeks or 4 months after tamoxifen treatment (n=3
 mice for each genotype and time point). After euthanasia, the prostates
 were dissected out and minced with scalpel, and then processed for 1h
 digestion with collagenase/hyaluronidase (#07912, STEMCELL Technologies)
 and 30min digestion with TrypLE (#12605010, Gibco). Live single prostate
 cells were sorted out by flow cytometry as DAPI-. For each mouse, 5,000
 cells were directly processed with 10X genomics Chromium Single Cell 3’
 GEM, Library & Gel Bead Kit v3 according to manufacturer’s
 specifications. For each sample, 200 million reads were acquired on
 NovaSeq platform S4 flow cell. Reads obtained from the
 10x Genomics scRNAseq platform were mapped to mouse genome (mm9) using the
 Cell Ranger package (10X Genomics). True cells are distinguished from
 empty droplets using scCB2 package. The levels of mitochondrial reads and
 numbers of unique molecular identifiers (UMIs) were similar among the
 samples, which indicates that there were no systematic biases in the
 libraries from mice with different genotypes. Cells were removed if they
 expressed fewer than 600 unique genes, less than 1,500 total counts, more
 than 50,000 total counts, or greater than 20% mitochondrial reads. Genes
 detected in less than 10 cells and all mitochondrial genes were removed
 for subsequent analyses. Putative doublets were removed using the Doublet
 Detection package. The average gene detection in each cell type was
 similar among the samples. Combining samples in the entire cohort yielded
 a filtered count matrix of 48,926 cells by 19,854 genes, with a median of
 6,944 counts and a median of 1,973 genes per cell, and a median of 2,039
 cells per sample. The count matrix was then normalized to CPM (counts per
 million), and log<sub>2</sub>(X+1) transformed for analysis of
 the combined dataset. The top 1000 highly variable genes were found using
 SCANPY (version 1.6.1) (77). Principal Component Analysis (PCA) was
 performed on the 1,000 most variable genes with the top 50 principal
 components (PCs) retained with 29% variance explained.
 To visualize single cells of the global atlas, we used UMAP
 projections (https://arxiv.org/abs/1802.03426). We then performed Leiden clustering. Marker genes for each cluster were found with scanpy.tl.rank_genes_groups. Cell types were determined using the SCSA package, an automatic tool, based on a score annotation model combining differentially expressed genes (DEGs) and confidence levels of cell markers from both known and user-defined information. Heat-map were performed for single cells based on log-normalized and scaled expression values of marker genes curated from literature or identified as highly differentially expressed. Differentially expressed genes between different clusters were found using MAST package, which were shown in heat-map. The logFC of MAST output was used for the ranked gene list in GSEA analysis (48). The custom gene sets used in GSEA analysis are shown in Table S2. Gene imputation was performed using MAGIC (Markov affinity-based graph imputation of cells) package, and imputated gene expression were used in the heatmap. <strong>Analysis of public human gene expression datasets</strong>: To analyze <em>TP53</em> RNA expression in human prostate cancer samples, we obtained normalized RNA-seq data from prostate cancer TCGA (www.firebrowse.org) (3). To assess the role of TP53 loss on progression free and overall survival of ETS-positive prostate cancer, we used cBioportal. We queried MSKCC (queried on 10/12/2022), TCGA, and Stand-Up-To-Cancer cohorts, we selected patients with fusions of ERG, ETV1, ETV4, FL1, or TMPRSS2 and compared the overall survival between TP53 wild-type and TP53 mutant subsets. <strong>Statistics</strong>: All statistical comparisons between the two groups were performed by Graphpad Prism 7.0 software using a two-tailed unpaired <em>t</em>-test.
description: Please see the README document
 ("README_Dataset-ETV4_mediates_dosage-dependent_prostate_tumor_initiation.txt") and the accompanying published article: Dan Li, Yu Zhan, Naitao Wang, Fanying Tang, Cindy J. Lee, Gabriella Bayshtok, Amanda R. Moore, Elissa W.P. Wong, Mohini R. Pachai, Yuanyuan Xie, Jessica Sher, Jimmy L. Zhao, Makhzuna Khudoynazarova, Anuradha Gopalan, Joseph Chan, Ekta Khurana, Peter Shepherd, Nora M. Navone, Ping Chi, Yu Chen, ETV4 mediates dosage-dependent prostate tumor initiation and cooperates with p53 loss to generate prostate cancer. Science Advances 2023. Accepted. DOI: 10.5061/dryad.v41ns1s0s