3-D vascularized breast cancer model to study the role of osteoblast in formation of a pre-metastatic niche

Breast cancer cells (BCCs) preferentially metastasize to bone. It is known that BCCs remotely primes the distant bone site prior to metastasis. However, the reciprocal influence of bone cells on the primary tumor is relatively overlooked. Here, to study the bone-tumor paracrine influence, a tri-cellular 3-D vascularized breast cancer tissue (VBCTs) model is engineered which comprised MDA-MB231, a triple-negative breast cancer cells (TNBC), fibroblasts, and endothelial cells. This is indirectly co-cultured with osteoblasts (OBs), thereby constituting a complex quad-cellular tumor progression model. VBCTs alone and in conjunction with OBs led to abnormal vasculature and reduced vessel density but enhanced VEGF production. A total of 1476 significantly upregulated and 775 downregulated genes are identified in the VBCTs exposed to OBs. HSP90N, CYCS, RPS27A, and EGFR are recognized as upregulated hub-genes. Kaplan Meier plot shows HSP90N to have a significant outcome in TNBC patient survivability. Furthermore, compared to cancer tissues without vessels, gene analysis recognized 1278 significantly upregulated and 566 downregulated genes in VBCTs. DKK1, CXCL13, C3 protein and BMP4 are identified to be downregulated hub genes in VBCTs. Together, a multi-cellular breast cancer model and culture protocols are established to study pre-metastatic events in the presence of OBs.

Breast cancer is the most frequently diagnosed cancer and the leading cause of cancer-related death in females worldwide 1 . Once the tumor metastasizes to a distant site, the 5-year survival chance decreases from ~ 93 to ~ 22% (National Cancer Institute SEER database). There are five molecular subtypes of breast cancer of which the triple negative subtype is particularly challenging 2 . Triple negative tumors represent around 15% of invasive ductal breast cancers. They display distinctive epidemiological, phenotypic and molecular features with distinctive patterns of relapse, and a poor prognosis despite relative chemosensitivity 2 . Increasing clinical evidence has suggested that the most common site of triple negative breast cancer (TNBC) metastasis is bone [3][4][5][6] , which disrupts the balance between bone formation and resorption 7 . Patients with this condition have a median survival of around 2-3 years following initial diagnosis of bone involvement.
To understand better the propensity of breast cancer cells (BCCs) to metastasize preferentially to bone, an investigation of the mechanistic pathways that drive organ-specific metastasis is crucial. Conventionally, researchers have relied on two-dimensional (2-D) assays and animal models in order to study the metastatic process 8 . However, 2-D models do not recapitulate the tissue architecture and the cellular interactions of the tumor microenvironment. On the other hand, animal models have failed to translate pre-clinical success to patient outcome, thereby questioning its resemblance in the cellular level to the human tissue. With animal models, it is difficult to study the interaction of specific cell types in a controlled manner thereby limiting the possibility to conduct parametric studies. Additionally, ethical issues in using animals as models for human disease exist. An alternative approach is the use of three-dimensional (3-D) cultures such as hydrogels composed of either ECM molecules [9][10][11] or synthetic polymers 12 , scaffold-free organoids culture 13,14 , and microfluidic chips 15,16 .
Although there are a large number of prior studies focused on the fabrication of sophisticated 3-D in-vitro models utilizing the above-mentioned techniques, most models are not developed for the investigation of cancer metastasis progression. Out of the limited number of research articles that deal with breast-to-bone metastasis, the majority of these studies recapitulate the late-bone metastatic stage where the BCCs have already colonized the bone microenvironment, leading to further metastasis or dormancy 8,17,18 . In other words, these studies focus on post-metastatic, direct interaction between BCCs and bone. However, to understand the metastasis process in its entirety, the study of early-stage pre-metastatic progression of BCCs within the primary site is equally important.
One overlooked pre-metastatic phenomenon is the involvement of the distant secondary site in remotely regulating the primary site. Bone could act as an active participant, releasing cytokines and growth factors to the primary breast tumor site, and play a functional role in the initiation of tumor metastasis 19 . The concept that primary tumors can prime the secondary site remotely prior to the initiation of metastasis is not controversial 20,21 . However, the reciprocal crosstalk from the secondary site, although understudied, could also influence changes in the primary site. It is intriguing to note that high bone mineral density (BMD) in postmenopausal women has been linked with an increased incidence rate of breast cancer 22,23 . Furthermore, it was shown that low levels of bone turnover proteins c-terminal telopeptide (CTX1, bone resorption marker) and osteocalcin (bone formation marker), independent of BMD, correlated with increased breast cancer risk in postmenopausal women 24 . Although these phenomena are not fully understood and require further investigation, it begs the question of whether the bone can indirectly, via release of soluble factors, lead to breast cancer metastatic initiation and progression 19 .
To mimic this early-stage pre-metastatic cross-talk scenario in-vitro, there is a need to develop a 3-D indirect culture protocol where the secondary and primary sites are not in direct contact but interact via indirect paracrine signaling. Due to the lack of such models and protocols, the study of organ specificity (organotropism) observed in metastasis remains a challenge and several open questions like the role of future secondary sites in conditioning the primary tumor persists. The goals of our study were therefore to engineer a 3-D vascularized breast cancer tissue (VBCTs) representing the primary tumor site indirectly co-cultured together with OBs, representing the secondary site and to investigate whether and how this indirect interaction with bone cells could lead to alterations in the primary breast tumor microenvironment. Therefore, a quad-culture system was developed consisting of TNBCs (MDA-MB231), fibroblasts, and endothelial cells (HUVECs) for indirect coculture with OBs, the bone-forming cells. We opted for a scaffold-free model system due to its several advantages over scaffold-based models 25 . To engineer the scaffold-free models, we utilized layer-by-layer ECM coating and accumulation technique 26 . This technique allowed the formation of thick vascularized tissues and provided a microenvironment to the cancer cells with enhanced direct cell-cell interactions 26,27 . Fabrication of the models in cell-culture inserts allowed easy handling and the indirect interaction of the formed VBCTs (within the insert) with the OBs (on the well plate).
We show the formation of a vascularized breast cancer model using MDA-MB231 and the influence of OBs on the VBCT morphology, vasculature and gene expression. Furthermore, the effect of blood vessels (BVs) on the model in presence of OBs was elucidated. We observed that the morphology of MDA-MB231 changes significantly in 3-D scaffold-free tissues as compared to 2-D cultures. Following this, we investigated the influence of OBs and MDA-MB231 on the tissue vasculature. Finally, utilizing gene microarray analysis, the indirect paracrine effect of OBs on the VBCTs was studied. Since vascular re-modeling is one of the first steps in the initiation of a pre-metastatic niche, we also studied the influence of BVs in the breast cancer tissues co-cultured with OBs. Important hub genes and highly interconnected protein-protein interaction (PPI) clusters were identified to understand the influence of OBs and BVs in the VBCTs.
In short, we fabricated a simple and realistic biological system that was utilized to understand the influence of the secondary bone site on the primary breast cancer microenvironment. This system, in the future, could be used to dissect the role of secondary tumors that drives breast cancer organotropism.

