Transcriptional regulators and regulatory pathways involved in prostate gland adaptation to a hypoandrogen environment

Abstract Anti-androgen therapies, including orchiectomy, are effective at promoting prostate cancer remission, but are followed by progression to the more aggressive castration-resistant prostate cancer (CRPC). Castration promotes gland and tumor shrinkage. However, prostate adaptation to androgen deprivation involves striking parallel events, all requiring changes in gene expression. We hypothesized that transcription factors (TF) and other transcription-related genes are needed to orchestrate those changes. In this work, downstream analysis using bioinformatic tools and published microarray data allowed us to identify sixty transcriptional regulators (including 10 TF) and to integrate their function in physiologically relevant networks. Functional associations revealed a connection between Arnt, Bhlhe41 and Dbp circadian rhythm genes with the Ar circuitry and a small gene network centered in Pex14, which might indicate a previously unanticipated metabolic shift. We have also identified human homologs and mapped the corresponding genes to human chromosome regions commonly affected in prostate cancer, with particular attention to the PTEN/HHEX/MXI1 cluster at 10q23-25 (frequently deleted in PCa) and to MAPK1 at 22q11.21 (delete in intermediate risk but not in high risk PCa). Twenty genes were found mutated or with copy number alterations in at least five percent of three cancer cohorts and six of them (PHOX2A, NFYC, EST2, EIF2S1, SSRP1 and PARP1) associated with impacted patient survival. These changes are specific to the adaptation to the hypoandrogen environment and seem important for the progression to CRPC when mutated.


