Comprehensive Analysis of Long Noncoding RNAs and Messenger RNAs Expression Profiles in Patients with Marjolin Ulcer.

BACKGROUND Marjolin ulcer (MU) is an aggressive cutaneous malignancy. Typically, MU occurs over a period of time in post-burn and/or post-traumatic lesions and scars. However, the pathogenesis of scar carcinogenesis and MU development remains to be elucidated. The present study aimed to investigate the long noncoding RNA (lncRNA) and messenger RNA (mRNA) expression profiling in MU, which could provide new information on the potential molecular mechanisms of MU development. MATERIAL AND METHODS The lncRNA microarray analysis was conducted in normal skin, scar, and MU tissue, and quantitative real-time PCR experiment was carried out to validate the reliability of the microarray data. Furthermore, a series of integrative bioinformatic approaches were applied to decipher the function of differentially expressed lncRNAs. RESULTS A total of 7130 lncRNAs and 9867 mRNAs were differentially expressed among normal skin, scar, and MU tissues. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis demonstrated that these aberrantly expressed transcripts were mainly involved in cell cycle, immune response, and the p53 signaling pathway. Series Test of Cluster analysis indicated certain dysregulated lncRNAs were expressed with a gradually increasing or decreasing trend and might participated in malignant transformation of scar tissue postburn. Co-expression analysis showed 5 selected lncRNAs might regulate cell proliferation through the p53 signaling pathway. Finally, the competing endogenous RNA (ceRNA) network indicated that lncRNA uc001oou.3 might be implicated in ceRNA mechanism during MU development. CONCLUSIONS Taken together, our study implied the aberrant expression of lncRNAs may play an important role in the pathogenesis and development of MU, and the exact mechanism warrants further investigation.

Background

Marjolin ulcer (MU) is an aggressive cutaneous malignancy that typically follows post-burn and/or post-traumatic lesions and scaring, as well as some chronic wounds (e.g., osteomyelitis, pressure sores, venous stasis ulcers, and dermatitis artefacta) [1-3]. The incidence of malignant transformation from burn scar tissue or other pathologic lesions is about 1% to 2% [4,5]. It takes 11 to 75 years from the initial injury to onset of malignant transformation, with an average of 30 to 36 years [1,6]. Although some epidemiological studies have been conducted, the exact etiology and pathogenesis of MU remains ambiguous and controversial. Previous theories attempting to explain the occurrence of MU include repeatedly ulcerated stimulation, prolonged proliferation due to wound infection and long-term chronic inflammation stimulation, irritation by ongoing exposure to toxins released from damaged tissues, weakened immunity resulting from poor blood circulation and lymphatic regeneration in scar tissue, genetic mutations, and interaction between environment and heredity [2,7,8]. MU makes up 1.2% of all cutaneous tumors; and about 2% of squamous cell carcinomas (SCCs) and 0.03% of basal cell carcinomas (BCCs) originate in burn scars [1,9]. Compared to other skin carcinomas of similar histotype, MU has worse prognosis and higher recurrence and metastasis rates [10]. Therefore, it is of great necessity to carefully explore the mechanism of MU.

Although more than 90% of the genomic DNA is transcribed in eukaryotic genomes, only about 2% of the human genome is protein-coding genes, whereas the overwhelming majority of gene transcripts are non-coding RNAs (ncRNAs), in which long non-coding RNAs (lncRNAs) account for 4–9% [11,12]. LncRNAs, a class of ncRNA molecules longer than 200 nucleotides with no or little protein-coding capacity, can function as a regulator of gene expression on a number of levels including epigenetic, transcriptional, and post-transcriptional regulations. Accumulating studies have indicated that lncRNAs are important to maintain homeostasis and participate in the physiological or pathological process, including tumors [11,13]. They are extensively involved in the regulation of proliferation, DNA damage response, apoptosis, and cell cycle in cancer cells, which are responsible for tumorigenesis and tumor progression [14]. Recently, the competing endogenous RNA (ceRNA) hypothesis, which states that lncRNA can negatively regulate the activity of miRNAs by ceRNA mechanism and thus influence expression of protein-coding transcript, has been verified in several tumors (e.g., liver, lung, and breast cancer). Just like the abnormal regulation of a signaling pathway is responsible for malignant transformation and tumor progression, disorder of the ceRNA network may also lead to malignant transformation and tumor progression.

