The present study was approved by the Ethics Committee Board of the Chinese Academy of Medical Sciences and Peking Union Medical College. All the patients involved in this program gave informed consent to the work. We obtained keloid and corresponding normal skin from 16 patients who received surgery in Peking Union Medical College Hospital. None of the patients had received previous treatment for keloid. The diagnosis of keloid was confirmed by the pathological results. During the surgery, keloid was excised completely, and normal skin tissue was obtained after trimming. Tissue was frozen by liquid nitrogen and stored at −80°C.

Total RNA was extracted from the frozen tissue using a TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Each step followed the manufacturer's instructions (Invitrogen). The RNA quality was assessed by NanoDrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was assessed by standard denaturing agarose gel electrophoresis.

Keloid and normal skin tissue from 3 randomly selected patients were tested in microarray profiling. We used the ArrayStar Human LncRNA/mRNA Expression Microarray version 3.0 (Arraystar, Inc., Rockville, MD, USA). In total, 30,586 lncRNAs and 26,109 coding transcripts can be detected by this microarray. Public transcriptome databases and lncRNA publications were used to construct lncRNAs.

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technologies) with minor modifications. First, mRNA was purified from total RNA after removal of rRNA (mRNA-ONLY™ Eukaryotic mRNA isolation kit; Epicentre, Madison, WI, USA). Then, each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3′ bias utilizing a random priming method (Quick Amp Labeling kit, One-Color; Agilent p/n 5190-0442). The labeled cRNAs were purified by RNeasy Mini kit (Qiagen p/n 74104). The concentration and specific activity of the labeled cRNAs (pmol Cy3/μg cRNA) were measured by NanoDrop ND-1000. Hybridization adopted Agilent Gene Expression Hybridization kit (Agilent p/n 5188-5242). A total of 1 μg of each labeled cRNA was fragmented by adding 5 μl 10X blocking agent (Agilent p/n 5188-5242) and 1 μl of 25X fragmentation buffer (Agilent p/n 5188-5242), then the mixture was heated at 60°C for 30 min, finally 25 μl 2X GE hybridization buffer (Agilent p/n 5188-5242) was added to dilute the labeled cRNA. A total of 50 μl of hybridization solution 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 (Agilent p/n G2545A). The hybridized arrays were washed, fixed and scanned using the Agilent DNA G2505C microarray.

Agilent Feature Extraction software (version 11.0.1.1) was used to examine acquired array images. Quantile normalization and subsequent data processing were executed using the GeneSpring GX v11.5.1 software package (Agilent Technologies). After quantile normalization of the raw data, lncRNAs and mRNAs, at least 3 out of 6 samples, were chosen for further data analysis. Differentially expressed lncRNAs/mRNAs were classified by a fold change cut-off of 2.0 (upregulated or downregulated) combined with P-value <0.05 as selection criterion. The microarray analysis was performed by Kangcheng Biology Engineering Co., Ltd., (Shanghai, China).