Formation of 3-D vascularized breast cancer tissues.
To fabricate the scaffold-free breast cancer model, we utilized the layer-by-layer ECM coating and accumulation technique. Normal fibroblasts (NFs) were coated with fibronectin (FN) and gelatin (G) proteins to form nano-layer of ECM moieties around the cells 28 . The coated individual NFs were mixed with HUVECs and MDA-MB231 BCCs followed by accumulation in a cell-culture insert. We observed that MDA-MB231 morphology was altered in 3-D scaffold-free tissues as compared to 2-D culture (Fig. 1a,b,e). We observed that the BCCs interacted with the BVs either by extending along the vascular length or by remaining in a circular morphology in direct contact with the BVs (Fig. 1c,d). Majority of the BCCs showed an extended morphology and elongated filopodia in the vascularized 3-D culture (Fig. 1b,c,e). Similar filopodia structures in MDA-MB231 BCCs were observed previously in a 3-D collagen matrix 29 . We quantified the BCCs aspect ratio and circularity and noted that the BCCs morphology in 3-D ranged from low (circular) to high (elongated) while BCCs in 2-D flat surface did not demonstrate elongated morphology (Fig. 1f,g).
Recently, it was shown that MDA-MB231 morphology fluctuated depending on the ECM composition of the substrate 30 . The presence of FN led to an increase in the formation of tunneling nanotubes, whereas the presence of matrigel led to highly elongated cell morphology 30 . Furthermore, MDA-MB231 morphology has been shown to change drastically in collagen substrates of varying concentrations suggesting that the morphology is highly dependent on its surrounding microenvironment 30 .
Since our model is comprised of cell-secreted ECM rather than preferentially pre-selected ECM hydrogels with designated concentration, we believe that the scaffold-free models depict an unbiased tumor-microenvironment and therefore predict better the in-vivo cellular morphology.
In order to validate the scaffold-free in-vitro model for a broad utilization in cancer research, we screened for additional breast cancer subtypes and observed their morphology within the 3-D microenvironment (Fig. 2a-f). Presence of MDA-MB231 and OBs modulates vasculature. Firstly, to confirm that the OBs seeded on the well plate ( Fig. S1f) can influence the BCCs within the cell-culture insert, we conducted a 2-D migration assay to observe the migration of BCCs due to OBs. We noted an enhanced migration of MDA-MB231 cells in the presence of OBs or 100 ng/ml C-X-C-chemokine ligand 12 (CXCL12) ( Fig. S1a-g). Previous studies have shown that the positive-feedback loop initiated by the C-X-C chemokine receptor type 4 (CXCR4), highly expressed in malignant BCCs, to the CXCL12 expressed in bone lead to breast cancer migration 33,34 . We noted that the pre-treatment of BCCs with CXCR4 antagonist AMD3100 (25 μg/ml) restricted this migration of BCCs towards both OBs and CXCL12 (Fig. S1d,e,g). Therefore, we confirmed that cell culture inserts could be utilized for indirect co-culture assays where the secreted factors of one cell could influence the other cell type. Following this, we fabricated the VBCTs within the insert and seeded OBs on the bottom of the well plate to study the influence of OBs on vascularization of the tissue. A schematic illustration of the setup is shown in Fig. 3a. For testing the effect of BCCs and OBs on vessel fate, four conditions were tested namely vascularized non-cancerous tissue (NH), vascularized tissue with MDA-MB231 (NHM), vascularized non-cancerous tissue indirectly co-cultured with OBs (NHO) and, vascularized tissues with MDA-MB231 indirectly co-cultured with OBs (NHMO). For clarity, the conditions and their cell constituents are summarized in Table 1.
We observed that the NH and NHO tissues formed an organized vascular network, whereas NHM and NHMO models showed abnormal vascularization with disorganized architecture (Fig. 3b). A pictorial representation is shown in Fig. 3c. We noted that the vessel area percentage and junction density was significantly lower in NHMO samples as compared to NH and NHO samples (Fig. 3d,f). Similarly, the average vessel length was significantly lower in NHMO as compared to NH and NHO (Fig. 3e). The average vessel length was also lower in NHM as compared to NH (Fig. 3e).
These readings indicate that in a scaffold-free system, the inclusion of BCCs leads to detrimental alterations in the surrounding vasculature and this process is accentuated when the cancerous tissue is exposed indirectly to OBs.
To further understand the influence of OBs and BCCs derived soluble factors on angiogenesis, we conducted the VEGF-ELISA on day 7 of culture time. Previous reports have indicated that the VEGF levels in the (d) Analysis shows a significant decrease in vessel area percentage in NHMO samples as compared to NH and NHO samples (n = 3, with 3 different tissue area per sample). (e) Significantly decreased average vessel length (mm) in NHMO as compared to NH and NHO, and significantly decreased vessel length in NHM as compared to NH, (n = 3, with 3 different tissue area per sample), (f) significantly decreased junction density in NHMO as compared to NH and NHO, (n = 3, with 3 different tissue area per sample), (g) VEGF analysis shows a significant difference in VEGF levels between NH and NHM, NHO and NHMO, and NH and NHMO, VEGF level in OBs monolayer alone was statistically different to NHM, NHO and NHMO but not significantly different to NH (p-values not shown for graph simplicity), (n = 3). Values in all the graphs are mean ± standard deviation; analysis was done using ordinary one-way Anova. www.nature.com/scientificreports/ intra-tumoral region of TNBC tissues were significantly higher than non-TNBC tissues 35 . Elevated VEGF level in breast cancer patients has also been linked to diminished progression-free survival and overall survival of the patient 36 . OBs have also been shown to be a source of VEGF, particularly during bone repair 37 . We observed the highest VEGF secretion in NHMO condition followed by NHM, NHO, NH and OBs respectively. VEGF levels from NHMO were significantly different to those from NHM, NH, NHO, NH and OBs (Fig. 3g). VEGF levels from NHM were significantly different to NH and OBs while NHO was significantly different to OBs (Fig. 3g).
We can therefore conclude that the presence of MDA-MB231 in the in-vitro tissue led to increased production of VEGF, and the indirect co-culture with OBs elevates the VEGF levels.
In the present cancer model, the vessel area percentage or junctional density in the vascularized tissue did not correlate with the increase in VEGF levels. Instead, we observed segmented, constricted and abnormal BVs in the presence of MDA-MB231 that led to the decrease in vascular density.
The current cancer model comprised BCCs distributed randomly within the 3-D tissue. To better deliberate the interaction of BCCs and the vasculature, we fabricated a separate model where cancer cells were seeded on top of the vascularized tissue construct instead of being randomly distributed (Fig. S3a). As the MDA-MB231 cells were grouped in a single layer, it was expected that the direct interaction of BCCs with the vascular bed beneath would better demonstrate the fate of vessels. At the intersection between cancer cells and vasculature, we observed bright PECAM-1 + spots that indicated broken or open vascular ends (Fig. S3c). Magnified images of the spots showed that the open vascular ends have thin irregular sprouts protruding out of them (Fig. S3b,d).
In order to further deliberate the paracrine interaction between the VBCTs and OBs, we conducted Interleukin-6 (IL-6) and CC-motif chemokine family (CCL2) ELISA after day 7 of culture time (Fig. S2a,b). CCL2, also known as MCP-1 (Monocyte Chemoattractant protein-1) is observed in prior studies to be overexpressed in tumor and the surrounding stroma 38 . CCL2 is known to have an influence in bone remodelling and angiogenesis 39,40 . We observed that the CCL2 levels were not significantly different between NHM and NHMO (Fig. S2a). Furthermore, there was no significant difference between NH and NHM, or NHM and NHO. It is to be noted that the presence of BCCs in the tissue slightly increases (non-significant) the soluble CCL2 levels as compared to conditions without BCCs. We could however observe a significant difference between NHMO and NH, and NHMO and NHO indicating that the presence of both MDA-MB231 and OBs in the in-vitro system leads to elevated levels of CCL-2 (Fig. S2a).
In prior research, IL-6 has shown to have clinical significance in breast-to-bone metastasis, and anti-IL6 therapies has been investigated 41 . Furthermore, IL-6 also influences the vasculature by regulating VEGF and disrupting pericytes coverage of blood vessels in the tumor microenvironment 42 . We observed a significant increase in IL-6 levels in NHMO condition as compared to NH, NHO, and NHM. This observation indicates that the presence of both OBs and MDA-MB231 lead to enhanced levels of IL-6. Together, we show via VEGF, CCL2 and IL-6 level analysis, the paracrine interaction between different cell types in the constructed in-vitro model system.