Accordingly, we speculated that lncRNAs may play important roles in scar carcinogenesis postburn and progression of MU. In this study, we conducted large scale analysis of lncRNA and mRNA expression profiling in MU. Furthermore, we delineated comprehensive functional landscapes of lncRNAs in MU by bioinformatics approaches. We hope that a comprehensive analysis of the differentially expressed (DE) lncRNAs may provide new information regarding the potential role of lncRNAs in MU or provide therapeutic strategies for the treatment of MU.

Material and Methods

Samples collection and preservation

Tissue samples including normal skin tissues, para-cancerous scar tissues, and MU tissues were obtained from surgical specimens of 11 patients with MU at Xiangya Hospital Central South University (clinical characteristics was shown in Supplementary Table 1). All samples were confirmed by hematoxylin-eosin (HE) staining after operation (Supplementary Figure 1). Informed consent was obtained from all patients. The study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Xiangya Hospital Central South University (ethical approval code: 201603079). The tissue samples were immediately put into the liquid nitrogen for preservation after excision.

RNA extraction

Total RNA of each sample was extracted using TRIzol® Reagent (Invitrogen, USA) according to the manufacturer’s instructions. RNA quantity and quality were measured by NanoDrop ND-1000 and Agilent 2100 Bioanalyzer, respectively. RNA integrity was assessed by standard denaturing agarose gel electrophoresis. Then, the RNA samples were stored at −80°C until further use.

LncRNA microarray

Arraystar Human lncRNA Microarray V3.0 detected 30586 lncRNAs and 26109 coding transcripts, and the lncRNAs were carefully selected from the authoritative public transcriptome database (including Refseq, UCSC known genes, Gencode, etc.) and highly-respected publications. Fold-change (MU vs. scar, MU vs. normal skin, and scar vs. normal skin) and P-values were calculated from the normalized expression levels. Fold-change ≥2.0 and P<0.05 were considered differential expression. The microarray work was completed by Kang Chen Bio-tech, Shanghai, P. R. China. (GEO accession number: GSE110138)

Quantitative real-time PCR (qRT-PCR) validation

A total of 1 μg purified RNA sample was reverse transcribed using the PrimeScript RT reagent kit (Takara, Japan), and quantitative real-time PCR (qRT-PCR) was performed using SYBR® Fast qPCR Mix (Takara, Japan) on ABI VII7 PCR System (Applied Biosystems, USA) following the manufacturer’s protocols. The PCR conditions included an initial step at 95°C for 30 seconds, followed by 40 cycles of amplification and quantification (95°C for 5 seconds, 60°C for 34 seconds). The sequences of the primers used are listed in Supplementary Table 2. The expression levels of RNA were normalized to internal control β-actin, and then calculated using the 2−ΔΔCt method. LncRNAs were selected for validation as follows: 1) differentially expressed among MU, scar and normal skin tissues simultaneously; 2) dysregulated significantly, with raw signal intensity >200; and 3) RNA length <2000, associated coding genes of which were related to tumor or inflammation. The mRNAs associated with tumor were selected.

Gene Ontology (GO) and KEGG pathway analysis

Gene Ontology (GO) has developed 3 structured ontologies that describe gene products in terms of their associated biological processes, cellular components, and molecular functions. Fisher’s exact test was used to find if there was true difference between groups. Pathway analysis was based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database () usually used to obtain insight into the underlying biology of differentially expressed gene.

Series Test of Cluster (STC) analysis

Series Test of Cluster (STC) analysis was performed based on previous research [15]. First, we set MU, scar, and normal skin as different points, and then chose DE lncRNAs among them according to the random variance model corrective ANOVA, which was used to determine genes that were expressed separately and differentially. Subsequently, a set of unique and representative temporal expression profiles were identified. After clustering algorithm, each lncRNA was assigned to the model profile which most closely matched the lncRNA’s expression profile as determined by the correlation coefficient.

Construction of lncRNA-mRNA co-expression network

