Samples from 5 colon adenocarcinoma patients including tumor and paired non-cancerous tissues were collected and immediately cryo-preserved in liquid nitrogen following surgery at The Sanming First Hospital Affiliated to Fujian Medical University from October 2014 to November 2014. All diagnoses of colon adenocarcinoma were confirmed by pathology, and the non-cancerous samples were obtained 5 cm away from the edge of the tumor and were free of tumor cells, as evaluated by a pathologist. Tumor stage was evaluated based on the tumor-node-metastasis (TNM) staging system of the International Union against Cancer (Table I). At the time of the study, no patient had received chemotherapy, radiotherapy or other preoperative treatments. The study procedures were reviewed and approved by the Ethics Committees of the Sanming First Hospital Affiliated to Fujian Medical University [permit number: 2013(96)]. Informed written consent was obtained from all of the subjects.

Total RNA from all 8 samples were extracted using TRIzol reagent (Invitrogen Life Technologies) following the manufacturer's protocol. RNA quantity and quality were assessed by a NanoDrop ND-1000 spectrometer (NanoDrop Technologies), and RNA integrity was evaluated by standard denaturing agarose gel electrophoresis methods.

Arraystar Human LncRNA Microarray V3.0, which is designed for the global profiling of human lncRNAs and protein-coding transcripts, was used in this study. Microarray V3.0 was updated from the previous V2.0, and approximately 30,586 lncRNAs and 26,109 coding transcripts can be detected by this third-generation lncRNA microarray. All the sequences were collected from authoritative data sources, such as National Center for Biotechnology Information (NCBI) RefSeq, University of California-Santa Cruz (UCSC), Ultra Conserved Regions (UCRs), and related literature. Each transcript is represented by a specific exon or splice junction probe which can accurately distinguish between individual transcripts. The microarray hybridization and bio-information analysis was performed by KangChen Bio-tech (Shanghai, China).

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technologies) with minor modifications. Briefly, mRNA was purified from total RNA using mRNA-ONLY Eukaryotic mRNA isolation kit (Epicentre) after removal of rRNA. Subsequently, each specimen was magnified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3′ bias using Arraystar Flash RNA labeling kit (Arraystar), which is a random priming method. The labeled cRNAs were depurated using the RNeasy Mini kit (Qiagen). The concentration and specific activity of the labeled cRNAs (pmol Cy3/µg cRNA) were measured by NanoDrop ND-1000 (NanoDrop Technologies). One microgram of each labeled cRNA was fragmented by adding 5 µl 10X blocking agents and 1 µl 25X fragmentation buffer, and the mixture was heated at 60°C for 30 min. Next, 25 µl 2X GE hybridization buffer was added to dilute the labeled cRNA. Fifty microliters of hybridization solutions was dispensed into the gasket slide and assembled to the lncRNA expression microarray slide. The slides were incubated for 17 h at 65°C in an Agilent hybridization oven. The hybridized arrays were washed, fixed, and scanned using an Agilent DNA Microarray scanner (part number G2505C).

Acquired array images were analyzed using Agilent Feature Extraction software (version 11.0.1.1). Quantile standardization and subsequent data processing were performed by GeneSpring GX v12.1 software package (Agilent Technologies). Differentially expressed lncRNAs and mRNAs with statistical significance between the two groups were identified through p-value/FDR filtering and fold change filtering (fold change ≥2.0). Hierarchical clustering and combined analysis were performed using homemade scripts.

Gene co-expression networks were built according to the normalized signal intensity of specifically expressed genes. The network construction procedures included the following: i) Preprocess data: the same mRNAs with different transcripts taking the median value represent the gene expression values, without special treatment of the lncRNA expression value; ii) screen data: remove the subset of data according to the lists showing the differential expression of lncRNAs and mRNAs; iii) calculate the Pearson's correlation coefficient and use R-value to calculate the correlation coefficient between lncRNAs and mRNAs; and iv) screen by Pearson's correlation coefficient: select the Pearson's correlation coefficient greater than 0.98 as the meaningful value and draw the lncRNA/mRNA co-expression network by using the Cytoscape program. To make a visual representation, only the strongest correlations (≥0.9994 for DKK4, MMP7, SFRP4 and WIF1; ≥0.994 for MYC, CTNNB1 and PLCB4; ≥0.9975 for SFRP2 and WNT11; ≥0.98 for LEF1) were included in these renderings.

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacture's protocol. The first strand cDNA was synthesized using SuperScript™ III First Strand kits (Invitrogen). GAPDH was used as an internal control. The difference in expression of lncRNA was calculated using the 2^−ΔCt method described by Schmittgen and Livak (11). Primers used in this study are listed in Table II.

Statistical analysis was performed using SPSS 17.0 software package. Significant differential expression levels of lncRNAs or mRNAs were analyzed by Student's t-test and FDR filtering was used for comparative analysis. The p-value ≤0.0.5 (two-tailed) and fold change ≥2.0 were considered statistically significant.

