Genome-wide analysis of Ollier disease: Is it all in the genes?

Background Ollier disease is a rare, non-hereditary disorder which is characterized by the presence of multiple enchondromas (ECs), benign cartilaginous neoplasms arising within the medulla of the bone, with an asymmetric distribution. The risk of malignant transformation towards central chondrosarcoma (CS) is increased up to 35%. The aetiology of Ollier disease is unknown. Methods We undertook genome-wide copy number and loss of heterozygosity (LOH) analysis using Affymetrix SNP 6.0 array on 37 tumours of 28 Ollier patients in combination with expression array using Illumina BeadArray v3.0 for 7 ECs of 6 patients. Results Non-recurrent EC specific copy number alterations were found at FAM86D, PRKG1 and ANKS1B. LOH with copy number loss of chromosome 6 was found in two ECs from two unrelated Ollier patients. One of these patients also had LOH at chromosome 3. However, no common genomic alterations were found for all ECs. Using an integration approach of SNP and expression array we identified loss as well as down regulation of POU5F1 and gain as well as up regulation of NIPBL. None of these candidate regions were affected in more than two Ollier patients suggesting these changes to be random secondary events in EC development. An increased number of genetic alterations and LOH were found in Ollier CS which mainly involves chromosomes 9p, 6q, 5q and 3p. Conclusions We present the first genome-wide analysis of the largest international series of Ollier ECs and CS reported so far and demonstrate that copy number alterations and LOH are rare and non-recurrent in Ollier ECs while secondary CS are genetically unstable. One could predict that instead small deletions, point mutations or epigenetic mechanisms play a role in the origin of ECs of Ollier disease.


Background
Enchondroma (EC), a benign cartilage forming tumor in the medulla of the bone, is thought to be a precursor of secondary central chondrosarcoma (CS). EC develops either as a single, solitary lesion or as multiple lesions in the context of Ollier disease [1]. Ollier disease is the most common subtype of enchondromatosis and shows multiple ECs with marked unilateral predominance [1,2]. The risk of malignant transformation towards central CS is up to 35% in Ollier disease [1,3]. There is no marker that would indicate progression towards malignancy, thus there is a vital need to understand the genetics of these tumors which may help to develop markers for early diagnosis. A comprehensive understanding of the molecular events in ECs and central CS also enables the identification of possible targets for treatment [4].
While the genetics of enchondroma is poorly understood, the involvement of the EXT genes is well established in the development of solitary as well as hereditary multiple osteochondromas (MO) (OMIM 133700), benign cartilage tumors at the surface of bone [5]. The lack of EXT function seems to disturb hedgehog signalling in MO, while activated hedgehog signalling in mice seems to underlie the development of the Ollier related phenotype [4]. Heterozygous mutations in PTH1R are found in a small subset of Ollier patients [6][7][8]. It is however unclear whether these mutations in PTH1R are causing or modifying the disease [7], and since~90% of Ollier patients lack PTH1R mutations, we aimed to study Ollier related ECs by mapping genetic changes using genomic arrays.
We hypothesized Ollier disease to be a germ-line mosaic condition due to the fact that it is non-hereditary and because of its unilateral predominance feature [3,9]. An early postzygotic mutation resulting in asymmetric involvement of skeletal structures can be expected, as was also shown for polyostotic fibrous dysplasia [10]. One could speculate that an inactivating mutation in a tumor suppressor gene, similar to EXTs in osteochondroma, would have occurred in the developing mesoderm early after gastrulation. In case of a tumor suppressor gene, later on, an additional hit may be required for the formation of ECs with subsequent genetic changes causing progression towards central CS. We tested this hypothesis using high-resolution SNP array combined with expression array on DNA derived from tumor tissue and, whenever available, normal DNA from Ollier patients.
SNP arrays provide an excellent possibility for largescale, genome-wide, high-resolution analysis of both DNA copy number alterations (CNA) and loss of heterozygosity (LOH) in cancer cells. It also provides a feasible means of detecting genotyping alterations in the tumors of individual patients and, in principle enables the identification of new areas with common allelic imbalance that could harbor potential tumor suppressor genes which helps in the identification of novel candidate genes affected by genomic abnormalities [11,12]. In the present study, we used Affymetrix Genome-Wide Human SNP Array 6.0 to obtain a comprehensive registry of genetic aberrations in 37 tumors of 28 patients with Ollier disease and correlate it with gene expression using Illumina Human-6 v3 Expression Array and qRT-PCR and protein expression using tissue microarray (TMA). Based on the obtained genomic profiles with limited and non-recurrent genetic changes in Ollier ECs, we conclude that they are genetically heterogeneous and that the reported CNA in this study are likely to be secondary random events in ECs.  [13]. Diagnoses were originally made in the multidisciplinary teams of the centers of origin. The histology was revised and representativity was assessed on the available paraffin or frozen tissue by one experienced bone tumor pathologist. For SNP array analysis, 14 Ollier ECs, and 23 Ollier CS (13 grade I, 8 grade II, 2 grade III) from 28 Ollier patients were used. As controls, normal DNA derived from fresh frozen muscle tissue

