Chromosome compartmentalization alterations in prostate cancer cell lines model disease progression

Prostate cancer is the predominant new cancer diagnosis in males in the United States. It is also the second most common cause of male cancer-related deaths, second only to lung cancer (Siegel et al., 2020). Patients with late-stage prostate cancer present with a higher incidence of metastases to trabecular bone (Jacobs, 1983; Bubendorf et al., 2000; Hernandez et al., 2018).

The specific mechanisms that promote metastasis to bone are not understood. However, disseminated tumor cells can be detected in the blood of ∼25% of prostate cancer patients with localized disease. The abundance of these circulating tumor cells positively correlates with metastatic occurrence (Moreno et al., 2005; Danila et al., 2007; Todenhofer et al., 2016).

The human genome's packaging into the nucleus's constrained space requires a systematic organization, from chromosomal territories (Meaburn and Misteli, 2007) to transcriptionally active and inactive chromatin compartments (Lieberman-Aiden et al., 2009). Within compartments, chromatin further organizes into topologically associating domains (TADs), which are segregated from each other by the insulator protein CTCF (Dixon et al., 2012, Nora et al., 2017), and finally into chromatin loops (Nuebler et al., 2018). With the advent of technologies such as genome-wide chromosome conformation capture (Hi-C), it has become clear that these structures are essential for proper gene regulation, DNA replication, and repair (Hnisz et al., 2016; Dekker and Mirny, 2016; Pope et al., 2014; McCord and Balajee, 2018). Further, rearrangement of these domains can impact both the nuclei's ability to squeeze through tight spaces during metastatic migration and the expression patterns of oncogenes (Gerlitz and Bustin, 2010; Stephens et al., 2018; Hnisz et al., 2016; Barutcu et al., 2015). Nuclear atypia is a common diagnostic tool in prostate cancer (Verdone et al., 2015; Diamond et al., 1982), and there is evidence that that the nuclear lamina content of prostate cells differs between normal epithelium, benign prostatic hyperplasia (BPH), and cancer (Partin et al., 1993). This suggests that genome architectural changes may occur in prostate cancer and might influence cancer-promoting gene expression profiles and the nuclear malleability necessary for cancer cells to metastasize.

Mutations, structural alterations, and gene regulatory changes at numerous genomic loci have been associated with a higher risk for prostate cancer (Ahmadiyeh et al., 2010; Helfand et al., 2015; Du et al., 2016). In particular, the TMPRSS2-ERG translocation in chromosome 21 has been found to associate with poor patient prognosis (Zhou et al., 2020; Demichelis et al., 2007; Tomlins et al., 2005). It has been shown that overexpression of ERG results in chromatin conformation changes (Rickman et al., 2012), and FISH experiments combined with irradiation have suggested that such induced chromosome rearrangements alongside DNA breaks can encourage the formation of the TMPRSS2-ERG translocation (Mani et al., 2009). Indeed, it is known that relative chromosome proximity can influence the pattern of translocations that occur (Zhang et al., 2012; Balajee et al., 2018), so this type of local rearrangement could result from chromosome compartmentalization changes that increase contact frequency among the different loci (Engreitz et al., 2012).

Recent studies of the 3D genome structure associated with prostate cancer (Rhie et al., 2019; Luo et al., 2017; Taberlay et al., 2016) have contributed some insight into regulatory chromatin loops, epigenetic alterations, and the influence that variable structures have in transcription. However, these studies did not address whether there are early changes in genome architecture that persist throughout progression. In this context, changes in compartment identity from transcriptionally repressed heterochromatin to euchromatin poised for transcription could identify genes required for early oncogenesis and those necessary for metastasis. Further, changes in TAD positioning or shifting of TAD boundaries could reveal altered interactions of neighboring promoter-enhancer regions.

In the present study, we use Hi-C to characterize the genome organization across a cohort of nine cell lines that model prostate cancer progression from the normal epithelium to bone metastasis, including two bone metastatic cell lines of African American lineage. We further assess the different hierarchical levels of genome organization: from large interchromosomal translocations to compartment identity to TAD location and TAD boundary shifting. We have identified a cohort of 386 genes that change compartment identity across prostate cancer progression. Interestingly, most of these genes switch compartments as proximal clusters. These compartment identity changes are accompanied by distinct structural features at higher resolution such as gain or loss of TAD structure, stalled transcriptional loops and structural deserts, and TAD boundary appearance, disappearance, or positional shifting. Finally, our results revealed several “genomic architecture hotspots” whose structural changes are persistent throughout the metastatic models; these include WNT5A, CDK14, androgen receptor (AR), and the TMPRSS2-ERG locus, among others. These results suggest that the 3D genome structure may serve as a prognostic marker for the progression of prostate cancer to bone metastasis.

To determine how the genome's organization is affected by disease stage, we selected a cohort of nine cell lines that model the progression of prostate cancer, as follows: RWPE1 (Bello et al., 1997) was used to represent normal epithelium. This cell type is epithelial and nontumorigenic, and although it is an artificial cell line model, we verified as described below that its genome structure properties are similar to primary prostate PrEC cells. LNCaP (Horoszewicz et al., 1983), originally isolated from lymph node metastasis, was used as an early adenocarcinoma model. VCaP (Korenchuk et al., 2001) and MDaPCa2a/b (Navone et al., 1997) were used as models for prototypical osteoblastic bone metastasis; these cells are of Caucasian and African American origins, respectively. 22RV1 (Sramkoski et al., 1999) and LNCaP C4-2B (Thalmann et al., 2000) were included in the study as cells that, although isolated from human sources, are models of murine bone metastasis. Atypical metastatic cell lines used in this study include PC3 (Kaighn et al., 1979; osteoclastic and androgen-independent) and Du145 (Stone et al., 1978; brain metastasis; Fig. 1).

