Comprehensive analysis of mRNA-lncRNA co-expression profile revealing crucial role of imprinted gene cluster DLK1-MEG3 in chordoma

Chordoma is a rare bone tumor with high recurrence rate, but the mechanism of its development is unclear. Long non-coding RNAs(lncRNAs) are recently revealed to be regulators in a variety of biological processed by targeting on mRNA transcription. Their expression profile and function in chordoma have not been investigated yet. In this study, we firstly performed the comprehensive analysis of the lncRNA and coding genes expression analysis with three chordoma samples and three fetal nucleus pulposus tissues. lncRNA and gene microarrays were used to determine the differentially expressed lncRNAs and protein coding genes. 2786 lncRNAs and 3286 coding genes were significantly up-regulated in chordoma, while 2042 lncRNAs and 1006 coding genes were down-regulated. Pearson correlation analysis was conducted to correlate differentially expressed lncRNAs with protein coding genes, indicating a comprehensive lncRNA-coding gene co-expression network in chordoma. Cis-correlation analysis showed that various transcripts of MEG3 and MEG8 were paired with the most differentially expressed gene DLK1. As located in the same locus, we further analyzed the miRNA clusters in this region, and identified that 61.22% of these miRNAs were significantly down-regulated, implying the silence of the imprinted gene cluster DLK1-MEG3. Overexpression of MEG3 suppressed the proliferation of chordoma cells. Our study pointed out the potential role of lncRNAs in chordoma, presented the lncRNA-coding genes co-expression profile, and revealed that imprinted gene cluster DLK1-MEG3 contributes to the pathogenesis of chordoma development.


INTRODUCTION
Chordoma is a rare bone tumor with an incidence rate of 0.08/100000. It is widely acknowledged that chordoma is arisen from the remnants of the notochord, with the most common site being the sacrum, skull base and mobile spine [1]. Chordomas are histologically classified as low-grade neoplasm, whereas, its high recurrence rate makes this tumor very similar to malignant neoplasms [2,3]. Currently, the available management for controlling chordoma includes surgery followed by radiation therapy, which is thought to be successful. However, the tendency to recur and metastasize makes this treatment still in a challenge [4]. Recent studies indicated that the growth of chordoma is possibly relying on multiple mechanisms, including genetic and epigenetic modifications. However, the precise mechanisms that promote the development of chordoma remain poorly understood [5,6]. Epigenetic is the study to understand the various mechanisms that result in mitotically heritable changes in gene expression that do not involve any changes in DNA sequences [7][8][9]. The major epigenetic changes include DNA methylation, histone modification and post-transcriptional regulation controlled by non-coding