The lncRNA-mRNA co-expression network based on calculating the Pearson’s correlation coefficient (PCC) between the expression levels of coding and noncoding genes was constructed. LncRNA were initially selected if they: 1) had an antisense property; 2) were from significant profiles No.3, 7, 12, and with significant change in MU tissues in corresponding profile; and 3) were validated by qRT-PCR. After integration, 5 lncRNAs were further analyzed. Additionally, mRNAs were screened from our microarray using the selection parameters PCC ≥0.97 and P<0.05.

Competing endogenous RNA (ceRNA) network analysis

The mRNAs sharing a meaningful correlation and similar tendency of expression with this 5 lncRNAs were subjected to further analysis. The potential miRNA response elements (MREs) were searched on the sequences of lncRNAs and mRNAs, and the overlapping of the same MRE both on lncRNA and mRNA predicted lncRNA-miRNA-mRNA interaction. To better enhance the reliability of the ceRNA network, miRNAs harboring less than 2 miRNA binding-sites with its targeted lncRNA were filtered. The miRNA binding sites were predicted by miRcode (), while the miRNA-mRNA interactions were predicted by Targetscan ().

Statistical analysis

The data were presented as the means ±SD. All statistical data were analyzed using GraphPad PRISM software (version 6, San Diego, CA, USA). The unpaired t-test was performed to analyze the statistical significance of the microarray and qRT-PCR results. P<0.05 is considered statistically significant.

Results

Dysregulated lncRNAs and mRNAs among normal skin, scar and MU

Based on our microarray result, a total of 7130 lncRNAs and 9867 mRNAs were identified to be differentially expressed among normal skin, scar, and MU tissues. There were 543 lncRNAs that had sense exonic overlap with protein-coding genes, 1044 lncRNAs that were intronic sequences, and 918 lncRNAs that were natural antisense sequences. Visualization of the aberrantly expressed genes between different tissues was illustrated by bar graph (Figure 1A). Among them, 490 lncRNAs and 1440 mRNAs were upregulated in MU compared with para-cancerous scar as well as normal skin, while 1563 lncRNAs and 711 mRNAs were downregulated (Figure 1B, 1C). In addition, there were 1211 lncRNAs and 1790 mRNAs dysregulated specially between MU and para-cancerous scar group, meanwhile, 3306 lncRNAs and 5049 mRNAs between MU and skin group (Figure 1B, 1C). These results suggested that those transcripts may be aberrantly expressed in different phase from the normal skin to scar to MU, and thus exert corresponding effects.

Figure 1

LncRNA and mRNA expression profiles based on microarray data. (A) Bar graph of the number of dysregulated lncRNAs and mRNAs among normal skin, scar, and MU tissues. Venn diagrams showed the numbers of overlapping and nonoverlapping differentially expressed lncRNAs (B) and mRNAs (C) in MU compared with scar and normal skin, respectively. MU – Marjolin ulcer; lncRNA – long noncoding RNA, mRNA – messenger RNA.

Validation of qRT-PCR

To validate the reliability of the microarray data, the qRT-PCR experiment was carried out to detect the expression levels of those selected lncRNAs and mRNAs (described in the method section) in normal skin, scar, and MU samples from 8 patients (Figure 2). Interestingly, the fold changes of qRT-PCR results were smaller than that in microarray results, such as for lncRNA TCONS00007922, uc001oou.3, and mRNA PML, suggesting that the microarray data may enlarge the difference to some degree. Overall, the qRT-PCR results were consistent with the expression trends of the microarray results.

Figure 2

Validation for the expression of candidate transcripts among normal skin, scar, and MU by qRT-PCR. Differential expression from microarray analysis of 4 upregulated lncRNAs (A), 4 downregulated lncRNAs (C), and 4 mRNAs (E) in MU compared with scar and normal skin. The raw intensity of those transcripts was quantile normalized. The relative expression levels of the same transcripts were determined by qRT-PCR (B, D, F). Gene expression was normalized to β-actin, data are presented as mean ±SD, n=8. * P<0.05, ** P<0.01. MU – Marjolin ulcer; qRT-PCR – quantitative real-time PCR; lncRNA – long noncoding RNA; mRNA – messenger RNA.

Functional enrichment analysis

