Bone metastasis unique plasma exosomal microRNA signature in non-small cell lung cancer

Background 20–40% of lung cancer patients develop bone metastasis (BM) with significantly decreased overall survival. Currently, BM is mainly diagnosed by CT scan or MRI when symptom develops. Novel biomarkers with higher prediction value of BM are needed. Methods Prospective analysis was undertaken on non-small cell lung cancer (NSCLC) patients with (BM+) and without BM (BM-). Plasma exosomal RNA was isolated and sequenced from peripheral blood of patients. Differential expression analysis and weighted gene co-expression networks analysis (WGCNA) of mi-RNA sequencing data were performed between two groups. Results Hierarchical clustering based on the total miRNA profile can clearly separate cancer patients and healthy individuals (H), but not patients BM + or BM-. WGCNA identified three consensus clusters (A, B, C) of highly correlated miRNAs, among which cluster B (144 miRNAs) showed significantly differential expression in lung cancer patients, especially in BM + group. Three differentially expressed miRNAs between BM + and BM- patients within cluster B were identified as miR-574-5p, a suppressor of Wnt/β-catenin pathway, was down-regulated, while miR-328-3p and miR-423-3p, two activators of the same pathway, were up-regulated in BM + patients. Pathway analysis of cluster B miRNAs revealed enrichment in metabolic pathways that may involve in preconditioning of the metastatic niche. Cluster A miRNAs (n = 49) also showed trend of upregulation in BM + patients. Interestingly, pathway analysis indicated that 43 of them are associated with chromosome14, which has been suggested to promote EMT and bone metastasis. These data indicated that a cluster of mi-RNAs showed significantly differential expression in BM + group, including miR-574-5p, miR-328-3p and miR-423-3p.

easy detection of bone metastasis.
MicroRNAs (miRNAs) which are small 18 to 24 nucleotides non-coding RNAs target the 3' untranslated region of messenger RNAs (mRNAs) and regulate gene expression, resulting in mRNA cleavage or suppression of protein translation [11]. MiRNAs often located in the fragile regions of the chromosome which have a high frequency of deletions, rearrangements and amplification and have been implicated in malignancy [12]. At present, miRNAs are being studied as diagnostic and prognostic biomarkers for many cancers including NSCLC [13,14], breast cancer [15], prostate cancer [16], hepatocellular cancer [17]. Many observations also strongly implicated the possibility of developing miRNAs as non-invasive circulating biomarker for the early detection of solid cancers [18,19]. MiRNAs are ideal for clinical detection for many reasons: miRNAs are very stable in body fluid. The structure of exosomes released by cells into the blood further protects miRNAs from degradation. Second, next generation sequencing has provided a highly robust and accuracy system for detection of miRNAs in body fluid on a genome-wide scale [20]. In this study, we focus on the miRNAs in plasma-derive exosomes to investigate the potential of using miRNAs as biomarkers for early detection of bone metastasis in NSCLC patients.

Methods
Patients and clinical sample collection A total of thirty EGFR/ALK positive NSCLC patients were enrolled in this study including sixteen phase IV patients with bone metastasis and fourteen phase IV patients without bone metastasis. Among them, twenty-five patients were diagnosed with adenocarcinoma (ADC). The median age at diagnosis is 55, ranging from 35 to 78. Peripheral blood was collected from each patient on a regular basis from routine clinical care, and plasma sample was prepared within two hours of blood drawn and then stored at -80 o C.
Plasma exosome isolation 1 mL of plasma sample was centrifuged at 10,000 g for 30 min at 4 o C to remove any cell debris. The collected supernatant was then subjected for ultra-high speed centrifugation at 150,000 g for 70 min at 4 o C. Pellet containing exosome was resuspended in 200 µL PBS for downstream applications.
Exosomal RNA isolation and small RNA sequencing Total RNA including miRNA was extracted from plasma-derived exosome using miRNeasy Serum/Plasma Kit (QIAGEN) following manufacturer's instructions. The quantification and size distribution of the extraction were analyzed by Qubit 4.0 and Agilent Bioanalyzer 2100 (Agilent), respectively. Quantified miRNA was subjected for sequencing library preparation using NEBNext® Small RNA Library Prep Set for Illumina® (NEB Biolabs) following manufacturer's instructions. Briefly, isolated miRNA was subjected for 3' and 5' adaptor ligation, followed by 17 cycles of PCR amplification. PCR products from library preparation were subjected for gel electrophoresis on 6% Novex® TBE PAGE gel (Thermo Fisher Scientific) and DNA fragments between 140-150 bp were recovered from the gel. Purified miRNA cDNA library was quantified by Qubit 4.0 and the size distribution was analyzed on Agilent Bioanalyzer 2100. miRNA cDNA libraries from different plasma samples were pooled and sequenced on Illumina HiSeq4000 platform.
miRNA-seq data analysis miRNA identification and reads counting in each miRNA were performed using miRDeep2 [21]. After trimming the 3' adaptor sequence, all sequences ranging in length from 18-26 nt were recorded in a non-redundant file along with reads count. To identify known miRNAs, the miRNA tags were aligned against miRNA precursor sequences reported in the miRNA database 'miRBase' (release 21) using the 'quantifier.pl' script within miRDeep2. Differential expression (DE) analysis of miRNA sequence data was performed with the Bioconductor package edgeR [22]. miRNAs with read counts per million mapped reads (CPM) > = 2 in at least 20% of all samples were identified as expressed miRNAs. DE between different groups was evaluated by fitting a negative binomial generalized linear model and then adjusting the P-value for multiple testing using the Benjamini-Hochberg correction with a false discovery rate of 0.1 and a minimum log 2 (CPM) of 4.