Research Paper
Oncotarget 112624 www.impactjournals.com/oncotarget RNAs. Based on previous studies, DNA hypermethlation and histone deacetylases were found in chordoma tissues, suggesting aberrant epigenetic modifications may attend in the development of chordoma [6].
Highly expressed non-coding RNAs compose most of the cellular RNA, and their relevance in cell functionality has been known for years. In chordoma, non-coding RNAs have been demonstrated to participate in various biological processes. For instance, microRNAs (miRNAs), one of the subtypes of non-coding RNAs, were found to be ectopically expressed in chordoma samples and cell lines. miRNA-1, one of the miRNAs, was observed to significantly under express in chordoma [10,11], and reexpression of miRNA-1 in chordoma cell line decreased cell growth and proliferation by targeting on MET [12]. Besides, recent findings suggest that another subgroup of non-coding RNAs, which have the transcripts longer than 200 nucleotides in length, are also observed to abnormally express in human cancers [13][14][15]. Accumulating evidence demonstrated that these long noncoding RNAs (lncRNAs) are able to regulate cancer development, metastasis and chemo/radio-resistance by governing a wide repertoire of molecular biological and genetic processes, including transcription, translation, splicing and chromatin modification, etc [16][17][18][19][20][21][22]. lncRNAs can be classified as cis and trans acting molecules. Cis-acting lncRNAs function at the site of transcription and influent their neighbouring genes expression. Trans-acting lncRNAs function away from the site where it synthesized. Several cis-acting lncRNAs are able to guide epigenetic regulator to the transcription site, in order to create an anchor to recruit chromatin remodeling related proteins [23][24][25][26]. Cis-acting lncRNAs are found to control the epigenetic regulation of some imprinted genes. Imprinting depends on the parental origin of the genes, which play crucial roles in development [27]. The expression of imprinted genes must be tightly controlled, and many imprinted gene locus express lncRNAs that were found to regulate the neighbouring imprinted protein-coding genes expression in cis, allele specifically [28]. For example, the lncRNA AIR can silence the nerighbouring imprinting genes including SLC22A3 and IGF2R [29].
In the current study, we performed microarrays to detect the lncRNAs and protein coding genes expression profile of chordoma and the control tissue fetal nucleus pulposus (FNP), in order to find out the lncRNAs involved in chordoma biological processes. In addition, we constructed the co-expression network of lncRNA-coding gene to predict the potential function of the differentially expressed lncRNAs. Among all the differentially expressed coding genes and lncRNAs screened, we identified that the imprinted gene cluster in the DLK1-MEG3 locus might contribute greatly to chordoma development. Our study could be useful to better understand the interplay between non-coding RNAs and coding gene transcripts in chordoma biology.

Demographic and clinical characteristics of three sacral chordoma
Three patients diagnosed with primary sacral chordoma under tumor resection were enrolled in our study. The images, including X-ray, computed tomography (CT) and magnetic resonance imaging (MRI), showed typical image manifestation of chordoma in the sacrum.
On the X-ray ( Figure 1A) and CT images ( Figure 1B), large lytic tumor lesion was observed in the midline and an associated soft-tissue mass, and sacrum bone destruction could be seen. Based on MRI ( Figure 1B), the tumor lesion showed hyper-intense on T2 weighted images. HE staining ( Figure 1C) and immunohistochemical staining (Cytokeratin and S100, data not shown) further indicated that all the three samples collected were chordomas belong to the classic category. The results confirmed that the tissue we harvested were classic chordomas.

Overview of the protein-coding genes and lncRNA expression profile in chordoma and FNP tissue
In order to determine differentially expressed protein-coding genes and lncRNAs, we performed microarray analysis on chordoma and fetal nucleus pulposus tissue. Based on the expression level of all the RNAs screened in the microarray, hierarchical clustering analysis was performed to group protein-coding genes and lncRNAs. The results suggested that the samples come from the same cluster within the case or control group ( Figure 2).
Differentially expressed protein-coding genes and lncRNAs were selected when they were altered by fold change no less than 2.0 (P value < 0.05). A total of 4292 protein-coding genes and 4899 lncRNAs were found to differentially express in chordoma. Concerning to the protein-coding genes expression with more than 100 fold change, there were 9 up-regulated (HLA-DQA1, HLA-DPB1 and C1QB, etc.) and 5 down-regulated genes (DLK1, UCMA and LECT1, etc.). Most of the differentially expressed genes belonged to the 2 ≤ FC < 5 group ( Figure 3A). Similar to the differentially expressed protein-coding genes, 16 up-regulated lncRNAs (NONHSAT060785, NONHSAT116806 and NONHSAT104550, etc.) and 6 down-regulated lncRNAs (NONHSAT113429, NONHSAT113428 and NONHSAG044080, etc.) were identified with expression change more than 100 fold, and the group with most differentially expressed lncRNAs was 2 ≤ FC < 5 ( Figure 2B). The list of differentially expressed proteincoding genes and lncRNAs were shown in Supplementary  Table 1 and Supplementary Table 2. To further exhibit the distribution of these differentially expressed protein-Oncotarget 112625 www.impactjournals.com/oncotarget coding genes and lncRNAs, we constructed the volcano plots. As shown in Figure 3C, The significance versus fold change was plotted into volcano plots, which reflected the distribution of the screened protein coding genes and lncRNAs between chordoma and FNP samples.
To validate the microarray results, we selected 5 coding genes and 4 lncRNAs for quantitative polymerase chain reaction (qPCR). The results showed that IGF-1R, EGFR, MET and PDGFRB expression was significantly higher in chodoma, and DLK1 expression was decreased as shown in microarray data ( Figure  3D). qPCR analysis of the randomly selected lncRNAs revealed that NR_033359.1, NONHSAT054600 and ENST00000553465 were significantly lower in chordoma, and NONHSAT097304 expression was much higher in chordoma as reflected by microarray analysis ( Figure 3E).