OBs modulates gene expression in VBCTs.
Our next focus was to elucidate the influence of OBs on the tumor tissue. VBCTs were cultured indirectly with OBs for 7 days and subsequently harvested. Gene microarray analysis was conducted to identify differentially expressed genes (DEGs) and crucial hub genes. We analyzed the DEGs between NHM and NHMO samples to elucidate the influence of OBs on the VBCTs. Upregulated genes (Fold change, FC ≥ 2, p < 0.05) and downregulated genes (FC ≤ − 2, p < 0.05) were selected as the significant DEGs.
We observed 1476 significantly upregulated genes and 775 significantly downregulated genes in models that were co-cultured with OBs. The volcano plot depicts the upregulated and downregulated genes (Fig. 4a). The top five upregulated DEGs and the top five downregulated DEGs are shown in Table 2. Each of the highly upregulated genes was loaded into TNM plotter (BRCA database) to observe the differences in their expression in normal, cancerous and metastatic tissues (Fig. S5a-e) 43 . We observed that all the top five upregulated genes in NHMO samples as compared to NHM samples were expressed more in cancerous and metastatic breast cancer tissues as compared to normal breast tissues (Fig. S5a-e).
The gene signature consisting of all the five upregulated genes was observed to be expressed more in both TNBC and broad BRCA tissue database (fed into GEPIA2 application) as compared to normal tissue database (Fig. S5f). The functional enrichment analysis was conducted on the upregulated and downregulated DEGs separately (Fig. 4b,c). Five top upregulated pathways were RNA metabolism, the establishment of protein localization to organelle, cellular response to stress, positive regulation of hydrolase activity and vesicle mediated transport (Fig. 4b). Five top downregulated pathways were sensory perception of smell, regulation of body fluid levels, GPCRs-Class A Rhodopsin like anion transport and cardiac conduction (Fig. 4c).
Next, the four highly interconnected hub genes from the upregulated PPI network were identified namely HSP90N, CYCS, EGFR, and RPS27A. The expression of the four hub genes was tested in tumor and normal tissue database (fed into GEPIA2 application) specific for TNBC patients (Fig. 5a-d). The fold change and p-values of the upregulated hub genes are provided in Table S1.The result from the TNBC database indicated that HSP90N and CYCS are significantly present in the tumor tissue as compared to normal breast tissue (Fig. 5a,c), and EGFR is expressed significantly more in normal breast tissue as compared to breast cancer tissues (Fig. 5b). Disease-free survival plots of the hub genes specific for TNBC patients are shown (Fig. 5e-h). Out of the four hub genes, the survival plots indicated that HSP90N had a significantly unfavorable impact on the survival of the patients (Fig. 5g).
Next, four clusters of highly interconnected upregulated genes were determined (Fig. S6). Reactome pathway enrichment analysis of each cluster was noted (Fig. 5i). The pathways included terms such as antigen processing, neddylation, metabolism of RNA, collagen biosynthesis and modifying enzymes, and ECM organization.  44 . The comparison between cancerous tissues with vessels, and cancerous tissues without vessels could elucidate the influence of BVs in cancer progression. Therefore, we compared the gene expression between day 7 NHMO and NMO samples. We observed 1278 upregulated genes and 566 downregulated genes in the presence of BVs. The volcano plot depicting the upregulated and downregulated genes is shown in Fig. 6a. Table 3 shows the list of the top five most highly upregulated and downregulated genes.
The upregulation of endothelial-related genes was expected as one tissue condition consisted of BVs and the other condition was non-vascularized. The ECs-related genes such as ANGPT2, CD34, LYVE1, VWF and PECAM1 were highly upregulated in the vascularized tumor tissue as compared to non-vascularized tissue due to the presence of BVs (Table 3).
Interestingly, several of the downregulated genes in NHMO vs NMO may have implications in OBs fate and breast cancer progression. The five highly downregulated genes were PCSK5, PTGDS, CXCL14, ADGRD1, and SMOC2 (Table 3). SMOC2, a calcium-binding protein is known to inhibit osteogenesis and prevent mineralization in OBs 45 . ADGRD1 also known as GPR133 could also have implications in bone diseases 46 . CXCL14 has shown to be expressed in the stromal cells in malignant breast cancer 47 and has been shown to induce lung cancer metastasis to bone 48 .
Next, the top functional enrichment pathways of the upregulated genes were determined namely blood vessel development, epithelial cell differentiation, regulation of cell adhesion, chemotaxis, regulation of GTPase activity,  www.nature.com/scientificreports/ leukocyte migration, positive regulation of cell migration and pathways in cancer (Fig. 6b). The genes involved in these upregulated pathways are shown in Table S3. Two highly interconnected gene clusters are determined from the upregulated PPI network in Fig. S7. The downregulated pathways were cell surface receptor signaling pathway, reproductive structure development, regulation of ossification, heart field specification and response to BMP (Fig. 6c). The genes involved in these downregulated pathways are shown in Table S4. Four hub genes from the upregulated PPI network were identified namely KDR, ITGAM, CXCR4, and PECAM1. Since the upregulated hub genes were mostly associated with endothelial cells, we speculated that the downregulated hub genes would provide more information about the influence of vessels in the VBCTs. The four downregulated hub genes were determined namely BMP4, C3 protein, CXCL13, and DKK1. BMP4 and C3 were significantly expressed in normal tissue as compared to cancerous breast cancer tissues (Fig. 6d,e) whereas CXCL13 and DKK1 is observed to be expressed in breast cancer tissue as compared to normal breast tissue (Fig. 6f,g).

