Single-Cell Transcriptomics Analysis Identifies Nuclear Protein 1 as a Regulator of Docetaxel Resistance in Prostate Cancer Cells

Prostate cancer is the most frequently diagnosed cancer and the second leading cause of cancer-related death in men (1). Advanced disease is typically treated with androgen deprivation therapy (ADT); however, patients develop resistance to ADT in what is termed castration-resistant prostate cancer (CRPC; ref. 2). In 2004, docetaxel was approved for the treatment of CRPC (3, 4). However, 50% of patients do not respond to docetaxel therapy and the majority of those that do respond relapse within 3 years of initiating the therapy (3). Thus, understanding mechanisms through which docetaxel resistance develops is critical to improving therapy of CRPC.

Docetaxel has been shown to have two main mechanisms of antitumor activity; namely, microtubule depolymerization and inhibition of Bcl-2 expression (5, 6). However, multiple mechanisms of resistance have previously been identified. Structural changes to β-tubulin block docetaxel from affecting microtubules (7, 8). Prostate cancer cells can also upregulate Bcl-2 to overcome docetaxel-induced apoptosis (9). However, these previous studies of docetaxel resistance were performed on bulk populations of prostate cancer cells. It is unknown whether the resistant cell population develops from a subset of resistant cells already present in the parental population or develops de novo subsequent to therapy. To study the process of acquisition of resistant characteristics, different methods of investigation are necessary.

Advances in next-generation sequencing have allowed for the investigation at the single-cell level. Single-cell RNA-sequencing (scRNA-seq) allows researchers to investigate the variability and complex gene expression across all the individual cells, instead of a more homogeneous expression profile from traditional bulk RNA-seq of tissues. In prostate cancer, circulating tumor cells have been shown to have heterogeneous androgen resistance gene profiles within individual patients (10). This single-cell heterogeneity may play a role in the clinical and therapeutic heterogeneous response (11). In addition, with scRNA-seq it is possible to trace the origin of cancer cells to the cell type of origin (12). Similarly, scRNA-seq of docetaxel-sensitive versus -resistant cells may uncover mechanisms leading toward development of docetaxel resistance in prostate cancer.

In this study, we conducted scRNA-seq of docetaxel-sensitive and -resistant variants of the DU145 and PC3 prostate cancer lines. We analyzed the heterogeneity within and between the cell lines. We identified similar gene expression changes and functional gene pathways across both cell lines that are important to the acquisition of resistance. We also identified potential regulators of the resistant gene expression profile shared between both cell lines and identified a specific mediator of chemoresistance. Our findings uncover the heterogeneity in prostate cancer, as well as identified signaling pathways important for the acquisition of docetaxel resistance.

Parental (sensitive) DU145 and PC3 prostate cancer cells were obtained from ATCC. Docetaxel-resistant variants of the DU145 and PC3 prostate cancer cells have been described previously (13). Cell identification was confirmed every 6 months using short tandem repeat analysis. Cells were tested every 3 months for Mycoplasma. Cells were used within three passages from the time of thawing. Cells were cultured in RPMI1640 (Invitrogen Co.) supplemented with 10% FBS and 1% penicillin–streptomycin (Life Technologies, Inc.). Resistant DU145 and PC3 cells were maintained with the addition of docetaxel in DMSO to a final concentration of 10 nmol/L, while sensitive DU145 and PC3 cells were maintained with the addition of DMSO to a final level of 0.1%.