Subgroups classification of differentially expressed lncRNAs
We continuously adapted specific probes for lncRNAs to further classify the differentially expressed lncRNAs in accordance with their different features. As shown in Figure 4, there were 141 up-regulated and 104 down-regulated sense lncRNAs, and 149 up-regulated and 85 down-regulated anti-sense lncRNAs were found as well. The profile data suggested 66 bidirectional lincRNAs were up-regulated, and 37 were down-regulated. As for intergenic lncRNAs, 709 were up-regulated and 973 were down-regulated. Meanwhile, we also identified 1721 upregulated and 843 down-regulated intronic lncRNAs in chordoma.

Construction of the lncRNAss and proteincoding genes co-expression network
To further explore the potential regulatory mechanism of the differentially expressed lncRNAs, the lncRNAs and protein-coding genes co-expression network was constructed based on correlation analysis. Pearson correlation coefficient analysis was performed to evaluate the correlation between lncRNAs and coding genes, and the lncRNAs and coding genes with Pearson correlation coefficient no less than 0.90 and the P value less than 0.01 were picked up for network construction. The results of lncRNA and coding genes correlation analysis were listed in Supplementary Table 3. As DLK1 was found to be the most differentially expressed protein-coding gene in chordoma (down regulated with FC 2771.79), we analyzed the lncRNAs that correlate with DLK1 specifically ( Figure  5A). Interestingly, we found that 5 transcripts of maternal expressed gene 3 (MEG3) and 3 transcripts of MEG8 were correlated with DLK1 ( Figure 5A), indicating the possible regulatory role of MEG3 and MEG8 for DLK1.

