Impact of extracellular matrix properties on neuroblastoma genomic heterogeneity

Increased tissue stiffness is a common feature of malignant solid tumors, often associated with metastasis and poor patient outcomes. Vitronectin, as an extracellular matrix anchorage glycoprotein related to a stiff matrix, is present in a particularly increased quantity and specic distribution in high-risk neuroblastoma. Furthermore, as cells can sense and transform the proprieties of the extracellular matrix into chemical signals through mechanotransduction, genotypic changes related to stiffness are possible.


Abstract Background
Increased tissue stiffness is a common feature of malignant solid tumors, often associated with metastasis and poor patient outcomes. Vitronectin, as an extracellular matrix anchorage glycoprotein related to a stiff matrix, is present in a particularly increased quantity and speci c distribution in high-risk neuroblastoma. Furthermore, as cells can sense and transform the proprieties of the extracellular matrix into chemical signals through mechanotransduction, genotypic changes related to stiffness are possible.

Methods
We have applied high density SNPa and NGS techniques to in vivo and in vitro models (orthotropic xenograft vitronectin knock-out mice and 3D bioprinted hydrogels with different stiffness) using two representative neuroblastoma cell lines (the MYCN ampli ed SK-N-BE(2) and the ALK mutated SH-SY5Y), to discern how tumor genomics patterns and clonal heterogeneity of both cell lines are affected.

Results
We describe a remarkable subclonal selection of some genomic aberrations in SK-N-BE(2) cells grown in knock-out vitronectin xenograft mice that also emerged when cultured for long times in stiff hydrogels.
Specially, we detected an enlarged subclonal cell population with chromosome 9 aberrations in both models. Similar abnormalities were found in human high-risk neuroblastoma with MYCN ampli cation. Genomics of the SH-SY5Y cell line remained stable when cultured in both models.

Conclusions
Focus on heterogeneous intratumor segmental chromosome aberrations and mutations, as a mirror image of tumor microenvironment, is a vital area of future research.

