Neuroendocrine prostate cancer has distinctive, non-prostatic HOX code that is represented by the loss of HOXB13 expression

HOX gene-encoded homeobox proteins control body patterning during embryonic development; the specific expression pattern of HOX genes may correspond to tissue identity. In this study, using RNAseq data of 1019 human cancer cell lines that originated from 24 different anatomic sites, we established HOX codes for various types of tissues. We applied these HOX codes to the transcriptomic profiles of prostate cancer (PCa) samples and found that the majority of prostate adenocarcinoma (AdPCa) samples sustained a prostate-specific HOX code whereas the majority of neuroendocrine prostate cancer (NEPCa) samples did not, which reflects the anaplastic nature of NEPCa. Also, our analysis showed that the NEPCa samples did not correlate well with the HOX codes of any other tissue types, indicating that NEPCa tumors lose their prostate identities but do not gain new tissue identities. Additionally, using immunohistochemical staining, we evaluated the prostatic expression of HOXB13, the most prominently changed HOX gene in NEPCa. We found that HOXB13 was expressed in both benign prostatic tissues and AdPCa but its expression was reduced or lost in NEPCa. Furthermore, we treated PCa cells with all trans retinoic acid (ATRA) and found that the reduced HOXB13 expression can be reverted. This suggests that ATRA is a potential therapeutic agent for the treatment of NEPCa tumors by reversing them to a more treatable AdPCa.

www.nature.com/scientificreports/ cells, and HOXB13 is considered a marker for prostatic cancer origins, to distinguish metastatic PCa from other cancers 14,15 . However, it remains unclear whether HOXB13 plays a tumor suppressive or pro-oncogenic role in PCa. On the one hand, HOXB13 expression is elevated in castrate-resistant PCa and the induced expression of HOXB13 promotes androgen-independent growth of LNCaP cells 16 . On the other hand, constitutive expression of HOXB13 suppresses colony formation in these cells 17 .
In this study, we established a list of HOX codes for various tissue sites and applied them to PCa samples. We found that NEPCa samples had different HOX codes from AdPCa samples. Also, using immunohistochemical staining, we evaluated the protein expression of HOXB13, the most prominently changed HOX gene in NEPCa. The results revealed that HOXB13 expression was maintained in AdPCa but was reduced or lost in NEPCa. Finally, we report that the decreased HOXB13 expression in PCa cells can be reverted.