Hi-C was performed on cells under normal cell culture conditions (see Materials and methods). Additionally, publicly available datasets for RWPE, PrEC, and LNCaP-C42B were used as comparisons (Rhie et al., 2019, Taberlay et al., 2016). Hi-C mapping and quality control statistics for all samples can be found in Table S1.

The frequency of chromosome contacts is represented in a heatmap where the XY axes are the chromosomal coordinates. The color intensity reflects the frequency with which two particular locations were found to be in contact. At a resolution of 1-Mb bins, the Hi-C heatmaps reveal chromosome territories. At a 250-kb resolution, it is possible to see a plaid pattern of interaction strength, which represents the spatial segregation of A and B compartments. We classify each genomic region as belonging to the A or B compartment using principal component analysis (PCA). Positive values of the first eigenvector (eigen1 or PC1) represent A compartment regions (typically open euchromatin) while negative values denote B compartment (heterochromatin). Finally, at 40-kb resolution, distinct TADs are evident, which are regions of enhanced contacts (which may promote interactions between promoters and enhancers) segregated by the insulator protein CTCF (Fig. 2 A).

Whole-genome contact maps for cells in our model reveal that the highest incidence of interactions occurs in cis: chromosomes primarily interacting within themselves (Fig. 2, B–E). Chromosomal translocations are evident as very strong interactions occurring in trans between different chromosomes. Comparing our Hi-C results with published spectral karyotyping (SKY) data (Pan et al., 1999; van Bokhoven et al., 2003; Roh et al., 2008) and karyotyping information from the American Type Culture Collection (ATCC), we found that our Hi-C data detects 84% of all previously reported translocations. Owing to the high resolution of Hi-C data, we also characterized several previously unreported, smaller translocation events (Fig. S1, Fig. S2, and Fig. S3). Interestingly, chromosomal territories remain well defined throughout progression, with subtle changes in intrachromosomal interactions, as shown in 250-kb resolution heatmaps of each chromosome (refer to Fig. S2 for examples). From these contact maps, it is evident that some chromosomes in some cell lines are broken into multiple pieces (Fig. S4). However, when we perform compartment analysis on each of these broken pieces separately, we find that their underlying A/B compartmentalization is largely similar to cell lines with unbroken chromosomes.

Examination of higher resolution (250-kb) z-score correlation heatmaps of cis interactions for each chromosome shows that the characteristic plaid contact pattern, associated with spatial compartment identity, changes among cell lines (Fig. 3 A). These patterns can be observed throughout all chromosomes, and the changes in compartmentalization are specific to each cell line (Fig. S5). Motivated by an apparent compartment “erosion” in VCaP, we quantified the overall strength of interactions within A compartment regions and within B compartment regions in the different cell types. We find that all prostate cancer cell lines show a loss of A compartment interaction strength relative to RWPE, reflecting an increased intermixing of the A compartment with B compartment regions. Meanwhile, B compartment strength experiences less change in LNCaP and VCaP and even increases in MDAPCa cell lines. Overall, this leads to an imbalance of B and A compartment strength in all prostate cancer cell lines (Fig. S6, A and B).

Systematic analysis of the A/B compartment tracks, per chromosome, for all the cell lines in our model, revealed regions where the compartment identity remained the same and where it was changed (Fig. S5). We note that there is a general level of similarity between epithelial-derived cell lines (Fig. S6 C). By comparing compartment tracks for each cancer cell line to the nontransformed control (RWPE), we find that the changes among cell lines are specific and localized. For example, eigenvector tracks for chromosome 21 show three distinct hotspots of compartment identity switches between the RWPE and cancer cell lines (Fig. 3 B). Interestingly, the rightmost region encompasses the TMPRSS2-ERG locus, a local translocation site that correlates with worse progression and metastasis in prostate cancer (Tomlins et al., 2005; Demichelis et al., 2007; Hägglöf et al., 2014).

A comparison of genome-wide compartment classifications for all cell lines by hierarchical clustering demonstrates that RWPE indeed clusters very closely with primary prostate epithelium PrEC cells (Fig. 3 C). We further demonstrate that RWPE and PrEC compartment profiles are consistent at important genomic locations of interest (Fig. S7). Because RWPE is a reproducible cell culture model, while PrEC cells vary by donor over time, we remain focused on RWPE for the remainder of the comparisons. When we examine the compartment patterns by nearest neighbor analysis (SPRING plot; Weinreb et al., 2018; see Materials and methods), we observe that LNCaP is a central node in a model axis to metastasis (Fig. 3 D, axis highlighted in yellow), connecting RWPE (nontransformed epithelium) to bone metastatic cell lines (VCaP and MDAPCa). Interestingly, the compartment patterns of MDAPCa2A and 2B, which were derived from an African American patient, is distinct from VCaP, which is of Caucasian lineage. This analysis also clusters atypical metastatic lines DU145 and PC3 together along a third axis radiating from the LNCaP central node.

