Stromal androgen-receptor (AR) action is essential for prostate development, morphogenesis and regeneration. However, mechanisms underlying how stromal AR maintains the cell niche in support of pubertal prostatic epithelial growth are unknown. Here, using advanced mouse genetic tools, we demonstrate that selective deletion of stromal AR expression in prepubescent Shh-responsive Gli1-expressing cells significantly impedes pubertal prostate epithelial growth and development. Single-cell transcriptomic analyses showed that AR loss in these prepubescent Gli1-expressing cells dysregulates androgen signaling-initiated stromal-epithelial paracrine interactions, leading to growth retardation of pubertal prostate epithelia and significant development defects. Specifically, AR loss elevates Shh-signaling activation in both prostatic stromal and adjacent epithelial cells, directly inhibiting prostatic epithelial growth. Single-cell trajectory analyses further identified aberrant differentiation fates of prostatic epithelial cells directly altered by stromal AR deletion. In vivo recombination of AR-deficient stromal Gli1-lineage cells with wild-type prostatic epithelial cells failed to develop normal prostatic epithelia. These data demonstrate previously unidentified mechanisms underlying how stromal AR-signaling facilitates Shh-mediated cell niches in pubertal prostatic epithelial growth and development.
Androgen- and androgen receptor (AR)-mediated signaling is essential for male sexual development (Cunha et al., 1987; Cooke et al., 1991). Murine prostatic development initiates at embryonic day (E) 17.5 from the endodermal urogenital sinus (UGS) in response to rising levels of testicular androgens (Cunha et al., 1987). Mutation of the Ar gene in testicular feminized (Tfm) mice results in the complete absence of prostate development (Cunha and Chung, 1981). Earlier tissue recombination studies using combinations of wild-type or AR-deficient Tfm urogenital sinus mesenchyme (UGM) with urogenital sinus epithelium (UGE), demonstrated that mesenchymal, rather than epithelial, AR plays a decisive role in early prostatic epithelial development (Cunha and Lung, 1978; Cunha, 1984). This provided the first evidence of stromal AR expression acting to maintain the cellular niche responsible for supporting prostate stem- and progenitor-mediated prostatic epithelial development. After birth, the AR continues to play a crucial role in prepubescent prostatic gland branching and morphogenesis, as well as in pubertal prostatic growth and maturation (Sugimura et al., 1986; Cunha et al., 1987; Hayashi et al., 1991). However, owing to the technical limitation of previous tissue recombination models only being able to use embryonic tissues, the role of stromal AR in the cell niche responsible for supporting prostatic stem/progenitor cells in pubertal prostatic epithelial morphogenesis and growth have not been appropriately characterized.
The Shh signaling pathway plays a crucial role in regulating both embryonic and adult stem cells in a variety of self-renewing organs (Ahn and Joyner, 2005; Kugler et al., 2015). In the prostate, the Shh ligand and its downstream effector, Gli1, are exclusively expressed in epithelial and mesenchymal cells, respectively (Doles et al., 2006; Peng et al., 2013). Upon activating Shh signaling, prostatic stromal Gli1-expressing cells are able to repopulate the stroma following involution and regeneration (Peng et al., 2013). Shh signaling also regulates prostatic epithelial development, morphogenesis and regeneration through mesenchymal-epithelial interactions (Peng and Joyner, 2015; Bushman, 2016). However, Gli1 expression has been shown not to be essential for mouse prostate development and morphogenesis (Doles et al., 2006; Peng et al., 2013). Interestingly, selective deletion of AR expression in mesenchymal Shh-responsive Gli1-expressing cells significantly impaired early prostate development, morphogenesis and growth (Le et al., 2020). Although this study indicates a crucial role of AR in Gli1-expressing cells and their lineage, the function and identity of AR- and Gli1-expressing cells in the prostate stem cell niche and related mechanisms remain poorly understood.
In this study, we directly investigated the cellular properties of stromal AR and Gli1-expressing cells and their roles acting as the cellular niche in regulating prepubescent and pubertal prostate morphogenesis and growth. Single-cell transcriptomic analyses showed that selective deletion of AR expression in prepubescent Gli1-expressing cells dysregulates androgen signaling-initiated stromal-epithelial paracrine interactions, resulting in robust defects in pubertal prostate growth and development. Specifically, stromal AR loss elevates Shh-signaling activation in both prostatic stromal cells and adjacent epithelial cells in a paracrine fashion, resulting in stark inhibition of prostatic epithelial growth and pubertal development. Single-cell pseudotime-trajectory analyses further identified aberrant differentiation fates of prostatic epithelial cells altered by stromal AR deletion. In vivo recombination of AR-deficient stromal Gli1-lineage cells with wild-type prostatic epithelial cells failed to develop normal prostatic epithelia. These lines of scientific evidence demonstrate novel mechanisms underlying stromal AR signaling in facilitating Shh-mediated cell niches in pubertal prostatic epithelial growth and development.
Deletion of the AR in prepubescent Gli1-expressing cells diminishes pubertal prostate growth
An indispensable role of AR in Gli1-expressing cells during embryonic prostate development, prepubescent morphogenesis and pubertal growth has been identified recently (Le et al., 2020). To directly address the unknown mechanisms underlying stromal AR in prepubescent Gli1-expressing cells acting as the cellular niche that supports prostatic pubertal epithelial development, we profiled the single-cell transcriptomes from prostate tissues of ROSA26RmTmG/+:ArLoxP/Y:Gli1CreERT2/+ (R26mTmG/+:ArL/Y:Gli1CreER/+) mutant (ARKO) mice and ROSA26RmTmG/+:Gli1CreERT2/+ (R26mTmG/+:Gli1CreER/+) control littermates. Tamoxifen (TM) was administered at postnatal day (P) 14 to activate CreER-mediated recombination of the floxed loci, switching from expression of membrane-bound tandem dimer Tomato (mT) to membrane-bound green fluorescent protein (mGFP) and/or deleting AR expression in Gli1-expressing cells in control and ARKO mice, respectively (Muzumdar et al., 2007; Le et al., 2020) (Fig. 1A). These mice were then grossly examined at P17, P35 and P56 during prostate prepubescent and pubertal growth and development (Fig. 1B). At P17, there was no significant difference between the prostates isolated from ARKO and control mice. At P35 and P56, the prostates of the control mice showed normal structure and a significant increase in the weight and size between these time points, reflecting normal pubertal prostate growth in response to elevated androgens (Fig. 1D1,D2; Fig. S1C, left panel). In contrast, the prostates from ARKO mice appeared to be significantly smaller in size and weight than their control counterparts at P35 and P56, and showed almost no difference between these two time points (Fig. 1D3,D4; Fig. S1C, left panel), demonstrating a significant androgen-dependent growth defect. Abnormally small seminal vesicles were also observed in ARKO mice at the above time points (Fig. S1A,B). Histological analyses of prostate tissues from control mice showed normal prostatic epithelia with mature glands filled with prostatic secretions by P56 (Fig. 1E1,E1′,E2,E2′). However, noticeable pathological changes featuring prostatic epithelial growth retardation and aberrant prostatic morphogenesis of all four prostatic lobes, including a reduction in gland size with more dense stromal tissue surrounding the glands and irregular glands with enlarged and bulbous tips, appear at the onset of puberty in P35 ARKO samples (Fig. 1E3,E3′). These abnormalities became even more striking at P56, highlighted by significantly undersized prostatic glands, devoid of secretions and with irregularly shaped vacuolated luminal epithelial cells (Fig. 1E4,E4′). Examining the levels of serum testosterone showed no significant differences between samples collected from ARKO and control littermates. Both genotypes had lower serum testosterone at P35 than at P56, demonstrating normal pubertal increases in testosterone (Fig. S1C, right panel), ruling out the potential effect of insufficient circulating androgens in the above prostatic defects. Analyses of Gli1-activated mGFP expression showed positive staining exclusively in stromal cells adjacent to the epithelium at all time points in samples of both genotypes, confirming the stromal cell properties of Gli1-expressing cells and their descendants, further referred to as Gli1-lineage cells (Fig. 1F1-F6). The above data demonstrate a crucial role for AR expression in Gli1-lineage cells in supporting androgen-dependent prostate growth during puberty, which are also consistent with our previous report (Le et al., 2020). In addition, these data also corroborate the use of P35 mouse prostate tissues to characterize stromal AR action and related mechanisms via single-cell RNA-sequencing (scRNA-seq) and other analyses (Fig. 1C).
Deletion of AR in prepubescent Gli1-expressing cells alters pubertal prostatic cell populations
We performed scRNA-seq analysis using mouse prostate tissues isolated at P35. In one of the representative experiments, ∼8000 and 9500 cells isolated from R26mTmG/+:ArL/Y:Gli1CreER/+ and R26mTmG/+:Gli1CreER/+ littermates, respectively, were prepared and sequenced. Following sequencing and alignment to the mm10 reference genome with an added mGFP sequence (Zhang et al., 1996), both the ARKO and control samples underwent multiple steps of filtering (see Materials and Methods). Post-filtering, 6871 and 7148 cells from ARKO and control samples, respectively, with an average of 2454 genes and 11,664 unique molecular identifier (UMI) counts per cell were used in the analyses (Fig. S2).
The scRNA-seq data from both ARKO and control samples were initially visualized individually using t-distributed stochastic neighbor embedding (tSNE) (Fig. 2A), and then merged and re-clustered based on their transcriptomic profiles using Seurat's integrated merging method (Stuart et al., 2019) (Fig. 2B). A total of seven epithelial and seven stromal (non-epithelial) cell subsets were identified (Fig. 2C), and aligned in separated ARKO and control tSNE plots (Fig. 2D). These data demonstrate comparable overall cell populations sequenced between these two samples, and support further analysis to assess the molecular changes induced by AR loss between ARKO and control samples at P35 using these identified cell clusters. The cellular properties of each cell subset were determined and then grouped as epithelial and stromal cell populations based on the cell-specific markers reported in previous human and mouse scRNA-seq and in-situ hybridization analyses (Wang et al., 2001; Toivanen et al., 2016; Henry et al., 2018; Kwon et al., 2019; Noda et al., 2019; Guo et al., 2020; Joseph et al., 2020) (Fig. 2E). Specifically, the epithelial cell subsets included prostate luminal (LE) and basal (BE), urethral luminal (UrLE), seminal vesicle (SV), myoepithelial (MyoE), proliferating epithelia (ProE) and a small epithelial cell cluster of unclear origin labeled as other epithelia (OE). Stromal cell subsets consisted primarily of fibroblasts (Fb) and smooth muscle cells (SMC) as well as endothelial (Endo), pericyte (PC), proliferating fibroblast (ProS) and two leukocyte (Leu1 and 2) clusters. In tSNE expression plots, the expression of Epcam, Vim and Myh11 was visualized within epithelial, stromal and SMC subsets as defined above (Fig. 2C-F). Ar expression was detected in both prostatic epithelial and stromal cell subsets (Fig. 2F). However Gli1 and mGFP expression was observed selectively within prostatic fibroblasts and SMCs, demonstrating the stromal cell properties of Gli1-lineage cells in pubertal prostate tissues (Fig. 2F). These high-resolution data validate the cellular properties of AR and Gli1-lineage cells in pubertal prostate tissues. Analyses of total cell distributions for each cell type in ARKO and control samples showed basal and luminal epithelial cells along with SMCs and fibroblasts being the four most abundant cell types (Fig. 2G), consistent with the cell composition of prostate tissues observed histologically (Fig. 1E3-4, 2H1-6). Most epithelial cell subsets, except basal cells, significantly decreased in ARKO samples in comparison with the controls. In contrast, a proportional increase in stromal SMCs and fibroblasts appeared in ARKO samples due to the large reduction in epithelial cells compared with control counterparts, resembling the reduced overall gland size and prostatic epithelia with relative dense stroma observed in P35 ARKO prostate tissue sections (see Fig. 1E3). Although observations of robust growth retardation in prostatic epithelial cells in ARKO samples demonstrate a dominant role of stromal AR in prostate epithelial growth and proliferation, reduced proliferation was also observed in stromal cell populations in ARKO samples (Fig. 2G).
Using co-immunofluorescence (Co-IF) analyses, we further determined stromal cellular properties of Gli1-lineage cells and the effect of AR deletion in prostate cell proliferation. In both ARKO and control prostate tissues, Gli1-CreER-activated mGFP-positive cells show negative E-cadherin staining (Fig. 2H1,H1′,H4,H4′) but positive staining for smooth muscle actin (SMA) or Vim (blue arrows, Fig. 2H2-H3′,H5-H6′). In addition, a stark decrease in Ki67- and PCNA-positive cells was revealed in ARKO samples in comparison with the controls (Fig. 2I1-6). These data support the above scRNA-seq data, further demonstrating the stromal cell properties of pubertal Gli1-lineage cells and the growth-promoting role of stromal AR in this cell lineage.
Stromal AR expression regulates Shh signaling in prostatic stromal cells
Given the significant prostate growth defects observed in pubertal prostate tissues in stromal AR-deficient R26mTmG/+:ArL/Y:Gli1CreER/+ mice, we re-clustered stromal cells after separating them from the epithelial cells to gain higher resolution insight into their cell properties. Eleven cell clusters were identified upon re-clustering (Fig. 3A,B), with increased resolution leading to the further separation of fibroblasts and SMCs into multiple clusters (Fig. 3B). Because Gli1-lineage cells are restricted within SMC and fibroblast cell subsets, we then further separated these two cell subsets from other stromal cells yielding two SMC and four fibroblast subclusters based on their transcriptional profiles (Fig. 3C; Fig. S3A,B; Table S1). Despite more in-depth cell clustering, the newly identified cell subclusters aligned between ARKO and control samples, with no populations specific to only one genotype (Fig. 3D). Analysis of cell distributions showed that fibroblast cluster 4 (Fb4) and SMC cluster 2 (SM2) were proportionally larger in controls, whereas Fb1-3 clusters were larger in ARKO samples (Fig. 3E). Gli1-CreER-mediated mGFP expression appeared to be consistent in fibroblast and SMC clusters in both ARKO and control samples; however, reduced expression of mGFP was observed in the Fb3 cluster, especially in ARKO cells, suggesting divergent cellular properties between these Gli1-lineage cells (Fig. 3F, left panel). In alignment with expectations for our model, a significant decrease in Ar expression appeared in all mGFP-expressing stromal cell clusters of ARKO samples (Fig. 3F, right panel). Analyses of tSNE expression plots further showed the expression of Ar and Gli1 in stromal cell clusters. Interestingly, we noted elevated Gli1 expression in ARKO cells (Fig. 3G), suggesting an inverse relationship in expression between Ar and Gli1 in those prostate stromal cells. We further validated the above observation using Co-IF and immunohistochemistry (IHC) analyses in both pre-pubertal and pubertal prostate samples. Uniform nuclear AR staining in mGFP-positive cells was detected within prostatic stromal areas in all control samples isolated at P17, P35 and P56 time points (Fig. 3H1-H3, blue arrows). In contrast, reduced, or nearly full loss of, AR expression at the three time points appears selectively in mGFP-positive cells of ARKO samples (Fig. 3J1-J3, yellow arrows). IHC analyses also showed Gli1 expression restricted to stromal cells of both ARKO and control samples at all time points (Fig. 3I1-I3,K1-K3, arrows). Importantly, increased Gli1 expression was detected in ARKO samples compared with control counterparts, particularly at P35 and P56 time points (Fig. 3K1-K3 versus I1-I3), further demonstrating the inverse relationship in expression of AR and Gli1 proteins. In agreement with the above observation, an increase in the number of Gli1-positive, Ar-negative cells was identified in the ARKO stromal cells (Fig. 3L), accompanied by an overall increase in Gli1 expression in ARKO stromal cell populations in comparison with the control samples. These multiple lines of evidence demonstrate an inverse relationship between AR and Gli1 expression in prostatic stromal cells, implicating a new potential mechanism by which stromal AR expression in Gli1-lineage cells regulates pubertal prostate growth and development.
Stromal AR signaling regulates developmental signaling pathways in Gli1-lineage cells
Identifying an inverse relationship between AR and Gli1 expression in prostatic Gli1-lineage cells is new and intriguing. Using higher read-depth bulk RNA-seq, we further assessed the transcriptomic profiles of Gli1-lineage cells to investigate the regulatory role of stromal AR in Gli1-lineage cells. Analyses were performed using RNA samples that were prepared from sorted mGFP-positive cells of R26mTmG/+:ArL/Y:Gli1CreER/+ and R26mTmG/+:Gli1CreER/+ mice at P35 with TM administration at P14 (Fig. 4A). A total of 731 differentially expressed genes (DEGs) (391 up- and 340 downregulated) were identified comparing ARKO with control samples (Fig. 4B; Table S2). Gene set enrichment analysis (GSEA) using pre-ranked gene lists identified an enrichment of a variety of up- and downregulated signaling pathways altered by AR deletion in sorted Gli1-lineage cells (Fig. 4C). Among them, upregulation of Hedgehog signaling, epithelial-to-mesenchymal transition (EMT), TGF-β, Wnt, focal adhesion, calcium and Erbb pathways, and downregulation of interferon α and γ response, oxidative phosphorylation, Kras, DNA repair and estrogen late response-related pathways were identified in ARKO cells. Quantitative reverse-transcriptase real-time PCR (RT-qPCR) analyses also showed reduced Ar expression alongside increased expression of Shh downstream targets, Gli1, Ptc1 (also known as Ptch1) and Smo in ARKO Gli1-lineage cells in comparison with the controls (Fig. 4D). These data are consistent with previous results (see Fig. 3G,I,K,L) and demonstrate an upregulation of Shh signaling in ARKO prostatic Gli1-lineage cells. In this study, we also performed both chromatin immunoprecipitation (ChIP)-qPCR and ChIP-seq analyses to examine the downstream targets of AR using sorted mGFP-positive cells isolated from ARKO and control prostate tissues. However, we did not identify a direct role of AR on Gli1 transcription (see the Discussion).
To gain in-depth mechanistic insight into stromal AR action, we characterized single-cell transcriptomes of Gli1-lineage fibroblasts and SMCs separately in ARKO and control samples. DEGs detected in more than 5% of cells with a P-value <0.05 and average log fold change >0.25 were detected using a Wilcoxon Rank Sum test by comparing either mGFP-positive fibroblasts or SMCs, referred to as Gli1-lineage cells as described earlier, between ARKO and control samples (Fig. S4A; Tables S3 and S4). GSEA using the above DEGs identified both positively and negatively enriched signaling pathways comparing ARKO with control Gli1-lineage fibroblasts or SMCs (Fig. 4E,F). Intriguingly, many significantly altered signaling pathways identified via bulk RNA-seq analyses were repeatedly detected in these analyses (Fig. 4C). Specifically, upregulation of Shh-related target signaling pathways as well as other developmental and metabolic stress-related pathways, including Wnt and hypoxia, was observed in both Gli1-lineage fibroblasts and SMCs from ARKO samples. In addition, oxidative phosphorylation, extracellular matrix receptor interactions and innate immune response-related pathways were consistently downregulated in those ARKO fibroblast and SMC cells (Fig. 4E,F). GSEA plots further showed comparable enrichment of Shh target genes in both ARKO Gli1-lineage SMCs and fibroblasts (Fig. 4G). The upregulation of Shh downstream targets, including Dcn, Cdkn1a and Mknk2, was observed in both ARKO Gli1-lineage fibroblasts and SMCs in comparison with controls from the above gene sets (Fig. 4H). In addition, a subset of Shh downstream targets were regulated selectively either in ARKO fibroblasts, such as Lum, Igfbp3 and Apod, or in SMCs, Filip1l and Cystm1 (Fig. 4H). Taken together, these various analyses demonstrate upregulated Shh signaling in both ARKO Gli1-lineage SMCs and fibroblasts at single-cell resolution.
Stromal AR signaling in Gli1-lineage cells regulates pubertal prostatic epithelial cell growth through Shh signaling
It has been implicated that Shh signaling regulates prostatic epithelial development and morphogenesis through mesenchymal-epithelial interactions (Peng and Joyner, 2015; Bushman, 2016). To directly assess the effect of AR loss in stromal Gli1-lineage cells on prostate epithelia, we re-clustered prostatic epithelial cell subsets after excluding non-prostatic epithelial cells from the initial scRNA-seq analyses (Fig. 5A,B). Seven epithelial cell clusters were identified based on their single-cell transcriptomes (Fig. S5A,B), which included two basal (BE), three luminal (LE), one urethral luminal (UrLE) and one proliferating epithelial cell cluster (ProE) (Wang et al., 2001; Toivanen et al., 2016; Noda et al., 2019; Joseph et al., 2020) (Fig. 5C). A significant decrease in cell numbers of LE3 and ProE clusters were identified in ARKO samples compared with the controls in tSNE plots (Fig. 5D, red arrows). Analysis of prostatic epithelial cell distributions consistently showed a reduction of luminal epithelial clusters in the ARKO sample, with almost a complete absence of LE3 and ProE (Fig. 5E). In contrast, a relative increase in basal cell populations was revealed in the above ARKO samples (Fig. 5E). Using tSNE expression plots, we further assessed the cellular properties of prostatic epithelial cells at single-cell resolution (Fig. 5F). The expression of Ar was identified in all epithelial clusters, but at lower levels in BE clusters (Fig. 5F). The expression of Pbsn, an AR downstream gene, was mainly detected in LE clusters, with elevated expression in LE2 and LE3, implying their highly differentiated and mature luminal cell properties. The proliferation marker Ki67 (Mki67) gene was mainly localized in ProE, and the expression of Shh was identified within BE clusters (Fig. 5F). Genes highly specific to each epithelial cluster were also assessed (Fig. S5A,B).
The above data provide a high-resolution depiction of cellular changes in pubertal prostatic epithelial cells altered by stromal AR deletion. Specifically, observation of restricted expression of Shh in prostatic basal cells is intriguing. Prostatic basal epithelial cells physically lie between the luminal epithelial and stromal cells, and possess prostatic progenitor properties capable of luminal epithelial differentiation and growth (Ousset et al., 2012; Toivanen et al., 2016). To determine the paracrine effects of stromal AR deletion on prostatic epithelial development defects, we assessed the transcriptomic changes of basal epithelial cells between ARKO and control samples. In total, 467 DEGs were identified in the analyses, of which 167 are upregulated and 300 downregulated (Fig. S4B; Table S5). GSEA based on the above DEG lists identified upregulated prostatic developmental pathways (Hedgehog, Tgf-β, Wnt and Notch pathways) in ARKO BE cells (Fig. 5G), corresponding to similar changes observed in ARKO stromal cells (Fig. 4F,G). In addition, oxidative phosphorylation, interferon α response, and epithelial morphogenesis and differentiation pathways appeared to be downregulated in the ARKO BE cells (Fig. 5G). Observations of downregulated metabolic pathways along with upregulated hypoxia signaling in the analysis suggest a potential regulatory mechanism by which stromal AR loss induces metabolic stress on prostatic basal epithelial cells during pubertal development. In addition, decreased expression of Hoxb13, Etv4, Mmp14 and Bmp7, crucial regulators of branching morphogenesis and prostatic differentiation, and elevated expression of Shh and Hes1, a Notch downstream effector (Su and Xin, 2016), were detected in ARKO prostatic basal cells (Fig. 5H).
It has been shown that stromal AR-expressing cells can produce a variety of growth factors and cytokines, termed andromedins, to regulate prostate epithelial cells in a concentration gradient-dependent manner (Yan et al., 1992; Brennen and Isaacs, 2018). Given our observations of elevated Shh transcripts in prostatic basal cells adjacent to stromal ARKO Gli1-lineage cells, we performed IHC analyses in prostate tissues isolated from both ARKO and control samples to assess corresponding protein levels. Decreased expression of AR was revealed in stromal cells of ARKO prostate tissues in comparison with the controls, confirming specific AR deletion in prostatic stromal cells (Fig. 5J1,J1′ versus I1,I1′, arrows). Again, more intense staining of Gli1 was observed in stromal cells from the above ARKO samples compared with controls (arrows in Fig. 5J2,J2′ versus I2,I2′). Furthermore, consistent with scRNA-seq data, increased staining of Shh was observed in basal epithelial cells adjacent to stromal Gli1-lineage cells in ARKO tissues in comparison with controls (Fig. 5J3,J3′ versus I3,I3′, arrows), demonstrating a corresponding increase in Shh signaling in stromal and epithelial cells in ARKO prostatic tissues. In this study, using ex vivo culture approaches, we further identified the inhibitory effect of Shh on prostate epithelial growth in prostate samples isolated from ARKO or control littermates. Of note, an apparent partial rescue in growth of ARKO samples was observed when treated with cyclopamine, an Shh inhibitor in the experiments (Fig. S6A-G). These results directly demonstrated that aberrant Shh activation inhibits prostatic epithelial growth, which is consistent with previous studies (Freestone et al., 2003; Wang et al., 2003; Berman et al., 2004; Pu et al., 2004). Taken together, these multiple lines of evidence implicate a regulatory mechanism by which stromal AR in Gli1-lineage cells governs pubertal prostate epithelial growth and development through Shh signaling pathways.
Stromal AR deletion in Gli1-lineage cells regulates the differentiation and fates of prostatic epithelial cell lineages
To determine the effect of stromal AR in Gli1-lineage cells in prostate epithelial differentiation, we applied algorithms to construct dynamic cellular trajectories, unveiling unsynchronized single-cell transcriptomic changes that relate to cell fate decisions using separated ARKO and control samples (Qiu et al., 2017) (Fig. 6A,B). Three well-defined cell branches were identified in control epithelial cells using single-cell trajectory analyses, one mainly comprised of basal cells and two of luminal cells (Fig. 6A, top panel). The basal cell branch was set as the point of origin in the trajectory plot, based on their prostatic stem/progenitor cell properties (Toivanen et al., 2016). Using expression plots, we further assessed the cellular properties of these cell branches. Strong expression of Krt5 was observed in the basal branch tip, whereas high expression of Krt19, Ar and its downstream targets, Pbsn and Nkx3.1 (Nkx3-1), was most localized within both luminal branch tips (Fig. 6C-E). Interestingly, Ki67-positive cells were focused in the central region between basal and luminal branch tips as well as the L1 branch, suggesting asymmetrical proliferation events leading to both basal and luminal differentiation along with propagation of mature luminal epithelial cells (Toivanen et al., 2016; Shafer et al., 2017). The expression of Etv4, a regulator of the Shh axis in branching morphogenesis (Lu et al., 2009; Herriges et al., 2015; Zhang et al., 2016), appeared in the transitional region between basal and luminal branches. Of note, only slight Shh expression was detected in basal and other cell branches in R26mTmG/+:Gli1CreER/+control samples. In contrast, the trajectory plot of epithelial cells from ARKO samples showed two basal branches and one luminal branch, varying significantly from the trajectory observed in controls (Fig. 6B). The smaller B2 branch falls between the main basal and luminal branch tips, suggesting the presence of an intermediate cell population, resulting from hindered prostatic epithelial differentiation from basal to luminal cells. Of note, the B2 branch mainly comprised cells from the BE2 cluster (Fig. 6B, top panel), and expressed squamous cell differentiation marker genes, including Krt6a, Sfn and Aqp3 (Fig. S5A,B) (Bollag et al., 2020; Gindele et al., 2020), confirming its aberrant differentiation fate. Intriguingly, elevated Shh expression was primarily observed in basal cell branches of ARKO samples (Fig. 6E). A lack of Ki67 expression was also observed in the above cell branches. Moreover, the expression of Etv4 and Hoxb13, two master regulators for prostatic branching morphogenesis and epithelial differentiation, were significantly reduced (Fig. 6D,E) (Economides and Capecchi, 2003). Consistently, the LE3 cluster, a well-differentiated and mature luminal cluster, only appeared in the control samples with high expression of Cd52, which encodes a protein with anti-adhesion properties to increase motility (Yamaguchi et al., 2008), typically detected on the cell surface of sperm (Fig. 6E). The data from single-cell trajectory analyses further demonstrate aberrant differentiation of prostatic epithelia in ARKO samples and provide novel insight into the mechanisms underlying such changes.
Stromal AR is indispensable in prostate epithelial morphogenesis and growth
Using in vivo tissue recombination assays, we directly assessed the role of stromal AR in Gli1-lineage cells acting as the cellular niche to support pubertal prostate epithelial cell growth. Both AR-deficient R26mTmG/+:Ar L/Y:Gli1CreER/+ mutant and R26mTmG/+:Gli1CreER/+ control mice were administered TM at P14 and sacrificed at P21, at which time Gli1-CreER-activated mGFP-positive cells were isolated, sorted and combined with sorted prostatic epithelial cells from 6-week-old wild-type mice (Lawson et al., 2007; Karthaus et al., 2014). The mixed cell suspensions were implanted under the kidney capsule of SCID mice and analyzed 8 weeks later (Fig. 7A). Gross examination showed that grafts combining control mGFP-positive cells with wild-type epithelial cells were transparent and heavier in weight than those with ARKO mGFP-positive cells and wild-type epithelial cells (Fig. 7B). Fluorescent microscopic analyses showed the presence of mGFP-positive cells surrounding epithelial glandular cells in both samples with control or ARKO stromal cells, confirming the presence of Gli1-CreER-activated mGFP-positive stroma (Fig. 7D1,F1). Histological analyses showed a glandular epithelium resembling secretory prostatic glands in normal control grafts (Fig. 7C1-C3). In contrast, grafts composed of ARKO mGFP-positive cells and wild-type prostatic epithelial cells revealed only a few undeveloped glands, lacking normal prostatic secretions (Fig. 7E1-E3). In the samples with control stromal cells, Krt5-positive basal cells appeared under luminal cells in normal prostatic glands (Fig. 7D6). Positive staining for Ki67 was also detected within luminal cells along with positive staining for Hoxb13 (Fig. 7D7,D8), but no clear staining for Gli1 was observed (Fig. 7D9). In contrast, glandular cells in the grafts of ARKO stromal cells and wild-type epithelial cells showed positive staining for AR (Fig. 7F2) but negative staining for prostatic differentiation markers, probasin and NKX3.1 (Fig. 7F2-F4). These glandular cells also showed positive staining for both Krt5 and Krt8, indicating their intermediate and immature cell characteristics (Fig. 7F5,F6). Unlike in control samples, there was almost no staining for Ki67 and Hoxb13 (Fig. 7F7,F8), suggesting their under-differentiated and less proliferative cellular properties. In support of our previous observations, noticeable staining for Gli1 was observed in surrounding stromal cells of these grafts (Fig. 7F9). The above data provide an additional line of scientific evidence reaffirming the decisive role of stromal AR in Gli1-lineage cells and provide novel insight into how they mediate prostatic epithelial differentiation and development.
Prostatic induction, morphogenesis and maturation are fully reliant on androgen-mediated paracrine interactions between prostatic mesenchyme and epithelium (Prins and Putz, 2008). Early tissue recombination studies demonstrated that mesenchymal AR but not epithelial AR is essential for prostatic embryonic induction and development (Cunha and Lung, 1978). These findings provided the first evidence that a cell niche in the UGM plays an essential role in supporting stem/progenitor-mediated prostate differentiation and expansion. In vivo, expression of AR in mesenchymal Gli1-expressing cells has been identified recently to be essential for early prostatic development, implicating their role in supporting prostate development and morphogenesis (Le et al., 2020). In the postnatal prostate, androgen signaling still remains essential for controlling prostate growth, homeostasis and regeneration through mesenchymal-epithelial interactions (Brennen and Isaacs, 2018). Mouse prostate morphogenesis primarily occurs starting at P15, and its pubertal growth and maturation continues during puberty (P25-P40) when circulating androgens levels rise (Sugimura et al., 1986; Hayashi et al., 1991; Staack et al., 2003). Currently, the roles and regulatory mechanisms for stromal androgen signaling in facilitating prostatic stem/progenitors during pubertal prostatic epithelial morphogenesis and growth are largely unknown. Specifically, the interaction between Shh and androgen signaling pathways in regulating pubertal prostate development is also unclear.
In this study, we showed that selective deletion of AR expression in prepubescent stromal Gli1-expressing cells significantly impairs pubertal prostate morphogenesis and growth in R26mTmG/+:Ar L/Y:Gli1CreER/+ mice. Examination of prostate tissues isolated from the above ARKO mice at P35 and P56 identified aberrant prostatic morphogenesis and robust epithelial growth retardation, specifically diminishing androgen-induced prostatic glandular expansion. Using scRNA-seq approaches, we identified a significant decrease in proliferative stromal and epithelial cell populations in ARKO samples and provided evidence for the growth retardation defects at single-cell resolution. In addition, a well-differentiated and mature luminal cell cluster was almost completely missing in ARKO samples. The data generated from the above scRNA-seq analyses corroborate prostate growth defects identified grossly and histologically and provide new, high-resolution assessments for cellular and molecular changes induced by AR deletion in Gli1-lineage cells.
To determine the decisive role of stromal AR in controlling the cellular niche necessary for pubertal prostate differentiation and growth, we pursued a series of experiments to assess the crucial role of stromal AR in Gli1-lineage cells. Analyses of the transcriptomic profiles of Gli1-lineage cells showed upregulation of Shh signaling and other developmental pathways such as Wnt, Notch and Tgf-β, appeared in ARKO mGFP-positive stromal cells compared with their control counterparts. In addition, an increase in Gli1, Ptc1 and Smo expression was identified in ARKO mGFP-positive cells in qRT-PCR assays. IHC analyses further revealed an inverse expression pattern between AR and Gli1 in prostate stromal cells. Aberrant transcriptomic changes were also identified in ARKO SMCs and fibroblasts at single-cell resolution through scRNA-seq analyses. Elevated hedgehog, as well as hypoxia, and cell cycle arrest-related signaling pathways were identified in both ARKO fibroblasts and SMCs, implying that loss of stromal AR results in elevated Shh signaling and places prostate cells under metabolic stress. These data implicate a new regulatory loop between AR and Shh signaling within prostatic stromal cells. However, using both ChIP-qPCR and ChIP-seq approaches, we were unable to identify the direct role of AR in controlling Gli1 transcription, suggesting that other indirect regulatory mechanisms may be involved. In support of this hypothesis, we observed elevated Shh expression in prostatic basal cells which, via canonical Shh signaling, would explain elevated Gli1 expression (Lamm et al., 2002). This would imply that basal Shh expression rises in response to an unidentified mesenchymal-to-epithelial signal downstream of stromal AR-deletion. Of note, a similar inverse expression between AR and Gli1 has also been reported in prostate cancer cells (Chen et al., 2009). Therefore, future investigation of this regulation should provide new mechanistic insight into the interaction between androgen and Shh signaling in both prostate development and tumorigenesis.
Specific deletion of stromal AR in Gli1-lineage cells resulting in robust prostatic epithelial growth defects in R26mTmG/+:ArL/Y:Gli1CreER/+ indicate a paracrine mechanism underlying AR-initiated reciprocal interactions between prostatic epithelial and mesenchymal cells during prostate pubertal development and growth. Because prostatic basal epithelial cells lay between stromal and luminal epithelial cells in prostatic glands, and possess prostatic progenitor/stem cell properties (Ousset et al., 2012), we specifically examined the single-cell transcriptomic changes between ARKO and control basal epithelial cells to determine aberrant signaling pathways altered by stromal AR deletion. Interestingly, the upregulation of Shh signaling was among those dysregulated signaling pathways, corresponding to the changes observed in ARKO stromal cells. Elevated expression of Shh mRNA transcripts appeared in ARKO basal cells at single-cell resolution. IHC analyses further demonstrated increased expression of Shh in prostatic basal cells adjacent to ARKO stromal cells that possess increased Gli1 expression. It has been shown that a variety of andromedins are produced in mesenchymal AR-expressing cells and affect prostate epithelial cells through a concentration gradient dependent manner (Yan et al., 1992; Brennen and Isaacs, 2018). Thus, identifying an increase in SHH expression in prostatic basal cells directly adjacent to AR-deleted stromal cells provides scientific evidence for stromal AR-initiated stromal-epithelial paracrine regulation during puberty. In this study, using ex vivo culture approaches, we further identified the inhibitory effect of Shh in prostate epithelial growth using prostate tissues isolated from ARKO and control mice (Fig. S6), which is consistent with earlier reports (Freestone et al., 2003; Wang et al., 2003; Berman et al., 2004). These data, along with the data from scRNA-seq analyses showing that upregulated hypoxia related signaling existed in ARKO stromal cells, implicate that loss of stromal AR in Gli1-lineage cells triggers Shh activation, resulting in metabolic stress and growth retardation in adjacent epithelial cells.
In this study, we also assessed the role of stromal AR in Gli1-lineage cells in governing epithelial cell fates using single-cell trajectory analyses. Aberrant differentiation paths were observed in both basal and luminal epithelial cells in ARKO samples with apparent defects in basal-to-luminal differentiation. Data from our in vivo tissue recombination assays functionally demonstrate the crucial role of stromal AR action in Gli1-lineage cells, acting as a cellular niche, stimulating prostatic epithelial growth and differentiation. Through these experiments, it appears that prostatic Gli1-lineage cells possess a diverse population of both fibroblasts and SMCs. Results from scRNA-seq analyses also showed very similar molecular changes in both ARKO fibroblasts and SMCs. Although these Gli1-lineage fibroblasts and SMCs appear to be a heterogeneous population, results from this current study along with previous reports (Lai et al., 2012; Yu et al., 2012) suggest that prostatic Gli1-lineage cells contain a population (or populations) that are functionally distinct from the general prostate stromal cell populations previously targeted in prostate development, growth and morphogenesis. For example, selective deletion of AR either in prostatic fibroblasts, alone or in combination with SMCs, showed only slight pubertal prostate growth defects in previous mouse models (Lai et al., 2012; Yu et al., 2012). Nevertheless, further characterization of cellular properties of Gli1-lineage cells will be extremely important to better understand the crucial role of stromal cell niches in prostate development and tumorigenesis.
MATERIALS AND METHODS
All experimental procedures and care of animals in this study were carried out according to the Institutional Animal Care and Use Committee (IACUC) at Beckman Research Institute at City of Hope and approved by the IACUC. Euthanasia was performed by CO2 inhalation followed by cervical dislocation.
Mouse genotypes and experiments
Gli1CreER and ROSAmTmG mice were obtained from Jackson Laboratories (stocks 007913 and 7676). ArLox/Y mice were obtained from Dr Guido Verhoeven (De Gendt et al., 2004). The wild-type and mutant littermates were generated and used in this study. To elicit genetic recombination, mice were intraperitoneally injected with 125 µg/g body weight of TM (Sigma-Aldrich) suspended in corn oil (Sigma-Aldrich) at P14 (Lee et al., 2015; He et al., 2018).
Histology and immunostaining
Prostate tissues with four different lobes were fixed in 10% neutral-buffered formalin (American MasterTech Scientific) and processed into paraffin. We cut 5 μm serial sections and processed them from Clearify (American MasterTech Scientific) to PBS through a decreasing ethanol gradient. For histological assessment, Hematoxylin and Eosin (H&E) staining was performed as previously described (Lee et al., 2015; He et al., 2018).
For IHC, slides were treated by boiling in 0.01 M citrate buffer (pH 6.0) for antigen retrieval, incubated in 0.3% H2O2 for 15 min, blocked in 5% normal goat serum (Gibco) for 1 h and incubated with appropriate antibodies (see Table S12) diluted in 1% normal goat serum at 4°C overnight. Slides were then incubated with biotinylated secondary antibodies for 1 h followed by horseradish peroxidase streptavidin (Vector Laboratories) for 30 min and visualized using a DAB kit (Vector Laboratories). Slides were counterstained with 5% (w/v) Harris Hematoxylin (Thermo Fisher Scientific) and coverslips were mounted.
For detecting mT and mGFP signals, whole prostate tissues containing four different lobes were fixed in 10% neutral-buffered formalin at 4°C overnight, cryoprotected in 30% sucrose at 4°C overnight and embedded in OCT (Tissue-Tek). We washed 5 μm sections on slides three times with PBS. To detect mT and mGFP signals, slides were directly mounted using Vectashield Mounting Medium with DAPI (Vector Laboratories). For Co-IF staining, slides were treated for antigen retrieval as described above, blocked in 5% normal goat serum for 1 h and incubated with primary antibodies diluted in 1% normal goat serum at 4°C overnight. Slides were washed in PBS then incubated with fluorescent-conjugated secondary antibodies for 1 h, and then mounted as described above. Detailed information regarding antibodies that were used in this study is provided in Table S12.
Microscope image acquisition
H&E and IHC slides were imaged using an Axio Lab A1 microscope with 10× and 40× Zeiss A-Plan objectives. Images were taken using a Canon EOS 1000D camera and analyzed using Axiovision software (Carl Zeiss). Images of immunofluorescent staining and mTmG signals were acquired on an Olympus Motorized Inverted Research Microscope Model IX81 using 20× and 40× Olympus Plan Fluor objectives, a QImaging RETIGA 2000R camera and Image-Pro 6.3 software (Media Cybernetics).
Serum testosterone measurement
Mouse serum testosterone levels were measured using a Mouse/Rat Testosterone ELISA kit (Alpco Diagnostic). Mouse blood samples were collected and tested following the manufacturer's protocol. The concentration of each sample, corresponding to mean of absorbance value, was then calculated from the calibration curve.
Preparation of single prostate cells
The wild-type or mutant mice (n=2) were injected with TM at P14 and prostate tissues were collected at P35 to prepare single cells for scRNA-seq or bulk RNA-seq. Two individual pairs of wild-type and mutant littermates were used for bulk RNA-seq and scRNA-seq. Single-cell suspensions from dissected whole prostate tissues containing four different lobes were prepared as previously described (Le et al., 2020). Briefly, mouse prostate tissues were collected, minced into small pieces and digested with 1 ml of collagenase (10 mg/ml, StemCell Technologies) in DMEM/F12 (Corning) with 10% fetal bovine serum (FBS), DHT (10 nM) and Y-27632 (10 µM) (StemCell Technologies) at 37°C for 90 min. Cells were then digested with TrypLE (1 ml, Gibco) at 37°C for 15 min and were centrifuged at 300 g for 5 min. Cells were passed through 40 µm nylon mesh (Fisherbrand) to get single cells, washed once with DMEM/F12 10% FBS and dissolved in 35 µl of PBS-0.05% bovine serum albumin for single-cell sequencing.
Single-cell suspensions from the above wild-type and mutant mice were used for library preparation with 10x Genomics Chromium Single Cell v3 Chemistry Solution following the manufacturer's protocols. The library purity and size were validated by capillary electrophoresis using the High Sensitivity DNA Kit (Agilent Technologies). The library quantity was measured fluorometrically using Qubit dsDNA HS Assay Kit from Invitrogen. Approximately 8000 and 9500 cells from the above ARKO and control mice, respectively, were captured on a 10x Chromium device using 10x Genomics Single Cell 3′ Solution V3 Kit, and the libraries were sequenced on Illumina Novaseq 6000 S4 flow cell (Illumina) to a depth of 60,000-70,000 reads per cell. Raw sequencing data were processed using the 10x Genomics Cell Ranger pipeline (version 3.1.0) to generate FASTQ files, and aligned to the mm10 reference genome with an added mGFP sequence (Wang et al., 2009) to generate gene expression count. Following alignment and initial quality control filtered feature, a total of 7480 and 8230 cells from ARKO and control samples, respectively, were uploaded as filtered feature matrices to R (3.6.1) using the Seurat package (18.104.22.16802) (Stuart et al., 2019). Post-filtering, 6871 and 7148 cells from ARKO and control samples, respectively, with an average of 2454 genes and 11,664 UMI counts per cell were used in the analyses (Fig. S2). Normalized and scaled data were clustered using the top principal components generated from 5000 highly variable genes and a resolution of 0.5 using Seurat. Cells then underwent further filtering to remove potential empty droplets and doublets (750<nFeature_RNA<7500) as well as low-quality cells with high percentages of mitochondrial RNA (percent.mt<15) (Fig. S2A). The Seurat package was used for tSNE visualization and gene expression analysis. Samples were then merged using Seurat's integrated merging method using 20 dimensions and a resolution of 0.5 (Stuart et al., 2019). Trajectory analysis and generation of pseudotime was performed by converting Seurat objects into CellDataSet format. Following conversion of the data to a CellDataSet, trajectory analysis was performed using the Monocle2 package (2.12.0) in R (Qiu et al., 2017). Pathway analysis was performed using GSEA 4.0.3.
RNA extraction and RNA-seq library preparation and sequencing
RNA samples were extracted with Trizol (500 µl, Zymo Research). RNA-seq libraries were prepared using SMARTer® Ultra™ Low Input RNA Kit for Sequencing v4 (TaKaRa Clontech) and KAPA Hyper Prep Kit (KAPA Biosystems) according to the manufacturers’ protocol. The resulting double-stranded cDNA was sheared using a Covaris LE220 Plus (Covaris) with a 200 bp peak in the 50 µl volume setting. The fragmented cDNA underwent end repair, 3′ end adenylation and ligation with barcoded adapters. The libraries were validated using the Agilent Bioanalyzer DNA High Sensitivity kit (Agilent) and quantified using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). The sequencing library templates were prepared for sequencing using the Illumina HiSeq SR Cluster V4 Kit. Sequencing runs were performed on an Illumina HiSeq 2500 using the single read mode of 51 cycles of read 1 and 7 cycles of index read with the SBS V4 Kit. Real-time analysis (RTA) 2.2.38 software was used to process the image analysis and base calling.
RNA-seq data analysis
For RNA-seq data analysis, the quality of sequence reads from the RNA-seq data were assessed and low-quality reads were filtered using the FastQC tool (Babraham Bioinformatics). Sequence alignment and quantification were performed using the STAR-RSEM pipeline (Li and Dewey, 2011; Dobin et al., 2013). To reduce systemic bias between samples, the Trimmed Mean Method (TMM) was applied to gene-level expression counts (Robinson et al., 2010). Genes were filtered out and excluded from downstream analysis if they failed to achieve raw read counts of at least two across all the libraries. Based on the normalized values, we identify DEGs between AR knock-out and wild-type samples using the previously reported statistical method (Shin et al., 2017). Briefly, T-value, rank-sum difference and log2-median-ratio were computed for each gene. Empirical distributions of the null hypothesis were estimated by calculating T-value, rank-sum difference and log2-median-ratio for the genes after randomly permuting the samples. For each gene, the adjusted P-values of the observed T-value, rank-sum difference and log2-median-ratio were computed using their corresponding empirical distribution by two-tailed test. Then, we combined these three P-values into an overall P-value using Stouffer's method (Hwang et al., 2005) and computed the false discovery rate (FDR) using Storey's method (Storey, 2002). DEGs were defined as the genes with an FDR less than 0.05 and absolute log2-median-ratio greater than 0.58 (1.5-fold). For the functional enrichment analysis, we used GSEA software (Subramanian et al., 2005). Briefly: the gene sets were obtained from MSigDB (Liberzon et al., 2011), genes in the RNA-seq dataset were sorted in descending order using the log2-median-ratio, we computed enrichment score (ES) using a Kolmogorov–Smirnov running sum statistic for the gene set, and significance of each gene set was computed using a distribution of null hypothesis which was generated from a random gene set by 1000 permutations.
Quantitative real-time PCR
qPCR reactions were performed in triplicate using an Applied Biosystems 7900 Fast sequence detector with SYBR Green PCR master mix (Thermo Fisher Scientific). Primers were designed using PrimerQuest (IDT) (Table S13). All reactions were normalized to expression of the housekeeping gene PP1A (Ppp1ca).
Tissue recombination assay
For tissue recombination assays, wild-type and mutant mice were injected with TM at P14 and prostate tissues were collected at P21. Dissociated prostate cells were prepared as described in the previous section. Cells were sorted for mGFP-positive and tdTomato-negative, or CD24- and CD49f-positive for prostatic epithelial cells (Karthaus et al., 2014). After sorting, cells were dissolved in DMEM/F12 (100 µl) with 10% FBS and counted using Trypan Blue (Gibco). The purity of mGFP-positive cells was confirmed by counting the number of mGFP-positive cells compared with the total number of cells that stained negative for Trypan Blue. All of the samples used in the study possessed >99% purity. Sorted GFP-positive and adult prostate epithelial cells were combined in collagen and grown in DMEM/F12 with 10% FBS supplementation overnight in 37°C incubator. Various combinations of GFP-positive cells and adult prostate epithelial cells were made as described in Fig. 7, including control samples with prostatic epithelial and stromal cells only. Approximately 40,000-80,000 epithelial cells and 160,000-200,000 stromal cells were used for each recombinant. The combined tissues were then implanted under the kidney capsule of 8-week-old male SCID mice the following day, and the grafts were analyzed 8 weeks later.
The wild-type and mutant mice were injected with TM at P14 and individual lobes of prostate whole mounts were dissected and plated on cell culture inserts (20 µm pore size; Millipore) individually in serum-free medium. Serum-free medium was made of DMEM/F12 plus ITS (Gibco), and 100 units/ml penicillin and 100 mg/ml streptomycin, modified as previously described (Lipschutz et al., 1997). Shh (2 µg/ml; R&D), cyclopamine (5 µM; Toronto Research Chemicals) and testosterone (10–9 M; Sigma-Aldrich) was added to the medium at the beginning of the cultures. The medium was changed every other day. The cultures were collected at day 10 and fixed in 10% formalin in phosphate buffer (0.1 M, pH 7.4) overnight at room temperature before they were processed for immunocytochemistry.
Quantification and statistical analysis
Statistical analyses were performed using GraphPad Prism 6. All data are presented as mean±s.d. Two group comparisons were analyzed with a two-tailed Student's t-test and a value of P<0.05 was taken as statistically significant.
Conceptualization: A.W.O., V.L., Z.S.; Methodology: A.W.O., V.L., J.W., A.H., W.K.K., D.-H.L., J.A., X.W., S.Y., Z.S.; Software: X.W., M.K., S.Y.; Validation: V.L., A.H., W.K.K., D.-H.L., G.R.C., S.Y., Z.S.; Formal analysis: A.W.O., V.L., J.W., W.K.K., X.W., M.K., S.Y., Z.S.; Investigation: A.W.O., V.L., J.W., W.K.K., G.R.C., Z.S.; Resources: Z.S.; Data curation: Z.S.; Writing - original draft: A.W.O., Z.S.; Writing - review & editing: A.W.O., J.A., Z.S.; Visualization: D.-H.L., G.R.C., Z.S.; Supervision: Z.S.; Project administration: Z.S.; Funding acquisition: Z.S.
This work was supported by National Institutes of Health grants R01CA070297, R01CA166894, R01DK104941 and P30CA033572. Deposited in PMC for release after 12 months.
The single cell mRNA-sequencing and bulk RNA-sequencing data have been submitted to GEO under accession number GSE156168.
The authors declare no competing or financial interests.