Results
Establishment of HOX codes of various tissue origins. It has been suggested that different tissues have unique HOX codes. Using the RNA-seq data collected by the Cancer Cell Line Encyclopedia (CCLE) 18 , we analyzed the expression of HOX genes in 1019 cell lines that were derived from various tissue origins and calculated the HOX code for each tissue type. The workflow is presented in Fig. 1A and the median correlation scores between cell lines of the same tissue origin and the HOX codes of various tissues are presented as a heatmap (Fig. 1B). This heatmap reflects how closely each tissue type is correlated with its HOX code. As shown in Fig. 1B, most tissue types had a rather strong correlation with their calculated HOX codes. Among these tissues, the most notable was the prostate group (consisting of 8 cell lines), followed by autonomic ganglia and skin, all displaying a strong correlation with their HOX code. To further validate the HOX codes, we applied the CCLEderived HOX codes to TCGA datasets of various cancer types and found that most HOX codes displayed high correlation with the tissues they represented (sFigure 1). Among these, prostatic HOX code is a prominent and exclusive one.
NEPCa tumors display a change in HOX code. The strong correlation between the tissue types and their HOX codes suggests that the HOX codes could be used to identify the tissue origin of a given sample. Additionally, HOX codes could be used to study how closely a sample is related to its tissue origin. A "lineagekeep" sample would display a high correlation with the HOX code of its tissue origin, whereas a "lineage-switch" sample would show a low correlation. Therefore, we applied the HOX codes to a RNA-seq dataset of human PCa including 34 CRPCa and 15 NEPCa samples 5 . The correlation scores between the expression of HOX genes in each sample (row) and the HOX codes for various tissues (column) were visualized in a heatmap (Fig. 1C). The vast majority (31/34) of CRPCa samples displayed a high correlation (Pearson correlation score > 0.3) with the HOX code of prostate but not of other tissues. More importantly, the majority of (13/15) NEPCa samples showed a low or negative correlation with the prostate HOX code. Similar results were also observed in other NEPCa datasets including SU2C dataset 19 (sFigure 2A), GSE32967 20 (sFigure 2B), GSE41192 21 (sFigure 2C), GSE66187 22 (sFigure 2D). While NEPCa tumors displayed a loss of the prostate specific HOX code, they did not show a consistent correlation with the HOX codes of any other tissues ( Fig. 1 and sFigure 2).
Taken together, these results indicate that HOX code, to a certain extent, could reflect the identity of tissues and that the loss of the prostate specific HOX code in NEPCa reflects a loss of prostatic identity in these tumors.
The expression of HOXB13 is decreased in NEPCa. Among all the HOX genes, HOXB13 displayed the highest mRNA level in prostate specimens (sFigure 3A). Compared with other posterior HOX genes (HOX13s), HOXB13 also displayed a more consistent expression pattern in TCGA prostatic samples (sFigure 3B). Therefore, we specifically examined HOXB13 expression in prostatic tissues.
Analysis of PCa TCGA dataset indicates that consistent with a previous report 23 , the mRNA level of HOXB13 is higher in PCa than in normal prostatic tissues (sFigure 3C). Additionally, analysis of NEPCa RNA-seq datasets 5,[19][20][21][22]24 indicates that HOXB13 is the most changed HOX gene in NEPCa when compared with AdPCa. As shown in Fig. 2, the mRNA level of HOXB13 was lower in NEPCa samples than AdPCa in all the NEPCa datasets examined. In contrast, the mRNA levels of several other members of the HOXB family including HOXB3, HOXB5, HOXB6 and HOXB7 were higher in NEPCa than AdPCa.
Further, immunohistochemical (IHC) staining was conducted to assess the protein expression of Hoxb13. As shown in Fig. 3A-C, in murine PCa derived from TRAMP mice (n = 10), Hoxb13 expression was detected in prostatic intraepithelial neoplasia (PIN) lesions (Fig. 3A). However, the expression of Hoxb13 was decreased in the NEPCa area (Fig. 3A), which was highlighted by the positive stain of NEPCa markers Foxa2 (Fig. 3B) and chromogranin A (Fig. 3C).
The expression of HOXB13 was also evaluated in human prostatic tissues including benign prostate, AdPCa and NEPCa. As shown in Fig. 3D-I, HOXB13 expression was detected in both benign prostatic tissues (Fig. 3D) and prostate adenocarcinomas (AdPCa, Fig. 3E). In benign prostate samples, HOXB13 was stained positive in 26 of 28 BPH. Compared to that of BPH, the intensity of HOXB13 staining was greater in adenocarcinomas with Gleason scores of 7, 8, 9, and 10 albeit not statistically significant, possibly due to the wide variations in H-score observed in both BPH and the adenocarcinomas (Fig. 3J). The largest variation in H-score, ranged from 15 to 280, was noted in adenocarcinoma with Gleason score of 9. The highest H-score of 280 was found in adenocarcinomas with Gleason scores of 9 and 10. There was a positive correlation between the expression of HOXB13 (H-score) and Gleason scores of AdPCa (Spearman correlation test, Rho = 0.348, p < 0.05).
In contrast to the positive stain in AdPCa, no HOXB13 immunohistochemical reactivity was observed in 8 of 9 NEPCa including 8 small cell carcinoma and 1 neuroendocrine carcinoma, not otherwise specified. Of Majority tissues correlated well with their HOX codes. It is noteworthy that prostate tissue displayed a high correlation with its HOX code (C) AdPCa and NEPCa tumors displayed high and low correlation with prostate HOX code, respectively.   www.nature.com/scientificreports/ the 9 NEPCa, there was only 1 small cell carcinoma stained weakly and focally (H-score of 2) with HOXB13 immunomarker. In summary, HOXB13 expression was detected in both benign prostatic tissues and prostate adenocarcinomas but its expression was lost or reduced in NEPCa tumors.

