A predictive biomarker panel for bone metastases: Liquid biopsy approach.

Bone metastases is one of the common metastatic site and leading cause of cancer-related mortality in progressive cancer patients. The purpose of the present study is to establish a liquid biopsy based multi-gene classifier and associated signalling pathways for early diagnosis of bone metastases. We used publically available microarray datasets and analysed them in a platform/chip-specific manner using GeneSpring software. Analyses of gene expression datasets identified 15 consistently over-expressed genes with statistical significance. Further, expression profile of same set of 15 genes were compared in breast and lung cancer exosome derived mRNA with (n = 10) and without (n = 10) bone metastases against healthy controls. ROC curve analysis performed individually for all the 15 genes shortlisted the 5 most relevant genes with significant sensitivity and specificity in both cancers. This liquid biopsy-based bone metastases predictor using multi-gene panel is a unique approach with potential clinical applications for effective management of aggressive cancers.

Introduction

Bone metastasis is the leading reason for mortality and morbidity globally and its occurrence is immensely growing [1]. The bone is the third most common site for metastasis, after the lung and liver [1], [2]. Cancer metastases to the bone are more widespread among patients with progressive cancer of the breast (73%), prostate (68%), or lung (36%) [3]. Bone metastases lead to skeletal-related events (SREs) in final-stage cancer patients which define as spinal cord compression, the essential requirement for radiation or surgery to bone, pathologic fracture, and hypercalcemia [4]. Unfortunately, current beneficial therapeutic options are in most cases the only palliative and, although not curative, surgery remains the noteworthy effective for bone metastasis treatment [5]. The growing global incidence of bone metastasis along with the distinct lack of noteworthy development in overall survival rate during the past four decades requires newer methodologies to recognize the molecular mechanisms of the disease and to classify potential markers for early detection, prognosis, and as targets for therapy [6]. Consequently, identifying impending pathological genes and signalling pathways inducing the bone metastasis process is very important.

It is a major target in primary cancer research to regulate the genetic mechanisms that support the bone metastatic processes, which include: tumor cell intravasation, cell survival during circulation, extravasation into new tissues, and successfully inhibited growth at a secondary site [7]. A recent study of exosomes has suggested that early stage tumor release extracellular vesicles carrying numerous types of tumor markers [8]. Exosomes are particles released from almost all kinds of cells that are enclosed by a lipid bilayer and cannot replicate [9]. The mechanism of exosome biogenesis and contained cargo within are still not yet totally understood, nevertheless studies have indicated that exosomes can transport cargo of DNA, RNA, proteins and modulate target cells [10], [11]. Therefore, exosomes provide a wide range platform for the discovery of new liquid biopsy based biomarker for cancer using combine and individual cargo molecules [12]. Nowadays, numerous studies have to specify that tumor-derived exosomes are responsible for cancer progression which, shows that TDEs may hold abundant potential for cancer diagnosis, prediction, and treatment response assessment. However, the significant function of exosomes assists in the pathological communication between primary tumor and bone cells within the bone microenvironment remains a developing field.

Investigative approaches that can develop existing databases and classify markers concordant through studies may be an advantageous methodology to narrow down to markers of increased assurance and clinical applicability. Moreover, relevant studies have been reported on cancer of numerous sites. Meta-analysis based approaches in head and neck cancer have led to the identification of novel targets and candidate biomarkers [13]. Parmigiani and the team excellently applied meta-analysis of gene expression for the molecular classification of lung cancer [14]. In prostate cancers, markers causal to the carcinogenic process were identified utilizing a similar approach [15]. Besides, markers highly related to diagnosis and treatment outcome prediction in breast cancer were also recognized through similar meta-analysis based approaches [16], [17]. In this study, we hypothesize that a meta-analysis of publicly available bone metastasis genomic expression datasets of numerous cancer types can classify a common metastatic signature of bone metastasis. We tested this hypothesis by applying a cross-platform and cross-study meta-analysis method on multiple microarray datasets and subsequently validate them in clinical samples to introduce the benefits of exosomes as liquid biopsy based approaches.

Materials and methods

Search criteria and data mining

