New insight into long non-coding RNAs associated with bone metastasis of breast cancer based on an integrated analysis

Bone is the most common site of metastatic breast cancer, and it is a leading cause of breast cancer-related death. This study aimed to explore bone metastasis-related long non-coding RNAs (lncRNAs) in breast cancer. Four mRNA datasets and two lncRNA datasets of bone metastasis, lung metastasis and liver metastasis of breast cancer were downloaded from Gene Expression Omnibus (GEO) database. The differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) in group of bone metastasis vs lung metastasis and bone metastasis vs liver metastasis, as well as the overlap of the two groups, were identified. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and protein–protein interaction (PPI) network construction of DEmRNAs were conducted. The cis nearby-targeted DEmRNAs of DElncRNAs were obtained. Quantitative real-time polymerase chain reactions (qRT-PCR) was used to detect the expression levels of selected DEmRNAs and DElncRNAs. LOC641518-lymphoid enhancer-binding factor 1 (LEF1) pair was selected to verify its role in migration and invasion capability of breast cancer cells by wounding healing assay and transwell invasion assay. A total of 237 DEmRNAs were obtained in bone metastasis compared with both lung metastasis and liver metastasis. A total of three DElncRNAs in bone metastasis compared with both lung metastasis and liver metastasis were obtained. A total of seven DElncRNA-nearby-targeted DEmRNA pairs and 15 DElncRNA-nearby-targeted DEmRNA pairs in group of bone metastasis vs lung metastasis and bone metastasis vs liver metastasis, were detected, respectively. Four cis LncRNA-mRNA interaction pairs were identified, which are LOC641518-LEF1, FLJ35024-Very Low Density Lipoprotein Receptor (VLDLR), LOC285972-Retinoic Acid Receptor Responder 2 (RARRES2) and LOC254896-TNF receptor superfamily member 10c (TNFRSF10C). qRT-PCR using clinical samples from our hospital confirms the bioinformatics prediction. siRNA knocking down LOC641518 down-regulates LEF1 mRNA expression, and reduces the migration and invasion capability of breast cancer cells. We concluded that four LncRNA-mRNA pairs, including LOC641518-LEF1, may play a central role in breast cancer bone metastasis.


Introduction
Globally, breast cancer is the most common malignancy in women [1]. Bone is the most common site of metastatic breast cancer, accounting for nearly three-quarters of patients with metastatic breast cancer [2]. Bone metastasis is often accompanied by bone pain and skeletal-related events, leading to increased mortality and seriously affecting the quality of life of patients [3]. The clinical complications of breast cancer bone metastasis bring heavy burden to the individual and society [4]. Bone metastases secondary to breast cancer are incurable, and the molecular mechanisms are complex and involved in multiple process [4].
Over the past few decades, a great deal of information about lncRNAs has been generated, which facilitated the studies on diverse cancer etiology, including breast cancer [5]. Increasing evidence demonstrated that lncRNAs participate in many cellular process, such as cell fate decision, immune response, cancer cells proliferation and metastasis [6]. Gooding et al. indicated that the lncRNA BORG drives breast cancer metastasis to lung and disease recurrence [7]. Metastasis associated lung adenocarcinoma transcript 1 (MALAT1), previously described as a metastasis-promoting lncRNA, was reported to suppress breast cancer metastasis [8]. HOXA11-AS was demonstrated to promote breast cancer invasion and lymph node metastasis through affecting epithelial-mesenchymal transition [9]. However, there are few studies on bone metastasis-related lncRNAs in breast cancer. Li et al. reported that lncRNA MAYA (MST1/2-Antagonizing for YAP Activation) promoted breast cancer bone metastasis by activating ROR1-HER3 and Hippo-YAP pathways, which indicated important role of lncRNA in breast cancer bone metastasis [10]. Currently, there is also lack of systemic LncRNA-mRNA network analysis based on human database, so it is necessary to explore important lncRNAs associated with breast cancer bone metastasis.
In this present study, a comprehensive integrated analysis of lncRNAs and mRNAs expression profiles of bone metastasis, lung metastasis and liver metastasis of breast cancer downloaded from GEO database were performed. The differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) were acquired. The cis nearbytargeted DEmRNAs of DElncRNAs were obtained as well. Finally, we established in vitro cell based assay to investigate and confirm the role of novel LncRNA-mRNA axis in tumor metastasis. This study aimed to provide additional knowledge on the molecular mechanisms of breast cancer bone metastasis, especially the role of probone metastasis LncRNAs.