For 1 week, cells were transferred to docetaxel-free media. Cells were trypsinized in 0.05% Trypsin EDTA for 5 to 10 minutes at 37°C and washed with media. The cell suspension was loaded into in the Fluidigm C1 Machine and processed into single-cell cDNA libraries according to the manufacturer's protocol (PN 101-4981). Briefly, full-length mRNA-seq libraries were generated from single cells captured using the Fluidigm C1 Single Cell mRNA Seq IFC, 10–17 μm (PN 100-5760) and Fluidigm C1 Single-Cell Reagent Kit for mRNA seq (PN 100-6201). Each chip was visually inspected to identify which wells contained cells. Wells containing one cell were included in library preparation ad sequencing. The capture rate was between 78% to 96% across all chips used in this study. Full-length cDNA was converted into sequence ready libraries using SMART-seq v4 Ultra Low Input RNA Kit for Fluidigm C1 System (Takara Bio, catalog no. 635025) and SeqAmp DNA Polymerase (Takara Bio, catalog no. 638504). Library preparation was completed using Nextera XT DNA Library Prep Kit (Illumina, catalog no. FC-131-1096) and Nextera XT DNA Library Prep Index Kit (Illumina, catalog no. FC-131-1002). Samples following PCR reactions as called for in each kits' manufacturer's protocol were purified using Agencourt AMPure XP (Beckman Coulter, item no A63880). Samples were sequenced on Illumina HiSeq 2500 Rapid for DU145 cell line variants and Illumina HiSeq4000 with single end option for PC3 cell line variants. Reads that were below the minimum quality controls were discarded. Each sample was aligned to the Human Genome hg38 (14) using bowtie alignment tool (15). We captured a total of 324 cells across all cell lines. Poor quality cells were removed on the basis of low number of reads, as determined using the Fluidigm Singular package (https://www.fluidigm.com/software/; Supplementary Fig. S1). A total of 12 cells were removed. Sixty-four DU145 sensitive (DU145-S) cells, 71 DU145 resistant (DU145-R) cells, 89 PC3 sensitive (PC3-S) cells, and 88 PC3 resistant (PC3-R) cells were included in all downstream analysis. To identify genes for downstream analysis, we used the Fluidigm Singular Package. In brief, genes that were expressed in at less than 10 cells in each cell line were excluded. For the remaining genes, the lowest 15% of expressed genes were excluded. Finally, genes need to be identified in all four cell line samples to be included in the final gene list. This resulted in 12,862 genes being included for all the subsequent downstream analyses.

Dimensional reduction of cells was conducted using the t-distributed stochastic neighbor embedding method (tSNE; ref. 16). All high-quality genes were included for clustering unless noted. Individual clusters were identified using the method proposed in the Seurat pipeline (17). In brief, the distance metric for identifying clusters was based on principal components (PC) from principal component analysis. We next embedded the cells in a K-nearest neighbor graph structure and used the Louvain algorithm to iteratively group cells together. For the DU145-only cells, the top 10 PCs were used with a resolution of 0.1. For the PC3-only cells, the top five PCs were used with a resolution of 0.3. For all cells, the top 10 PCs were used with a resolution of 0.1.

Differential expression analysis was conducted on cells using the scDE R package (18). All genes were included in differential expression analysis. Genes with a FDR-adjusted P < 0.05 and a log fold change > 1 or <−1 were considered differentially expressed. To test for significance in overlap in Venn diagrams, we used a hypergeometric probability.

Functional enrichment analysis for cellular location and molecular function ontologies were conducted using PANTHER protein classifications (19, 20). Gene pathway analysis and upstream regulator analysis was conducted using Ingenuity Pathway Analysis (IPA) package (Qiagen; ref. 21). Potential regulators must have a P value less than 0.05 and an activation score above 2 or below −2 to be included in both analyses of DU145 and PC3 cell lines.

We used the web-based database called ONCOMINE (www.oncomine.org) to analyze public microarray cancer data provided from over 700 independent datasets (22). We analyzed a prostate cancer dataset containing 10 primary and 21 hormone refractory metastatic samples (23). Normalized intensities for nuclear protein 1 (NUPR1) for each sample were analyzed. Student t test was used to determine statistical significance.

Cells (2 × 10⁵) were plated with 1 mL of 10% RPMI medium without antibiotics in a 12-well plate. PC3-R and DU145-R cells were transfected with human NUPR1 siRNAs (catalog no. AM16708, Thermo Fisher Scientific) or Negative Control #1 siRNA (catalog no. AM4635, Thermo Fisher Scientific) by using Lipofectamine RNAiMAX (catalog no. 13778030, Thermo Fisher Scientific). PC3-S and DU145-S cells were transfected with p8 (NUPR1) (NM_001042483) Human Tagged ORF Clone (catalog no. RG222237, OriGene) or pCMV6-AC-GFP Vector (catalog no. PS100010, OriGene) by using Lipofectamine 2000 Reagent (catalog no. 11668030, Thermo Fisher Scientific) while the cells grow to 80% confluency. Twenty-four to 72 hours later, the validation of expected change in NUPR1 protein expression was followed by Western blot analysis.