To consistently mathematically classify regions that significantly changed compartments between cell lines, PC1 values per bin per cell line were normalized by subtracting the corresponding value derived from the RWPE cell line. Significant changes in compartment identity $(|X¯+1.5σ|)$ were identified and the corresponding genomic areas annotated. Interestingly, genomic regions that switch compartment in the LNCaP-RWPE comparison tend to have a switched identity in bone metastatic cell lines as well (Fig. 3 E). These regions thus change compartment status early in progression and then maintain their new spatial localization rather than switching back or experiencing progressively more shift in compartment. Overall, we have identified 181 genomic bins (250 kb in size) whose compartmentalization changes from the B to the A compartment in either (1) six or more cancer cell lines compared with RWPE or (2) all cells in the progression axis, compared with RWPE. These genomic bins contain 300 genes that are therefore moving from the B to the A compartment (Table S2), suggesting that these genes become more likely to be transcribed. These include AR, WNT5A, CDK14, and genes located close to TMPRSS2, such as BACE2. The majority of these genes (256) are grouped in 48 proximal clusters (Table S3). We have also identified 86 genes whose genomic loci switch compartments from the A to B compartment, suggesting a genome structure rearrangement more likely to result in gene silencing (Table S2). Such is the case for certain cadherins, annexins, and mediators of inflammation. Of these genes, 64 are grouped in 16 clusters (Table S3).

A comparison between the bone metastatic cell line (VCaP) and the adenocarcinoma cell line (LNCaP) revealed that 345 more bins switch from the B to the A compartment, comprising 313 genes. Interestingly, these events appear to be VCaP-specific, suggesting that the MDA bone metastatic lineages undertake a different pathway to oncogenesis.

Since a switch from the B to A compartment could make a region permissive for transcriptional activation, we used the Clariom-S microarray to profile expression levels in the targets identified. Using the expression level of RWPE as a baseline, we have found that 47% of genes that switch from B to A compartment show a significant transcriptional induction in LNCaP (twofold or higher). In turn, 49% of those genes are even further up-regulated in VCaP compared with LNCaP, suggesting an exacerbation of this expression pattern with progression. These genes include, among others, AR, TMPRSS2, CDK14, and WNT5a (Fig. 4 A). In contrast, only 18% of the genes that show induction in LNCaP are induced further in the MDAPCa cell lines (Fig. 4 B).

Finally, we found 25 genes underwent a B to A compartment switch, were up-regulated in LNCaP versus RWPE, and were further up-regulated in both MDAPCa cell lines and VCaP relative to LNCaP. This set includes CDK14, WNT5a, and BACE2, which is a close neighbor of TMPRSS2 (Fig. 4 C). These data suggest that these transcriptional hotspots are shared in osteoblastic metastatic cells and might be necessary for colonization and survival in a secondary bone site. Interestingly, while amyloid precursor protein (APP), the proteolytic target of β secretase 2 (BACE2), does not change compartment between RWPE and LNCaP, it does switch both compartments from B to A and shows an expression increase when we compare VCaP to LNCaP.

TADs contribute to the 3D architecture of the genome by sequestering enhancers and promoters with their target genes, with a low likelihood of interaction across boundaries (Lupiáñez et al., 2015; Nora et al., 2012; Dixon et al., 2012). The disappearance of TAD boundaries has been identified as a potential activator of oncogenes (Hnisz et al., 2016). To analyze whether there are changes in 3D genomic structure at the local level surrounding genes associated with compartment identity switches, we used higher-resolution Hi-C heatmaps (40 kb).

The WNT5a-ERC2 locus (Fig. 5, Cluster 1) is in the B compartment in nontransformed control, switching to the A compartment in all four members of the model axis to progression: LNCaP, VCaP, and both MDAPCa cell lines. While the right TAD boundary location is relatively consistent, the location of the left boundary shifts positions: in LNCaP, DU145, MDAPCa2B, and PC3, WNT5A localizes in a TAD with LRTM1 instead of ERC2. Interestingly, in VCaP, these three genes are clustered in the same TAD. FZD1 and CDK14, a receptor and activator cyclin of noncanonical WNT signaling, and which have been associated with the function of WNT5A, also switch to the A compartment (Fig. 5, Cluster 2). In this case, the genes are separated by a TAD boundary that rests atop the gene body of CDK14, but increased loops are evident on the separate TADs in the cancer cell lines, compared with RWPE. Remarkably, CDK14 is one of the genes whose transcription is consistently up-regulated in all metastatic cell lines in the progression axis. TAD boundary shifting, appearance, and disappearance were also observed in other clusters that switch compartment identity (Fig. S8).

A detailed survey of all compartment-switched loci revealed five possible categories of changes in local structure. We classified all compartment switch regions into these categories. Strikingly, the AR locus displays all five types of structure change across the different cell lines. First, a highly disorganized area (no TADs, uniformly distributed interactions across a region of the heatmap) becomes highly organized, or vice versa (Fig. 6 A). In AR, the gene is in the B compartment in RWPE, and the area around it is highly disorganized. With progression, the area becomes organized into TADs and sub-TADs for both LNCaP C4-2b and VCaP. The third type of structural conformation is the appearance of a high incidence of interactions along the edge the TAD (Fig. 6 B). We call these interactions “loops” or “stalled loops” because these have previously been associated with the phenomenon of cohesin becoming stalled as it extrudes loops but encounters RNA polymerase, CTCF, or other barriers. This increased loop formation is evident at the AR locus in LNCaP, MDAPCa2a, and MDAPCa2b. The fourth structural feature is a complete absence of structure (loss of contacts in the entire region) associated with the gene of interest. This “structural desert” is observed at the AR locus in the 22RV1 cell line (Fig. 6 C) and is reminiscent of previously observed structural features at highly transcribed long genes (Heinz et al., 2018; Leidescher and Nübler, 2020 Preprint). Finally, the local structure can remain unchanged. For AR, this is the case in PC3 and DU145, where the gene remains in a disorganized region of the B compartment (Fig. 6 D).