Weight Co-Expression Networks
Weight co-expression network of miRNAs was performed in accordance to the protocol of WGCNA package in the R language [23]. MiRNAs were aggregated into modules by hierarchical clustering and refined by the dynamic tree cut algorithm. Thereafter, module eigenvalues were calculated. The  Table 1). Using exact test, 91 miRNAs were identified as differentially expressed miRNAs (DE miRNAs) between H and C including 71 up-regulated miRNAs and 20 down-regulated miRNAs in patient group (Supplementary Table 2). Based on the DE exosomal miRNA profile, samples of healthy individuals (H) could be separated from samples of NSCLC patients (C) using supervised hierarchy clustering (Supplementary Figure S1). All these data proved that NSCLC patients exhibited a unique miRNA profile compared to healthy population.

Detection of co-expression clusters in exosomal miRNAs
To characterize the correlation pattern and predict the function of miRNAs, we applied weighed gene co-expression network analysis (WGCNA) to the 1287 miRNAs detected in the plasma exosomes.
WGCNA is widely used in genomic data analysis which can detect clusters of highly correlated genes based on pairwise correlations [23]. As shown in Supplementary Figure S2, we identified three clusters of co-expressed miRNAs which were represented by different color codes (Brown: cluster A; Turquoise: cluster B; Blue: cluster C). Cluster B had 144 miRNAs which is the largest cluster among three. 49 miRNAs were assorted into Cluster A while 95 miRNAs were in Cluster C (Supplementary Table 3).
Furthermore, BM + group showed significant up-regulation in cluster B eigengene value compared to healthy population and BM-group which suggested miRNAs in cluster B might related to the initial of bone metastasis (Fig. 1B). The 144 miRNAs from cluster B basically differentiated BM + group from BM-group using unsupervised clustering which also shed light on the function of the miRNA cluster in bone metastasis ( Fig. 2A). We performed pathway analysis of the 144 miRNAs in cluster B using miRNA enrichment analysis and annotation (MiEAA) [24] and top 20 miRNA pathways were shown in Table 1. Interestingly, cluster B was enriched in metabolism processes such as pyruvate metabolism (p = 0.010575), glycolysis and gluconeogenesis (p = 0.014983), purine metabolism (p = 0.014983), propanote metabolism (p = 0.014983) and pyrimidine metabolism (p = 0.027086) (Supplementary Table 4).
In cluster A, BM + group exhibited a trend of increase in eigengene value compared to healthy population and BM-group (Fig. 1A). miRNA enrichment analysis was done with 49 miRNAs in cluster A, however, not many related pathway was found. Nevertheless, we discovered that 43 out of 49 cluster A miRNAs were expressed by chromosome 14 (Supplementary Table 5).
In cluster C, the eigengene values were comparable within healthy population group, BM + and BMgroup (Fig. 1C) which indicated that this cluster was not related to bone metastasis. The cluster A miRNAs failed to differ BM + from BM-with unsupervised clustering in Fig. 2B also suggested that cluster C miRNAs were not involve in the bone metastasis. Related pathways of cluster C were shown in Supplementary Table 6. miRNAs in cluster C seemed to related to all kinds of signaling pathways including cell cycle (p = 0.000322), proteasome and lysosome (p = 0.000322), p53 signaling pathway (p = 0.00037), insulin signaling pathway (p = 0.000322) and Ras pathway(p = 0.00037).
These DE miRNAs might regulate cancer metastasis through Wnt/β-catenin signaling pathway and might be potential biomarker for bone metastasis in NSCLC.