Cis-regulation of lncRNAs
As lncRNAs are physically linked to the locus from which they are encoded, they may have the ability to function during or immediately following the transcription process. To understand the cis-regulation of lncRNAs in chordoma, the lncRNAs and the associated protein coding genes that screened by Pearson correlation coefficient (Correlation ≥ 0.90. q value < 0.01) were further mapped to their genomic locus. lncRNAs which locate within 300K windows upstream or downstream of the given protein coding genes were further analyzed. The most 200 over or under expressed lncRNAs were selected in this analysis, in which, 129 were located. The lncRNA-coding gene pairs for cis-regulation were shown in Figure 5B. We noticed that several of the down-regulated lncRNAs were all targeting on DLK1, the most under expressed coding gene in chordoma.
There are several isoforms of MEG3 and MEG8, among which, five of them were thought to cis-correlate with DKL1, including ENST00000398474 (MEG3-005), ENST00000455531 (MEG3-008), ENST00000452514 (MEG3-023), ENST00000521256 (MEG3-024), ENST00 000556407 (MEG3-025), ENST00000556475 (MEG8-006), ENST00000553584 (MEG8-015) and ENST00000 553465 (MEG8-016) ( Figure 5B). Concomitant downregulation of the DLK1 and MEGs imply that imprinted genes may contribute to the development of chordoma. Concomitant down-regulation of the imprinted genes DLK1 and MEG3 potentially silenced the corresponding miRNA clusters expression As reported, a cluster of imprinted genes at 14q32.2, the DLK1-MEG3 locus, were found silence in several types of cancers [30], and loss of imprinting genes impairs various biological processes [31]. Our lncRNA-coding genes profile data showed that the most differentially expressed protein coding gene in chordoma was DLK1, and Cis-regulation analysis revealed that this gene was targeted by a bunch of lncRNAs, most of which were identified to be different transcripts of MEG3 and MEG8. They were concomitantly down-regulated in chordoma. As miRNAs cluster consists the major part of the DLK1-MEGs locus, and DLK1-MEG3 silencing induced miRNA suppression were found to appear in a few pathological processes [32]. Therefore, we performed miRNA microarray analysis on the chordoma samples used in the lncRNA-coding gene expression analysis. The clusters of miRNAs that locate in the 14q32.2 locus were analyzed, for which, the miRNAs expression with no less than two-fold change were considered as different. The relative expression value of the miRNAs clusters in the DLK1-MEG3 locus was presented in Figure 6A (miRNA cluster A) and 6B (miRNA cluster B), showing that most of the miRNAs were down-regulated. Totally, 30 (61.22%) miRNAs were identified to be down-regulated, and 3 (6.12%) miRNAs were up-regulated in the cluster of 14q32.2 in chordoma, and 16 miRNAs remain showed the distribution of differentially expressed protein coding genes and lncRNAs identified in chordoma. The protein coding genes or lncRNAs with more than 2 fold change were noted as red dots(positive) and green dots(negative). Randomly selected protein coding genes (D) and lncRNAs (E) were validated by qPCR. ( * represents for P value < 0.05, and ** represents for P value < 0.01).
Oncotarget 112628 www.impactjournals.com/oncotarget unchanged ( Figure 6C). Furthermore, we did KEGG analysis of the genes that correlated with MEG3 and MEG8, revealing that these genes were enriched in p53, MAPK and TNF et al. signaling pathways ( Figure 6D), most of which are active in cancer biology. The results above validated that DLK1-MGE3 locus was silenced in chordoma, and it might contribute to the development of this neoplasm.

Over-expression of MEG3 suppresses chordoma cells proliferation
As determined by qPCR, MEG3 expression increased significantly in chordoma cells after transfection with the pcDNA-MEG3 ( Figure 7A). The expression of its down stream target genes including P53 and Bcl2 was also changed significantly ( Figure 7B), indicating successful  Oncotarget 112629 www.impactjournals.com/oncotarget over-expression of MEG3 in chordoma cells. Increased P53 and decreased Bcl2 expression may also imply the enhanced apoptosis in chordoma cells with MEG3 overexpression. Cell cycle analysis revealed that MEG3 increased the percentage of cells in the G1/G0 phase compared with control ( Figure 7C). Meanwhile, EdU staining further demonstrated that the proliferation rate was significantly lowered in MEG3 over-expression group comparing with the control group ( Figure 7D and 7E). The results suggest that MEG3 can effectively suppress chordoma cells proliferation.