To better understand the potential functions and mechanism of dysregulated lncRNAs in MU, the GO analysis and KEGG pathway analysis were used. Our data showed that the upregulated lncRNAs (MU vs. scar and MU vs. normal skin) were primarily involved in cell cycle, mitotic cell cycle, protein metabolic process, and immune system process (Figure 3A, 3B), while the downregulated lncRNAs were involved in immune response, immune system process, growth, and regulation of growth (Figure 3C, 3D). The KEGG pathway analysis revealed that the upregulated transcripts were mainly involved in the p53 signaling pathway, steroid biosynthesis, and protein processing in the endoplasmic reticulum (Figure 3E, 3F). Meanwhile, the highly-enriched pathways that corresponded to the downregulated transcripts were chemokine signaling pathway, B cell receptor signaling pathway, and melanogenesis (Figure 3G, 3H).

Figure 3

Functional enrichment analysis of the differentially expressed lncRNAs. Top 10 GO terms of biological processes for upregulated lncRNAs between MU and scar (A) and between MU and normal skin (B). Top 10 GO terms of biological processes for downregulated lncRNAs between MU and scar (C) and between MU and normal skin (D). Top 10 KEGG pathways for upregulated lncRNAs (E) and downregulated lncRNAs (G) between MU and scar. Top 10 KEGG pathways for upregulated lncRNAs (F) and downregulated lncRNAs (H) between MU and normal skin. GO – Gene Ontology; lncRNA – long noncoding RNA; MU – Marjolin ulcer; KEGG – Kyoto Encyclopedia of Genes and Genomes.

STC analysis of deregulated lncRNAs

STC analysis was implemented to reveal gene expression dynamics, which explicitly took into account the temporal changes in certain important lncRNA categories that were dysregulated during the phase of scar formation and scar hyperplasia, as well as scar carcinogenesis. Sixteen model profiles were used to summarize the expression pattern of 8112 lncRNAs (Figure 4A). Eight patterns (profiles No.3, 0, 2, 4, 7, 15, 12, 11) were statistically significant (P<0.05), and profile No.3 was the most significant pattern. LncRNAs in profile No.3 (Figure 4B) showed a gradually upregulated trend from the normal skin to scar to MU, while profile No.12 (Figure 4C) showed a progressively upregulated trend, and profile No.7 (Figure 4D) was basically downregulated in MU. The STC results suggested that scar tissue may be the intermediate step from normal skin to MU.

Figure 4

STC analysis of the dysregulated lncRNAs among normal skin, scar, and MU. (A) The expression patterns of 8112 lncRNAs were analyzed and 16 model profiles were summarized. Each box represented a model expression profile, the upper number in the profile box was the model profile number and the lower one was the P-value. In total, 8 expression patterns of genes showed significant P-values (P<0.05; colored boxes). The profile No.3 (B) contained 684 genes and they were progressively downregulated from the normal skin to scar to MU tissues. The profile No.12 (C) contained 334 genes and they were progressively upregulated. The profile No.7 (D) contained 768 lncRNAs which were basically downregulated in MU. Letter A, B and C in the horizontal axis denoted normal skin, scar, and MU, respectively. The vertical axis displayed the time series of gene expression levels for the gene. MU – Marjolin ulcer; lncRNA – long noncoding RNA.

LncRNA-mRNA co-expression and function prediction

The co-expression network was composed of 5 lncRNAs and 334 coding genes (Figure 5). GO and pathway analysis for the co-expressed mRNAs revealed that they were highly enriched in cellular protein metabolic process (ontology: biological process), cytoplasm (ontology: cellular component), and RNA binding (ontology: molecular function), and that the p53 signaling pathway was the most highly enriched pathway (see Supplementary Figure 2). LncRNA NR_026995 was negatively correlated with YWHAZ, SESN3, PRKDC, CDKN1A, and CDK1 levels, while lncRNA TCONS_00007922 was positively correlated with SKP1, RRM2, RB1, CCND2, and CDC23 levels, and those mRNAs are mainly involved in the p53 signaling pathway and cell cycle. Intriguingly, TNFAIP1, SERPINB5, and RB1, which are well known as tumor suppressor, were all co-expressed with lncRNA ENST00000432694, TCONS_00007922, uc001oou.3, and ENST00000457834.

Figure 5