The Gene Ontology (GO) project offers a controlled vocabulary to label gene and gene product attributes in any organism (http://www.geneontology.org). Three fields are covered in the ontology, which are biological process (BP), cellular component (CC) and molecular function (MF). To avoid more than accidental overlap between the DE list and the GO annotation list exists, Fisher's exact test is applied. P-value represents the significance of GO terms enrichment in the DE genes. Pathway analysis is a functional analysis mapping genes to KEGG pathways. The P-value (EASE-score, Fisher-P value or Hypergeometric-P-value) indicates the significance of the pathway correlated to the conditions. The lower the P-value, the more significant is the pathway. The analysis was performed by Kangcheng Biology Engineering.

We chose 4 lncRNAs (uc002lfu.1, ENST00000522743, NR_038439 and ENST00000521141) to explore the connection between these lncRNAs and corresponding coding genes. First, the median value of all the transcripts expressed from the same coding gene was selected, while there was no special treatment with the expression value of lncRNA. Secondly, differentially expressed lncRNAs and mRNAs data were screened and taken away from dataset. Thirdly, Pearson correlation coefficient (PCC) between coding gene and lncRNA was calculated using R-value. Fourthly, PCC ≥0.999 was chosen as meaningful related pair. Finally, the network was drawn through Cytoscape (v2.8.1). In the network, we used blue nodes representing coding genes, and red/pink nodes representing lncRNAs. A solid line indicates a positive correlation, and a dashed line a negative correlation.

We used quantitative real-time PCR (qRT-PCR) to confirm the expression levels. Total RNA was extracted from frozen tissue using a TRIzol reagent (Invitrogen Life Technologies) and reverse-transcribed using RT reagent kit (Thermo Fisher Scientific) according to the manufacturer's instructions. Real-time PCR was performed using the SYBR-Green method with ViiA 7 real-time PCR system (Applied Biosystems) following the manufacturer's protocols. Three significantly upregulated lncRNAs (uc002lfu.1, ENST00000522743 and NR_038439) and one downregulated lncRNAs (ENST00000521141) were chosen to test the reproducibility of the data. After further research, NR_038439 drew our attention and was tested in 13 other pairs of tissues. Glyceraldehyde three-phosphate dehydrogenase (GAPDH) acted as endogenous control. Primer sequences were as follows: GAPDH: forward, 5′-GGGAAACTGTGGCGTGAT-3′ and reverse, 5′-GAGTGGGTGTCGCTGTTGA-3′; uc002lfu.1: forward, 5′-TGCTTGATCCAAATAATGCC-3′ and reverse, 5′-TTCAGTCCAGAGATGTGCC-3′; ENST00000522743: forward, 5′-AGACCCAAAGCTGACACCA-3′ and reverse, 5′-ATCTCCCCTCATCCAAACC-3′; NR_038439: forward, 5′-CTGGCAAGGTTGAGTAGGCT-3′ and reverse, 5′-TGCTCCCTTCACACGGTCA-3′; ENST00000521141: forward, 5′-AAGGATGTGGGAGTTTGAGAC-3′ and reverse, 5′-CAGAGGGGAGAGCGTGTTC-3′. The PCR conditions were as follows: PCR reaction under 95°C for 10 min, then 40 cycles of 95°C for 10 sec and 60°C for 60 sec. After setting the concentration of PCR product as 1, all the PCR products were diluted as follows: 1×10⁻¹, 1×10⁻², 1×10⁻³, 1×10⁻⁴, 1×10⁻⁵, 1×10⁻⁶, 1×10⁻⁷, 1×10⁻⁸ and 1×10⁻⁹. As the target genes and housekeeping gene went through PCR separately, the standard curves were drawn based on the dilution. The relative expression (RQ) was obtained through dividing the concentration of target gene by the one of housekeeping gene. The lncRNA expression differences between keloid and normal skin were analyzed using Student's t-test with SPSS (version 16.0; SPSS, Inc., Chicago, IL, USA). A value of P<0.05 was considered significant.

The microarray we used in the present study can detect 30,586 lncRNAs and 2,6108 annotated mRNAs from authoritative data sources including ‘RefSeq’, ‘UCSC known-gene’ and ‘Genecode’. We found a total of 16,710 lncRNAs expressed in keloid from lncRNA expression profiling data using microarray analysis (Fig. 1A). Through these data, we compared lncRNA expression levels among 3 pairs of human keloid and their adjacent normal skin tissue (Fig. 1C). There were an average of 2,988 upregulated lncRNAs (range from 2,885 to 3,071) and 1,152 downregulated lncRNAs (range from 1,036 to 1,316), which were significantly differentially expressed (≥2.0-fold) (Fig. 2A). To make the results more precise and reduce the possibility of false positive, we set the screen standard of lncRNA as follows: The abnormal lncRNA should exist in all 3 pairs of tissues. After screening, there were 1,731 constantly upregulated lncRNAs and 782 down-regulated lncRNAs in keloid compared with normal skin. Top 25 differentially expressed lncRNAs are listed in Table I. In the upregulated lncRNAs, the maximum fold change was 24.68, which belonged to uc002lfu.1. Moreover, among the downregulated lncRNAs, ENST00000449151 possessed the maximum fold change, which was 84.18. The dendrogram from hierarchical clustering analysis exhibits relationships among lncRNA expression patterns (Fig. 2C).

The mRNA expression profiling data showed a total of 18,788 mRNAs in keloid using microarray analysis (Fig. 1B). After comparing the three keloids with their paired adjacent normal skin tissue (Fig. 1D), we found an average of 1,799 upregulated mRNAs (range from 1,643 to 1,975) and 3,784 downregulated mRNAs (range 3,516 from 4,035) which were significantly differentially expressed (fold change ≥2.0) (Fig. 2B). In the three paired samples, 1,079 mRNAs were consistently upregulated and 3,282 were consistently downregulated. Top 25 differentially expressed mRNAs are listed in Table II. In the upregulated mRNAs, NM_003014 possessed the maximum fold change, which was 18.22. In the downregulated mRNAs, NM_053283 had the maximum fold change, which was 6366.20. The clustering analysis showed relationships among the mRNA expression patterns, which were present in the samples (Fig. 2D).

GO analysis was used to analyze the main function of the closest coding genes according to the GO database. We found that in upregulated transcripts from keloid, the highest enriched GOs were biological regulation (biological process) (Fig. 3A), membrane (cellular component) (Fig. 3B), and receptor activity (molecular function) (Fig. 3C). Also, in the downregulated transcripts, the highest enriched GOs were cellular process (biological process) (Fig. 3D), cell part/cell (cellular component) (Fig. 3E), and binding (molecular function) (Fig. 3F). In pathway analysis, there were a total of 11 pathways related with upregulated transcripts in keloid (Fig. 4A). The most enriched network was Neuroactive ligand-receptor interaction, Homo sapiens (human). Also, 44 pathways were found related with downregulated transcripts in keloid (Fig. 4B). Among them, the most enriched network was valine, leucine and isoleucine degradation, Homo sapiens (human). It is interesting that calcium pathway signaling was one of the items above. Considering calcium played a vital role in wound healing, a common process in keloid formation and calcium ion influx inhibitor was used for keloid treatment (17), calcium pathway signaling was important in pathological behavior.

We built a co-expression network using correlation analysis between differentially expressed lncRNAs (uc002lfu.1, ENST00000522743, NR_038439, ENST00000521141) and mRNAs (Fig. 5). LncRNAs and mRNAs with Pearson correlation coefficients not less than 0.999 were chosen to draw the network by Cytoscape (v2.8.1). There were 4 lncRNAs and 298 mRNAs in this network, 302 nodes possessed 304 pairs of connections between lncRNAs and mRNAs. The network was complex and indicated that one lncRNA was associated with multiple mRNAs, and vice versa.

In order to validate the results from microarray, we selected 3 upregulated (uc002lfu.1, ENST00000522743 and NR_038439) and 1 downregulated lncRNA (ENST00000521141) from the differentially expressed lncRNA profile, then performed quantitative RT-PCR in original tissues. All of the chosen lncRNAs exhibited abundant expression and had significant changes between the two groups in microarray. The results of qRT-PCR were consistent with microarray data (Fig. 6A). Because calcium signaling pathway may be related with keloid, we chose NR_038439 (lncRNA CACNA1G-AS1), an antisense RNA to CACNA1G which encodes a subtype of T-type Ca²⁺ channels, and performed qRT-PCR in another 13 pairs of examples. The result was consistent (P<0.05) (Fig. 6B).