DISCUSSION
Chordoma is a rare bone tumor with the characterization of chemo-and radio-resistance. The extremely high recurrence rate of chordoma poses this tumor on the borderline of low-grade and malignant neoplasm [1]. Although chordoma has been characterized for more than a century, little of the fundamental mechanisms and structures in the initiation and progression of chordoma is known. Growing evidence suggests the emerging role of lncRNAs in the various biological process influencing cellular metabolism and development [33]. The newly characterized non-coding RNAs have drawn the scientific attention since their indispensable roles in cancer development and progression [34]. Benefitting from the advance of new technology achieving unprecedented depths in RNA sequencing, a large pool of lncRNAs have been identified to spread over the mammalian genome. As the integration analysis of lncRNA-coding gene is one of the most common approaches to predict lncRNAs functions in different biological processes, therefore, we comprehensively analyzed the lncRNAs and mRNAs in order to predict their potential functions in chordoma.
Our micro-array data revealed that a large number of genes were differentially expressed in chordoma tissues, including EGFR [35,36], PDGFRβ [37], IGF-1R [38,39] and MET [40,41] etc., which have been identified to the potential therapeutic targets in chordoma with both clinical and molecular evidence. Unlike gene coding RNAs, lncRNAs are observed to express at a low level, and the expression profile in cancer development is specific [42]. The results also indicated that abundant lncRNAs have differential expression patterns in chordoma. Based on the relationship between lncRNAs and their associated Oncotarget 112630 www.impactjournals.com/oncotarget protein-coding genes, we classified the lncRNAs into different groups, including sense, anti-sense, bidirectional, intergenic and intronic lncRNAs [43]. Particularly, there are a large number of intergenic lncRNAs (lincRNA), interspersed in the regions between coding transcripts, found to differentially express in chordomas. lincRNAs were known to exert important functions in various cellular processes, such as maintaining stem cells pluripotency, promoting cell proliferation and accelerating cancer progression [44,45].
Genomic imprinting is an epigenetic modification of a specific parental chromosome imparted during gametogenesis that induces differential expression of the two alleles [46]. Imprinted genes are known to regulate various physiological processes. Several maternally imprinted genes possess tumor-suppressive potential and were found suppressed in several types of cancer. A large number of lncRNAs, which have been found at all imprinted clusters, are always expressed from the parental chromosome with the unmethylated imprint control region (ICR) [46]. Dysregulation of the lncRNAs in imprinted clusters can silence the imprinted genes, and loss of imprinted genes has been associated with cancer, imprinting-related disease and psychiatric disorders  changed (B, C). Cell cycle analysis showed the chordoma cells were retained in the G1/G0 phase with MEG3 over-expression. (D, E) EdU staining of the chordoma cells showed the decreased proliferation rate in MEG3 over-expression group.( * represents for P value < 0.05, and ** represents for P value < 0.01). www.impactjournals.com/oncotarget [32]. Comprehensive analysis of lncRNAs and mRNAs expression profiles gives us a better understanding of biological functions of the differentially expressed lncRNAs. Therefore, we performed the correlation analysis of the differentially expressed lncRNAs and protein coding genes, and identified abundant correlations between the transcripts of the two groups. Our lncRNAcoding gene profile revealed concomitant down-regulation of DLK1 and MEG3 in chordoma samples. Maternally expressed gene MEG3 is reciprocally imprinted with the paternally expressed gene DLK1, constituting an imprinting domain on human chromosome 14q32 [47]. MEG3, a large non-coding RNAs as its transcript lacks open reading frame, is positioned about 100 kb from the coding gene DLK1 [47,48]. MEG3 is involved in various physiological and pathological processes of cell biology, and act as a tumor suppressor through interaction with p53 and MDM2 [32]. MEG3 is highly expressed in normal tissues, including pituitary, brain, adrenal gland and placenta [49]. However, loss of MEG3 expression has been observed in neuroblastomas, hepatocellular cancers and gliomas [50][51][52]. The inhibition effect of MEG3 on tumor cells proliferation was tested in several cancer cell lines with a number of independent assays, including colony formation, growth curve and BrdU incorporation. Cancer cells with MEG3 overexpression grew significantly slower than the control cells [49,50]. It is reported that apoptosis induced by MEG3 has also been observed in cancer cells [52]. The collected evidence demonstrated that MEG3 restricts cancer cells biological behavior by inhibiting their proliferation and inducing apoptosis. Further in vivo observations confirmed that re-expression of MEG3 in PDFS cells can reduce tumor growth in nude mice [53]. The role of DLK1 in tumor is not clear, and it acts as either an oncogene or tumor suppressor. The DLK1 expression is observed to elevate in gliomas, and overexpression of DKL1 in glioma cell line simulated cells proliferation [54]. In contrast, Kawakami et al. [55] demonstrated that DLK1 expression is lost in human renal cell carcinomas, and its re-expression induced cell death and suppressed tumor growth in nude mice. In several previous reports on other cancer types, either MEG3 or DLK1 was reported to become deregulated, but not both. Under expression of both DLK1 and MEG3 are not always expected, as it is difficult to envision either allelic loss or epitope switching could lead to the concomitant down-regulation of these two imprinted genes. Notably, previous studies also found that miRNAs cluster, which is another major part of DLK1-MEG3 imprinted domain, has also been found to be suppressed or silenced significantly [56,57]. Therefore, we evaluated the expression of miRNAs in the DLK1-MEG3 locus with microarray, and the results showed that 30 of the 49 miRNAs expression was decreased or lost in chordoma samples, suggesting that silencing may extend across a large part or the entire imprinted gene cluster at 14q32. One of the possible reasons that explain the actual mechanism is novel repressed epigenetic state established during chordoma formation at the 14q32 gene cluster, and this has been demonstrated to occur in urothelial carcinoma [58]. However, further studies should be done to explore the mechanism that inactivates the DLK1-MEG3 imprinting locus in chordoma.
In conclusion, our study mapped the expression profile of lncRNAs and protein coding genes in chordoma, and constructed the lncRNA-coding gene co-expression network. A large number of lncRNAs and protein coding genes expression were found to be dysregulated. Cisregulation analysis and functional assays revealed that the suppression of imprinted gene locus DLK1-MEG3 is crucial to the development in chordoma. Meanwhile, further research work needs to be carried out to explain the function of other lncRNAs and their regulatory effect in chordoma.