DNA and RNA isolation
Tumor samples were selected that contained more than 80% of tumor cells, as estimated on haematoxylin and eosin-stained frozen sections. Most of the samples were macro dissected and L2099 sample was micro dissected to enrich the tumor percentage [14]. DNA from fresh frozen tissue was isolated using the wizard genomic DNA purification kit (Promega, Madison, WI), according to the manufacturer instructions. Blood DNA was isolated as described by Miller et al [15]. Saliva DNA was isolated using the Oragene DNA kit (DNA Genotek Inc Ontario, Canada) according to the protocol provided by the supplier and DNA was precipitated using sodium acetate precipitation. DNA concentration was quantified spectrophotometrically using Nanodrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA) and the fragment sizes were determined on 1% agarose gel.
RNA isolation from the fresh frozen tissue was performed as described previously [16]. RNA concentration was measured spectrophotometrically and the fragment sizes were determined by RNA 6000 Nano LabChip kit using Agilent 2100 Bioanalyzer (Santa Clara, CA, USA).
DNA and RNA from all the samples were good enough to continue with the experiment.

SNP Array and Data Analysis
We used the Affymetrix Genome-Wide Human SNP Array 6.0. Genomic DNA preparation, labeling and hybridization were performed according to Affymetrix's recommended protocols (Affymetrix, Santa Clara, CA, USA). Then, arrays were scanned with GeneChip™ GSC3000 7G Whole-Genome Association System (Affymetrix). Overall hybridization quality was estimated by the genotype call rate using the Birdseed genotype calling algorithm in Genotyping Console (version 3.0.2, Affymetrix). Average call rate was 97.83%. To analyze the data we used statistical language R version 2.8 and Nexus software version 4.1 (BioDiscovery, CA, USA). We did not use HapMap samples as baseline in this study to avoid the bias for the experiment performance at different labs, batch effect and hybridization quality. The analysis was performed on a subset (30 controls, 14 ECs and 23 CS of Ollier patients, Table 1) of a larger experiment of 92 samples including samples unrelated to Ollier disease to achieve a larger set of common controls. Than we performed copy number analysis using 92 samples as a baseline in R software and only 29 control samples of high quality as a baseline in Nexus software. Results that we obtained using different softwares and different baselines were comparable. We used CEL files in R software. For the genomic analysis using R, we did genotyping using the CRLMM algorithm in the Oligo package [17], copy number analysis using the aroma affymetrix package [18], and we constructed LAIR plots to visualize regions of LOH and allelic imbalance [19]. Chromosome-X was not analyzed to avoid gender-related issues [20]. In Nexus, we performed copy number analysis using CNCHP log-ratio files generated by genotyping console using 29 controls as a baseline. Hidden Markov model (HMM) based SNP-FASST segmentation was used to identify aberrant genomic regions. Here we considered at least 5 probes for each segment. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/ accession number GSE22965).