The cells were lysed with RIPA Lysis Buffer (Millipore) containing a complete protease inhibitor tablet (Sigma). The protein concentration of tumor extracts was determined using BCA Protein Assay Kit (Pierce). Lysates were cleared by centrifugation at 14,000 × g for 20 minutes. The supernatant fractions were used for Western blot analysis. Protein extracts were resolved by 15% SDS-PAGE and probed with NUPR1(p8) Polyclonal Antibody (catalog no.PA1-4177, Thermo Fisher Scientific). Antibody to GAPDH (Millipore) was used as loading control. The antibody binding was revealed by using a horseradish peroxidase–conjugated anti-rabbit IgG (1:3,000, Cell Signaling Technology) or anti-mouse IgG (1:3,000, Amersham Pharmacia Biotech). Antibody complexes were detected by SuperSignal West Pico Chemiluminescent Substrate or SuperSignal West Dura Extended Duration Substrate (Thermo Fisher Scientific) and ChemiDoc Imaging System (Bio-Rad).

A total of 2 × 10⁴ PC3-S, PC3-R, DU145-S, and DU145-R cells were seeded in a 96-well plate. PC3-R and DU145-R cells were transfected with human NUPR1 siRNAs or negative control #1 by using Lipofectamine RNAiMAX. PC3-S and DU145-S cells were transfected with p8 (NUPR1) (NM_001042483) Human Tagged ORF Clone or pCMV6-AC-GFP vector by using Lipofectamine 2000 reagent. After 24 hours, the transfected cells were treated with 0, 10, 20, 40, 80, and 160 nmol/L docetaxel. The cells were cultured for 2 days, and then the number of viable cells was measured by proliferation reagent MTS (Promega) as directed by the manufacturer. Docetaxel IC₅₀ was determined by GraphPad Software.

Total RNA from prostate cancer cells was extracted using RNeasy Mini Kit (Qiagen), 2 μg of total RNA was reversed-transcribed in a final volume of 20 μL. qRT-PCR was performed in triplicate using SYBR Green qPCR MasterMix (Qiagen) in a 10 μL reaction volume on Applied Biosystems QuantStudio 5 Real-Time PCR Systems (Thermo Fisher Scientific). Human NUPR1 PKM2 primers (catalog no. PPH19840F) were purchased from Qiagen. All analyses were performed in triplicate in two independent experiments. Measurements from triplicate C_t values were normalized to GAPDH (Invitrogen).

The accession number for the scRNA-seq dataset reported in this article is GSE140440.