LncRNA-mRNA co-expression network in MU. The network is based on Pearson’s correlation coefficient. In this network, squares in yellow represent upregulated lncRNAs, squares in green represent downregulated lncRNAs, nodes represent mRNAs, solid lines mean positive correlations, while dashed lines mean negative correlations. MU – Marjolin ulcer; lncRNA – long noncoding RNA; mRNA – messenger RNA.

Construction of a ceRNA network

We predicted a ceRNA network consisting of 4 lncRNAs (lncRNA NR_026995 was filtered for reasons described in the methods section), 75 mRNAs, and 97 miRNAs (Figure 6). Some of the involved mRNAs have been reported to be associated with cancers such as MDM4, NF2, TP53I11, and IGF1. Particularly, downregulated lncRNA uc001oou.3 is predicted to be involved in ceRNA of miR-6849-3p targeting MDM4, JAK3, and STAT2 (the latter 2 were excluded from the network for their expression trend was contrary to lncRNA uc001oou.3), and both lncRNA uc001oou.3 and MDM4 were validated for a downregulation in MU. The ceRNA network might bring novel insight for the tumorigenesis mechanisms of MU. Of note, further validation of candidate ceRNA will be indispensable.

Figure 6

Competing endogenous RNA (ceRNA) network in MU. The ceRNA network was based on lncRNA-miRNA and miRNA-mRNA interactions. In this network, nodes in cyan represent lncRNAs, nodes in blue represent mRNA, nodes in red represent miRNA, and the edges represent sequence matching. MU – Marjolin ulcer; lncRNA – long noncoding RNA; micRINA – microRNA; mRNA – messenger RNA.

Discussion

MU is an aggressive malignancy, and the malignant degeneration of post-burn scars and lesions is an ineluctable eventuality. However, the exact pathogenesis of MU is not yet fully understood. Nevertheless, a tumor-promoting microenvironment caused by chronic inflammation should be largely incriminated [16]. Hitherto, lncRNA expression profiles have been widely used to reveal the underlying molecular mechanisms contributing to pathogenesis of many human diseases. Wan et al. have reported that lncRNAs play an important role in regulating the skin physiology function and diseases development [17]. Furthermore, they are also implicated in skin wound healing and hypertrophic scar formation [18,19]. As lncRNAs are expressed in a disease-, tissue- or developmental stage-specific manner, they are closely associated with malignant processes and so are attractive diagnostic and therapeutic targets [20]. A number of lncRNAs have been found to be responsible for cancer initiation and progression [21].

In the present study, for the first time, lncRNA and mRNA expression profiling in MU was determined by microarray analysis. A total of 7130 lncRNAs and 9867 mRNAs were differentially expressed. Among them, 2053 lncRNAs were dysregulated between the MU and scar group as well as the MU and skin group. Meanwhile, 1211 lncRNAs were found to be dysregulated specially between the MU and scar group while 3306 lncRNAs between the MU and skin group. These expression pattern suggested that those transcripts may be aberrantly expressed in different phase from normal skin to scar to MU, thus exert corresponding effects.

Cytoskeleton-associated protein 4 (CKAP4), the associated gene of lncRNA ENST00000552486 which was upregulated in MU, has received great attention recently for its association with proliferation and prognosis in some cancers [22]. In accordance with previous studies, our data also showed an upregulation of lncRNA TUG1 in MU. It is well known that lncRNA TUG1 is commonly overexpressed in numerous tumors, and associated with proliferation, migration, chemotherapy resistance, and poor prognosis in SCC [23,24]. Sand et al. revealed that lncRNAs were differentially expressed in cutaneous SCC, a number of them (e.g., lncRNA ENST00000524045, NR_036580, and TCONS_00010751) were also abnormally expressed in MU and showed the same expression trend [25], which verified our data. Surprisingly, lncRNA HOTAIR, overexpressed in many kinds of SCC and other tumors, is downregulated in MU, which suggests that the regulation mechanism of some lncRNAs in MU might differ from other kinds of SCC, considering that SCC accounts for the overwhelming majority of pathological types of MU.