Background
Tumors are not homogeneous structures, but rather very complex tissues involving many cell types, such as tumor, stromal and in ammatory cells [1][2][3]. Cells are primarily supported in solid tumors by the extracellular matrix (ECM), a three-dimensional (3D) dynamic network composed mainly of brous proteins, glycoproteins and proteoglycans [1,4]. In tumors, cells respond to the biochemical signals of their ECM, but also to physical forces such as tension (traction and compression forces that maintain their stability), known as biotensegrity [2]. In fact, increased stromal stiffness is a classic hallmark of cancer [1,5], as is transformation of this stiffness into chemical signals through mechanotransduction for cell advantage [5][6][7]. Data suggest that an aberrant ECM may promote genetic instability and can even compromise DNA repair pathways necessary to prevent malignant transformation [6,[8][9][10]. Among the cell population trying to respond to ECM stiffness, only the most genotypically and phenotypically adapted ones will succeed in invading new metastatic niches [11,12]. The most advantageous cells will also segregate anchor proteins and transform the ECM into a more biochemically and tensegrally favorable microenvironment, which will also promote growth and migration [6,13]. A continuous interplay therefore exists between tumor ECM and cell genomics and behavior [13].
ECM stiffness and composition are currently of such importance that, together with mechanotransduction, are considered promising therapeutic targets in many cancer types [14][15][16]. It is also the case of neuroblastoma (NB), a peripheral sympathetic tumor that causes 15% of childhood deaths from cancer [17]. We previously de ned an aggressive pattern of rigid ECMs in high-risk NB (HR-NB) [18]. NB stiff ECM is rich in cross-linked collagen III bers, poor in glycosaminoglycans, supporting sinusoidal vascular structures (blood and lymphatic) and with a high quantity of territorial vitronectin (VN, located in the cytoplasmic compartment and in a thin layer around the tumor cells) [19][20][21][22]. NB stiff ECM is also rich in collagen I bers and bronectin but in less amount and with a more physiological-like topology than collagen III and VN [18,23] Moreover, VN is a glycoprotein containing multiple cell receptor binding sites including integrins, uPAR and PAI-1 [24], and which can therefore form transitory ECM elements-cell junctions. It seems to anchorage to ECM bers and proteoglycans, as well as disrupted cell adhesion, promoting spreading and metastasis in NB [22,25]. Some studies have shown an important role for VN in initiation and progression of hepatic [26], breast [27] and lung [28] cancers, and its particular involvement in migration in ovarian cancer [29,30]. In a previous report we used inoculation of the NB cell lines studied here into the adrenal gland in VN knockout mice (VN-KO) to study the systemic role of VN in tumor growth. Interestingly, we found VN-dense cytoplasmic expression in tumor cells from both VN-KO and wild type (VN-WT) mice, with no differences in VN staining pattern or in tumor growth rate in early passages [22]. However, in this study we sought to resolve the particular question of whether host lack of VN has any impact on tumor genomics.
Besides MYCN ampli cation (MNA), segmental chromosomal aberrations (SCAs), ALK mutation and ampli cation, TERT rearrangement and ATRX mutation and deletion are also genetic characteristics of HR-NB [31,32]. In HR-NB patients these genomic aberrations usually show high intratumor heterogeneity, resulting in resistance to treatment and metastasis [33][34][35]. We have reported a link between a SCA (1p chromosome arm deletion, 1p-) and reduction of glycosaminoglycans in the ECM of high-risk NB (HR-NB) [36]. Nevertheless, the association between other ECM features and HR-NB genomic aberrations, and clonal selection, is a neglected area of study.
3D cultures are emerging as good models for ECM-related biological studies, also reducing use of animal models. These scaffolds can reproduce complex ECM features absent from traditional, 2D or monolayer cultures, and allow us to add degrees of complexity to delve deeper into the physiopathological tumor microenvironment [37,38]. We recently published a study using 3D bioprinted hydrogels composed of methacrylated gelatin (GelMA) and increasing concentrations of methacrylated alginate (AlgMA) to recreate different ECM stiffness. We described a pattern of aggressive behavior in NB cells cultivated for long times in rigid hydrogels [39]; however, this study did not analyze how scaffolding stiffness and growing time affected genomic heterogeneity.
In this work the described models (xenograft VN-KO mice and 3D bioprinted hydrogels) are used to discern how tumor genomics patterns and clonal heterogeneity in two NB cell lines are affected by the VN host status, and tumor microenvironment stiffness and also to extrapolate the results in primary human HR-NB. This study shows that tumor genomics and clonal genetic heterogeneity are directly related to tumor surroundings.

In vivo models
In vivo mice models were obtained as previously reported [22]. Four-to six-week-old genotyped female or
Genomic analysis by High Density Single Nucleotide Polymorphism arrays (HD-SNPa) DNA was extracted from 2D cultures, xenografted tumors and 3D bioprinted hydrogels using the salting out method for either frozen and para n-embedded tissue, as previously described [33]. Liquid biopsy study (circulating tumor DNA, ctDNA) was carried out on plasma samples with QIAamp Circulating Nucleic Acid (Qiagen) following the provided protocol. CytoScan HD was used for DNA from frozen tumors and OncoScan CNV gene chips (Affymetrix, Santa Clara, CA, USA) for DNA from para nembedded tissue and plasma. DNA ampli cation, tagging, and hybridization to gene chips were performed according to the manufacturer's protocol. Data were analyzed using Chromosome Analysis Suite 3.2 (ChAS) software (Affymetrix, ThermoFisher Scienti c) and Nexus 10.0 Copy Number Discovery (BioDiscovery). We estimated the percentage of cells affected by each alteration (Fig. 1) with the average of the smooth signal of probes distributed throughout the whole alteration [41]. Log 2 ratio of each aberration and aberrant cell percentage of samples grown with SK-N-BE(2) cells was estimated by TuScan tool and summarized in Table S1. To analyze which genetic pathways were most affected by the detected SCAs, we used the PANTHER tool of Gene Ontology.

Mathematical analysis of detected genomic aberrations
All genomic aberrations detected with ChAS for SK-N-BE(2) cells were dichotomized by presence or absence in each sample. Multivariate analysis was performed using PAST software. We used the Jaccard coe cient to obtain neighbour-joining tree clustering and a matrix of similarity and distance indices.

New Generation Sequencing NB-mechanopanel
We employed the SeqCap EZ HyperPlus kit from Roche Diagnostics, and the NextSeq 550 sequencer (Illumina) to obtain at least 100x coverage. Next Generation Sequencing (NGS) raw data were ltered for quality with PrintSeq and mapped against human genome GRCh38 with BWA. Mutations of mice genome were removed with two different software tools: Xeno lteR and Disambiguate. Table S2 shows the NB-mechanopanel and mutations based on an allele frequency of at least 0.1 and 100x coverage.

Results
Chromosomal aberrations detected in 2D cell culture Genomic pro ling of trypsinized SK-N-BE(2) cells grown traditionally in monolayer culture revealed the presence of ten SCAs detected by HD-SNPa, as previously reported [42]: three SCAs typically observed in NB and seven SCAs classed as atypical for being present with lower frequency in NB. This cell line also showed a numeric chromosome alteration (NCA). As already described, SK-N-BE(2) presented genome ampli cation in 2p24.3 involving the MYCN, MYCNUT, MYCNOS and GACAT3 genes [42,43]. A chromothripsis-like phenomenon was also identi ed at chromosome 21 involving sixteen fragments.
HD-SNPa detected no differences in SCAs between SH-SY5Y cells in 2D cultures or when xenografted in VN-KO mice. Interestingly, in VN-WT mice without + 2p(ter-16.3, 49 Mb) the P1 tumor failed to grow in any of the eight tentative tumor xenograft attempts (Fig. S1).

SCAs of SK-N-BE
Hierarchical cluster analysis is a good tool for identifying sample similarities using genomic data Neighbor-joining tree clustering and an array of similarity and distance indexes based on the Jaccard coe cient (Fig. 3A) were obtained, according to presence or absence of the SCAs detected with ChAS software in each sample derived from SK-N-BE(2) cells (Fig. 3B). Mathematical analyses showed the highest similarity indexes between the genomics of 2D cultured cells and hydrogels with less stiffness and/or culture time. Separated from them were tumors from the control mice, with few similarities to stiff hydrogels and experimental tumors, especially with the lower indices in P3 to P5. In general, VN-WT tumors also had low similarity between their own passages.
Grouped together were the genomic pro les of cells derived from the stiffer and/or longer-cultivated hydrogels, and VN-KO tumor passages, allowing us to clearly observe the differences with the softer and/or shorter time 3D-cultures and with VN-WT tumors. Note that tumors of different VN-KO tumor passages had better similarity rates.

Biotensegrity had impact in single nucleotide mutations
We sequenced a small customized panel (NB-mechanopanel) of speci c genes, with mutations previously described in NB related to cytoskeleton remodeling and extracellular matrix changes, and other genes related to NB and cancer in general (Table S2). In SK-N-BE(2) cell line and in all samples derived from it, we identi ed two mutations described as p.C135F in TP53 [45] and p.N755K in ATRX [46], de ned as pathogenic in the COSMIC database. The p.I1590 = and p. SH-SY5Y cell lines grown in 2D, as well as in 3D hydrogels and mice, showed the pathogenic mutation described as p.F1174L in ALK [47], as well as p.G12V in KRAS not previously reported in this cell line [48]. We did not nd any new signi cant mutations (Table S2).

Chromosome 9 could be a hot spot for structural aberrations in HR-NB patients
We next sought to review the chromosome 9 pro les of homogeneous MNA primary tumor in HR-NB patients, due to the fact that the FSCAs and SCA of chromosome 9 were the main unusual genetic changes of the SK-N-BE(2) cells in the VN-KO tumor passages and in rigid and long-culture 3D-scaffolds. Surprisingly, we observed greater than expected presence of genomic aberrations on chromosome 9 in the 43 tumor pro les checked. Speci cally, 39.53% of cases had aberrations affecting the chromosomal positions of the DOCK8 and/or KANK1 genes (Fig. 4).

Discussion
HR-NB is a disease with risky spatial-temporal genetic heterogeneity and high chromosome instability related to poor prognosis and therapy resistance [49,50]. It is not yet clear whether this genomic variability could be the cause or consequence of selective evolutionary pressure in the tumor microenvironment [11]. We assumed that clonal genetic selection and/or evolution of aggressive NB cells could be promoted in part by the stiffness and composition of their surrounding matrix. To get closer to con rming this hypothesis, we selected two representative NB cell lines characterized by genomic alterations typically observed in HR-NB (MNA in SK-N-BE(2) [42] and ALK mutation in SH-SY5Y) [47]. Both cell lines were grown separately in hydrogels with ECM rigidity similar to the ECM of human HR-NB and xenograft tumors [39]. To make the 3D models even more biomimetic, malignant neuroblasts were cocultured with 10% Schwann cells in some cases, since even in poor stroma NB a low proportion can be found accompanying the undifferentiated or poorly differentiated neuroblasts [51]. To analyze the impact of ECM properties on genomic heterogeneity, we selected the HD-SNPa approach, complemented by the NGS technique, as one of the most commonly used procedures for detecting both MNA and SCAs in routine genetic diagnosis of NB [52]. Unlike this study, most research related to ECM tumor mechanotransduction has focused on uncovering changes in gene expression, which are di cult to translate to the clinic [53][54][55]. We looked for typical and/or atypical HR-NB SCA changes that could be associated with NB physics, and rapidly used for clinical translation to diagnosis and therapy. We also explored whether MNA and/or ALK-mutated cells respond similarly to ECM properties.
In previous studies, we demonstrated that elevated VN secretion by tumor cells was founded in rigid ECMs, and was related to poor patient outcomes and to growth of orthotopically inoculated NB cell lines [22]. Moreover, we observed that the high amount of VN is secreted forming tracks and we postulated that VN participates in stiffness, mediating biotensegrity and promoting tumor migration in HR-NB [22,56]. We also reported that SK-N-BE(2) cells grown in stiff 3D hydrogels show more e cient adaptation, increased cell proliferation, high expression of the Bcl2 antiapoptosis marker, and high mRNA processing rate (24). In the present study, genomic analysis of SK-N-BE(2) cell line revealed a number of equivalent aberrations in VN-KO tumors and in cells grown in stiffest and/or longer cultured times hydrogels. In these samples, we observed a positive selection of cells that in addition to MNA, contained: i) stable SCAs (1p-, 3p-and 17pq-); ii) SRO generated from previous SCAs (+ 2p and chromotripsis-like of chromosome 21) and iii) new SCAs and FSCAs (1p-, +1q, 9p-, +9p, 9q-, +9q). We also observed negative pressure which triggered progressive loss of the remaining SCAs (+ 7q, + 11q, 13q-, 19q-, 20p-). Schwann cells seem to enhance clonal selection in NB in a more biomimetic context; however, the genomic differences found between hydrogels with vs without Schwann cells are minor and affect only the percentage of cells with some SCAs. The evolutionary adaptation we have detected may be due only to deterministic events (clonal selection) or also to interaction with stochastic processes (generation of mutations and genetic drift). However, the fact that genomes in two independent models, both with tumor spatial structures related to cellular aggressiveness, show similar genomic aberrations suggests a predominance of the deterministic evolutionary force. The increase in selection pressure, driven by adaptation to niches deliberately altered (by their lack of systemic VN, and by high stiffness and/or long 3D growth time) in both models, could prompt this similar clonal selection towards genetic aberrations detected by probably advantageous cellular phenotypes [57]. Both niches, which initially lacked VN, showed a large amount of VN secreted by neuroblats [22]. Our results show new evidence for the probably role of VN in biotensegrity mechanotransduction by mediating the cellular response to ECM stiffness [22,56]. The adaptation process was faster in VN-KO tumor cells but less manageable than in hydrogels. Furthermore, cells grown in rigid 3D hydrogels showed faster selection for SCAs and less presence of intratumoral heterogeneity than those cells with long and soft 3D growth, which could re ect a more e cient adaptation to niche [39] as occurs in aggressive patients tumors with rigid ECMs [58].
Conversely, in 3D hydrogels with short culture times, cells may had not yet been subjected to the selective process, due to their incipient adaptation to ECM. As in the early stages of patient tumors, cells in 3D bioprinted hydrogels need to grow and survive, and then adjust genetically and epigenetically to acquire the features they need to take on new roles such as migration. For this among other characteristics, 3D models are increasingly important, reducing the use of animal models and controlling the parameters focused on the study, adding degrees of complexity to gain deeper insight into the tumor microenvironment.
Interestingly, some genes of the new SCAs undergoing positive selection in VN-KO tumors and stiff and/or long-cultured hydrogels are involved in ECM composition and architecture, mechanotransduction and cell migration. In chromosome 1 there is special focus on the collagen gene COL11A1 (in the nearer centrometic fragment of 1p-), the genes related to actin ARPC5 and ACTN2, some laminins and the integrin α10 (in + 1q), while the deleted FSCA of chromosome 9p contains the genes KANK1 and DOCK8, both involved in actin polymerization [59][60][61]. In our customized NB-mechanopanel, COL11A1 and DOCK8 genes speci cally showed certain gene variants that have disappeared in the tumors from VN-KO mice with the mentioned heterozygote deletions. The extent to which this occurred, even when the cell population with these mutations was not inconsiderable in the other growing conditions (allelic frequency: 0.4), raises the question of whether these variants might play a role in adaptation. VN of ECM could have a swift impact on the mutational pro le of these migration-related genes. Chromosome 9 aberrations have not been a focus of interest so far in HR-NB, unlike those of chromosome 1 (1p-and + 1q) [32,62,63]. Therefore, we decided to delve on 9 chromosome abnormalities in a cohort of primary human HR-NB. Chromosome 9 aberrations involving KANK1 and/or DOCK8 genes were detected in 39.53% of the MNA primary HR-NB analyzed. Consequently, we recommend a spotlight on chromosome 9 aberrations [64][65][66] affecting DOCK8 and KANK1 in NB, and especially on the sometimes denigrated FSCAs, as they could involve invasion in MNA HR-NB, as occurs in other malignant tumors [67]. The relevance of these new genetic alterations found in human HR-NB, especially the frequency and prognostic impact of DOCK8 and KANK1, should be determined in collaborative studies.
The SK-N-BE(2) cell line has been characterized as particularly heterogeneous [42]. Likewise, high genetic instability of xenograft mice, showing gains and losses of SCAs with passages, has already been described [57]. In all our in vitro and in vivo samples only two SCAs remained stable (with the same length and break points) in most cells: partial deletion of 3p and 17pq, adding MNA. 3p loss is a common event in HR-NB [68]; the region contains certain important suppressor genes, including RHOA. RhoA protein plays an signi cant role in mechanotransduction, as it mediates focal adhesions and stress ber formation [69]. Its reduced signal (associated with 3p) has been linked to growth and metastasis in colorectal tumors [70,71]. Moreover, in these tumors RHOA downregulation is used to match the activation of Wnt signal pathway and inactivation of TGFβ. All together they participate in initiation and progression of tumors [70]. According to Gene Onthology, the Wnt pathway is one of the most affected by genomic aberrations of SK-N-BE(2) cells, and most of the altered genes of these pathways are in SCAs that changed (were positive or negative selected) in VN-KO tumors and 3D bioprinted models (+ 7q, +9q, + 11q, 13q-). TGFβ is located in the 19q deleted region which also had reduced presence in the samples mentioned. This gene has dual action: progression of tumors when underexpressed, and degradation of the ECM when its levels increase [72]. Meanwhile, 17pq-contains the tumor suppressor gene TP53, which plays an important role in tumorogenesis and aggressiveness [73]. This gene was also mutated in all samples from SK-N-BE(2) cell line with an alleleic frequency close to 1. The double mechanism involving a mutation in one of the alleles and a deletion or CNLOH of the other one has been highlighted for its combined negative impact [73]. In addition, the reported imbalance of chromosome 17, in which a loss of 17pq has been detected, is usually found with the gain of the remaining fragment of the chromosome, that is, a typical + 17q [74]. 17q gain is precisely the most common genomic aberration of NB [31]. In fact, the hidden + 17q became evident in some of our hydrogels and in VN-WT tumors. A probable explanation for this is that in addition to overexpression of tumor progression-promoting genes within the 17q fragments, the loss of function of genes localized on 17p-or CNLOH of 17p (such as TP53) may play a signi cant role in neuroblastoma development [74].
The impact of ECM properties on genomic heterogeneity of the SH-SY5Y cell line was not evident when using our models. These cells may be less sensitive, or would adapt to ECM stiffness/composition differently, without SCA changes or longer-term SCA changes. Two research groups have described morphological and gene expression changes in this cell line in stiff environments [53,54]. It is possible that ALK-mutated cells respond to ECM changes through epigenetic alterations or genetic mutations undetected by the techniques we employed [53,[75][76][77].

Conclusions
Our study is further proof of the impact of biotensegrity on cancer evolution, in which dominant clones could be selected from modi ed sub-clones of the primary tumor by its adaptation to tumor microenvironment changes. We have described a similar genomic response of a MNA NB cell line grown in two different models (xenograft VN-KO mice and stiff hydrogels) with a deliberately altered niche by its composition and stiffness. Whether the emergence of tumor cells with genomic pro les distinct from the parental cell lines was the result of clonal selection and/or stochastic evolutional processes should be de nitive addressed by single cell DNA sequencing of the tumors. As MNA human HR-NB also showed chromosome 9 aberrations affecting interesting genes (DOCK8, KANK1), functional validation as well as cell proliferation and migration studies in similar models with a large number of HR-NB cell lines are needed for future targeted treatments. Although 3D bioprinted and in vivo models could have differences in clonal selection pressures, more accurate and biomimetic 3D scaffolds are useful to study genomics and clonal evolution of tumor cells under different stiffness conditions. A better understanding of mechanotransduction pathways related to ECM stiffness and VN involment will expand therapeutic possibilities, and 3D bioprinted scaffolding represents an excellent tool to test the e ciency, toxicity and impact of treatment on genomics and cell heterogeneity. Moreover, strengthened studies in mouse model system that give rise to spontaneous NB tumors that recapitulated the full range of NB progression would also be adequate. Finally, large collaborative studies with cells derived from primary and metastatic biopsies or tumoroids grown under different conditions, which progressed in the context of a syngeneic interaction between tumor and stromal cells, are necessary to associate SCAs with the biotensegrity response of MNA and ALK-mutated NB cells.

Consent for publication
Not applicable

Availability of data and materials
The datasets used and analyzed in the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.

Funding
This study was supported by ISCIII (FIS) and FEDER (European Regional Development   (25)). When the percentage is less than 30% the background color of the aberration is lighter. The striped background and the asterisk (*) after cell percentage indicates aberrations only found in the HD-SNPa of ctDNA.   (A) Neighbor-joining tree clustering and matrix of similarity and distance indices based on the Jaccard coe cient. The highest Jaccard coe cient between samples (equal or near to 1) are represented in blue in the matrix, and with nearby branches of the same color in the tree. Lower Jaccard coe cients show a progressive fading towards yellow (equal or close to 0) and a greater distance between tree clusters. Note that softer and/or shorter cultivated hydrogels had higher similarity indices with VN-WT P0-P2 (remarked in black boxes) than with VN-KO tumors (remarked in red boxes) and are also nearer in the tree.
Comparably, stiffer and/or longer cultivated hydrogels correlated better with VN-KO tumors (remarked in purple boxes) than with VN-WT ones (remarked in pink boxes), showing even smaller indices with P3 to P5. Stiffer and/or longer-time cultured hydrogels and VN-KO have been clustered closer in the tree. (B) Genomic aberrations and their relationships with the models according to their presence and length, and related to the Jaccard proximity at the tree. Multivariate analysis was performed with PAST software in function to the chromosomal aberrations detected in each sample derived from the SK-N-BE(2) cell line.
Presence of the SCAs in each sample is represented by light grey squares, SCAs as CNLOH by dark grey squares, and absence of genomic aberration is denoted by the remaining colors.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.