We have previously created and characterized docetaxel-resistant DU145 and PC3 prostate cancer cells (13). However, the underlying cellular heterogeneity associated with the development of drug resistance response was not investigated. To investigate the cellular heterogeneity of docetaxel resistance in prostate cancer, we conducted scRNA-seq of both the resistant variant (resistant) and parental (sensitive) prostate cancer cell lines (Fig. 1A). Using the Fluidigm C1 system, we isolated 312 total high-quality cells across all the cell lines. Across all four cell lines, we observed genes with higher expression also detected in a higher number of cells (Supplementary Fig. S2). The DU145 and PC3 variants, the sensitive and resistant variants, formed isolated clusters (Fig. 1B). In the DU145 variants, four sensitive cells clustered with the resistant population (Fig. 1B, right). In the few sensitive cells that clustered with the resistant cells in DU145, we observed similar number of total and mapped reads (Supplementary Fig. S3A–S3C). We do note the few resistant transcriptome-like–sensitive DU145 cells did have a lower number of unique genes and a higher number of unique transcripts compared with the remaining sensitive and resistant cells (Supplementary Fig. S3D and S3E). However, in the PC3 lines, none of the sensitive cells clustered with the resistant cells (Fig. 1B, left). We did observe gene expression differences between the sensitive and resistant cell lines in both DU145 and PC3 (Fig. 1C; Supplementary Table S1). To focus on specific candidates that could promote resistance, we examined the top 12 differentially expressed genes in either the DU145 or PC3 sensitive and resistant cells. In the DU145 cells, we observed upregulation of APOBEC3C expression in the resistant cells compared with sensitive cells (Fig. 1D, left). In addition, KRT8, KYNU, and ICAM1 had downregulation in the resistant cells compared with the sensitive cells (Fig. 1D, left). For the four DU145-S cells that clustered with the DU145-R cells (DU145 resistant transcriptome-like–sensitive cells), we investigated which genes drove the clustering results. We identified 1,788 genes expressed in all four DU145 resistant transcriptome-like–sensitive cells, of which 436 were also upregulated in the DU145-R cells in comparison with the DU145-S cells indicating approximately 25% of the upregulated genes in the resistant transcriptome-like–sensitive cells were genes also upregulated in the DU145-R cells. This large concordance may account, in part, for the reason these four sensitive DU145 cells clustered with the resistant cells. We additionally investigated whether the 436 overlapping genes were highly expressed in the DU145-R cells. We ranked all 12,862 genes in the DU145-R cells based on overall single-cell gene expression across all the DU145-R cells (Supplementary Table S2). We then noted where each of the 436 overlapping genes were located on the ranked list (Supplementary Fig. S3F; Supplementary Table S2). The majority of genes were found in the top 30% of the ranked genes, indicating they would play major role in the clustering algorithm. These results suggest resistant transcriptome-like–sensitive cells are present in the sensitive population but that they do not fully recapitulate the gene expression profile of the drug-resistant cells.

The two cell lines, although both being prostate cancer cell lines, expressed a different set of genes following development of drug resistance suggesting different genetic pathways of resistance (Fig. 1). To determine whether there were common changes in gene expression among the two different cell types at the single-cell level, we combined all 312 high-quality cells and subjected them to clustering analysis. In the resulting tSNE plot, each of the four cell lines clustered separately from each other (Fig. 2A). As observed previously, several sensitive cells clustered with their respective resistant counterpart in the DU145 cell lines. We additionally visualized the gene expression of the top 25 expressed genes in each cluster (Fig. 2B). To identify common changes in gene expression among all resistant cells compared with all the sensitive cells, we compared the differentially expressed genes in the DU145 and PC3 resistant cells versus their sensitive counter parts (Fig. 2C). Among all the differentially expressed genes between the two cell types, 182 of the same genes were differentially expressed in both comparisons. Of these, 139 were changed in the same direction; namely, four were upregulated and 135 were downregulated in the resistant versus sensitive cells; whereas, 43 were altered in different directions between the DU145 and PC3 resistant cells and their sensitive counterparts (Fig. 2C; P < 4.8e-9). Using only those 139 genes, we clustered all the cells again to see whether the cells would cluster based on their drug sensitivity as opposed to their cell line origin. The sensitive cells of both lines clustered together (Fig. 2D). However, the resistant cells of DU145 and PC3 still clustered separately except for three PC3-R cells that clustered with the sensitive cells (Fig. 2D). Because the changes in differential expression were in the same direction for both resistant cells lines, this observation suggests that there are still underlying differences in the gene expression other than regulatory changes that result in the two resistant cell lines clustering separate from each other while the two sensitive cell lines clustering together. To investigate this further, we investigated the gene expression levels of the top 12 expressed genes. In the first six genes in the violin plot, both resistant variants had higher expression compared with their sensitive counterparts (Fig. 2E). However, in the last six genes in the violin plots, the resistant variants had lower expression compared with their sensitive counterparts (Fig. 2E). Together, this suggests that while similar gene directional changes are induced by docetaxel resistance in both cell lines, there are still variation in the cells lines due to the differences in gene expression magnitude.