Overall, for those genes that switch from the B to the A compartment, all types of changes happen with fairly even probability: 12.64% lose structure at the local level, 14.37% gain structure, 28.16% present stalled loops, 21.26% associate with structural deserts, and 23.56% do not change. However, genes located in the A to the B compartment switches predominantly lose structure (39.06%) or do not change (21.88%). Only ∼3% of A to B switch loci gain structure. The remaining loci are distributed evenly between stalled loop areas and structural desert change (Fig. 6 E). From this survey of prostate cancer cell lines, therefore, we also gain basic insight about the types of local structure change that most often accompany compartment level changes.

As previously mentioned, the incidence of the TMPRS2-ERG translocation has a positive correlation with prostate cancer progression and a poor prognosis. For the TMPRSS2-ERG locus, we observe that the region that encompasses TMPRSS2, BACE2, PLAC4, MX1, MX2, and FAM3B switches from the B compartment in RWPE to the A compartment in all cancer cell lines. In the model cell line for nontransformed epithelium RWPE, this genomic area is enclosed in a large TAD, including genes downstream of TMPRSS2 (Fig. 7 A), as is also the case for LNCaP, 22RV1, MDAPCa2a, and PC3, with a slight shifting of the TAD boundary location. In LNCaP-C42B, an osteoblastic cell line derived from LNCaP, a sub-TAD appears within this locus, sequestering the MX1 gene. Sub-TAD fragmentation also occurs in DU145, MDAPCa2b, and VCaP. In contrast, ERG is in the A compartment in RWPE, switching to the B compartment only in LNCaP and MDAPCa2b. The TAD boundary location around ERG does not change. Of note, the whole genomic area surrounding these two clusters is located in the A compartment in VCaP, a known carrier of the TMPRSS2-ERG translocation, evident as a high-interaction location in the Hi-C heatmap (Fig. 7 A, VCaP, rectangle).

In the bone metastatic cell line VCaP, the TMPRSS2-ERG translocation is apparent in the Hi-C heatmap of the long arm of chromosome 21 at a 250-kb resolution (Fig. 7 B). While a translocation is not evident in the MDAPCa cell lines, there is a distinct higher incidence of interactions close to that area: a log2 ratio comparison against control (RWPE) shows that that chromosomal region is enriched for interactions in adenocarcinoma (LNCaP; Fig. 7 C), and that this phenotype is aggravated in both bone metastatic MDAPCa cell lines.

To evaluate how the findings of our cell line–based model relate to patient data, we examined the expression data from the Prostate Adenocarcinoma cohort (PRAD) of The Cancer Genome Atlas (TCGA). A total of 549 patients were taken into consideration, comprising 52 normal and 497 tumor samples, stratified by Gleason score, an indication of cancer severity. First, we examined all genes in regions that switch from the A to the B compartment across our cell lines (as in Table S2). We observe a significant decrease in the pooled set of genes in tumor samples versus normal samples overall (Fig. 8 A). We find that 55% of genes in this set show significant decreases in tumor versus normal samples, many of which are relevant to prostate cancer (Fig. 8, B and C). This is a higher proportion of down-regulated genes than would be expected at random, as in 200 random gene sets of the same size, we never found any sets in which 55% or more of the genes were down-regulated in tumors (Fig. 8 C).

In contrast, when we consider the set of genes that switch from B to A (as in Table S2), we observe that they are actually not up-regulated more often than would be expected at random in PRAD tumor samples (Fig. 8 D). This is consistent with our model that movement into the A compartment is not synonymous with gene activation, but instead only provides an opportunity for genes to be activated. But when we looked at the subset of genes which both change from the B to A compartment and change expression in our cell line models, we found that a more significant proportion of genes (47%) is up-regulated in patient tumor samples (Fig. 8, E and F). Only 6.5% of 200 random gene sets of the same size show this proportion of up-regulation. But some key genes in this set, such as WNT5A and HUNK, clearly do not show an increase in expression in patient tumor samples versus normal overall (Fig. 8 E). Instead, these genes show increasing expression with increasing Gleason score, a measure of more aggressive disease (Fig. 8 G). Indeed, most of the genes (75%) in this B to A switched subset that did not change expression in overall tumor versus normal instead show an increase in expression over cancer progression (Fig. 8 H). Thus, overall, a significant proportion of genes that switch from the B to the A compartment and changed expression in our cell line models either increase expression in patient tumors in general or increase their expression gradually with worsening disease in patients.

The deliberate, hierarchical organization of chromatin within the eukaryotic nucleus’s constrained space is necessary for adequate DNA maintenance, repair, and gene transcription or silencing, all of which contribute to the cell's homeostasis. In this study, we have characterized the genomic architecture in a cohort of cells that model prostate cancer progression. At the genome-wide level, given the inherent high resolution of Hi-C data, we have identified small translocations that, to our knowledge, have not previously been reported in the literature.