HOXB13 is expressed in prostatic luminal epithelial cells. Dual immunofluorescence staining was
conducted to examine the expression pattern of HOXB13 in prostatic glands. As shown in Fig. 4, the expressions of HOXB13 and basal epithelial cell marker (Cytokeratin 5 or Cytokeratin 14) were mutually exclusive. However, HOXB13 was co-expressed with prostate luminal epithelial cell marker NKX3-1. The luminal expression of HOXB13 was further supported by an analysis of a single cell RNA-seq dataset of normal human prostate (GSE120716), where HOXB13 positive cells were enriched in luminal epithelial subpopulation (sFigures 4A,B) 25,26 . A similar result was seen in a single cell RNA-seq analysis of mouse prostate (sFigure 4C) 27 . DU145, a commonly used PCa cell line, is more similar to NEPCa cell line H660. HOX code analysis of the eight prostatic cell lines that are collected in the CCLE dataset 18 indicated that H660 (a NEPCa cell line), DU145 (a widely used AdPCa cell line), and PRECLH cell lines displayed a low association with the prostate HOX code whereas the other five PCa cell lines displayed a strong correlation with prostate HOX code (sFigure 5). This is consistent with the results of clustering analysis of HOX genes in these cell lines, which revealed that the cell lines that have lost prostate HOX code (DU145, H660, and PRECLH) were clustered together, distinct from the other five cell lines (Fig. 5A). Compared with that of the other five cell lines, the mRNA levels of HOXB13 were lower in the three cell lines that have lost prostate HOX code. The reduced HOXB13 expression in H660 and DU145 cells was confirmed by Western Blot analysis ( Fig. 5B and C). These results indicate www.nature.com/scientificreports/ that DU145, labeled as an AdPCa cell line, is more similar to a NECPa cell line than other AdPCa cell lines with regard to its HOX code.

The reduced expression of HOXB13 in PCa cells can be reverted. Accumulating evidence sup-
ports that NEPCa tumors emerge via the trans-differentiation of AdPCa and that transcriptome reprogramming plays a key role in this process. Our HOX code analysis suggests the trans-differentiation is accompanied with a change in HOX code. Here using HOXB13, the most prominently changed HOX gene, as an indicator of epigenetic reprogramming in NEPCa, we tested whether the reduced HOXB13 expression in NEPCa cells could be reverted. We treated DU145 and NCI-H660 cells with all-trans retinoic acid (ATRA), a commonly used agent to induce the differentiation of stem cells 28 . The established NEPCa NCI-H660 cell line was used in this experiment. Because the DU145 cell line was shown to display a similar HOX code to H660 and previous studies have shown that DU145 cells have some NEPCa features 29 , this cell line was also included in the ATRA experiments. The AdPCa LNCaP cells, which express high level of HOXB13, were used as a control cell line. As shown in Fig. 5D and E, the mRNA levels of HOXB13 were induced in both H660 and DU145 but not LNCaP cells that were treated with ATRA. However, ATRA did not induce the protein expression of HOXB13 in these cells (data not shown). To further evaluate whether ATRA has any functional roles in altering the proliferation of PCa cells, we conducted cell proliferation assays and found that ATRA at 10 μM did not significantly change the proliferation of DU145 cells (Fig. 5F) or LNCaP cells (data not shown). Taken together, these data indicate that ATRA can revert the mRNA expression of HOXB13 in NEPCa but not AdPCa cells. However, ATRA could not induce the protein levels of HOXB13 nor affect the proliferation of PCa cells.