Having identified differentially expressed genes following docetaxel resistance, we investigated the functional implications of the differential expression. To do this, we analyzed the cellular location and molecular functions of the differentially expressed genes in either DU145 or PC3 using the PANTHER database (19, 20). The PANTHER database combines gene function classifications, pathways, and statistical analysis tools to enable researchers to explore large-scale sequencing experimental data. There were little differences in the cellular location of the proteins encoded for by the differentially expressed genes across the DU145 and PC3 cell lines (Fig. 3A; Supplementary Table S3). This suggests cellular location is similarly altered during drug resistance. In context of the molecular functions of the differentially expressed genes, we observed a higher percentage of genes related to catalytic activity in DU145 variant comparison than in the PC3 variant comparison (38% vs. 32%; Fig. 3B; Supplementary Table S4). However, a higher percentage of genes related to binding function were identified in the PC3 variant comparison than in the DU145 variant comparison (37% vs. 41%; Fig. 3B; Supplementary Table S4). In addition, when we examined the subclasses of the two largest molecular functions, we observed similar percentages (Fig. 3B; Supplementary Table S4). However, there were large differences in the percentage of genes upregulated or downregulated in the subclasses between the two cell lines. DU145 genes (1.3%) and 1.4% of PC3 genes were identified as having lyase activity. However, two of the 10 genes in DU145 were upregulated in the resistant cells; whereas, 10 of the 12 PC3 genes were also upregulated in the resistant cells (Fig. 3B; Supplementary Table S4). This suggests there are some underlying gene function differences and similarities important to the acquisition of drug resistance in either cell lines.

To investigate this further, we conducted causal network analysis using the IPA advanced analytics package (21) on the differentially expressed genes. A total of 162 pathways were identified across both cell lines (P < 0.05 in at least one comparison) as being differentially regulated in resistant versus sensitive cells lines, of which 38 were common to both cell lines (Fig. 3C). Of the top 10 pathways enriched in either cell line (based on P value), eight were common to both cell lines (Fig. 3D, shared pathways marked in red). The top differentially regulated pathway in both cell lines was protein ubiquitination (Fig. 3D). A total of 265 genes were associated with this pathway, in which 157 were identified in the DU145 comparison and 169 were identified in the PC3 comparison (Supplementary Table S5). Ubiquitination is a posttranslational modification, which can mark a protein as target for degradation. In addition, it is also associated with other functions such as cell signaling, cell cycle, or DNA repair (24). Thus our finding provides strong evidence that our analytic method is capable of detecting pathways of resistance. However, this analysis is limited in determining the enrichment direction of each pathway. This suggests the poorly studied pathways we have identified are good candidates for further study in relation to docetaxel resistance in prostate cancer.