Discussion
In this retrospective study, we analyzed the plasma-derived exosomal miRNAs from stage IV NSCLC patients with or without bone metastasis. A total of 1287 miRNAs were identified and assorted into three major clusters A, B, C by WGCNA. Cluster A had a trend of increase in BM + group compared to BM-and majority of cluster A miRNAs was transcripted by chromosome14. Cluster B showed significant difference between BM + and BM-and the three DE miRNAs (hsa-miR-574-5p, hsa-miR-328-3p and hsa-miR-423-3p) identified between BM + and BM-belonged to cluster B as well. Pathway analysis revealed that cluster B miRNAs were mostly related to metabolism processes. Hsa-miR-574-5p, hsa-miR-328-3p and hsa-miR-423-3p were reported actively involved in the Wnt/β-catenin signaling pathway and might be candidates as biomarkers for bone metastasis in NSCLC patients.
Cluster C showed no difference among healthy population, BM + and BM-, thus was not relate to bone metastasis in lung cancer.
The metabolic processes we found associated with cluster B miRNAs were highly associated with tumor metastasis. In general, metabolism changes often accompany tumor progression and were known to be associated with establishment of metastatic niches which are preconditioned for the arrival of metastatic disseminated cancer cells [31]. On the other hand, cancer cells undergo metabolic alteration and acquire metastatic traits to adapt to multiple environments. Upregulation of glycolysis [32], pyruvate kinase [33], pyrimidine phosphorylase [34] was observed in many cancers and played important role in tumor metastasis and aggressiveness. Recently, purine signaling pathway has been reported to contribute to the bone marrow metastasis in neuroblastoma [35]. All these evidence supported our finding that cluster B was metastatically relevant.
Wnt/β-catenin signaling pathway plays an important role in the epithelial-to-mesenchymal transition (EMT) and contribute to cancer progression and metastasis in different types of malignancies. MiRNAs are the major regulators of Wnt/β-catenin signaling pathway which make them ideal for therapeutic targets against metastatic tumor [36]. Here, we reported three differential expressed miRNAs which might be signatures for bone metastasis in NSCLC. Intriguingly, hsa-miRNA-574-5p, a suppressor of Wnt/β-catenin signaling pathway, has been reported to promote metastasis of NSCLC [37]. As for the two activators of Wnt/β-catenin signaling pathway, hsa-miR-328 was implicated in the high glucoseinduced EMT [38] and hsa-miR-423-3p was considered as potential biomarker for lung cancer diagnosis [39]. These three miRNAs might act together through Wnt/β-catenin pathway to promote EMT and could be unique makers for bone metastasis in lung cancer. The detection of these miRNA biomarkers in plasma exosomes has the potential to become a specific, sensitive and non-invasive method for monitoring bone metastasis in clinical setting for it only require blood sample which is easy to obtain. In future, we can validate these markers by monitoring miRNA levels in NSCLC patients before and after they develop bone metastasis. More study on the function of these miRNAs need to be carried out to further address their role as biomarkers or therapeutic targets for advanced NSCLC patients.
Overall, the 43 miRNAs from cluster A which were transcripted by chromosome 14 showed increasing expression in BM + group. One possibility of this increase in BM + might due to amplification of chromosome 14 in stage IV patients with bone metastasis. However, ctDNA or tumor samples from patients were not available for copy number variation detection to confirm this hypothesis. So far, chromosome 14 amplification was not reported in NSCLC or other cancers. But there was study on chromosome 14 allelic loss which is common in nasopharyngeal carcinoma and essential tumor suppressor gene loss in tumorigenesis [40].
To sum up, by comparing the plasma-derived exosomal miRNA features of NSCLC patients with or without bone metastasis, we identified a cluster of bone-metastasis related miRNAs and three DE miRNAs which might be applied to prediction of bone metastasis in NSCLC patients in future.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. QingZhouSupplementarytable.xlsx QingZhouSupplementaryInformation.docx