Discussion
The theory that there is already a prior interaction between primary tumor and secondary site before metastasis is not new 49 . For instance, previous reports have shown that lung tumors can remotely activate the OBs, and the OBs, in turn, supply the lung tumors with tumor-promoting neutrophils even before the inception of metastasis 50 .  www.nature.com/scientificreports/ Another in-vivo study showed that OBs-derived hypoxia-inducible factor led to distant breast cancer growth and dissemination, partly due to the enhanced production of CXCL12 19 . The study also showed that osteoclasts did not possess a similar ability to regulate distant BCCs, which highlights the importance of OBs in influencing primary tumor site 19 . A recent article showed that colorectal cancers release extracellular vesicles that activate fibroblasts in distant organs 51 . The fibroblasts then produces inflammatory cytokines that promote primary tumor growth 51 . Therefore, it is evident that the pre-metastatic indirect interaction between secondary and primary sites could illuminate several hidden yet important pathways that initiate the metastasis process. In human bone, OBs account for 4-6% of the total bone cells and deposit bone matrix 52 . OBs maintain a homeostatic bone environment together with osteoclasts, the bone-degrading cells and resident osteocytes 52 . Breast cancer metastasis to the bone disrupts this homeostasis and leads either to enhanced bone deposition (osteoblastic) or enhanced bone resorption (osteolytic). Several studies in the past have already elucidated the role of BCCs in the fate of OBs 53,54 . These studies illuminate how the BCCs remotely prime the secondary bone sites before or after homing. However, the role of the secondary bone site on the primary tumor prior to metastasis hasn't been studied in its entirety.
Cancer cells not only possess the intrinsic ability to proliferate and metastasize but also rely on their surrounding tumor-microenvironment for their progression. The TME is a crowded milieu comprising of multiple cell types such as fibroblasts, endothelial cells, and ECM moieties that constantly interact with one another to modulate tumor fate. The scaffold-free VBCTs that we fabricated are composed of the cancer cells at an interface of vasculature and connective tissue and provide a realistic alternative to conventional 3-D models. Our results indicate that the presence of MDA-MB231 in VBCT resulted in sectioned and abnormal vessels and this abnormality was accentuated in the presence of OBs.
Solid tumors require vessels for growth and metastasis and therefore enforce the production of more BVs. However, cancer cells not only rely on angiogenesis to initiate the metastatic process but can also hijack preexisting vessels in the tumor microenvironment in a non-angiogenic process called vessel co-option 55,56 . Cancer cells, during rapid proliferation in the tissue, also compress the surrounding vasculature in the microenvironment thereby leading to vessel collapse 57 . Recently, it was shown in an in-vitro 3-D model that cancer cells' interaction with the vasculature led to compression and destruction of the vessel 58 . In another study, it was shown that pancreatic ductal carcinoma cell (PDAC) ablates endothelial cells both in in-vitro and in-vivo settings that lead to hypo-vascularization in the PDAC tumor tissue 59 .
In a confined tumor-microenvironment, cells compete for space and nutrients. Cancer cells are highly proliferative in comparison to endothelial cells, which could lead to the elimination of the 'weaker' cell population (endothelial cells) by the 'fitter' cell population (MDA-MB231) 59 . The scaffold-free model fabricated in a confined cell-culture insert area includes three different cell types in large numbers, each competing for the available nutrients and space, which could explain the elimination of endothelial cells. Furthermore, it could be hypothesized that the presence of OBs and their secreted factors such as CXCL12 could lead to enhanced cellular proliferation and therefore further lead to reduced availability of the limited nutrient to individual cells. We observed via cell proliferation XTT assay that CXCL12 leads to a slight increase in MDA-MB231 and HUVECs proliferation, and a significant increase in proliferation of fibroblasts (Fig. S1h).
Several other factors could affect angiogenesis and vascular fate. One such factor is the progression stage of the tumor. For instance, it has been shown previously in gliomas that the vessels undergo three stages of development: (a) tumor-vessel interaction phase (b) transition phase where vessel apoptosis takes place followed by (c) angiogenic phase 60 . The increase in VEGF and the presence of thin vascular sprouts in VBCT suggest a transition phase preceding angiogenesis. For these phases to occur, tumor cells have to be in close proximity to the vessels since these phases require direct cell-cell interactions. The advantage of our model is the ability of cancer cells to directly interact with BVs (Fig. S4). Lastly, the influence of the highly proliferating TNBC on vasculogenesis (apart from angiogenesis) should also be taken into consideration since MDA-MB231 cells are mixed together with endothelial cells and not introduced to fully-formed vessels. Future experiments where MDA-MB231 cells are introduced after the formation of BVs in the model could elucidate different angiongenic mechanisms.
We show a number of interesting upregulated hub genes from the gene array data between NHMO and NHM samples where we studied the effect of OBs on VBCTs. We observed HSP90N as an upregulated hub gene that has been previously linked to aggressiveness in breast tumors 61,62 . HSP90N inhibitors have been used in the past for cancer therapy and could also have an implication in bone metabolism 63 . Interestingly, HSP90N expression is also associated with enhanced bone metastasis from prostate and breast tumors 64,65 .
The other hub genes EGFR, RPS27A and CYCS have also been previously identified as hub genes in TNBC extracted from several gene expression datasets 66 . The overexpression of EGFR, another upregulated hub gene, leads to poor prognosis in breast cancer patients 67,68 . EGFR results in cancer cell proliferation, aggressiveness and metastasis 69 . Interestingly, expression of EGFR has also be linked to increased metastasis of cancer cells to bone 69 . For instance, previous results have shown that EGFR signaling can lead to osteoclastogenesis by influencing OPG and MCP1 signals in osteoblasts 70 . However, when the GEPIA2 database was utilized for gene analysis a lower expression of EGFR in malignant TNBC samples as compared to normal human breast tissue was seen. This inconsistency of the expression pattern of EGFR has also been discussed previously in reference to ovarian cancer 71 where it is pointed that the inconsistencies might arise as the datasets such as GEPIA or ONCOMINE comprise of solely mRNA data, which correlates to 40% of the total protein change. We agree with the authors, however, we also believe that the expression of EGFR in breast cancer is still not deliberated in its entirety. For instance, feeding EGFR protein expression into the Human protein atlas database which also accounts for protein expression indicates that EGFR is (a) not prognostic in breast cancer (b) is less frequently detected in the immunohistochemistry samples of breast cancer patients 72 . New research have now indicated that a combined screening of HER3-EGFR is more suitable as it has been shown that one protein may downregulate and be compensated by the upregulation of the other protein. This could explain why EGFR as a prognostic biomarker is still uncertain 73  The changes in gene expression between NHMO and NMO were studied to deliberate the impact of BVs. We hypothesized the presence or absence of BVs would lead to modified gene expression. DKK1, an antagonist of the Wnt/β-catenin signaling pathway was detected as a downregulated hub gene. Previously, it has been shown that DKK1 functions as a tumor suppressor 74 . Interestingly, upregulation of DKK1 has shown to be involved in specific breast cancer metastasis to bone 75,76 . Furthermore, tumor-derived DKK1 also dictates the OBs fate by inhibiting OBs differentiation and facilitates osteolysis 77 . It is intriguing that in the presence of BVs in the cancerous tissue, DKK1 was observed as a downregulated hub gene. In malignant gliomas, it was previously reported that hypoxia led to an upregulation of DKK1 as compared to normoxia 78 . It could be theorized that in the absence of BVs and their secreted growth factors in the non-vascularized tissue, DKK1 was upregulated which in turn could result in the enhanced preference of breast cancer metastasis to bone. A more thorough future investigation into the role of BVs in DKK1 expression could highlight many possible mechanisms in bone metastasis. Another downregulated hub gene BMP4 has previously shown to inhibit proliferation but accentuate cell migration in breast cancer 79 . BMP4 can act as both tumor suppressor and promoter 80 . One study showed that BMP4 upregulation is associated with reduced metastasis and the downregulation of BMP4 has been associated with more aggressive breast cancer 81 . Another study showed that mice treated with BMP4 showed a positive correlation with bone metastasis 82 . Finally, downregulated hub genes complement protein C3 and CXCL13 both could have implications in breast cancer progression and modulation in bone microenvironment 83,84 .
Together, the developed in-vitro model is suitable to study several pathways involved in the progression of metastasis and organotropism. In this study, we show that indirect culture with OBs leads to a change in the gene expression profile of VBCTs. This emphasizes the point that reciprocal OBs secreted factors could be an important parameter in early metastasis events in breast cancer. In the future we will be investigating deeper into the protein expression pattern of the VBCT not only due to the interaction with osteoblasts but also with the interaction with other bone cells such as osteoclasts to elucidate the role of different cell types in cancer progression or dormancy.