Many regions across the genome are unchanged in their spatial compartmentalization across cell lines. This suggests that there are inherent 3D genome structure features of prostate epithelium that arise during initial differentiation and tissue patterning (Flyamer et al., 2017; Ke et al., 2017) and that these are persistent, regardless of malignancy status. These results suggest that genomic loci that switch compartment identity between the normal epithelium and cancer cells are associated with an oncogenic genomic architecture profile and that those features result from concerted biological events.

It is noteworthy that the majority of the compartment changes we identify involve a switch from the B to the A compartment. Correspondingly, we find in our compartment strength analyses that the A compartment becomes more intermixed, interacting more broadly, while the strongest B compartment regions remain more spatially segregated in the cancer cell lines. Both of these results point to a general shift in the prostate cancer lines toward a more open/poised-for-activation chromatin environment, which could lead to misactivation of oncogenes.

Our genome-wide compartment analysis revealed that compartment identity alone is enough to stratify prostate epithelium in a continuum of progression. LNCaP is a central node that connects the nontransformed epithelium (RWPE) to bone metastatic cell lines (MDAPCa2a/b and VCaP). These results provide further insight into the importance of cell line selection in progression studies: the atypical metastatic cell lines DU145 and PC3 cluster together away from the primary axis. This emphasizes that as researchers select cell lines for study, it is important to consider the differences in these cell lines and that not all will capture the most common pathways of metastasis.

It is important to consider that the compartment switching events observed here involve clusters of genes, which in some cases expand through several TADs. This observation echoes previous results in a plant system, where clusters of genes often changed their spatial compartmentalization and expression together (Nützmann et al., 2020). Interestingly, the clusters we observe include both genes previously identified in prostate cancer progression and others seemingly unrelated. Such is the case of the WNT5a locus, which is a known target of prostate tissue development, patterning, and cancer aggressiveness (Allgeier et al., 2008; Dai et al., 2008; Yamamoto et al., 2010; Huang et al., 2009). Our results show that in RWPE, WNT5a is sequestered within a TAD with the ELKS/RAB6-interacting/CAST family member 2 (ERC2). In contrast, TAD boundary shifting or eviction in the proximity of WNT5a in the bone metastatic cell lines represented in our primary axis for progression (VCaP and MDAPCa2a/b) results in TAD-limited interactions with leucine-rich repeats and transmembrane domains-containing protein 1 (LRTM1) instead. While this shift does not result in transcriptional induction of LRTM1, it could potentially lead to aberrant interactions between the promoters for both genes, ultimately resulting in the overexpression of WNT5a, as is also observed with increased Gleason score in patients. Since TAD structure is essential for proper gene regulation (Lupiáñez et al., 2015; Rhie et al., 2019; Guo et al., 2018), this phenomenon requires further exploration. Still, it is attractive to consider TAD-targeted therapies that hold the potential to reverse the deleterious effects of TAD shifting. Interestingly, another target of noncanonical WNT signaling, CDK14 (reviewed by Davidson and Niehrs, 2010), also switches from the B to the A compartment and is transcriptionally activated in our metastatic axis.

The clustering of different genes in the described compartment switches raises the interesting question of whether certain genes could act as a driver of compartment identity switches while neighboring genes act as “passengers.” For example, the transcriptional activation of one gene could influence a whole genomic region to switch to the A compartment. Indeed, previous work has shown that binding of transcriptional activators can prefigure spatial compartment alterations (Stadhouders et al., 2018; Therizols et al., 2014). This compartment switch would then result in the switching of neighboring genes, which may increase their probability of later becoming activated as well. Recent work has shown that spatial reorganization of a chromosome region can make it more permissive for de-repression, even if the structural switch does not immediately change its expression level (Manjón et al., 2021,Preprint). Is it possible that earlier transcriptional events required for the cell to survive a particular insult trigger a full compartment shift? Is this, in turn, a potential trigger for oncogenic transcriptional activity? Such seems to be the case of the observed APP–BACE2 axis along chromosome 21. APP presents with enhanced expression in the LNCaP–MDAPCa–VCaP axis compared with RWPE (Fig. S9). APP is also proximal to one of the persistent compartment switches in chromosome 21, but it does not change compartment itself: the A compartment identity atop the gene gets stronger with progression. Mounting evidence from the Alzheimer's field, where abnormal amyloid processing results in aggregation and neurodegeneration, has shown that this protein and derived peptides serve a crucial role as antimicrobials and are necessary for mounting an appropriate host response to infection (thoroughly reviewed by Moir et al., 2018). If an environmental signal such as infection results in the prostate epithelium being exposed to excessive or chronic APP, could this trigger expression of its processing enzyme BACE2? Our evidence of the BACE2 gene switching compartmentalization and experiencing higher levels of transcription in all three prototypical bone metastatic lines suggests this may be so. Critically, these events would imply the need for the compartment switch around the TMPRSS2 locus.