Multiplex Ligation-Dependent Probe Amplification (MLPA)
MLPA was used to confirm copy number gains and losses found within the selected candidate genes by SNP array. MLPA probes were designed using NCBI Build 36.1. We used two probes for TCRA and one probe for POU5F1, ANKS1B, FAM86D and PRKG1. Probe sequences can be obtained upon request. MLPA was performed as described previously [21]. Sample series of SNP array (Table 1) and in addition to that DNA from ten HapMap samples was used. Data analysis was done using SoftGenetics GeneMarker version 1.70. We set 1.2 and 0.8 as a threshold to detect the gains and losses respectively.

Expression Array and Data Analysis
In total, we analyzed 7 ECs while 2 growth plate and 4 articular cartilage samples served as a control (Table 1). Expression array was performed using Illumina Human-6 v3.0 Expression BeadChips (Illumina Inc., San Diego, CA). For Illumina BeadArray assay, cRNA was synthesized with an Illumina RNA Amplification Kit (Ambion, Austin, TX, USA), purified, labeled and hybridized as per the manufacturer's instructions. Then, arrays were scanned with Illumina BeadArray reader (500GX, Illumina) and scanned images were imported in BeadStudio software version 3.1.3.0 (Illumina). For the expression array, unprocessed data were collected from BeadStudio and processed using VST [22] and quantile normalization with the R bead array package [23]. The processed data was analyzed using Linear Models for Microarray Data (LIMMA) analysis which uses moderated t-test to detect differentially expressed genes between two groups by taking into account natural variance within these groups and correcting for multiple testing using false discovery rate [24].

Real-time quantitative Reverse Transcriptase PCR (qRT-PCR)
qRT-PCR was used to confirm the results obtained by the integration of SNP and expression array. Here we have used 5 controls (2 growth plates and 3 articular cartilage) and 8 ECs (Table 1). RNA used for expression array was also used for the qRT-PCR. cDNA was synthesized using 1 μg of total RNA with AMV reverse transcriptase (Roche Applied Science, Almere, The Netherlands) in combination with oligo-dT and random hexamer priming and qPCR was performed as described elsewhere [25]. Primer details can be provided upon request. Expression of the genes of interest (FAM86D, POU5F1 and ANKS1B) was normalized by geometric averaging of multiple internal control genes using the geNorm program [26]. Out of four normalization genes the best three were selected within this program: GPR108, SRPR and TBP. Relative quantification was performed using standard curves, followed by adjustment with the normalization factor calculated by geNorm program [25]. The average relative expression of gene of interest in ECs was compared to controls. To see the EC specific gene expression changes, relative expression of ECs with change in copy number was compared with ECs without change in copy number.

Tissue Microarray (TMA)
Five TMA blocks were constructed using TMA master (Zeiss, 3D Histech, Hungary) and each block contain maximum 70 cores including seven control tissues (growth plate, articular cartilage, breast carcinoma, prostate, colon, skin and tonsil) for orientation purpose. In total 5 TMA blocks contain 86 tumors, of which 65 Ollier related and 21 solitary ECs and CS (Table S1, Additional file 1). The TMA blocks contain 2-mm cores of each tumor in triplicate.