Dataset collection
To acquire the mRNA and lncRNA expression profiles of bone metastasis, lung metastasis and liver metastasis of breast cancer tissues, datasets in GEO database with the following criteria were retrieved: datasets should be whole-genome mRNA/lncRNA expression profile by array; these data were derived from bone metastasis, lung metastasis and liver metastasis of breast cancer tissues; datasets were normalized or original. After screening, four mRNA datasets and two lncRNA datasets were enrolled in this study and shown in Table 1.

Differential expression analysis
MetaMA [11] and limma [12] package in R were utilized to acquire the DEmRNAs and DElncRNAs. P-values were calculated using a significance threshold of p-value < 0.05. With R package "pheatmap", hierarchical clustering analysis of DEmRNAs and DElncRNAs were performed.

GO and pathway analysis of DEmRNAs
GeneCodis3 (http:// genec odis. cnb. csic. es/ analy sis) was applied to perform GO and KEGG pathway enrichment analysis of DEmRNAs. The threshold was set as false discover rate (FDR) < 0.05.

Cis nearby-targeted DEmRNAs of the DElncRNAs
DEmRNAs transcribed within a 100-kb window upstream or downstream of DElncRNAs were searched, which were defined as cis nearby-targeted DEmRNAs of DElncRNAs, to obtain the targeted DEmRNAs of DEl-ncRNAs with cis-regulatory effects.

Quantitative real-time polymerase chain reaction (qRT-PCR) validation
Twenty-five patients with bone metastasis of breast cancer, 25 patients with lung metastasis of breast cancer, and 25 patients with liver metastasis of breast cancer were incorporated in our study. Meanwhile, 25 primary breast cancer patients without metastasis were enrolled, and their adjacent non-cancer breast tissues was considered as normal tissue to be used as control in the present study. Their tumor tissues from metastatic sites were used for RNA isolation. Ethical approval was obtained from the ethics committee of the First People's Hospital of Chengdu and informed written consent was obtained from all of subjects. All the patients in the current study were same as described in the previous publication [13]. Total RNA of the tumor tissues from metastatic sites was extracted using the RNA liquid Reagent (TIAN-GEN Bio, Beijing, China) according to the manufacturer's protocols. One microgram RNA was applied to synthesize cDNA by Fast Quant Reverse Transcriptase (TIAN-GEN Bio). Then real-time PCR was performed in an ABI 7300 Real-time PCR system with SYBR R Green PCR MasterMix (TIANGEN). All reactions were carried out in triplicate. GAPDH was used for the internal reference. Relative gene expression was analyzed by fold change method. And all the primer sequence was shown in Additional file 1: Table S1.

Cell culture and transfection with plasmid and siRNA
Human breast cancer cell line MCF-7 and MDA-MB-231, and non-carcinoma human breast epithelial cell line MCF-10A were purchased from ATCC (Manassas, VA, USA) and were cultured in RPMI 1640 medium (Invitrogen, CA, USA) supplemented with 10% fetal bovine serum (Invitrogen) and incubated at 37 °C and 5% CO 2 . Commercial RNA interfering sets targeting LOC641518 were purchased from Ribobio technology (Shenzhen, Guangdong, China), and scramble siRNA was also transfected into breast cancer cells to regard it as control. Exogenous human full-length LEF1 plasmid vectors were purchased from Sino biological Inc. (Beijing, China); and empty vector was also transfected as negative control. Lipofectamine 3000 (Invitrogen) was used for all transfection studies by following the manufacturer's protocol. Transfection efficiency was examined using qRT-PCR. Briefly, total RNA was extracted from cells using Trizol (Invitrogen, CA, USA) and further purified by RNAeasy mini kit plus DNase I treatment (Qiagen, Germany). The relative mRNA level for each gene was quantified by real-time RT-PCR with SYBR Green (Applied Biosystems, CA, USA), using GAPDH as a control.

Wound healing assay
The migration ability of cells was measured by woundhealing assay. Briefly, 2 × 10 4 cells were inoculated in 6-cm tissue culture dishes, cultured overnight, scratch was performed when cell growth fusion reached to 90%. Migration images were captured 24 and 48 h after scratching, the migration rate (width) of cells were calculated.

Transwell migration
Corning ™ BioCoat ™ Matrigel ™ Invasion Chamber with Corning ™ Matrigel Matrix (Thermo Fisher Scientific, Waltham, MA, USA) was used for cell invasion assay. One hundred microliters of cells with concentration of 5 × 10 5 /mL were added in the upper chamber with (8 μm pore), 600 μL medium containing 10% serum was added in the lower chamber. After cultured for 6 h, medium in the chamber were removed and unmigrated cells were swabbed. Cells was fixed by 4% polyformaldehyde for 10 min, stained with crystal violet for another 10 min. The filter membrane was photographed under the inverted microscope (200×) after sealed with the neutral gum. Cells were counted by Image Pro Plus Version 6, three wells in each group and five vision fields of each well were randomly selected to calculate the average number of cells.

Identification of DEmRNAs in tumor tissues from bone metastatic sites compared with other metastatic sites of breast cancer
Compared with lung metastatic sites, a total of 1280 DEmRNAs (645 up-and 635 down-regulated mRNAs) were detected in tumor tissues from bone metastatic sites of breast cancer. Compared with liver metastatic sites, a total of 1482 DEmRNAs (851 up-and 631 down-regulated mRNAs) were detected in tumor tissues from bone metastatic sites of breast cancer. Hierarchical clustering analysis of top 100 up-and down-regulated DEmRNAs was exhibited in Fig. 1A, B, respectively. Top 10 up-and down-regulated DEmRNAs were displayed in Table 2. After overlapped these 1280 DEmRNAs and 1482 DEm-RNAs, a total of 237 DEmRNAs (149 up-and 88 downregulated mRNAs), which were differentially expressed in bone metastasis compared with both lung metastasis and liver metastasis, were obtained ( Fig. 1C, D, Additional file 2: Table S2, Additional file 3: Table S3).

Identification of DElncRNAs in tumor tissues from bone metastatic sites compared with other metastatic sites of breast cancer
Compared with lung metastatic sites, a total of 71 DEl-ncRNAs (38 up-and 33 down-regulated lncRNAs) were detected in tumor tissues from bone metastatic sites of breast cancer (Additional file 5: Table S5). Compared with liver metastatic sites, a total of 91 DElncRNAs (40 up-and 51 down-regulated lncRNAs) were detected in tumor tissues from bone metastatic sites of breast cancer (Additional file 6: Table S6). Hierarchical clustering analysis of DElncRNAs was exhibited in Fig. 5A, B,   respectively. Top 10 up-and down-regulated DElncRNAs were displayed in Table 2. After overlapped these 71 DEl-ncRNAs and 91 DElncRNAs, a total of three DElncRNAs (two up-and one down-regulated lncRNAs), which were differentially expressed in bone metastasis compared with both lung metastasis and liver metastasis, were obtained (Fig. 5C, D, Table 3).

Correlation study of DElncRNAs and DEmRNAs
The correlation study between DELncRNA and DEm-RNA expression was analyzed by Person's test, the correlation coefficient and p-value were shown in Table 4 (bone metastasis vs lung metastasis and bone metastasis vs liver metastasis), only pairs whose P value was less than 0.05 were shown.

Cis nearby-targeted DEmRNAs of the DElncRNAs
A total of seven DElncRNA-nearby-targeted DEmRNA pairs in group of bone metastasis vs lung metastasis, involving in seven DElncRNAs and seven DEmRNAs, were detected (Table 5). A total of 15 DElncRNA-nearbytargeted DEmRNA pairs in group of bone metastasis vs liver metastasis, involving in 12 DElncRNAs and 15 DEmRNAs, were detected (Table 5). Combined with correlation study, total 8 pairs were selected with significant correlation, and four pairs had correlation coefficient higher than 0.7, which were LOC641518-LEF1, FLJ35024-VLDLR, LOC285972-RARRES2 and LOC254896-TNFRSF10C respectively.

qRT-PCR validation
These four pairs whose correlation coefficient was above 0.7 were selected for verification by qRT-PCR using clinical samples from our hospital. As shown in Fig. 6A, results showed that the relative expression of LEF1, VLDLR, RARRES2 and TNFRSF10C were significantly up-regulated compared to adjacent non-tumor tissues and other metastatic sites (P < 0.05). As shown in Fig. 6B, RNA levels of LncRNA LOC641518 and FLJ35024 showed same trend with bioinformatics data. Although RNA levels of LOC285972 and LOC254896 in the bone metastasis were also higher than adjacent non-tumor tissues, they were not different compared with lung and liver metastasis. All the mRNA and LncRNA levels of tested genes from adjacent non-tumor breast cancer was significantly lower than metastasis tumor tissues (Fig. 6).
After that, we used the RNA relative levels of each pair LncRNA-mRNA from metastatic bone tumor tissues to perform a correlation analysis, which was shown in Fig. 6C. Results demonstrated all the four LncRNA pairs had significantly positive correlation, although the coefficient is not high as bioinformatics prediction. The validation test by qRCR proved that most of bioinformatics predictions were reliable.

LOC641518 promotes breast cancer metastasis via activation LEF1
Due to LEF1 has been long considered as a metastasis mediator [14,15], in the current study, the Loc641518-LEF1 axis was selected to investigate whether Loc641518 regulates LEF1 expression to promote breast cancer metastasis by in vitro study, which can help to confirm the speculation from bioinformatics data mining. The RNA levels of LOC641518 and LEF1 in non-carcinoma human breast epithelial cell line MCF-10A and two breast cancer cell lines, MCF-7 and MDA-MB-231, were examined by qRT-PCR (Fig. 7A, B), which could find that RNA levels of LOC641518 and LEF1 is higher in cancer cells than non-carcinoma cells, and highest in MDA-MB-231.
siRNA was used to knock down the expression of LOC641518 in both MCF7 and MDA-MB-231 cells. As shown in Fig. 7C, the expression of LOC641518 was significantly reduced by siRNA compared with scramble siRNA control. Meanwhile, the LEF1 mRNA level was also reduced in siRNA-LOC641518 groups compared with scramble groups (Fig. 7D). The wound healing assay and migration assay both demonstrated that knocking down LOC641518 significantly reduced capability of migration and invasion (Fig. 7E, F), and exogenous overexpression of LEF1 could reverse these phenomena.

Discussion
Bone is a primary site for breast cancer metastases, and patients often suffer from osteolytic bone metastases, which is incurable [16]. Previously, we conducted integrated bioinformatics analysis to find several possible key genes which are related with bone metastasis [13]. Currently, we used same gene dataset, expanding to find differential LncRNAs which are associated with bone metastasis, and further elucidate more possible LncRNA-mRNA interaction pairs. The purpose of this study is to introduce a broaden vision that more complicated gene networks are involved into the bone metastasis of breast cancer.  The most merit of the current study is that 71 DElncR-NAs (bone metastatic vs lung metastasis) and 91 DEl-ncRNAs (bone metastasis vs liver metastasis) have been identified. By correlation study between these DElncR-NAs and DEmRNAs, combined with cis nearby-targeted analysis, total nine LncRNA-mRNA interaction pairs were identified. By qRT-PCR study using clinical samples from our hospital, we confirm that most predictions are accurate and reliable. All these enriched data provide a landscape for researcher to clarify the role of LncRNA in breast cancer bone metastasis.
By knocking down the LOC641518, we observed reduced migration and invasion ability which was consistent with lower mRNA levels of LEF1, and overexpression of LOC641518 reversed its effect. LEF1, a member of the T cell Factor (TCF)/LEF1 family of high-mobility group transcription factors, is a downstream component of the Wnt signaling pathway [17]. LEF1 and HOXB9, two Wnt target genes, are identified as promoters of lung adenocarcinoma metastasis, and WNT/TCF pathway is identified as a determinant of metastasis to bone in lung adenocarcinoma [14]. In addition, it was reported that LEF1 was implicated in cell invasion and metastasis of breast cancer [18,19]. LEF1 antisense RNA 1 (LEF1-AS1/LOC641518) is an antisense lncRNA encoded in the LEF1 locus, which is related to the metastasis in various tumors, including colorectal cancer and esophageal squamous cell carcinoma [20,21]. Current study identified that LOC641518 cis-regulated the expression of LEF1 in group of bone metastasis vs lung metastasis, which might indicate that LOC641518 involved in breast cancer bone metastasis through regulating LEF1.
LOC285972-RARRES2, is another interaction pair identified in the current study. RARRES2, also known as chemerin, was considered as a chemoattractant and an adipokine [22]. The correlation between RARRES2 and cancers has been reported. Reduced expression of RARRES2 is observed in human hepatocellular carcinoma [23]. Lu et al. indicate RARRES2 that promotes tumorigenesis and metastasis in oral squamous cell carcinoma [24]. In addition, RARRES2 is identified as a potential biomarker for chordoma which is a rare, low-malignant bone tumor [25]. Low level of LINC00996 was report to be associated with colorectal carcinogenesis and metastasis [26]. In this study, LOC285972 was found to be down-regulated in both two groups and cis-regulate RARRES2 in the group of bone metastasis vs lung metastasis. Hence, we hypothesize that LOC285972-RARRES2 interaction pair is associated with the development of breast cancer bone metastasis. Very low density lipoprotein receptor (VLDLR) has been considered as a multiple function receptor due to binding numerous ligands, causing endocytosis and regulating cellular signaling [27,28]. Previous studies have confirmed that high VLDLR expression is correlated with lymph node and distant metastasis in gastric and breast cancer patients, that VLDLR may be a clinical marker in cancers, and has a potential link with β-catenin signaling pathway [27]. In the current study, we predict that FLJ35024 cis-regulate VLDLR, but very rare report about FLJ35024 so far. Therefore, their interaction deserves further study.
TNFRSF10C copy number variation has been reported that is associated with metastatic colorectal cancer [29], and its hypermethylation was associated with NSCLC [30]. One publication reports that TNFRSF10C plays a role in triple-negative breast cancer cell survival and metastasis regulated by SOX9 [31]. LOC254896 was expected to cis-regulate its expression by present bioinformatics analysis, but LOC254896 is still mysterious and lack of detailed study. One of limitation is that we only verify the role of LOC641518-LEF1 interaction pair in migration and invasion capability of breast cancer cells, and the remaining interaction pairs also deserve further investigation. The other limitation is that there is no in vivo data to be provided in the current study, further, the experimental bone metastasis model should be established to investigate their role in mice. These findings may make contribution to understand the process of bone metastasis through lncRNA regulation, and further be helpful to develop new strategies that effectively predict occurrence of bone metastasis, also provide novel targets for drug design in the future clinical practice. Fig. 6 The qRT-PCR results of DEmRNAs and DElncRNAs in tumor tissues from bone metastatic sites compared with other metastatic sites of breast cancer. A Expression of LEF1, VLDLR, RARRES2 and TNFRSF10C in the tissue of patients with breast cancer bone metastasis, lung metastasis and liver metastasis by qRT-PCR. The x axis and y axis presented gene name and relative expression, respectively. B Expression of LOC641518, FLJ35024, LOC285972 and LOC254896 in the tissue of patients with breast cancer bone metastasis, lung metastasis, and liver metastasis by qRT-PCR. The x axis and y axis presented gene name and relative expression, respectively. C Dot-plot graphic of correlation between LOC641518-LEF1, FLJ35024-VLDLR, LOC200772-AGXT, LOC285972-RARRES2 and LOC254896-TNFRSF10C. #P < 0.05, adjacent non-tumor breast tissue compared with three site metastases; *p < 0.05, bone metastasis compared with lung and liver metastasis over-expression LEF reverses its effect on wound closure. Cell migration was assessed by recover of the scratch. The area of the wound was measured at the two time points in every group, and % reduction of initial scratch area was compared with scramble; and the results were expressed as percentage relative to the corresponding scramble control. C Transwell invasion assay. MCF-7 and MDA-MB-231 cells penetrating the membrane were fixed and 0.1% crystal violet stained after 24 h as described in experimental procedures, and cell number penetrating the membrane were counted for each group. *P < 0.05, **P < 0.01. Data were the means of three measurements and the bars represented SD of the mean