With common gene pathways being dysregulated during drug resistance, we hypothesized these alterations were being driven by similar regulators. To test this, we conducted upstream regulator analysis using IPA to identify potential regulators (21). We identified 560 potential regulators with P < 0.05 in either cell line, of which 471 were common to both cell lines (Fig. 4A; Supplementary Table S6). A total of 121 were activated and 200 were inactivated in both comparisons (Fig. 4A). To further investigate these putative regulators of docetaxel resistance, we focused on those with an absolute z-score greater than 2. A total of 42 regulators met this criterion (Supplementary Table S6). Of the top 10 regulators (based on P value), 5 were common to both cell lines (Fig. 4B). Of these, four had activation in the same direction across both cell lines and one had different activation patterns in each cell line (Fig. 4B). TP53 was the top regulator with the most significant P value in both cell lines. This provides validation of our analysis method because TP53 has been shown to drive docetaxel resistance in prostate cancer (25). The second most significant upregulated regulator (based on P value) was NUPR1. NUPR1 has not previously been shown to play a role in docetaxel resistance in prostate cancer; whereas, it has been shown to promote resistance in breast cancer (26). To determine whether single-cell analysis offered an advantage to nominating NUPR1 as a putative candidate over RNA-seq of bulk cell, we analyzed bulk samples of each cell line variant. While we found NUPR1 to be upregulated in resistant cells, it was ranked as the 90th and 418th candidate gene in DU145 and PC3 cells, respectively, as opposed to third and second as found in single-cell analysis (Supplementary Fig. S4). This suggests, even though NUPR1 was identified in both bulk and single-cell samples, the single-cell samples provided stronger evidence for nominating NUPR1 as a candidate in the development of docetaxel resistance. We investigated further why NUPR1 ranked significantly higher as a mediator of resistance using the scRNA-seq compared with the bulk sample RNA-seq. This analysis was based on IPA regulator analysis, which does not look at the expression of the transcription regulators, but rather the expression of the regulator's targets to provide an approximation of the regulator's functional state. In many instances, transcriptional regulation occurs independently of the regulator's expression, but rather regulator activity is based on the protein conformation (e.g., HSP) or its localization (e.g., androgen receptor), which provides a rationale for IPA regulator analysis methodology. Thus, we looked at the average expression of NUPR1 pathway mediators in addition to the number of dropout events in both the single-cell and bulk-sequenced samples (Supplementary Fig. S5A and S5B). In the single-cell data, while not all gene targets were expressed in all of the cells, we observed higher average expression in many of the resistant cells compared with the sensitive cells (Supplementary Fig. S5A). However, in the bulk samples, although there was expression in all the samples, there was minimal differential expression between the resistant and sensitive variants (Supplementary Fig. S5B). This finding of higher overall differential expression of the NUPR1 pathway mediators in resistant versus sensitive cells in the single-cell analysis compared with bulk analysis accounts, in part, for the higher ranking of NUPR1 in the single-cell versus bulk analysis. The observations of inconsistency between pseudobulk constructed from the average of all of the cells from a scRNA-seq analysis versus true bulk is a common finding. In many cases, there is approximately 80% coefficient of determination (R²) between pseudobulk and true bulk in static cell samples (27, 28); however, in instances where a functional study is performed, such as a drug treatment, correlation may be as low as 2% (29). In our study, the R² between single-cell and bulk fold change was 44% for DU145 and 52% for PC3 (Supplementary Fig. S5C and S5D). This level of fit may be due to the relatively small sample size of single cells used. Thus, the observation that the fit between bulk and single cell is only between 44% and 55% could contribute to differences in clustering or ranking of target genes. Finally, to provide one more strategy to account for differences between single cell and bulk, we evaluated all of the genes that were statistically significantly differentially expressed in the single cells and determined whether they were also significantly differentially expressed in the bulk cells (Supplementary Fig. S5E). We found that of 12 genes significantly differentially expressed in both cell lines based on the single-cell data, three were not significant in both cell lines based on bulk; and that five were only significant in five of the cell lines; thus, overall there were only three (or 25%) of the genes significantly expressed in both bulk cell lines compared with that observed in 100% of the genes when identified in the cell lines in which both genes were differentially expressed on the basis of single-cell analysis. These data provide another rationale as to why NUPR1 was shown to rank higher in the single-cell versus bulk analysis.

To validate the prediction of the upstream regulator analysis, we analyzed mRNA and protein expression of NUPR1 in all cell lines. NUPR1 mRNA and protein expression were upregulated in both resistant cell lines compared with their sensitive counterparts (Fig. 4C and D). These findings were consistent with the increased NUPR1 activity predicted by the upstream regulator analysis. To provide some clinical correlation to these findings, we examined for NUPR1 expression in prostate cancers in the Oncomine database. We found that NUPR1 was upregulated in prostate cancer metastases compared with primary tumors (Fig. 4E).