Immunohistochemistry
Tissue sections were cut from TMA blocks and dried overnight at 60˚C. Slides were kept in Xylol for 20 minutes. Immunohistochemistry was performed as described earlier [27]. Detailed information on the antibodies used to check protein expression of POU5F1 and NIPBL are given in table S2, additional file 1. We used power vision (poly-HRP-GAM/R/R, Immunologic) as a secondary antibody. Visualization was carried out with liquid Dab + substrate chromogen system (DAKO, Glostrup, Denmark). As a negative control primary antibody was omitted. Immunostained TMA slides were scanned using a high resolution Mirax Desk instrument (Zeiss, Mirax 3D Histech, Hungary) and scored independently by two observers (JVMGB and TCP) with the Mirax viewer TMA module software version 1.1.12 (Zeiss) and discrepancies were discussed. In brief, the intensity (0 = no staining, 1 = weak, 2 = moderate, 3 = strong) and percentage of positive tumors cells (0 = 0%, 1 = 1-24%, 2 = 25-49%, 3 = 50-74%, 4 = 75-100%) were assessed. Cores with a negative internal control and loss of tissue were excluded from the analysis. A sum score ≥ 4 was considered positive. Statistical analysis was done using Pearson Chi-Square test in SPSS (version 16.0, Chicago, Illinois, USA).

Results
Genome-wide detection of copy number alterations and loss of heterozygosity using SNP array Samples L1975, L1974, L1980, L813, L1083 and L286N were excluded from further analysis because of quality issues and low call rates. All samples including 29 controls and 32 tumors (14 ECs, 18CS) showed DNA copy number aberrations mainly restricted to known variable regions ( Figure 1) [28,29]. CNA were more frequently found in CS II and CS III as compared to ECs and CS I (Figure 1). The number of copy number changes in saliva and blood DNA were comparable.

Genetic alterations in Ollier Enchondromas
We used paired analysis which is based on the comparison of tumor and corresponding normal DNA (available for 5 ECs from 4 patients) to study LOH and CNA. Although sequences of homozygous SNPs were identified that could indicate LOH, these same sequences were also observed in the corresponding normal sample. We could not find any LOH in these 5 ECs using both R and Nexus softwares. We have identified 7 EC specific CNA by paired approach ( Table 2). Selection of candidate genes TCRA, ANKS1B and PRKG1 for further validation is based on copy number change in minimum 10 probes within the gene.
Unpaired analysis revealed absence of LOH in the majority of ECs. We confirmed the loss of chromosome 6 with LOH in L206 ( Figure 1B) which was in agreement with the results published previously by our group using array CGH [30] proving validity of the assay. Also loss of one copy of chromosome 3 and 6 with LOH was found in L1683. These results were confirmed using R and Nexus software. Furthermore, an unpaired approach (29 controls as a baseline) was used to find most common copy number gains and losses in at least 5 out of 14 ECs (Table S3, Additional file 1) using Nexus. None of these were confirmed in paired comparison which suggests that they are not tumor specific changes. When ECs are located in the phalangeal bones, cellularity is increased [1]; we could not find differences between ECs or CSI of the hand versus those of long bones at the genomic level although the sample sizes are small.
Homozygous deletion of FAM86D at chromosome 3p12.3 was found in two ECs in one patient (L206 and L910) and selected for further validation (Figure 2). Interestingly both ECs were located in different digits of the hand. It was not possible to get normal DNA from the same patient to investigate tumor specific loss of FAM86D.

Genetic alterations in Ollier chondrosarcomas
Ollier CS showed large genomic alterations and LOH at various locations in the genome (Table 3, Figure 1C). Frequent recurrent CNA involves chromosomes 3p, 5q, 6q and 9p.

Verification of gene copy number alterations by MLPA
To validate regions identified by SNP arrays, we performed MLPA on 37 Ollier tumors (Table 1) and 38 controls for the candidate genes FAM86D, TCRA, ANKS1B and PRKG1 (Table 4). In addition, a probe was designed at POUF51 which is at chromosome 6p21.31 in order to confirm loss of chromosome 6 in L206 and L1683. We have used 38 controls as a baseline. Results of validation experiments are summarized in Table 4. In short we confirmed three (FAM86D, PRKG1 and monosomy 6) out of five candidate regions. We confirmed the homozygous loss of FAM86D in two ECs of the same patient ( Figure 2). Gain at PRKG1 in L1490 is confirmed which seems to be EC specific in this case but was not found in other ECs and CS. Loss at chromosome 6 using probe at POU5F1 was confirmed in L206 while loss was not confirmed in L1683 at given thresholds. TCRA region at chromosome 14q11.2 is known as a highly polymorphic region [31] and its frequent rearrangements are observed in blood lymphocytes [32]. We used majority of blood samples as a baseline and found copy number gain in all ECs at TCRA as a result of SNP array unpaired analysis. Paired analysis with SNP array showed copy number gain at TCRA in L1978 and was not confirmed with MLPA. Indeed this gain was due to loss in the corresponding normal DNA. Loss at ANKS1B in L2220 was not confirmed by MLPA at given thresholds however the peak was lower in L2220 compared to corresponding normal.