Sample preparation and RNA extraction
This study was approved by the Ethical Committee of the First Affiliated Hospital of Soochow University. Written consent was obtained before all subjects were enrolled. The chordoma samples were collected from 3 patients who were identified and treated by tumor resection surgeries at the Department of Orthopedics, the First Affiliated Hospital of Soochow University. The detailed information of the three patients was provided in Supplementary Table 4. FNP specimens were obtained from aborted fetuses with a gestational age of 8-24 weeks in the Department of Gynecology and Obstetrics. As shown in Figure 1, imaging and histopathological analysis confirmed the properties of chordoma and fetal nucleus pulposus tissues. Following informed consent at the time of acquisition, the samples were collected and stored in liquid nitrogen. Tissue samples were homogenized in TRIzol reagent (Ambion). In our study, 3 chordoma samples and 3 FNP samples were subjected to the microarray analysis. Total RNA was quantified by the NanoDrop ND-2000 (Thermo Scientific), and the RNA integrity was assessed using Agilent Bioanalyzer 2100 (Agilent Technologies).

Microarray assays
The chordoma and fetal nucleus pulposus samples RNAs were employed to synthesize double-stranded complementary DNA (cDNA), and then labeled with the Agilent Human lncRNA (4 * 180K, Design ID:062918) and Agilent Human miRNA (8 * 60K, Design ID: 046064). Microarray hybridization and washing were performed based on the manufacturer's standard protocols. Briefly, total RNA were transcribed to double strand cDNA, then synthesized into cRNA and labeled with Cyanine-3-CTP. www.impactjournals.com/oncotarget The labeled cDNAs were hybridized onto the microarray. After washing, the arrays were scanned by the Agilent Scanner G2505C (Agilent Technologies).

Microarray data analysis
Feature Extraction software (version10.7.1.1, Agilent Technologies) was used to analyze array images in order to get the raw data. Genespring was employed to perform the basic analysis with the raw data. Firstly, the raw data was normalized with the quantile algorithm. Differentially expressed genes or lncRNAs were then identified through fold change as well as P value assessed by student t-test. The threshold set for differentially expressed genes were fold change more than 2.0, and the P value less than 0.05. Afterward, Hierarchical Clustering was performed to display the distinguishable genes' expression pattern among samples. Treeview software (Stanford university, CA, USA) composed by Java was utilized to visualize the microarray data.