Conclusion
In summary, we present a 3-D co-culture system utilized to understand the interaction between several cell types involved in breast cancer metastasis to the bone. The fabricated VBCTs provide a realistic microenvironment to the TNBC cells with the presence of fully formed vasculature, surrounding fibroblasts and cell-secreted ECM. The simplicity of the culturing and harvesting of the 3-D model allowed a gene-level study of the influence of OBs and BVs on the engineered cancer tissue. These results provide a proof-of-concept for the utilization of such models for the understanding of time-dependent metastatic progression. Such models in the future can be used in predicting essential biomarkers for diagnosis and for pre-clinical testing of novel therapeutic agents.

Methods and materials
Cell culture and media condition. GFP + MDA-MB231 cells were kindly provided by Dr. Julia Steitz, Institute of Laboratory Animal Science, RWTH University Hospital, Aachen, Germany. The cells were cultured in DMEM (11965092, Thermofisher Scientific) with 10% fetal bovine serum (FBS). Osteoblasts were procured from Sigma Aldrich (406-05F), and cultured in Osteoblast growth media (417-500, Sigma Aldrich). Fibroblasts were obtained by isolation of skin biopsies that were surgically extracted from healthy volunteers (Department of Dermatology, RWTH Aachen University Hospital, Germany). The extraction and isolation were conducted in accordance with the Declaration of Helsinki principles and was approved by the ethical committee of RWTH Aachen University Hospital, Germany, Project Number EK188/4. Fibroblasts were cultured in DMEM with 10% FBS. HUVECs (Lonza, Passage < 4) were cultured in an endothelial growth medium (C-22111, Promocell). For the culture of VBCTs, DMEM: EGM (1:1) with 10% fetal bovine serum (FBS, Biowest) and 50 µg/ml Ascorbic acid (A92902, Sigma) was used.
Migration assay. OBs (passage < 5) were cultured in 24 well plates (7 × 10 4 cells) and allowed to cultivate for a day. MDA-MB231 cells were cultivated in serum-free DMEM media for 24-h. Required number of cells was treated with 25 μg/ml AMB3100 (EMD Millipore) for 1 h. Cell culture inserts (8 μm pore size, Corning) were pre-treated with DMEM medium for 2 h prior to cell addition. 1 × 10 4 treated or untreated MDA-MB231 cells were added to the upper compartment of the 8 μm pore size cell culture inserts in 250 µl of serum-free DMEM (with 0.1%BSA) whereas 600 μl of serum-free DMEM (with 0.1% BSA in media) was added to the lower compartment. After 1 h of attachment time, inserts were transferred into the well plate containing OBs and fresh media was added. CXCL12 (Sigma, SRP3276) in serum-free DMEM media was added to appropriate chambers depending on the experimental requirement. After 18 h of the assay, the lower compartment was stained with DAPI and images were taken in a fluorescent microscope (10× magnification). Three separate regions were selected per insert. Experiments were done in triplicates.
XTT assay. 5 × 10 4 fibroblasts, HUVECs or MDA-MB231 cells were seeded on 96-well plates. In each cell conditions, 0 ng/ml (vehicle: DMSO), 50 ng/ml and 100 ng/ml of CXCL12 were added for 24 h in respective serum-free media (0.1% BSA). XTT (11465015001, Sigma) was utilized as per the manufacturer's protocol. The specific absorbance at 460 nm was observed and the unspecific absorbance at 630 nm was subtracted.
ECM-coating and accumulation to fabricate cancer models. Fibroblasts were coated layer-by-layer using ECM coating and accumulation technique as previously described 26 . Briefly, fibroblasts were exposed to 0.04 mg/ml FN (sc-29011A, SCBT) and 0.04 mg/ml G (G9391, Sigma) solutions with an intermediary PBS washing step. After 9 steps of alternate exposure to FN and G, the coated cells were mixed with HUVECs in 15 www.nature.com/scientificreports/ ratios and 2.5 × 10 4 GFP + MDA-MB231 were mixed with the HUVECs and coated fibroblasts. The models were cultured for 24 h in 0.4 μm pore sized inserts (3470, Corning) and the inserts were transferred to well plates containing OBs. Prior to the fabrication of models, 7 × 10 4 OBs were seeded on 24-well plates and cultured for 2 days. To fabricate 3-D models comprising of MCF-7 and MDA-MB468, we mixed coated fibroblasts with HUVECs (15:2 ratio) and 2.5 × 10 4 MCF-7 or MDA-MB468 depending on the required model. Both the breast cancer cell lines were first live stained with CellTracker™ Red CMTPX Dye (C34552, Invitrogen) for 20 min according to the manufacturer's protocol and then subsequently mixed with fibroblasts and HUVECs. The models were cultured for 72 h with media changes every day.
Immunostaining of tissues. Tissues were extracted out of the cell culture insert, washed twice with PBS.
The harvested tissues were fixed in 4% PFA for 30 min at room temperature (RT). The tissues were washed and suspended in 0.1% Triton-X 100 for 10 min. The unspecific antibody binding was blocked by suspending the tissue in 3% BSA (Sigma, Germany) in PBS for 1.5 h in RT. After PBS wash (2×), mouse anti-human CD-31 primary antibody (MA5-13188, Thermofischer Scientific, 1:100 dilution in PBS) was added to the tissue and incubated overnight in 4 °C. After 24 h, the tissues were washed with PBS (3×, 5 min interval in the shaker) and goat anti-mouse Alexa 594 (A11005, Invitrogen) 1:200 dilutions was added for 3 h in RT. Tissues were washed with PBS (3×, 5 min interval in shaker). The tissues or 2-D MDA-MB231 were further stained with DAPI (1:1000, 5 min incubation in RT) if required. Finally, the samples were prepared for confocal microscopy (Leica TCS SP8).
Microscopy and image analysis. Stained VBCTs were extracted from the cell culture inserts, fixed on a glass slide and investigated under the confocal microscope. Samples were observed for vessels (CD31) using an excitation wavelength of 594 nm (red), GFP + MDA-MB231 using an excitation wavelength of 488 nm (green) and DAPI (nucleus) with an excitation wavelength of 405 nm. Images were taken at a resolution of 1024 × 1024 pixels using either 10× dry objective, 20× Oil immersion objectives or 63× Oil immersion objectives depending on the experimental requirements. Vessels were analyzed for vascular density, area percentage and junction point using AngioTool as per software instruction 85 . All illustrations in the paper were made using Biorender. com.