Expression array and its integration with copy number alterations
Expression array was performed using 7 ECs and 6 controls with Illumina Human-6 v3.0 (Table 1). We performed function based analysis by integrating the gene expression results with SNP array results. In total 1044 genes were differentially expressed in ECs compared to controls (adj. p-value < 0.01, Table S4, Additional file 1). We considered all up regulated genes (881 genes, adj. p-value < 0.01) with presence of gain in at least one EC and all down regulated genes (163 genes, adj. p-value < 0.01) with presence of loss in at least two ECs. We found NIPBL which was gained as well as up regulated while POU5F1 that was lost as well as down regulated in ECs compared to controls. The same pattern was found in CS.

Validation of integration approach of SNP and expression array using qRT-PCR
Using qRT-PCR, we confirmed lower expression of ANKS1B, FAM86D and POU5F1 (Table 4). L2220 showed one copy loss of ANKS1B and its expression was decreased in this EC as compared to the average of relative expression of other 7 ECs however there is no  difference in the expression comparing all ECs versus controls. SNP array results revealed homozygous loss in L206, L910 and one copy loss in L2103 for FAM86D. FAM86D expression was decreased in ECs compared to controls. There was no expression of POU5F1.

Protein expression
Using a TMA we demonstrated that POU5F1 protein expression was absent in all Ollier tumors which is in line with the expression array and qRT-PCR results. For NIPBL, 30% of Ollier tumors were positive ( Figure 3, Table 4). Only 34/65 Ollier and 11/21 solitary tumors could be analyzed since cores from the rest of the tumors were lost during the immunohistochemistry. There was no significant difference in the expression of NIPBL within the tumors of different grades in Ollier disease (Pearson Chi-Square, p-value = 0.1).