To uncover differential expression patterns of lncRNAs among normal skin, para-cancerous scar, and malignant tissues, we conducted a STC analysis. Based on the results, the lncRNAs in profile No.3 and No.12 exhibited a gradual decreasing or increasing trend from normal skin to scar to MU, and those in profile No.7 were basically downregulated in MU. Such an expression change in MU is somewhat consistent with that found in oral SCC and hepatocellular carcinoma [26,27]. Coussens et al. concluded that inflammation is a critical component of tumor progression [28]; and increasing evidence also has revealed chronic inflammation can predispose individuals to cancer. Of note, long-term chronic inflammation also persists in pathological scars and chronic wounds which may lead to MU. Just like hepatocellular carcinoma usually develops in patients with chronic hepatitis or cirrhosis over several decades, the STC results suggest that scar tissue may be the intermediate step from normal skin to MU, i.e., MU is a long-term complication of burn. Furthermore, studies have indicated that some genes and proteins play a crucial role in benign to malignant transformation in some kinds of human tumors [29,30]. In agreement with the study of Yarmishyn et al., our results also showed certain lncRNAs might take part in malignant transformation of scar and the development of MU [31].

Subsequently, GO and pathway analysis were applied to infer the potential function of DE lncRNAs in MU. GO annotation revealed that the top 10 decreased biological processes between MU and para-cancerous scar tissue mainly belonged to the immune system, which substantiated that the decline of immune function was closely related to scar cancerization, and supported a previous theory explaining the occurrence of MU as well [2,32]. Furthermore, KEGG pathway analysis indicated the p53 signaling pathway was the most highly enriched in upregulated lncRNAs between MU and para-cancerous scar tissue followed by some other classic cancer pathways, suggesting that deregulated lncRNAs might play a pivotal role in scar hyperplasia and carcinogenesis through regulating those pathways. Previous studies have verified that lncRNAs regulate tumor cell cycle and affect tumor cell proliferation by corresponding mechanisms [24]. Our results indicated that those DE lncRNAs, especially the upregulated lncRNAs, perhaps contribute to scar hyperplasia and MU progression by primarily regulating the cell cycle.

To further decipher the function of some key lncRNAs in MU, the lncRNA-mRNA co-expression network was constructed. In this network, some co-expressed mRNAs, such as TGFBI, ITGA6, HIF1A, and CHRM3, have been proven to play an important role in the field of tumor biology in several kinds of SCC and other cancers [33-36]. IN our study, pathway analysis of the co-expressed mRNAs showed that the p53 signaling pathway and cell cycle were highly enriched, indicating that these 5 lncRNAs might play important roles in cell proliferation. The network also showed that lncRNA ENST00000432694, TCONS_00007922, ENST00000457834, and uc001oou.3, positively or negatively, interacted with protein-coding transcript TNFAIP1, SERPINB5, and RB1, all of which are known as tumor suppressors. Taking together, it is reasonable to infer that lncRNAs likely contribute to carcinogenesis of MU by regulating their co-expressed genes that are mainly involved in the p53 signaling pathway. Yet, since the underlying mechanisms of the interactions between lncRNAs and mRNAs remain unclear, more investigation is indispensable for validation.

The ceRNA is widely involved in many biological processes, and aberrations of ceRNA networks can lead to human diseases, including cancer [37]. LncRNAs can serve as natural miRNA sponges, causing inhibition of miRNA function. In fact, recent evidence has illustrated that the functionality of lncRNAs in tumorigenesis is partially mediated by ceRNA cross-talk [38]. Therefore, it may be helpful in further elucidating the function of DE lncRNAs in MU through identifying well-established miRNAs binding to the lncRNAs. By means of bioinformatic analysis, we conducted the ceRNA predication. MDM4, a well-studied tumor-related gene directly binds to p53 thus inhibiting its tumor suppressor activity, has been identified as a contributor to carcinogenesis of numerous tumors including SCC [39]. Previous studies confirmed that miR-191-5p, miR-887, and miR-661 could target MDM4 then augment p53 activity [40,41], and our analysis showed that both lncRNA uc001oou.3 and MDM4 harbored more than 2 MREs of miR-6849-3p, which was supposed to be a potential ceRNA mechanism involved in MU. Undoubtedly, the ceRNA network brings to light an unknown regulatory network in MU, and the network warrants further investigation.

Conclusions