Discussion
In this paper, we report that HOXB13 is expressed in benign prostatic tissues and prostate adenocarcinomas, but its expression is decreased or lost in both human and mouse NEPCa. These data indicate that the loss of HOXB13 expression is a common feature of NEPCa. During embryo development, a nested and partially overlapping HOX expression pattern is established. A change of HOX expression results in deficiencies in axial formation and the loss of tissue identity 1 . In this study, using transcriptomic data of 1019 cell lines, we established a method to calculate the HOX codes for 24 tissues and applied the HOX codes to PCa RNA-seq datasets. Using this method, we identified samples that uphold or lose prostate identity. We showed that prostate HOX code was sustained in the majority of AdPCa samples. This result is consistent with previous reports that HOXB13 is a PCa (adenocarcinoma) marker 14 . In contrast, the majority of NEPCa tumors and the NEPCa cell line NCI-H660 displayed low correlation with the HOX code of prostates, indicating a loss of prostate identity in NEPCa cells. DU145, an AdPCa cell line that has certain NEPCa features 29 , also displayed a loss of HOX code. Taken together, these data suggest that during NEPCa development, prostate identity is gradually lost in PCa cells.
Though NEPCa tumors seemed to lose prostate identity, they did not show consistent correlation with the HOX codes of any other tissues, suggesting that they have not acquired a clear-cut new identity. This is in line with the consensus that NEPCa cells lose prostate differentiation, but this notion is not in complete agreement with the proposed acquisition of neuronal features in NEPCa. Our observation of the lack of clear-cut identity of NEPCa adds to this line of knowledge. Although the acquired neuronal "features" in some NEPCa tumors resemble neuronal characteristics in certain aspects, NEPCa cells have not acquired the neuron "identity".
Additionally, our HOX code method provides a unique technique for studying lineage switching during cancer progression. By applying the HOX codes to the transcriptome data of various prostate cell lines, we identified three cell lines that lose the prostate-specific HOX code. Among these, H660 is a NEPCa cell line and DU145 is a commonly used androgen-independent AdPCa cell line. Analysis of the HOX expression profile indicates that DU145 is more similar to NEPCa (H660) than other AdPCa cell lines. Given the NEPCa features in DU145 cells, these cells may represent a unique state of PCa-they have started losing prostatic differentiation but have yet to gain a full NE phenotype. This cell line could be a good model system for inducing NE differentiation.
Currently, there are no effective treatment options against NEPCa. A possible therapeutic approach is to revert the lineage switch of NEPCa. By inducing prostatic differentiation in NEPCa cells, we might be able to convert them back to the AdPCa state. For this purpose, we treated HOXB13-low NEPCa (NCI-H660) and NEPCa-like (DU145) cells with retinoic acid, a commonly used agent to induce differentiation 30 . We found that retinoic acid induced the mRNA expression of HOXB13 in both H660 and DU145 cell lines, but not in AdPCa LNCaP cells. This indicates that the loss of HOXB13 expression during the development of NEPCa is reversible. This notion also supports the trans-differentiation theory of NEPCa. However, ATRA did not induce the protein expression of HOXB13 in DU145 cells, suggesting the involvement of post-transcriptional regulatory mechanisms in HOXB13 expression. This may also explain the discrepancy between the mRNA and protein levels of HOXB13 in PC3 cells.
Our observation of the effects of ATRA is consistent with previous reports that ATRA induces the expression of HOXB13 in DU145 cells and inhibits NEPCa in TRAMP mice, PCa mouse models 31,32 . We showed here that ATRA reverts the lost expression of HOXB13 in NEPCa cells, providing evidence supporting the possibility of inducing prostate differentiation in NEPCa cells. Further research is warranted to study the therapeutic effects of retinoic acids and other differentiation-inducing agents in the treatment of NEPCa. www.nature.com/scientificreports/ HOX code for each tissue origin. We then applied the HOX codes (Supplementary Table 1) to each of the 1019 cell lines by calculating the Pearson correlation scores between the expression of HOX genes in that cell line and the HOX codes of various tissue origins. A heatmap was constructed to visualize the median correlation score in each group of samples (X-axis) and their correlation with the HOX codes of various tissue origins (Y-axis). To validate the HOX codes, we applied the HOX codes to TCGA datasets of various types of cancer. We also applied the HOX codes to publicly available NEPCa datasets and calculated the Pearson correlation scores to reflect the correlation between the HOX expressions in each prostatic sample and the HOX codes of various tissue origins.

Materials and methods
Sample collection. De-identified human prostate tissue specimens were obtained from LSU Health-Shreveport Biorepository Core, Overton Brooks VA Medical Center, Ochsner Health System Biorepository and Tissue for Research as described previously 33 . The specimens used in this research include 28 benign prostate hyperplasia, 44 AdPCa, and 9 NEPCa tissues. All the tissues were used in accordance with LSU Health-Shreveport IRB protocols. Archived tissue sections of TRAMP tumors were used for this study.
Immunohistochemical and immunofluorescence staining. Immunostaining was performed using Vectastain elite ABC peroxidase kit (Vector Laboratories, Burlingame, CA) as described previously 33  Cell proliferation assay. IncuCyte S3 live cell imaging system (Essen BioScience, Ann Arbor, MI) was used to evaluate the effects of ATRA on the cell proliferation. PCa cells (1000 cells/well) were seeded in 96-well plate and cultured in media containing DMSO or 10 μM ATRA. Cell culture media were refreshed every two days. Images were taken every four hours and the cell proliferation was monitored based on the cell confluency.
RNA extraction and RT-qPCR. RNA was extracted by using RNeasy Mini Kit (Germantown, MD).
Reverse transcription and qPCR were conducted by using iScript Reverse Transcription Supermix and SYBR green PCR Supermix (BioRad, Hercules, CA). GAPDH was used to normalize the gene expression. Primer sequences were listed in sTable 2.
Western blot. Cells were collected in PBS and lysed in passive lysis buffer (Promega, Madison, WI). Equal amounts of protein were loaded for Western blot analyses. ProSignal Dura ECL Reagent (Genesee Scientific, San Diego, CA) and Chemidoc (Bio-Rad, Hercules, CA) were used to visualize the proteins. Beta-tubulin was used as loading control.
Statistical analyses. T-test was used to compare the gene expression in AdPCa and NEPCa samples.
Spearman Correlation Coefficient test was used to analyze the correlation between tumor grade and protein expression (H-score). p < 0.05 was considered statistically significant. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.