Introduction
Prostate diseases in general and prostate cancer (PC) in particular are major concerns of public health care. One eighth to one sixth of males will develop PC and experience the risk of prostate cancer progression if not properly diagnosed, monitored and treated. Molecular markers for diagnosis and disease progression risk assessment are extremely necessary. One of the palliative treatments for advanced prostate cancer is androgen blockade, achieved by chemical or surgical castration. Besides the psychological and physiological side effects, the risk of androgen blockade resides in the common progression to the highly malignant and life-threatening form of castration-resistant prostate cancer (CRPC). Several mutations and chromosomal rearrangements have been associated with PC and CRPC. Recently, a series of chromosomal translocations including frequent bridging and rearrangements was de-scribed and demonstrated to occur in a few steps and to affect a series of important tumor suppressors (Baca et al., 2013). Nonetheless, the linear progression among progressive stages has been questioned, as metastatic lesions are clonally derived, show signatures particular to each affected individual and might result from less advanced local primary lesions (Holcomb et al., 2009, Haffner et al., 2013. Androgens, acting through the androgen receptor (AR), are required for prostate development and normal function (Roy et al., 1998). Hence, surgical or chemical castration defines a hypoandrogenic state in which a series of events sum up to promote organ (and tumor) shrinkage; androgen deprivation/blocking is the first line therapy for advanced prostate cancer. Gain-of-function mutations enabling the AR to recover activity in the hypoandrogen environment have been associated with the progression to CRPC (Feldman and Feldman, 2001), and include mutations, deletions and inversions at the ligand binding site leading to ligand-independent activation (Nyquist et al., 2013), and gene amplifications. In spite of these AR-centered modifications, a series of genes has been implicated in prostate cancer progression, including overall changes in gene expression (Abate- Shen and Shen, 2000, Shen and Abate-Shen, 2010 and complex chromosomal rearrangements frequently involving PTEN, NKX3.1,SPOP,CHD1,TP53,MAP3K7,FOXP1 and the T2-ERG fusion (Baca et al., 2013).
Furthermore, we consider that other physiological aspects of the adaptation to the hypoandrogen environment might be corrupted in cancer cells, and cooperate in establishing the selective pressure that contributes to the clonality of cells harboring chromosomal changes in general and AR modifications in particular.
Although epithelial cell apoptosis is a major event in prostate regression occurring in response to castration, it is not the sole event. For instance, remarkable reprogramming of immune system cells (Desai et al., 2004) and smooth muscle cells (Antonioli et al., 2004(Antonioli et al., , 2007 as well as reorganization of the extracellular matrix (Vilamaior et al., 2000) have been described and associated with a redefined functional state and immune barrier system. Additionally, we have reported the occurrence of desquamation as an additional phenomenon contributing to epithelial cell deletion (Rosa-Ribeiro et al., 2014a), and a relevant role for two macrophage subpopulations in both (a) the induction of epithelial cell death (Barbosa et al., 2019) and (b) the clearance of cell corpses and maintenance of the noninflammatory status (Silva et al., 2018).
However, little has been studied beyond the induction of apoptosis in epithelial cells. Progress has been made in terms of showing that epithelial expression of the AR is not necessary for epithelial cell death (Kurita et al., 2001), while the remaining cells develop resistance to androgen deprivation and preserve a differentiation-immature signature (Rosa-Ribeiro et al., 2014b).
The rat ventral prostate has been a valuable and robust in vivo model system to explore androgen regulation of gene expression (Wang et al., 1997, Kwong et al., 1999, Desai et al., 2004. The ventral prostate responds to androgen withdrawal with increased epithelial cell apoptosis, whereas the dorsolateral lobes show negligible cell death (Kwong et al., 1999). The responsiveness of the ventral prostate to androgens is also characterized by the number of differentially expressed genes after castration as compared to the dorsolateral lobe (1496 vs. 256 genes, respectively) (Desai et al., 2004). Most of these changes take place in the hypoandrogen environment, and are triggered by disengaging AR signaling.
We hypothesized that other transcription factors (TF) and transcriptional regulators (TR) are co-opted to coordinate the sequential changes observed in the gland after castration. In this work, we explore this idea in an attempt to find new genes, unforeseen metabolic process and chromosomal hotspots that might reveal possible connections between gland physiology under androgen deprivation promoted by castration and progression to CRPC after hormone therapies.
We have used bioinformatics to (a) select classspecific genes from a published list of genes and ESTs differentially expressed in response to castration and androgen supplementation after DNA microarray analysis, (b) to identify the regulatory networks in which the selected genes are involved, (c) to map the homologs of rat genes to human chromosomes, and (d) to find mutations and/or copy number alterations and changes in patient survival.
Accordingly, this study has unveiled a list of TF and TR genes and a series of unexplored physiological pathways, such as circadian rhythms (genes Arnt/Bhlhe41/Dpb) and peroxisome biogenesis (Pex14), hitherto neglected pathways in prostate biology. We also correlated the selected genes with chromosomal regions commonly deleted in prostate cancer, such as 10q23, which contains the PTEN/HHEX/MXI1 gene cluster, and 22q11.21, harboring the MAPK1 gene, found 20 genes mutated in at least 5% of three patient cohorts and six genes affecting patient survival when mutated.

Material and Methods
The microarray data from Desai et al. (2004) reported 1496 genes/ESTs differentially expressed in response to castration and testosterone supplementation. The list of gene bank accession IDs for all genes and ESTs was loaded into DAVID v6.7. Gene IDs and biological annotations are highly redundant within the vast array of public databases. The DAVID knowledge base collects and integrates various gene identifiers as well as more than 40 well-known publicly annotation categories, which are then centralized by the internal DAVID identifier in a non-redundant manner. A significant portion of input gene IDs failed to be mapped and were then processed using the gene ID conversion tool. All the identified IDs/gene names were listed by the gene name batch viewer. We further processed the identified IDs for the identification of functional annotations centered on TFs and TRs, and the identified genes were further studied to find their functional annotation clustering and possible integration in known biological functions.
The TFs/TRs were also studied for possible functional associations using the Ingenuity Pathway Analysis (IPA) software, with filtering for information in the rat, and choosing only direct interactions.
The human homologs to the rat genes were searched manually using the NCBI database, and their chromosomal location was used to map them to the human ideogram.
Finally, we assessed the cBioPortal (cbioportal.org) and checked three cohorts of prostate adenocarcinomas for the existence of mutations and/or copy number alterations and possible effect on patient survival (Armenia et al., 2018, Liu et al., 2018, Abida et al., 2019.
A limit of 5% mutations and a log rank test p-value smaller than 0.05% were set for each analysis, respectively.

Data processing
Using the gene accession conversion tool of DAVID v6.7, the program managed to convert 468 IDs from the list of 1477 total unique user IDs. The number of genes identified was similar to that obtained in the original work (Desai et al., 2004). Out of the 468 IDs, DAVID identified 60 TFs/TRs. Table 1 lists the detailed annotation and func-tional enrichment information that was retrieved using the terms TFs/TRs. The chromosomal location for each gene was determined using the NCBI databank. Twenty-two genes were identified as transcription factors (bold in Table  1).

Roles of transcription factors and functional associations among the selected genes
The selected TFs and TRs function in 17 important cellular pathways identified by DAVID ( Table 2). Some of Gene expression after castration 3 Table 1 -Genes with accession numbers, names and functions and chromosomal location in the rat chromosomes as well as the chromosomal location of the human homologs. Known and putative transcription factors are bold-faced.

S.No Gene Accession Number
Gene  the ontogenies were very general (such as "Pathways in cancer" or "Prostate cancer" or "Type II diabetes mellitus"); all but 3 (14/17) included Mapk1, and about half (9/17) contained Mapk9, related to areas of strong research or hubs in central signaling pathways. The ontogenies also pointed to TGF-b, Toll-like receptors and T-cell receptor signaling pathways. Novel ontogenies implicated the genes Arnt and Bhlhe41 in "Circadian rhythms" and Ets-1 (E26 oncogene homolog 1) in "Dorso-ventral axis formation".
Gene expression after castration 7  Table S1. The genes shown in gray are those identified in this work.
Upon further inspection, using IPA to set the interactions among the 60 genes, we found eight networks corresponding to known pathways (Figure 1; Table S1). They vary in terms of the number of individual nodes, but reveal interesting aspects of the yet-to-be proven physiology of the prostate gland in the hypoandrogen environment. Perhaps not surprisingly, they are ascribed to gene expression regulation, cell death and survival, and also to nucleic acid and carbohydrate metabolism and cancer. They also implicate particular pathways such as estrogen receptor, retinoic acid receptor, thyroid hormone receptor, NFkB signaling, TGF-b and establishing connections with the newly identified genes. It is interesting to note that Arntl and Bhlhe41, both involved in circadian rhythms (Pathway 12, in Table  2), appeared together in network number 2 (Figure 1). Arntl connected to the AR via either p300/EP300 or CREBBP acetyl transferases, and directly do Dbp, another circadian rhythm gene. Additionally, IPA retrieved one particularly interesting pathway, peroxisomal biogenesis and function, referenced by network number 8, which is centered on the gene Pex14.

Mutation rates and effect on patient survival
PTEN is a tumor suppressor frequently associated with prostate cancer (Abate- Shen and Shen, 2000, Shen and Abate-Shen, 2010, Baca et al., 2013. We found PTEN mutated in 16%, 21% and 33% of the patients for the three cohorts studied. Beside PTEN, only IKBKB was found mutated in at least 5% of the patients in the three cohorts. CEBPD and DMD showed mutations in more than 10% of the patients in one cohort. BTG2, SERP1, CKNH2, RXRG, TXN1P, UBE2I, ZFP37, BARHL1, RALGAPA1, CDK7, PARP1, PTEFGRN, TXN1, TMPO, TNFRSF1A and AGRN were mutated in at least 5% of the patients in at least one cohort. However, more than 50% of the studied genes showed deletions in the three cohorts. Survival curves existed for two of the three cohorts studied. PHOX2a and NFYC were associated with significant impact on patient survival (P < 0.05) in the first cohort and EST2, EIF251, SSRP1 and PARP1, in the second cohort.

Discussion
Sixty differentially expressed TF and TR were retrieved from the microarray data published by Desai et al. (Desai et al., 2004). These genes were assorted into 17 pathways by the DAVID knowledge base and into eight functional networks by Ingenuity Pathway Analysis (IPA). Though most these pathways were too general, and led to central metabolic hubs, such as PTEN (which was further validated in the parent study) and nuclear receptors pathways, some revealed particularly interesting and unforeseen aspects of prostate biology, such as circadian rhythms and peroxisome biogenesis. The selected genes were also mapped to chromosome regions frequently affected in prostate cancers, as their identification might serve as risk factors or therapeutic targets relevant to progression to CRPC. Additionally, 20 genes were found mutated in at least 5% of prostate adenocarcinoma patients.
We used the DAVID module Gene ID Conversion Tool (Huang et al., 2007a(Huang et al., , 2007b to identify gene IDs from the initial gene list. The number of genes was roughly the same as those uncovered in the parent study, perhaps indicating a potential limit for retrieving information from the early commercially available microarray chips. In order to further advance our understanding of the physiology and endocrinology of the VP gland and to facilitate the biological interpretation of prostate biology in a broad range of biological processes, the 468 genes were further processed to track down their functional classification and resulted in 60 different TF or TR, which were then assigned to known regulatory pathways using DAVID and to functional interaction networks using IPA, resulting in the identification of 22 transcription factors. Prostate cancer will affect one eighth to one sixth of men worldwide. In spite of enormous progress in the understanding of many facets of the disease, new concepts are Gene expression after castration 9 still emerging. One of these is the clonal origin of metastases (Haffner et al., 2013), punctuated rather than gradual progression (Baca et al., 2013) and the non-linear relationship of on-site and distant metastases, meaning that metastases might be generated from less advanced local tumor foci (Haffner et al., 2013). Accordingly, it is not completely understood as to which foci in the advanced stage of the disease will progress to CRPC. We followed the idea that the physiological adaptation of the gland to the hypoandrogenic castration-induced environment involves regulatory pathways to maintain the gland in a regressed, low proliferative and less functional (meaning less differentiated) state, and that these pathways in cancer cells might be corrupted and thereby contribute to the progression to CRPC. In this scenario, molecules with enhanced expression after castration might be found to be "tumor suppressors", particularly if they function as hubs in regulatory networks that are defective in cancer cells. Smad4, Ikbkb, Rara, ets1, Bhlhe41, Id2, Tnfrsf1a, Mxi1 and Dbp are new candidates, together with the well-known tumor suppressor Pten. As a matter of fact, these genes (except ID2 and Mxi1) were found deleted in prostate adenocarcinomas and advanced metastatic prostate cancers. This also raises the question of whether Pten-defective tumors should be submitted to androgen deprivation or blockade, as a major aspect of the prostate response to falling androgen levels relies on this phosphatase.
The functional characterization of the 60 genes revealed interesting attributes. DAVID retrieved 17 pathways, most of them centered on either Mapk1 or Mapk9, which is perhaps too general to indicate new physiological functions. It is worth noticing that MAPK pathway has been implicated in increased survival of castrate-resistant prostate cancer patients (Mukherjee et al., 2011).
Next, we uncovered circadian rhythms as a relevant pathway, centered on the genes Arntl and Bhlhe41. These genes were included in network number 2 retrieved by IPA, which also included the Dbp gene, also implicated in circadian regulation. Arntl is directly linked to Dbp and indirectly to Ar via the p300 and CREBBP acetyl transferases. Dbp has been reported to have peak expression at 8 h within the light portion of the 12h:12h light/dark cycle (Zeitegeber, ZT 8) in the rat prostate gland (Qi et al., 2009, Sunkel and in a similar fashion to other core clock genes in the mouse prostate (Bebas et al., 2009). Moreover, ARNTL polymorphisms have been significantly associated with susceptibility to prostate cancer (Zhu et al., 2009) . This evidence raises the possibility that AR-dependent and AR-independent circadian functions contribute to the prostate gland physiology, by opening a new connection to environmental factors, knowingly significant in prostate cancer risk and incidence.
The peroxisome biogenesis pathway, represented by the sole gene Pex14, is also another connection to the environment, as peroxisome proliferation and activity are related to several environmental (and dietary) factors, adding further complexity to the peculiar metabolic adaptations of the gland given its function in accumulating citrate in secretions (Singh et al., 2006).
We also found that some of the genes identified in the present investigation map to regions commonly deleted in prostate cancer. In particular, we refer to the Pten/Hhex/Mxi1 cluster at 10q23, which was characterized in detail before (Hermans et al., 2004). It will be interesting to investigate whether the differently sized deletions in this region might affect the behavior of prostate cancer cells, as MXI1 is usually lacking or inactive in prostate cancer (Eagle et al., 1995, Prochownik et al., 1998; its function is to suppress proliferation by antagonizing Myc (Taj et al., 2001), which in turn is commonly amplified in prostate cancer (Ishkanian et al., 2009). It is important to mention that Mxi1 expression is suggested to decay after castration, according to the parent study (Desai et al., 2004). In contrast, Hhex expression is increased, which functions as a coordinator of hematopoiesis and the development of endoderm-derived organs such as the liver and thyroid (Martinez Barbera et al., 2000).
Additionaly, we found the Mapk1 gene, whose homolog MAPK1 maps to 22q11.21, in a region between the 22q11.21 and 22q12.1 segments deleted in 29% and 33% of intermediate risk tumors, respectively, but not frequently observed in high risk cancers (Ishkanian et al., 2009). The MAPK1 gene product is better known as ERK-2 (or p42 MAPK) and funnels down a variety of extracellular signals to control several functions, particularly the G1-S transition within the cell cycle (Meloche and Pouysségur, 2007). The importance of MAPK1 is highlighted by the fact that it was enlisted as a node in 14 of the 17 pathways identified by DAVID, including prostate cancer among others, and the recent demonstration of the existence of identified mutations in members of the MAPK signaling pathway in the serum of 96% tested individuals harboring different tumors (Bettegowda et al., 2014). Given the particular association of deletions in this region and the intermediate but not high risk of cancer, MAPK1 might be a protooncogene contributing to PCa progression, metastasis and/or transition to CRPC, and its deletion might represent a lower risk of disease progression.
Finally, we found 20 genes mutated in at least 5% of patients. In contrast to PTEN, mutations in these genes are secondary. However, it has been noted that prostate cancer is commonly associated with diverse low frequency mutations (Armenia et al., 2018). Nonetheless, six of the identified genes (PHOX2a, NFYC, EST2, EIF251, SSRP1 and PARP1) were found associated with significant impact on patient survival.
The present analysis cannot distinguish between epithelial and stromal contributions to gene expression. Therefore, it is possible that some of the genes studied are expressed in the stroma. As a matter of fact, a previous approach from our laboratory has identified stromal and epithelial subsets of transcription factors .
In conclusion, this work provides insights into the vastness of physiological pathways involving multiple regulatory interactions among genes needed to adjust prostate biology to the reduced androgen levels achieved by surgical or chemical castration. These results are expected to help us understand the idiosyncrasies of prostate cancer.