To determine the potential functional role for NUPR1 in docetaxel resistance, we initially determined the impact of knocking down NUPR1 in the DU145-R cells. We observed a significant decrease in NUPR1 protein and mRNA expression following knockdown compared with control cells (Fig. 5A and B). We then evaluated the impact of NUPR1 knockdown on docetaxel resistance in the DU145-R cells. Knockdown of NUPR1 reduced the resistance to docetaxel of the DU145-R cells; albeit not to the level of DU145-S cells (Fig. 5C). This was reflected as a decreased IC₅₀ in the NUPR1-knockdown DU145-R cells compared with the control DU145-R cells, but not an IC₅₀ as low as that of the DU145-S cells (Fig. 5D). We also evaluated the impact of knockdown of NUPR1 on docetaxel resistance in the PC3-R cells and observed that similar to DU145 cells, knockdown of NUPR1 expression decreased the resistance and associated IC₅₀ to docetaxel in the PC3-R cells (Supplementary Fig. S6A–S6D). To provide further evidence for the ability of NUPR1 to promote docetaxel resistance, we overexpressed NUPR1 in the DU145-S cells. Both NUPR1 protein and mRNA expression were increased in the DU145-S cells transfected with NUPR1 compared with the cells transfected with empty vector (Fig. 5E and F). Overexpression of NUPR1 increased resistance to docetaxel in the DU145 cells; albeit, not to the level of DU145-R cells (Fig. 5G). This was reflected as a higher IC₅₀ for the DU145-S cells transfected with the NUPR1 expression vector versus the empty vector, but not as high as the IC₅₀ of the DU145-R cells (Fig. 5H). We also evaluated the impact of overexpressing NUPR1 on docetaxel resistance in the PC3-S cells and observed that similar to DU145 cells, overexpression of NUPR1 expression increased the resistance and associated IC₅₀ to docetaxel in the PC3-S cells (Supplementary Fig. S6E–S6H). Together, these data demonstrate that NUPR1 promotes docetaxel resistance in prostate cancer.

Defining mechanisms of chemotherapeutic resistance provides novel opportunities to overcome this resistance and provide a significant therapeutic impact. In this study, we investigated docetaxel resistance in prostate cancer at the single-cell level. We conducted scRNA-seq on two prostate cancer cell lines, DU145 and PC3, and their docetaxel-resistant variants. Our analysis revealed both independent and common gene changes in the two cell lines upon development of resistance. Furthermore, our strategy of using upstream regulator analysis using scRNA-seq data identified both a previously recognized mediator of prostate cancer chemoresistance (i.e., TP53) and a novel mediator of prostate cancer chemoresistance (NUPR1). We also identified that single-cell analysis outperformed bulk analysis for nominating NUPR1 as a putative mediator of chemoresistance.

Clinically observable drug resistance may occur through selection of preexisting resistant cells in the chemo-naïve tumor (30) or through development of de novo resistance through chemotherapy-induced genetic or epigenetic mechanisms (31, 32). Evidence exists for both theories across multiple cancer types (33, 34) and the mechanism may be due to the type of cancer and the type of chemotherapy. Our results suggest that a sensitive prostate cancer cell population contains a docetaxel-resistant subpopulation. However, we only observed this in one of the two cell lines we investigated. In addition, we did not see a continuum of cell transitions between the sensitive and resistant cells. Instead, the two subpopulations were discretely separated. This may be due to the fact that the resistant cells were selected in vitro over time and only represent the final resistant cells. Furthermore, the number of cells analyzed were relatively small compared with tumors in vivo, thus limiting our ability to evaluate the full spectrum of cells that could be present in a tumor, including limiting the ability to further determine whether these few cells were functionally resistant. To better understand how similar the resistant transcriptome-like–sensitive cells are to either treatment-naïve or drug-resistant cells, more single cells would need to be sequenced (35).

By subjecting the gene expression alterations common to both cell lines to upstream regulator analysis, we identified potential global regulators driving docetaxel resistance in prostate cancer. Some of the top regulators, based on P value in our analysis, have previously been shown to play a role in prostate cancer docetaxel resistance. For example, TP53 has been shown to drive docetaxel resistance in prostate cancer (25). In addition, MYC, the only regulator identified in both cell lines with negative activation (Fig. 4B), has been shown to play a role in docetaxel resistance in prostate cancer (36). The identification of a previously known regulator involved in docetaxel resistance suggests that the analytic method we applied to the scRNA-seq data is effective and may be useful to nominate additional putative regulators for further analysis. Furthermore, that NUPR1 was a higher candidate on the list based on the scRNA-seq compared with the bulk analysis indicates that there may be complementary value to both analytic methods.