Since its discovery, the TMPRSS2-ERG translocation in chromosome 21 has been recognized as an important indicator of poor prognosis and a higher risk of prostate cancer–related death (Tomlins et al., 2005; Perner et al., 2006; Demichelis et al., 2007; Hägglöf et al., 2014; Deplus et al., 2017). While a significant body of work has explored the role of ERG overexpression, its influence on AR signaling, and its association with chromatin structure (Li et al., 2020, Rickman et al., 2012), the genesis of this translocation and fusion remains poorly understood. Here, we implicate a 3D genome organization change of compartmentalization in this event. As mentioned before, TMPRSS2 switches from the B to the A compartment consistently in all prostate cancer cell lines we queried. Meanwhile, ERG switches from the A to the B compartment in LNCaP and its nearest neighbor in our experimental metastatic axis, MDAPCa2b. Remarkably, ERG remains in the A compartment for VCaP and MDAPCa2a, suggesting a potential dual switching event during progression to more aggressive phenotypes. All these spatial rearrangements in such proximity result in the enhanced contact incidence observed in all of these cell lines, even the early adenocarcinoma model LNCaP (as described in Fig. 7), and ultimately could lead to the translocation event, as it has been shown in other systems (McCord and Balajee, 2018, Zhang et al., 2012, Baca et al., 2013).

Recent efforts (Hawley et al., 2021 Preprint) have characterized 3D genomic profiles in prostate tumor cohorts. These studies recapitulate our findings that the 3D genome organization between malignant and benign prostate tissues remains largely consistent. We propose that prostate cancer progression is associated with specific changes in the 3D genome structure that arise early in the disease and facilitate an oncogenic expression phenotype. Based on these results, we can hypothesize that analyzing the 3D genome structure of patient-derived samples could be a prognostic marker for progression and bone metastasis.

RWPE-1, LNCaP, DU145, 22RV1, VCaP, and PC3 cell lines were obtained from the Physical Sciences Oncology Network Bioresource Core Facility, supported by ATCC (Manassas, VA).

All cell lines were cultured according to standard protocols, subculturing cells as they reached 80% confluency with the following media formulations: RWPE media was comprised of keratinocyte-specific media supplemented with EGF and bovine pituitary extract (Gibco; 17005042). RPMI (Gibco; 11835030) was supplemented to match the suggested formula by ATCC (4.5 g/liter glucose, 2.383 g/liter Hepes, and 0.11 g/liter sodium pyruvate) and 10% FBS (Corning; 35–010-CV). These media was used for LNCaP, 22RV1, and PC3. DU145 cells were cultured in DMEM supplemented with 10% FBS, and DMEM F12: Ham 1:1 (Gibco; 11–320-033) with 10% FBS was used for VCaP.

Cell lines MDAPCa2a and MDAPCa2b were a kind gift of Dr. Nora Navone (MD Anderson Cancer Center, Houston, TX). The identity of the received cell lines was verified by ATCC Human Cell STR testing services. These cell lines were grown in HPC1 media (Athena Enzyme Systems) supplemented with 10% FBS on standard cell culture T75 flasks coated with FCN Coating Mix (Athena Enzyme Systems) as per the manufacturer’s instructions.

All media formulations were supplemented with 100 µg/ml penicillin-streptomycin (Gibco; 15–140-122).

Cell pellets for Hi-C were prepared as previously described (Golloshi et al., 2018). Briefly, cells growing in monolayer in standard T75 flasks were quickly washed with 10 ml of HBSS (Gibco; 14–025-134) at room temperature and cross-linked with 10 ml 1% formaldehyde (Fisher Bioreagents; BP531-25, in HBSS) for 10 min on a shaking platform. The cross-linking reaction was quenched by adding glycine to a final concentration of 0.14 M (MP Biomedicals; ICN19482591), followed by a 5-min incubation at room temperature, with shaking. After cooling down the plates on ice for 15 min, the formaldehyde solution was aspirated from the plate and substituted with 10 ml of ice-cold HBSS containing 1× Halt Protease Inhibitor cocktail (Thermo Fisher Scientific; PI78438). 5 million cell aliquots were collected by centrifugation and snap-frozen in liquid nitrogen.

VCaP, MDAPCa2a, and PC3 Hi-C was conducted as previously described (Golloshi et al., 2018). Briefly, cross-linked cells were lysed in cold hypotonic buffer, followed by dounce homogenization, and then nuclei were permeabilized by incubation in 0.1% SDS. After the permeabilization reaction was terminated by a final concentration of 1% Triton X-100, chromatin was digested overnight using 400 U DpnII (New England Biolabs) at 37 C. Digested ends were filled in with biotin-dATP by Klenow DNA polymerase, and the blunt ends of interacting fragments were ligated together using T4 DNA ligase. DNA was then purified by phenol-chloroform extraction and ethanol precipitation. The ligated DNA was sonicated with a Covaris sonicator and then size-selected using Ampure XP beads to an average size of 200–400 bp. Biotinylated interacting fragments were pulled down using MyOne Streptavidin C1 beads (Invitrogen). For library preparation, the NEBNext Ultra II DNA Library prep kit (NEB) was used to carry out End Prep, Adaptor Ligation, and PCR amplification on streptavidin bead–bound DNA libraries.

Hi-C was also performed on LNCaP, DU145, VCaP, MDAPCa2a, MDAPCa2b, and PC3 using the Arima Hi-C (Arima Genomics) kit, following the manufacturer’s protocol A160141 v01 for library amplification using the NEBNext Ultra II kit (NEB; E7645S). Sequencing was performed by Genewiz on either an Illumina NovaSeq or HiSeq platform with 50 or 150 bp paired-end reads. Sequenced reads were mapped to a reference human genome (hg19), binned, and iteratively corrected according to established pipelines (Imakaev et al., 2012) using the dekkerlab-cMapping tool available at https://github.com/dekkerlab/cMapping.