Discussion
The origin of both solitary and Ollier related ECs is largely unknown. To address this, we performed genomewide analysis of ECs occurring in non-hereditary Ollier disease. Since these tumors are polyostotic, with a unilateral predominance, manifesting early in life, we postulated that there may be a germ-line mosaicism. We attempted to find causative genes for ECs of the Ollier disease using a high-resolution SNP array containing 1.8 million markers combined with expression array and obtained comprehensive genetic profiles of 37 Ollier disease related tumors. This is the first and largest genomewide molecular study on Ollier disease reported so far, which was possible through the collaboration of many different institutes within the EuroBoNeT Network and the European Musculo-Skeletal Oncology Society (EMSOS).  In general, the obtained genomic profiles showed absence of large genetic aberrations in Ollier ECs except loss of chromosome 6 in two ECs from two unrelated Ollier patients. Small non-recurrent genetic changes combining the SNP and expression array at FAM86D, PRKG1, ANKS1B, NIPBL and POUF51 were found in ECs. Most of these genes are not known to play an important role in cartilage formation. We found homozygous loss of FAM86D in two ECs of the same patient. Function of FAM86D is still unknown. We confirmed intronic gain at PRKG1 in one EC, which is involved in fatty acid metabolism [33]. We found loss at ANKS1B, while overexpression of ANKS1B is reported in pre-B cell acute lymphocytic leukaemia [34]. Gain at NIPBL was found in L2205 while at the protein level only 30% of Ollier ECs and CS expressed NIPBL. Inactivating mutations in NIPBL are associated with Cornelia de Lange syndrome and one of the characteristic features of this syndrome is reduction in limb growth (OMIM 122470). Loss of POU5F1 was found in two ECs with monosomy at chromosome 6 and its mRNA and protein expression was absent in all Ollier and solitary ECs and CS. The transcription factor POU5F1 (OCT3/4) is involved in regulating pluripotency and is normally expressed during early embryogenesis in embryonic stem and germ cells [35].
Here we studied extensively candidate genes obtained from paired analysis of enchondromas. Normal DNA was available from 4 Ollier patients enabling paired analysis. Despite the low number of paired samples our data suggest that no common CNA are associated with EC development. Extending the analysis with the unpaired samples we could not see any common CNA in all ECs. All aberrations we obtained in at least 5 out of 14 ECs are reported as common copy number variants in DGV database (http://projects.tcag.ca/variation/). Loss of chromosome 6 was the only recurrent change in ECs of two unrelated Ollier patients. Therefore, relatively small numbers of copy number alterations that we found per ECs are more likely to be secondary random genetic changes. Although SNP array technology is a powerful analysis tool, it can not detect balanced chromosomal translocations, inversions and whole-genome ploidy changes. However, balanced translocations and inversions have not been reported for the Ollier tumors in the literature so far [2].
Previously, PTH1R was reported to be the gene causing Ollier disease [6]. However, in subsequent studies it was shown that PTH1R is only mutated in a sub group of patients (~10%) decreasing receptor function tõ 70%, suggesting that it may contribute to the disease but is probably not causative [7]. PTH1R is a key player within the IHH pathway which is crucial for endochondral ossification. The presence of known point mutations (R150C, R255H, G121E and A122T) in PTH1R was excluded in the present series (unpublished data). Also, we could not find a deletion or LOH at the 3p21.31 region harboring the PTH1R gene. Recently, inactivating mutations in PTPN11 are reported in another enchondromatosis subtype called metachondromatosis in which multiple ECs are combined with osteochondroma-like lesions [2,36]. In our series we could not detect any CNA or LOH at PTPN11 by SNP array. Also, expression of PTPN11 was not decreased in ECs as compared to controls in expression array.
Large gains, losses and LOH were seen more often in Ollier CS as compared to ECs, which is in line with increased genetic instability and aneuploidy seen in solitary central chondrosarcoma progression. Most common CNA involve 3p, 5q, 6q and 9p in Ollier CS. However, we could not detect recurrent CNA in all Ollier CS that could have been responsible for malignant transformation of ECs. For Ollier CS, a deletion at the short arm of chromosome 1 [37], LOH at RB1 at chromosome 13 and p15/p16 loci at the short arm of chromosome 9 [38] were reported in single cases each. Our results show very few copy number alterations in ECs and an increased number of variable genomic alterations in CS. This is in support of the multistep model for CS development [39].
In conclusion, we present genome-wide copy number and expression profiles of the largest international series of Ollier ECs and CS reported so far. We show absence of recurrent CNA and LOH in majority of ECs suggesting that instead point mutations or other copy number neutral structural changes (inversions, insertions, balanced translocations) or deletions below the resolution of this platform involving a single or a few exons only [40] might have an important role in EC pathogenesis. This opens the possibility to study these tumors further using a next generation sequencing approach. An increased number of genetic alterations are found in Ollier CS, supporting the multistep genetic progression model.

Additional material
Additional file 1: Table S1, Table S2, Table S3, Table S4. Table S1 -Clinicopathological data of 86 tumors used in TMA -* Gender information was not available for four Ollier patients. Table S2 -Antibody  information Table S3 -Unpaired copy number changes in 35% of EC (min. 5 out of 14) -All the copy number gains and losses that we found in more than four ECs are an indicative of known copy number variation in DGV database of genomic variants. Table S4