Gene chip microarray.
To generate amplified sense-strand cDNA, 300 µg purified mRNA of each model were processed using the WT Expression Kit (Ambion, Austin, TX, USA). Affymetrix® GeneChip ® WT Terminal Labeling Kit Assay (Affymetrix, Inc., Santa Clara, CA, USA) was used for fragmentation and labeling. Prepared hybridization cocktails were applied to Clariom™ S assay (Thermo Fisher Scientific). The assays were washed and stained using the GeneChip ® Fluidics Station 450 (Fluidics Protocol FS450_0001). Arrays were scanned using Affymetrix GeneChip ® Scanner 3000 7G controlled by GeneChip ® Operating Software (GCOS) version 1.4 to produce CEL intensity files. www.nature.com/scientificreports/ (DFS) was analyzed for the hub genes using median cut-off (cut-off high 50%, cut-off low 50%) in GEPIA2 application 92 . DFS were plotted specifically by utilizing a subtype filter that only selects TNBC patients for the analysis. We omitted subtypes HER2 + luminal, luminal A and luminal B for the construction of the survival plots.
Statistics. Statistical analyses were performed using GraphPad Prism 9.0.0 software. All experiments were performed in triplicates. Data are presented as mean ± standard deviation unless otherwise stated. Data are analyzed using one-way ANOVA with Tukey's multiple comparison test. Significance was considered when p < 0.05. Individual p-values are stated in the graph.