Distinctive lncRNA and mRNA expression profiles were identified by high throughput screening in normal skin, scar, and MU tissues. Furthermore, a series of integrative bioinformatic approaches were applied to analyze the function of the lncRNAs. The results suggested that certain lncRNAs probably participated in malignant transformation of scar tissue and development of MU; the upregulated lncRNAs (including 5 selected lncRNAs) were highly involved in the p53 signaling pathway; and lncRNA uc001oou.3 might be implicated in ceRNA mechanism. Thus, our data not only provide new information regarding the potential role of lncRNAs in MU but also may lay a foundation for its diagnosis and therapy. However, the precise mechanism by which the lncRNAs operate in carcinogenesis of MU warrants further research.

Representative histological examination for patients with Marjolin ulcer: (A) normal skin tissue; (B) para-cancerous scar tissue; (C) Marjolin ulcer tissue.

GO and KEGG pathway analysis of the co-expressed mRNAs. (A) GO annotation for the co-expressed mRNAs, with the top 10 enrichment score-covering domains of biological processes, cellular components and molecular functions. (B) The top 10 enriched pathways for the co-expressed mRNAs. Enrichment score values were calculated as -log 10 (P values). GO – Gene Ontology; KEGG – Kyoto Encyclopedia of Genes and Genomes; mRNA – messenger RNA.

Supplementary Table 1

Clinical characteristics of patients with Marjolin ulcer.

Patient	Gender	Age*	Etiology	Localization	Latency period*	Ulceration period*	Histopathological type	Sample use
1	M	42	Flam burn	Right limb	40	1	Moderately differentiated SCC	Array
2	M	62	Scald	Right leg	30	10	Highly differentiated SCC	Array
3	F	49	Flame burn	Left leg	46	20	Highly differentiated SCC	Array
4	M	62	Snake bite	Right leg	30	1	Moderately differentiated SCC	qPCR
5	M	51	Scald	Left leg	40	3	Highly differentiated SCC	qPCR
6	F	33	Scald	Right leg	30	1	Highly differentiated SCC	qPCR
7	M	49	Snake bite	Left leg	40	1	Highly differentiated SCC	qPCR
8	M	47	Flam burn	Scalp	45	2	Highly differentiated SCC	qPCR
9	M	51	Scald	Left leg	40	3	Highly differentiated SCC	qPCR
10	M	54	Flame burn	Scalp	53	12	Highly differentiated SCC	qPCR
11	F	50	Flam burn	Scalp	50	3	Highly differentiated SCC	qPCR

Years;

SCC – squamous cell carcinoma.

Supplementary Table 2

Primers used for RT-qPCR validation in this study.

Gene	Forword (5′-3′)	Reverse (5′-3′)	Amplicon size (bp)
TCONS_00007922	TGTCAGCTCAGAAGACCAGG	ACCAGGACAGTTGTAGCACA	148
ENST00000457834	GGCCTTGAACTTTTCCACATC	GGACAGCACCAGCTCCTTTAT	165
ENST00000552486	GGAGCTGGCATTCTTACTCATT	GACAAGGGACCTAAAGATAGATGG	65
ENST00000414600	CATTTTGCCTTCTAGGTGCTTC	TTTGCTGGCACATAAGGGAGT	138
uc001oou.3	AAGTGCTGGGTGAATTCTTACT	TCACTGTCAGGGCCTTTACA	192
uc021uyg.1	CCTGACCTCCATCTATTTCTTC	TTTTGCGGATTCAGTTACATT	127
NR_026995	TCAGGGAGCTGCAGAATTCA	TAGGGAGTAGCTGAGTGGGT	130
ENST00000432694	CTGTAACACTTCCTGGCCTATTG	CCTGGTCACTTTGGTTTCCTC	55
PML	CCGACTTCTGGTGCTTTGAG	TTGGGGTTGGAGCAGAAGAT	172
SHC1	TGGGATGAGGAGGAGGAAGA	AATGTAGCTCCCAAGTGGCT	173
CHRM3	GACTACGTAGCCAGCAATGC	AAAGGAGATGACCCAAGCCA	150
MDM4	TTTCGTTGTTGGGCCTTGAG	TCTAACGACAGGAGCCGAAA	156
ACTB	CTGGAACGGTGAAGGTGACA	CGGCCACATTGTGAACTTTG	65