On the basis of the upstream regulatory analysis, we investigated a potential role for NUPR1 in docetaxel prostate cancer resistance. Previously, NUPR1 was shown to overcome docetaxel treatment in breast cancer cells (37). However, while NUPR1 has been identified as a tumor suppressor in prostate cancer (26), its role in development of docetaxel resistance in prostate cancer was not previously described. Our observation from the Oncomine database that NUPR1 is upregulated in advanced prostate cancer provides support for the potential importance of NUPR1 in the development of docetaxel resistance in clinical prostate cancer. However, a limitation to these data is the treatment status of these patients is unknown. We hypothesized that the primary cancer is not treated with docetaxel; whereas, most hormone refractory metastatic cancers are treated with docetaxel. Thus, our finding that NUPR1 is higher in the metastatic versus primary is supportive but not conclusive. Regardless, the increase in patients with advanced prostate cancer indicates a potential role for NUPR1 in prostate cancer. This observation combined with our functional studies, demonstrating that NUPR1 promotes docetaxel resistance, provides a degree of confidence that NUPR1 could contribute to docetaxel resistance in patients with prostate cancer. Further research work with clinical samples of pre- and postdocetaxel treatment could provide additional evidence of the role of NUPR1 in prostate cancer.

Our in vitro studies demonstrated that NUPR1 contributed to docetaxel chemoresistance in two different cell lines. However, upon knockdown of NUPR1 activity, the resistant cells, which had diminished resistance, still had an IC₅₀ greater than sensitive cells. Conversely, upon overexpression of NUPR1 activity, the sensitive cells, which gained resistance, still had an IC₅₀ lower than the resistant cells. Taken together, these results demonstrate that NUPR1 contributes to docetaxel resistance, but indicates other factors also contribute to docetaxel resistance. NUPR1 has been previously implicated in multiple pathways that could aid in drug resistance. In breast cells, NUPR1 regulates p21 through formation of a complex of p53 to impart a cell growth and survival advantage (37). NUPR1 also plays a role in TFGβ signaling, which can affect drug resistance (26). However, further research is necessary to determine the signaling pathways through which NUPR1 contributes to the development of resistance.

In conclusion, our study provides a view of the development of docetaxel resistance in prostate cancer at the single-cell level. Furthermore, it provides a strategy to exploit scRNA-seq data to nominate candidates for drug resistance. While each individual cancer may acquire resistance through independent means, it is likely there is a common subset of alterations that occurs across tumors from different patients that are important to the acquisition of resistance. Through focused analysis of these common genetic changes and their regulatory factors, we identified a novel mediator of docetaxel resistance in prostate cancer which offers a new therapeutic candidate to extend the effectiveness of docetaxel therapy for patients with prostate cancer.

P.M. Schnepp reports grants from the National Center for Advancing Translational Sciences during the conduct of the study. G. Shelley reports grants from NIH during the conduct of the study. H. Jiang reports grants from NIH during the conduct of the study. E.T. Keller reports grants from NIH during the conduct of the study. No potential conflicts of interest were disclosed by the other authors.

P.M. Schnepp: Formal analysis, writing-original draft, writing-review and editing. G. Shelley: Formal analysis, investigation, methodology, writing-original draft. J. Dai: Investigation. N. Wakim: Formal analysis. H. Jiang: Formal analysis, supervision. A. Mizokami: Conceptualization, resources, writing-review and editing. E.T. Keller: Conceptualization, supervision, funding acquisition, writing-review and editing.

This work was supported by NIH grants P01093900 (to E.T. Keller); P30CA046592 (to E.T. Keller); NCATS grant UL1TR002240 (to P.M. Schnepp); NIH NHGRI Genome Science Training Grant T32HG00040 (to N. Wakim); and P30CA046592 (to H. Jiang). We would like to thank the University of Michigan's Advanced Research Computing Technology Services for providing the computing resources. Research reported in this article was supported by the NCI of the NIH under Award Number P30CA046592 by the use of the following Cancer Center Shared Resource that supported this study: Rogel Cancer Center Single Cell Analysis Shared Resource (P30 CA046592). This work was supported by the NCI of the NIH under Award Numbers P01CA093900 to E.T. Keller.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.