In addition, the same analysis above was performed on fastq files from previously published Hi-C datasets for RWPE, 22RV1, and LNCaP C4-2B: GSE118629 and GSE73782 (Rhie et al., 2019; Luo et al., 2017). See Table S1 for all data sources and statistics. Newly generated LNCaP results were checked for consistency with previously published LNCaP Hi-C results (ENCSR346DCU; Taberlay et al., 2016). Since our LNCaP data had a dramatically higher cis/trans ratio than these previously published datasets, comparisons were difficult, and so we proceeded with only our newly generated data.

All Hi-C data analysis was performed using the existing cworld-dekker pipeline, available on github (https://github.com/dekkerlab/cworld-dekker) as follows.

Hi-C heatmaps were generated for genome-wide datasets at a 2.5-Mb resolution, and per chromosome at a 250-kb resolution, using the heatmap script.

Compartment analysis was performed via PCA on 250-kb binned matrices using the matrix2compartment script on each chromosome individually. Positive and negative PC1 values are initially arbitrary, so this script uses gene density to determine which set of bins (positive or negative) should be assigned as A (higher gene density) or B (lower gene density), respectively. PC1 profiles were visually inspected to ensure that they reflected a compartment-like pattern rather than a separation of chromosome arms or structural variant breaks. On some highly fragmented chromosomes, breaks between regions hindered proper compartment detection. In these cases, we detected breakpoints by finding extreme minima in insulation score profiles (Crane et al., 2015), creating separate matrices for each unbroken chunk of chromosome separately, and then running compartment analysis on each sub-region of the chromosome (Fig. S4).

The PC1 values per bin for the control cell line (RWPE) were subtracted from the values from each cancer cell line, resulting in a normalized distribution (referred to as ΔPC1analysis). Significant changes in compartment identity were defined as those bins whose subtracted value fell either under the mean minus 1.5 the SD or the mean plus 1.5 the SD, for at least six cell lines or all four cell lines in the primary axis, as defined by nearest neighbor analysis as described below (SPRING plot).

The genes contained in regions of interest determined from the ΔPC1 analysis were annotated using the knownGene primary table as referenced in the University of California Santa Cruz (UCSC) Genome Browser (https://genome.ucsc.edu/).

To facilitate the visualization of chromatin compartmentalization, heatmaps for cis interactions were generated by first calculating the z-score of the interactions at 250-kb resolution compared with an expected interaction at each distance from the diagonal and then taking the Pearson correlation of each row and column of the heatmap (z-score correlation matrices). Multi-compartment track figures and overlays were constructed using the visualization tool Sushi (Phanstiel et al., 2014)

Genome-wide TAD boundaries were determined using the matrix2insulation script on 40-kb binned matrices, following the insulation score approach with an insulation square size of 500 kb (Crane et al., 2015).

To calculate A-A and B-B compartment interaction strengths for each chromosome, distance corrected Hi-C intra-chromosomal interaction frequencies at 250-kb resolution were reordered according to their corresponding PC1 values (from strongest B to strongest A). Then, the reordered intra-chromosomal interaction matrix was smoothed at 500-kb resolution. Interactions were classified as A-A, B-B, and A-B and thresholded to include only the top 20% of interactions. The median value of each A-A, B-B, and A-B interaction is calculated.

Finally, the relative A-A and B-B compartment interaction strengths were obtained by subtracting the absolute A-B compartment interaction strength from the absolute A-A and B-B interaction compartment strengths, respectively.

For the average compartment interaction strength, the mean of the relative A-A and B-B interaction compartment strengths was calculated. A stronger A-A or B-B compartmentalization level compared with the A-B compartment intermixing would produce a higher positive value. On the other hand, a value close to zero suggests a weaker level of compartmentalization.

Saddle plots of compartment interaction strength were created by reordering a given chromosome matrix according to the PC1 values, from most negative (B compartment) to most positive (A compartment). Then, the interactions between each pair of bins were represented by the z-score, which compares the interaction frequency observed to the expected interactions between any two bins separated by that genomic distance. We subtract saddle plots from different cell types to observe the change in compartment interaction strength between them.

To construct the spring plot (Weinreb et al., 2018; Miura et al., 2020) of the PCa cells based on the Hi-C compartmental data, the compartment profile of the cells at the 250-kb resolution was binarized. For example, the genomic regions where the compartment strength are >0 were converted to 1 and the negative strengths to −1. The reason behind using that discretization step was to only consider the A/B compartment signature irrespective of the compartment strength. Once binarized, the data were then organized in a matrix format, where the rows are the genomic regions and columns are cells, and PCA was performed on that data. Since we were dealing with a small set of samples (cells) for our analysis, we kept all the principal components from the PCA transformation for further analysis. Then, a k-nearest neighbor graph with two nearest neighbors (includes the node itself) was constructed from the PCA-transformed data, and the network was visualized with a force-directed layout.

RNA was purified from 5 million cells at three different passages, per cell line, using the RNEasy mini kit (Qiagen; 74104) using QIAshredder (Qiagen; 79654) for homogenization. Purification was followed up by cleanup, and concentration using the RNase-Free DNase Set (Qiagen; 79254) and QIAquick PCR Purification Kit (Qiagen; 28104), respectively. RNA concentration was determined using a NanoDrop One (Thermo Fisher Scientific).

Clariom S microarray was performed through the Transcriptome Analysis Services (Transcriptome Profiling) from Thermo Fisher Scientific.

In addition, for cross-validation purposes, data for HG-U133 Plus2 microarray was also collected from the Encyclopedia of DNA Elements (ENCODE), as follows: RWPE (GSM966512, GSM966513, GSM966514), LNCaP (GSM2571978, GSM2571979, GSM2571980), DU145 (GSM1374469), PC3 (GSM1517530, GSM1517531, GSM1517532), LNCaP-C4-2B (GSM1565257, GSM1565258), and 22RV1 (GSM2571966, GSM2571967, GSM2571968).

Analysis of microarray data was performed using the Transcription Analysis Console from Applied Biosystems (Thermo Fisher Scientific), curating the log2 fold up-regulated/down-regulated genes (P < 0.05) with targets identified in the ΔPC1 analysis.

Gene expression of all the genes identified as switching compartments in our model cell lines (Table S2) was analyzed on the TCGA-PRAD using the USCS Xena Browser interface (Goldman et al., 2020). Data for normalized gene counts for all samples were downloaded from https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.PRAD.sampleMap%2FHiSeqV2.gz. Samples were annotated with Gleason score labels provided by TCGA. With custom R scripts, we statistically evaluated the expression change between tumor and normal samples for each individual gene (one-sided t test) and calculated the proportion of a given gene set with changes (P < 0.05) in the direction we expected (that is, up-regulation for B to A switched gene sets and down-regulation for A to B switched gene sets). We then selected random gene sets of the same size as our query set for comparison. 200 independent sets of random genes were selected for each comparison (that is, if we were considering 87 genes that switched A to B, we would choose 200 sets of 87 randomly selected genes). The proportion of genes called significantly up or down-regulated (P < 0.05 in one-sided t test) was calculated in each of the 200 random gene sets. The fraction of random gene sets containing the same or greater proportion of up or down-regulated genes was used to evaluate whether the amount of up- or down-regulation observed in the compartment switch gene sets was likely to be seen by chance. Note that by comparing proportions of genes that met a certain statistical threshold between real and random sets, we account for problems of multiple hypothesis testing. On the log count scale, most of the gene expression distributions are not highly skewed by outliers, so the mean expression was a representative metric for evaluation by the t test. We did not test for normality of the gene expression data distributions before using the t test because the t test does not require the underlying data to be normally distributed (Boneau, 1960; Norman, 2010). Instead, the t test assumes normality in the distribution of the sample means, a condition that is met, according to the central limit theorem, for sample sizes larger than 5–10, which we exceed in this case (Norman, 2010). Slopes of gene expression across Gleason score were calculated in R with a linear fit. The proportion of genes in a set with positive slope was evaluated and compared with the proportion of genes with positive slope in 200 sets of a matched number of genes chosen at random.

Fig. S1, Fig. S2, and Fig. S3 show that Hi-C data at 2.5-Mb resolution recapitulate previously reported and novel translocations in the cell lines in this study. Fig. S4 shows how we modified our compartment analysis to consider heavily fragmented chromosomes. Fig. S5 shows a comparison of Hi-C heatmaps, cis interaction heatmaps, and compartment analysis for chromosome X for RWPE, LNCaP, and VCaP. Fig. S6 characterizes differences in genome-wide chromosome compartmentalization for all cells in this study. Fig. S7 validates the consistency of the genome-wide compartmentalization between the RWPE cell line (used as control for this study) and PrEC, a human primary prostate epithelial cell line. Fig. S8 includes examples of gene clusters that switch compartments in advanced prostate cancer. Fig. S9 highlights the APP-BACE2 region, located on chromosome 21, and its potential role as a driver in compartment switching of neighboring genes. Table S1 includes the database reference and quality control metrics for this study. Table S2 lists the genes identified as switching compartment identity. Table S3 denotes the gene clusters that switch compartment identity, as related to prostate cancer progression.

All Hi-C and microarray data are available on GEO at accession no. GSE172099. Other processed data figures and a UCSC Genome Browser Track Hub containing all data are available for browsing at https://3dgenome.utk.edu/3d-genome-architecture-in-prostate-cancer-progression/.

We thank Jeremy Hughes (Web Communications Manager, University of Tennessee, Knoxville. Department of Arts & Sciences) for his help in constructing the accompanying website. We thank Nora Navone, Ph.D. (MD Anderson Cancer Center, Houston, TX) for providing the MDAPCa cell lines.

This work was supported, in part, by National Institutes of Health National Institute of General Medical Sciences grant R35GM133557 to R.P. McCord. R. San Martin was supported by a postdoctoral fellowship from the American Cancer Society (134060-PF-19-183-01-CSM).

The authors declare no competing financial interests.

Author contributions: Conceptualization: R. San Martin and R.P. McCord; investigation: R. San Martin and R. Golloshi; formal analysis: R. San Martin, R.P. McCord, P. Das, J.T. Sanders, and R. Dos Reis Marques; visualization: R. San Martin, R.P. McCord, P. Das, Y. Xu, and R. Dos Reis Marques; resources and validation: R. San Martin and J.M. Roberts; writing – original draft: R. San Martin and R.P. McCord; software: P. Das, Y. Xu, and R.P. McCord; writing – review and editing: all authors; supervision: R.P. McCord.