qPCR analysis
The total RNA was extracted from the chordoma samples and FNP tissues collected for microarray assay by using TRIzol reagent (Invitrogen, CA,USA). Then, the extracted RNA was reversely transcribed using PrimeScript RT reagent kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China) according to the manufacture's recommendations. The primers sequences for the selected lncRNAs and protein coding genes were listed in Supplementary Table 5. The real-time PCR was performed by using SYBRGreen (TaKaRa, Dalian, China), and GAPDH was employed as the internal control. The expression of each lncRNA was represented as fold change using 2 -ΔΔCt [59]. Student t-test was utilized to analyze the significance of the expression between chordoma samples and their control samples. The value of P < 0.05 was considered significant.

Correlation analysis between lncRNA and mRNA
The network between lncRNA and mRNA was constructed based on the correlation analysis among differentially expressed lncRNAs and protein coding genes. For each pair, Pearson correlation was performed to assess the correlation. For the pairs with absolute value of Pearson correlation coefficients not less than 0.90 were selected to draw the network using Cytoscape (National Resource for Network Biology).

KEGG pathway analysis
To identify the pathways of intersect, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg/) enrichment analysis. The Fisher's exact test, χ2 test and the threshold of significance were conducted for the analysis. The P value and FDR were calculated, and P < 0.05 was used as the screening criteria.

Cis-correlation analysis
The cis-regulation regions were identified by the following procedures. For each lncRNAs, we identified the mRNAs as "cis-regulated mRNAs" when: (1) the mRNAs locus are within 300k windows up-and downstream of the given lncRNA, (2) the Pearson correlation of lncRNA-mRNA expression is significant (P-value of correlation less than 0.05).

Cell culture and transfection
The human chordoma cell line U-CH1 was obtained from ATCC (USA). Chordoma cells were cultured in IMDM and RPMI1640 (1:4) supplemented with 10% fetal bovine serum (FBS), and penicillin/streptomycin (100 U/ ml) (Invitrogen, CA,USA). The cells were cultured in an incubator at 37°C and an atmosphere of 5% CO 2 , and the medium was changed every 3 days. The cells were subcultured once in confluence. The MEG3 overexpression vector (pcDNA-MEG3) and the corresponding control (empty pcDNA) were constructed. MEG3 sequence was synthesized and subcloned into pcDNA3.1. The gene of lncRNA MEG3 was confirmed by restriction enzyme digestion (Supplementary Figure 1). The constructed pcDNA-MEG3 and the empty pcDNA were respectively transfected into chordoma cells using Lipofectamine 2000 reagent (Invitrogen, CA,USA). RNA extraction and EdU staining assay were performed 3 days after transfection.

Cell cycle analysis
For cell cycle analysis, the chordoma cells were fixed with ice-cold 70% ethanol for 12 hours at −20°C. Then the cells were washed and re-suspended and stained with propidium iodide solution (Invitrogen, CA,USA) at the concentration of 5ug/ml. The cells were analyzed by Flow Cytometry (BD Biosciences, Franklin Lakes, NJ, USA).

EdU staining assay
EdU staining for chordoma cell proliferation was assessed using the Cell-Light™ EdU detection kit (RiboBio, Guangzhou, China) according to the manufacturer's instructions. The chordoma cells were incubated with 10 µM EdU for 2 hours before fixed with 4% paraformaldehyde. After fixation, the chordoma cells were permeabilizated using 0.5% Triton X-100 for 10 min. Then, the cells were stained with Apollo643 for 30 min, www.impactjournals.com/oncotarget and subsequently stained with Hoechst. Fluorescent cells were counted using confocal microscope.

Statistical analysis
Statistical analysis was performed with SPSS 17.0 software. Student t-test was used to compare the significance between control and MEG3 overexpression groups. All the data were expressed as mean ± standard deviation (SD). P value < 0.05 was considered as significant, and P value < 0.01 was considered as highly significant.

Author contributions
CKW ang YHL conceived and designed the experiments. CH and ZK performed the experiments. CH, ZK, LJ and WGZ analyzed the data. CH wrote the manuscript. All the other authors edited the manuscript.