Investigation of the generously accessible microarray datasets was carried out as per the PRISMA guidelines [18]. An electronic database, explored was performed using the Gene Expression Omnibus (GEO) (NCBI, http://www.ncbi.nlm.nih.gov/geo/) and Array Express (EBI) (http://www.ebi.ac.uk/arrayexpress) for the occurrence of raw data of microarray experiments carried out in bone metastases from the different primary site. The sequence of workflow for the analysis involved the following steps, i) data mining and retrieval of numerous databases to identify microarray studies that compared the expression of primary tumors of various origins (prostate, colorectal, renal, and breast) against distant bone metastases. ii) Studies containing only human tumors and iii) studies, including universal profiling of transcriptome using microarrays. Furthermore, Affymetrix [Affymetrix Inc., California, USA], Agilent [Agilent Technologies, California, USA] and Illumina [Illumina Technologies, San Diego, USA] three platforms were used for the meta-analysis. The common basic workflow followed the stepwise protocol recommended by Ramasamy et al for carrying out a meta-analysis [19]. The general workflow is detailed in Fig. 1.

Fig. 1

Bone metastases meta-analysis workflow. The freely and publically existing raw microarray data of bone metastases series were downloaded and grouped and analyzed in genespring statistical software.

Data analysis using GeneSpring

The raw data files were used for the analyses included .CEL (Affymetrix platform), .TXT (Agilent) and .CEL (Illumina) files that downloaded from the GEO and further analyzed using genespring software (http://genespring-support.com). Furthermore, the raw data were uploaded onto the genespring software then baseline transformed and normalized by Robust Multi-array Analysis (RMA) in Affymetrix or 75th percentile in Agilent platforms (single color) or Illumina platforms. The isolated sample files were then classified into 'Primary' and 'Metastasis' and re-analyzed as a single experiment. The experimental data at the gene level (arithmetic mean of all probes mapping to the same probe ID) was generated and quality control was approved using Principal Component Analysis (PCA) in GeneSpring. Moreover, in the PCA analysis outlier sample was removed and clustering carried out afterward to confirm a clear stratification between the two categories of primary and metastatic samples. Fold change analysis was then performed on the samples following which an unpaired t-test (unequal variance) were performed to obtain significant gene entities. The p-value computation (asymptotic) and multiple testing correction (Benjamini Hochberg FDR) were further performed to obtain gene entities with p-value < 0.05 and fold change (FC) of >2.0. Additionally, identified gene pathway analysis was done using the inbuilt pathway analysis tool. The individual gene entity list from each technology was extracted from the GeneSpring software and exported to excel files. A similar approach was adopted to identify all the common genes, and pathways across the three platforms of Affymetrix, Agilent, and Illumina (Fig. 1) [13].

Meta-gene signature prediction using a Venn diagram

The different site of primary to bone metastases sample was compared using a Venn diagram through to find out the gene signature for bone metastases. The comparison between up/down-regulated genes in bone metastases from primary breast, prostate, colon, and renal.

Functional annotation using multiple web source

The concordant gene list across the various platforms was analyzed in the diverse web resources to assess functional classes and protein–protein interactions. Gene ontology analysis and heat map were generated using FunRich and Gprofiler classification system to assess the functional classes of the genes [20]. The STRING and NetworkAnalyst database was used to predict and catalogue the protein–protein interactions between the concordant genes [21].

Patient-based validation

Sampling was carried out at Vedant and Civil hospital during routine Fine Needle Aspiration Cytology (FNAC) procedure as a part of the diagnostic workup. The validation of selected genes predicted from meta-analysis was done in exosomal RNA derived from serum sample using quantitative real-time PCR. 10 blood sample of healthy persons, 10 blood samples of primary breast and lung and 10 blood sample from patients with bone metastases were collected with prior consent. The study was approved by the Gujarat university ethics committee and all samples were collected after obtaining written informed consent from each patient. The median age of the patients was 60 years at diagnosis, ranging from 30 to 85. Clinical-Pathological details include tumor location, histopathology, age, gender, habit, the stage was noted in each case.

Isolation of exosome from primary and bone metastases patient serum

Exosomes were isolated with the miRCURY exosomes kit for Serum/Plasma (Qiagen Cat No. 76603) of patients with breast and lung patients with and without bone metastases. Briefly, samples were centrifuged at 3000g for 10 min to remove cells and cell debris. The 0.5 ml supernatant was mixed with 200 µl of precipitation buffer A and incubated at 4 °C for 60 min. After the incubation, the tubes were centrifuged at 1500g for 30 min and the exosome pellet was reconstituted in resuspension buffer and stored at −20 °C until further analysis [22].

Nanoparticle tracking analysis (NTA)

The isolated exosome’s size distribution and absolute quantification were estimated using Nanosight NS300 (NanoSight Ltd., Amesbury, UK) equipped with a 405 nm laser. A video of 60-sec duration was taken with a frame rate of 30 frames/sec, and particle movement was analyzed using NTA software (version 2.3; Nano-Sight Ltd.).

Flow cytometry analysis

Isolated exosome’s were labelled according to Els J van der Vlist., et al protocol. In brief, isolated exosome’s were incubated for 60 min in 10 µl of mouse anti-CD63-PE, CD-9, CD-81 antibodies. Before data acquisition, the samples were diluted with 480 µl of 0.22 µm pre-filtered PBS (final dilution 1:25). BD FACSCalibur flow cytometer was used for data acquisition and BD FACStation was used for calculation. Data were obtained using 2 individual apogee A50/Micro flow FCs equipped with 50mW 405-nm (violet), 488-nm (blue), and 638-nm (red) lasers [23].

Gene expression analysis by qRT-PCR

Total RNA was extracted using Aridia viral DNA and RNA extraction kit (Cat No./ID: AME0003) according to the manufacturer’s instructions. The integrity and concentration of the isolated RNA were determined using Bioanalyzer and 260/280 ratio for purity by Nanodrop (Epoch BioTek system). A total of 1.0 μg of RNA (260/280: 1.8–2.0) was reverse transcribed to cDNA using the high capacity cDNA synthesis kit (Qiagen; Cat no: 205411) as per the manufacturer's instructions. Real-Time PCR was performed in 20 μl volume that included 10 μl SYBR Green QPCR Master Mix (QuantiNova SYBR Green, Cat No./ID: 208052) containing 0.5 μl (200 nM) each of specific forward & reverse primer and 2 μl cDNA as template. 18 s rRNA was used as a housekeeping gene in each set of experiments. The list of all the primers used for the study (HSP90AA1, PTK2, SHC1, YWHAZ, MATR3, HSPD1, MMP9, VEGFA, IL3, NKTR, TUBGCP6, ACTG2, MYH11, CTTN, and SPP1) is shown in Table 1. Quantitative PCR using Sybr Green chemistry was carried out in QuantStudio® 5 ((QuantStudio® 5, Applied Biosystems, USA, CA)) in a 96-well reaction plate format with at the following thermal cycling conditions: 1 cycle of 5 min at 95 °C for the initial denaturation step and 40 cycles of 10 s at 95 °C for the denaturation step, 40 s at 56 °C / 60 °C for the annealing and extension step followed by melting curve detection for ensuring positive amplification of the target gene rather than non-specific products or primer dimmers. The fold change expression was evaluated using the ΔΔCT or 2ΔΔCT method. All experiments were performed in triplicate independently and the average CT value was calculated for the quantification of fold change analysis [24].

Table 1

Primer list.

Sr No	Gene name	Sequence	No of Base	Accession No
1	HSP90AA1	FP-5′-CCACTTGGCGGTCAAGCATT-3′	20	NM_005348
		RP-5′-AAGGAGCTCGTCTTGGGACAA-3′	21
2	PTK2	FP-5′-TATATGAGTCCAGAGAATCCAG-3′	22	NM_005607
		RP-5′-GCTTCACAATATGAGGATGGT-3′	21
3	SHC1	FP-5′-CACTTGGGAGCTACATTGCCTG-3′	22	NM_001130040
		RP-5′-GTGGTGGAGGTGGCATCTGTT-3′	21
4	YWHAZ	FP-5′-AGCCATTGCTGAACTTGATACA-3′	22	NM_145690
		RP-5′-AATTTTCCCCTCCTTCTCCTG-3′	21
5	MATR3	FP-5′-CAGCAGTCTACAAATCCAGCACC-3′	23	NM_018834
		RP-5′-CTGCATGTGTCTAGGTCCTTGC-3′	22
6	HSPD1	FP-5′-GCAAAGTTCCTCAGAAGTTGGT-3′	22	NM_002156
		RP-5′-GCAGCATCCAATAAAGCAGTT-3′	21
7	MMP9	FP-5′-GAGTGGCAGGGGGAAGATGC-3	20	NM_004994
		RP-5′-CCTCAGGGCACTGCAGGATG-3′	20
8	VEGFA	FP-5′-CTTGCCTTGCTGCTCTACC-3′	19	NM_003376
		RP-5′-CACACAGGATGGCTTGAAG-3′	19
9	ILF3	FP-5′-GTGTCCAATCACCAGTCCTG-3′	20	NM_012218
		RP-5′-GCTGAAGAAGTGGGAGTGTAGC-3′	22
10	NKTR	FP-5′-GCAAGCAGTTCAGAAGAGCCAAG-3′	23	NM_005385
		RP-5′-TCTCAGGCACTGGAGGAATCTC-3′	22
11	TUBGCP6	FP-5′-GGTGTTCAGAGACGCTTATGGC-3′	22	NM_020461
		RP-5′-CCACCTCTTTGGAGATGAGCAC-3′	22
12	ACTG2	FP-5′-CTGCCATGTACGTCGCCATTCA-3′	22	NM_001615
		RP-5′-GACATTGTGGGTGACGCCATCA-3′	22
13	MYH11	FP-5′-GTCCAGGAGATGAGGCAGAAAC-3′	22	NM_002474
		RP-5′-GTCTGCGTTCTCTTTCTCCAGC-3′	22
14	CTTN	FP-5′-TAATCCAATGAGGAATTTCCAG-3′	22	NM_005231
		RP-5′-TAGAGCCTGGTGCCTGGG-3′	18
15	18srRNA	FP-5′-GGAGTATGGTTGCAAAGCTGA-3′	21	GU198749
		RP-5′-ATCTGTCAATCCTGTCCGTGT-3′	20
16	SPP1	FP-5′-ACTCGTCTCAGGCCAGTTG-3′	19	NM_001040058
		RP-5′-CGTTGGACTTGGAAGG-3′	16

Hierarchical clustering

To study whether the expression profile in the primary cancer can exactly classify metastatic status, hierarchical unverified clustering was achieved using the hcluster method of R package “amap” and the plot was created using the heatmap.2 function of “gplots” package. Absolute Pearson and Pearson distances were used to calculate gene and sample distances respectively and gene linkages were done using the Ward algorithm. Inter-study normalization was completed with the Bioconductor package “inSilicoMerging” using an Empirical Bayes method [25].

Receiver operating characteristics curve analysis

The receiver operating characteristic (ROC) curve was generated to evaluate the predictive power such as accuracy and discriminating the specificity and sensitivity of each of the biomarkers. The optimal cut-point that yielded the maximum sensitivity and specificity was determined for each biomarker as cut off from ROC curves and area under the curve (AUC) was computed using MedCals (Belgium, Europe). The biomarker that has the largest area under the ROC curve was identified as having the strongest association with the presence of bone metastasis [13].

Results

DATA mining

The data mining of the different database documented a total 47 of series based on the search condition [Affymetrix n = 25, Agilent n = 10, Illumina n = 12]. After further filtration of data based on the inclusion and exclusion criteria, 18 data sets from Affymetrix, 1 from Agilent, and 4 from Illumina were included for further analysis in the study (Table 2). The chip and technology that mismatched with the analysis pipeline were excluded from the study. The total number of samples analyzed included 542 primary and 81 patients with bone metastases. The Affymetrix platform included a total of 380 primary and 63 metastases done with U133 plus 2.0: P = 241, M = 33, U133A: P = 116, M = 28, U95A: P = 23, M = 2 whereas the Agilent platform using 4X44K G4112F included 5 primary (P) and 5 samples with metastasis (M) samples. Similarly, the Illumina platform on HT-12 V 4.0 bead chip included 157 primary (P) and 13 samples with bone metastases. The representative site examined in the selected series were breast, prostate, colon, and renal cell carcinoma tumours metastasizing to the bone. The data set from each chip within each platform were analyzed individually using GeneSpring analysis software and as per the analytical pipeline to classify the concordant gene entities lists.

Table 2

Characteristic of individual studies retrieved from Gene Expression Omnibus for bone meta-analysis.

NO.	PUBLIC DATASET	ARRAY PLATFORM	SITE DETAIL	GENES VALIDATED	VALIDATION METHOD	PUBMED ID
1	GSE54323	HG-U133_Plus_2	Breast cancer metastasis to bone	–	–	25,888,067
2	GSE39494	Agilent-014850/4x44K G4112F	Breast cancer metastasis to bone	ABCC5	qPCR, IHC, WESTERN BLOT	23,174,366
3	GSE14017	HG-U133_Plus_2	Breast cancer metastasis to bone	–	–	19,573,813
4	GSE18549	HG-U133_Plus_2	Prostate cancer metastasis to bone	–	–	–
5	GSE32269	HG-U133A	Prostate cancer metastasis to bone	SOX9	qPCR, IHC, WESTERN BLOT	23,426,182
6	GSE68882	HG_U95A	Prostate cancer metastasis to bone	TAGLN, MSMB	qPCR, IHC	12,154,061
7	GSE101607	Illumina HumanHT-12 V4.0	Prostate cancer metastasis to bone	PSMB9, TAP1, HLA-A	qPCR, IHC	27,497,761
8	GSE101607	Illumina HumanHT-12 V4.0	Colon cancer metastasis to bone	–	–	27,497,761
9	GSE101607	Illumina HumanHT-12 V4.0	Renal cancer metastasis to bone	–	–	27,497,761

The data from the different series were analyzed as a single experiment in GeneSpring software and the PCA plot indicated the various groups showed different gene expression patterned. In breast cancer bone metastasis 12 series were analyzed in Affymetrix and Agilent platform that involved studies carried out on U133 plus 2.0 chips and 4X44K G4112F chip. In the analysis, few series were excluded since the samples were outliers during PCA and clustering. The PCA plot was created to analyze the behavior and clustering of samples of diverse subgroups and a total of 254 samples was incorporated in the final analysis (P = 241, M = 33). The statically noteworthy gene list having fold change value >2.0 and p-value <0.05 were considered for further analysis. In addition, breast cancer bone metastases samples included different breast cancer types such as basal, non-basal, invasive, lobular, ductal, invasive ductal, and BRCA1. In breast cancer bone metastases a total number of 20,680 down and 25,392 up genes were identified after the analysis of U133 plus 2.0 whereas 2387 down and 1837 up genes were identified from 4X44K G4112F. Similar study was carried out in prostate cancer bone metastases in that 8 series were analyzed in U133 plus 2.0 and U95A, a total number of 7350 down and 6456 up genes were recognized after the analysis on U133 plus 2.0 while 208 down and 329 up genes were identified using U95A chip. Parallel, a total number of 298 down and 4406 up genes were identified from Illumina platform HT-12V 4.0 bead chip in colon cancer bone metastases using 3 series. Comparative analysis in the 2 series of renal cancer to bone metastases using Illumina platform HT-12V 4.0 bead chip showed a total number of 10,568 down and 2278 up-regulated genes.

Meta-gene signature prediction

Intensive analysis of all the genes of the retrieved microarray studies was carried out and separate lists of differentially expressed genes according to the primary with bone metastases tumor were prepared. The meta-gene signature was identified by comparing primary tumors such as breast, prostate, colon, and renal cancers. The Venn diagram created using FunRich software to identify the meta-gene signature (Fig. 2A). The list of genes commonly identified between all primary bone metastasis tumors showed a total 985 gene list (p < 0.05 and FC > 2.0) were 478 down and 507 up-regulated significantly. From this a highly up and down-regulated metastatic 15 genes (HSP90AA1, PTK2, SHC1, YWHAZ, MATR3, HSPD1, MMP9, VEGFA, IL3, NKTR, TUBGCP6, ACTG2, MYH11, CTTN, and SPP1) panel was obtained based on which lead to the formation of bone metastases (Fig. 2B). The concordant list of 15 genes is provided in Table 3.

Fig. 2

(2A) Venn diagram. Venn diagram representing the overlapping genes. Each of the circles represents a dataset of bone metastases with various primary tumors breast cancer, colon cancer, prostate cancer and renal cancer. The numerals are the number of genes differentially expressed in the datasets represented by that area of overlap of the circles. (2B) Heat maps of identified 15 gene expression.

Table 3

The concordant list of top 15 genes description.

Entrez ID	Gene symbol	Description	Gene expression	Chromosome	Map location
3320	HSP90AA1	heat shock protein 90 kDa alpha (cytosolic), class A member 1	Up	14	14q32.33
5747	PTK2	Protein tyrosine kinase 2	Up	8	8q24.3
6464	SHC1	SHC (Src homology 2 domain containing) transforming protein 1	Up	1	1q21
7534	YWHAZ	Tyrosine monooxygenase / tryptophan 5-monooxygenase activation protein, zeta	Up	8	8q23.1
9782	MATR3	Matrin 3	Up	5	5q31.2
3329	HSPD1	heat shock 60 kDa protein 1 (chaperonin)	Up	2	2q33.1
4318	MMP9	matrix metallopeptidase 9	Up	20	20q13.12
6696	SPP1	secreted phosphoprotein 1	Up	4	4q22.1
7422	VEGFA	vascular endothelial growth factor A	Down	6	6p12
3609	ILF3	interleukin enhancer binding factor 3, 90 kDa	Down	19	19p13.2
4820	NKTR	natural killer cell triggering receptor	Down	3	3p22.1
85,378	TUBGCP6	tubulin, gamma complex associated protein 6	Down	22	22q13.31-q13.33
72	ACTG2	actin, gamma 2, smooth muscle, enteric	Down	2	2p13.1
4629	MYH11	myosin, heavy chain 11, smooth muscle	Down	16	16p13.11
2017	CTTN	Cortactin	Down	11	11q13

Functional annotation

Gene ontology delivers a collective descriptive background and functional annotation and classification for evaluate the gene sets data. GO groupings are structured into three groups: biological process, cellular component, and molecular function. We found GO terms for biological function significantly enriched in the metabolic process, cell adhesion and cell growth while, cellular component, the enriched GO term were cytoplasm, exosome and nucleus and for molecular function, the enriched GO term were structural molecular activity and the binding category significantly showed the major component with protein binding molecules (Fig. 3A).

Fig. 3

(3A) The top 15 enriched GO terms of differentially expressed genes. A. biological process B. molecular functions C. cellular component (3B) Molecular mechanisms connected with metastases progression. Expending the commercial pathway knowledge we established pathways that were developed or over-represented in the common metastatic signature. The figure represents the cross – talk amongst the various signal transduction pathways with main active genes such as HSP90, SPP1, VEGF, MMP-9, IL6, YWHAZ, PTK2 that direct developments such as cell proliferation, invasion, apoptosis and lead to the formation of bone metastasis. (3C) String protein interaction analysis. String output showing interaction of the common 15 genes that are specific for bone metastases.

To further assess the biological noteworthy pathway for the genes, we also done the KEGG pathway enrichment analysis. Hypergeometric test with P value < 0.05 was used as the standards for pathway detection. The most important pathway in our KEGG analysis was

mTOR signaling pathway, VEGF and VEGFR signaling pathway, and E cadherin and N cadherin signaling event. Osteopontin-mediated events, BMP signaling pathway, and p38 MAPK signaling pathway other significant pathways that are also dysregulated (Fig. 3B).

Protein-protein network and Hub-gene identification

Hub node significantly have a vital role in signaling mechanism. Therefore, the concordant list when further analyzed to STRING database to the investigation of network interaction between 15 genes. The major communication network consisted of two major interconnected groups with HSP90AA1 and PTK2 being the nodes of connection (Fig. 3C). Furthermore, NetworkAnalyst database delivered a protein network with 520 nodes, 637 edges, and 14 seeds. Twelve nodes with degrees above 10 were identified as hub genes, including HSP90AA1, YWHAZ, PTK2, CTTN, HSPD1, TUBGCP, MYH11, ACTG2, IL3, VEGFA as shown in Fig. 3C and table 4.

Table 4

The most 12 hub-gene identified based on PPT network analysis.

Hub-Gene	Degree	Betweness
HSP90AA1	238	93690.81
YWHAZ	85	31882.7
SHC1	85	30958.23
PTK2	68	24258.69
CTTN	39	13616.52
VEGFA	26	15830.24
HSPD1	21	6802.99
TUBGCP6	16	5648
MYH11	16	4747.28
MMP9	15	6688.74
ACTG2	12	2532.23
ILF3	10	4136.25

Characterization of the exosomes derived from breast and lung cancer bone metastatic cells

To explore the role of exosomes in bone metastases, exosomes were isolated from the serum of lung and breast patients with and without bone metastases, and their characteristic was confirmed by flow cytometer and nanoparticle tracking analysis using Nanosight based on the principle of particle size distribution. However, serum contains the heterogeneous population of extracellular vesicles comprising both round-shaped 30–100 nm diameter vesicles, consistent with exosomes. The exosomes isolated were approximately 32 nm in size which confirmed that these vesicles are exosomes (Fig. 4A). Furthermore, flow cytometry data also revealed that extracted exosome was express three specific exosomal marker CD63, CD81, and CD9 (Fig. 4B).

Fig. 4

Characterization of extracellular vesicles. (4A) Size and concentration evaluated by nanosight, indicates that sizes are compatible with exosomes. (4B) Exosomes markers (CD9, CD63, and CD81) were analyzed using Flow-cytometer. The data demonstrated that extracts were enriched with exosomal marker protein CD81 and CD63.

Differential gene expression analysis of exosomes: identifying a marker for liquid biopsy

In order to evaluate a gene expression pattern that can differentiate bone metastatic tumors from primary breast and lung tumors were identified as depicted in Fig. 5A. Quantitative gene expression patterns of HSP90AA1, PTK2, SHC1, YWHAZ, MATR3, HSPD1, MMP9, VEGFA, IL3, NKTR, TUBGCP6, ACTG2, MYH11, CTTN, and SPP1 were analyzed from exosomes mRNA of lung and breast cancer patients with and without bone metastases. Among these genes, in breast cancer patient HSP90AA1, IL3, VEGFA, and PTK2 were dramatically upregulated (≥2 fold differently expressed) in primary as well as bone metastases, whereas SHC1, YWHAZ, MATR3, and HSPD1 genes were downregulated (<2 fold differently expressed) in primary breast cancer as well as in bone metastases. On the other hand, the expression of NKTR, TUBGCP6, ACTG2, MYH11, MMP9, SPP1 and CTTN showed significant downregulation in primary breast and was upregulated in bone metastases (Fig. 5A).

Fig. 5

(5A) Gene expression. Primary breast and lung exosomes mRNA, advance stage breast and lung cancer with bone metastases, as well normal exosomes mRNA analysed with qRT-PCR. (5B) Unsupervised hierarchical clustering of primary breast (A) and lung (B) cancer with advanced stage bone metastatic samples. Clustering was based on 15 differentially expressed genes at a false discovery ratio level of 0.05. Breast and lung tumor identification looks at the top of the figure and each column represents gene expression of a single tumor. The colored bar specifies the variation in gene expression in target samples as compared to reference cells i.e., red, more expressed and cream, less expressed in target samples. Further, the black lines of the dendrogram stand for the support for each clustering. The metric performed was Euclidean distance, with complete linkage for distance between clusters. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In lung cancer patient HSP90, PTK2, MMP9, VEGFA, IL3, NKTR, TUBGCP6, ACTG2, and PTK2 were noteworthy upregulated in primary lung cancer with and without metastases. Curiously, it was seen that MYH11, SHC1, HSPD1, and CTTN genes were significantly upregulated in primary lung cancer compared to the bone metastases but contradictorily SPP1 was notable downregulated in primary lung cancer whereas upregulated in bone metastases. Furthermore, the expression of the MATR2 and YWHAZ genes were significantly downregulated in primary lung cancer as well as in bone metastases (Fig. 5A).

Further we performed the hierarchical cluster to checked whether the selected 15 genes would be beneficial in categorizing primary breast and lung cancer into sets that have diverse possible to develop bone metastases or not. This could only be predicted to occur if the gene expression profile related with metastases is already existing in a subset of cells in the primary tumor. The expression profile of the 15 genes was therefore used in hierarchical clustering to classify a group of 10 non-metastatic breast and lung tumor as depicted in Fig. 5B. The primary cancer breast and lung were clustered into two individual groups, based on their expression profile in primary and metastatic mRNA as highly correlating or not correlating with each other. We predicted that the tumours with a larger fold change in gene expression profile would have a worse diagnosis may be due to disease progression to metastases.

ROC analysis was carried out to determine the prediction of these identified genes, thus individual ROC was also studied with the sensitivity, specificity, and area under the curve are stated in table 5 (Fig. 6A & 6B). In detail, every 15 genes were able to classify between primary breast and lung cancer with and without bone metastases with an average sensitivity and specificity 90% and 100%, respectively. After the ROC curve analysis prediction of breast bone metastases and lung bone metastases could be formed with the 5 genes such as HAP90AA1, SPP1, IL3, SPP1 and PTK2 which shows the sensitivity and specificity results touched the 100%, respectively.

Table 5

ROC curve analysis of the individual genes. (A) Breast cancer with and without bone metastasis (B) Lung cancer with and without bone metastasis.

(A) Model Name	Associated criteria	Sensitivity	Specificity	Significance P (area = 0.5)	Youden Index J	AUC	95% CI
ACTG2	>0.4819	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
CTTN	>0.3103	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
HSP90AA1	>17.7448	100.00	90.00	<0.0001	0.9000	0.980	0.797 to 1.000
HSPD1	≤0.0039	80.00	70.00	0.5000	0.0920	0.710	0.467 to 0.888
MYH11	>0.9874	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
ILF3	>5.9969	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
MATR3	>0.1374	70.00	100.00	0.0815	0.7000	0.740	0.498 to 0.907
MMP9	>3.5146	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
NKTR	>0.6462	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
PTK2	>10.5421	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
SHC1	>0.0421	70.00	70.00	0.5417	0.4000	0.590	0.351 to 0.801
SPP1	>1.4261	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
TUBGCP6	>0.4268	100.00	100.00	<0.0001	1.0000	1.000	0.832 to 1.000
VEGFA	≤22.1618	90.00	100.00	<0.0001	0.9000	0.970	0.781 to 1.000
YWHAZ	≤0.0304	90.00	60.00	0.1036	0.5000	0.710	0.467 to 0.888

Fig. 6

ROC curve. Receiver operating characteristic curve analysis of the individual 15 genes in patients with primary breast (6A) and lung (6B) cancer with and without bone metastases.

Discussion

Metastasis is the single major cause of cancer related mortality [26]. The pathophysiology of metastasis progression includes the complex interplay between tumour cells and its microenvironment [27]. In the last few decades molecular landscaping of various metastatic cancers have found role of signal transduction molecules in cancer invasion, homing, and progression eventually leading to metastatic spread from primary site of origin. Therefore, there is dire need of emerging technologies and platforms to predict risk of developing site specific metastasis using liquid biopsy approach to avoid repeated tissue biopsy [28]. Bone metastasis is one of the most common metastatic site in solid tumours originating from breast, lung, prostate, colon, and kidney [3]. The biggest challenge in treating metastatic disease is due to its systemic spread and lack of appropriate treatment approaches [29]. Solitary metastatic site can be tackled through metastatectomy to show some survival improvement, while patients with multiple lesions are not always suitable for surgery [30]. This prompted us to establish a multi-gene expression panel for the early prediction of bone metastasis in patients with breast and lung cancer using a liquid biopsy approach. An early indicator of bone metastasis may provide survival advantage to patients and thereby provide them long disease free survival with timely intervention.

In the post human genome era high throughput transcriptome and miRNome profiling using NGS or microarray is a routine practice for establishing a diagnostic or prognostic marker(s). This approach is cost intensive and require access or collaboration with good clinical cancer care setup. Many such data sets are uploaded in publicly available databases for further research by academic institutions. Meta-analysis of such publicly available datasets can be utilized judiciously to derive some meaningful conclusion with clearly defined hypothesis and logically designed methodology. In the present study meta-analysis was carried out on microarray gene expression dataset derived from patients with and without bone metastasis for identification of possible biomarkers. All qualifying datasets containing microarray gene expression studies from various primary tumours e.g. colon, prostate, breast and renal were included for meta-analysis using Genespring software. The intention was to identify common gene expressed consistently in all bone metastasis patients irrespective of primary site of cancer. We found 15 genes which were most commonly over expressed in patients with bone metastasis as compared to no metastatic tumours. We speculated that since exosomes are routinely shaded off from primary tumours and enters into circulation if we can identify any of these 15 shortlisted gene transcript from patient’s blood sample it will be useful to establish a metastasis predictor. We therefore conducted a prospective study of realtime PCR based gene expression using exosomal mRNA from 10 lung/breast cancer patients harboring bone metastasis as compared to equal number of non-metastatic patients as control after normalized against healthy male and female as control.

Several studies based on meta-analysis of gene expression datasets derived from various cancer of different sites e.g. breast, osteosarcoma, prostate and pancreatic with or without bone metastasis have been carried out to establish biomarkers for diagnosis and/or prognosis [31], [32], [33], [34]. Meta-analysis is an efficient approach that enables the researchers to analyse data derived from various studies using similar or different platforms to identify probable markers of clinical significance [19]. These hypothesized biomarkers derived from global profiling data can later be validated prospectively using clinical samples for establishing a functionally applicable biomarker in specific malignancies. In the present study validation of multi-gene expression panel derived from meta-analysis result was carried out on a relatively small sample population consisting of only breast and lung cancer with and without bone metastasis.

Earlier study have reported that exosomes, small extracellular membrane derived vesicles released by the cells, are the key facilitators of tumor metastasis and early predictor of tumor invasion and progression [35]. Additionally, exosomes are secreted into various body fluids which can deliver tumour content to various target organs. [36]. Exosome can be isolated from various body fluids to reveal a patient's clinical direction and tumor signature [37]. Furthermore, Exosome based liquid biopsies have various benefits such as its minimally invasive nature, the main advantage of multiple samples can be collected at different time points during treatment. Analysis of exosomal cargo, therefore, has a significant advantage in establishing diagnostic, prognostic, and treatment monitoring markers [38]. We report here that all the 15 genes identified by meta-analysis are not significantly overexpressed in the exosomal cargo showing possibility of differential packaging of tumour-derived transcript during exosome biogenesis. However, a common gene signature consisting of HSP90AA1, SPP1, IL3, VEGFA, PTK2, and YWHAZ were found in both breast and lung cancer patients with bone metastasis.

HSP90, the 'cancer chaperone', plays a crucial role in cell proliferation and differentiation and is known to be associated with enhanced aggressiveness of the tumor and metastatic activity [39]. Further, IL3 contributes to carcinogenesis and metastasis by suppressing osteoblast differentiation and bone deposition [40], [41]. VEGF is also known to be a major player in the regulation of angiogenesis which stimulates cell migration and proliferation and has been studied widely in bone metastasis [42]. Additionally, PTK2 also known as focal adhesion kinase stimulate metastasis by controlling the process of cancer cell motility, invasion, and matrix metalloproteinase surface expression. Similarly, Mitra et al observed that FAK activity significantly increases MMP9 expression and impulsive breast cancer metastasis in genetically identical and orthotopic mouse models [43]. SPP1 also known as Osteopontin (OPN) shows multifunction characteristic in cancer progression and osteogenic differentiation [44], [45], [46]. Collectively, serum exosomal markers HSP90AA1, SPP1, IL3, VEGFA, and PTK2 found in the present study might be useful in detecting the early spread of bone metastasis leading to better clinical outcomes. Besides, these closely interacting gene interaction networks need to be further explored to understand their molecular mechanism and clinical relevance.

In this study, we calculated the prediction accuracy of the gene signature generated by ROC curve analysis to be more than 90%. Likewise, the clustering data also revealed a noteworthy correlation (p < 0.0001) between the gene expression patterns of the various genes and metastatic tumors. After further large scale validation these five gene panel (HSP90AA1, SPP1, IL3, VEGFA, and PTK2) can be established as a metastatic predictor panel based on realtime based gene expression in a cost effective and less invasive way.

In conclusion, meta-analysis of publicly available datasets using suitable statistical tool offers a cost effective way of generating hypothesis which can later be validated with more sensitive realtime based gene expression profiling. Using this approach we could identify a five gene panel with more than 90% accuracy for predicting bone metastasis in patients with breast and lung cancers using liquid biopsy.

Funding: The study was supported by Indian Council of Medical Research (Ortho/2018/NCD-I). The funding agency had no role in study design, data collection and analysis, preparation of the manuscript or decision to publish.

CRediT authorship contribution statement

Kinjal P. Bhadresha: Methodology, Software, Validation, Visualization, Writing - original draft. Maulikkumar Patel: Software. Nayan K. Jain: Supervision, Formal analysis, Investigation, Writing - review & editing. Rakesh M. Rawal: Conceptualization, Formal analysis, Funding acquisition, Investigation, Project administration, Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.