To distinguish the abnormality of lncRNAs in CRC, we contrasted the expression profile between CRC tissues and adjacent non-cancerous tissues from 4 colon adenocarcinoma patients (Fig. 1). Microarray data showed that a total of 3,523 lncRNAs and 2,515 mRNAs were consistently differentially expressed in the CRC tissues in comparison to the non-cancerous tissues. Among them, 1,989 lncRNAs were upregulated with 1,534 lncRNAs downregulated and 1,471 mRNAs were upregulated with 1,044 mRNAs downregulated (see additional file 1 and 2; http addresses are shown in Fig. 1 legend). The top 15 upregulated and top 15 downregulated lncRNAs and mRNAs are listed in Tables III and IV. Furthermore, 156 lncRNAs (108 upregulated and 48 downregulated) and 99 mRNAs (60 upregulated and 39 downregulated) were highly differentially expressed with an absolute fold change >10 [see additional file 1 and 2].

To deduce the function of differentially expressed lncRNAs and mRNA transcripts, GO analysis was used. GO enrichment analysis annotated the biological processes, cellular components and molecular functions of the transcripts (12). P-value assessment was performed to estimate the significance of GO term enrichment in the differentially expressed lncRNAs and mRNAs; a lower p-value indicated that the GO term was more significant than a higher value (p-value ≤0.05 was considered to be statistically significant). Our data showed that the most enriched GO terms were DNA metabolic processes (Fig. 2A), non-membrane-bound organelles (Fig. 2B) and protein binding (Fig. 2C) in upregulated transcripts (GO: 0006259 under biological process, p=4.71932E-11; GO: 0043228 under cellular components, p=7.26282E-14 and GO: 0005515 under molecular function, p=1.36123E-10, respectively). In downregulated transcripts, the most enriched GO terms were organic hydroxyl compound metabolic processes (Fig. 2D), extracellular regions (Fig. 2E) and steroid binding (Fig. 2F) (GO: 1901615 under biological processes, p=2.57819E-10; GO: 0005576 under cellular components, p=8.92426E-08 and GO: 0005496 under molecular function, p=1.19953E-05, respectively).

KEGG pathway analysis indicated that 22 pathways were related to the upregulated transcripts, with the most enriched pathway identified as hsa05322 (Systemic lupus erythematosus-Homo sapiens) (Fig. 3A), and 22 transcripts were associated with this pathway [see additional file 3, http address shown in Fig. 3A legend]. Moreover, 24 pathways correlated with the downregulated transcripts, with the most enriched pathway being hsa00982 (Drug metabolism-cytochrome P450-Homo sapiens) (Fig. 3B), and 13 transcripts were associated with this pathway [see additional file 4, http address shown in Fig. 3B legend). Among the top 10 upregulated pathways, the deregulation of the 'Wnt signaling pathway' was previously identified as a crucial initiation step of colorectal cancer (13), and was reported to be correlated with prognosis in CRC patients (14) (Fig. 3C). The gene category 'lupus erythematosus' has been proven to be associated with several types of cancer, including non-Hodgkin's lymphoma, leukemia, vulva, thyroid, lung, and possibly liver cancers (15). Moreover, the gene categories 'chemokine signaling pathway' (14) and 'PI3K-Akt signaling pathway' (16) have been shown to be associated with CRC. Among the top 10 downregulated pathways, the gene category 'Retinol metabolism' has been reported to be involved in colorectal cancer metastatic multiplicity and the gene category 'Steroid hormone biosynthesis' was found to play important roles in reducing the risk of CRC in many research studies (17) (Fig. 3D).

As mentioned above, the 'Wnt signaling pathway' was upregulated in CRC tissues compared to non-cancerous tissues. Therefore, we aimed to ascertain whether these differentially expressed Wnt signaling pathway components were correlated with these differentially expressed lncRNAs. To explore this problem, coding-non-coding gene co-expression network (CNC network) was constructed between the top 10 differentially expressed 'Wnt signaling pathway' genes and all differentially expressed lncRNAs. We found that the co-expression network was composed of 1,468 nodes and 9,940 connections between 1,458 lncRNAs and 10 coding genes, including 6,004 positively correlated and 3,936 negatively correlated lncRNA-mRNA pairs in the cancer/para-cancer groups [see additional file 5, http address is shown in Fig. 4 legend). These results suggested that the expression of Wnt signaling pathway genes correlated closely with lncRNAs in CRC, and that a possible functional correlation existed between these lncRNAs and Wnt signaling pathway genes in CRC. The representative co-expression networks are shown in Fig. 4.

To verify the microarray results, we randomly selected 4 upregulated lncRNAs and 3 downregulated lncRNAs from the differentially expressed lncRNAs in the microarray (Table II). qRT-PCR method was used to confirm the expression levels of the selected lncRNAs in the CRC tissues and paired adjacent normal tissues used in the microarray test. As shown in Fig. 5, despite the fact that the fold change of each selected lncRNA in the qRT-PCR test was not in full accord with the microarray data, the variation tendency of each lncRNA was consistent with the result of the micro-array. H19, an lncRNA which was reported to be upregulated in multiple solid tumors such as breast, bladder, esophageal and lung was also found to be increased by an average fold of 3.6 in CRC samples compared to normal tissues.