Bioinformatics analysis reveals the changes of peroxisome proliferator-activated receptor (PPAR) pathway in the development of Marjolin ulcers

Abstract

Objective: This study aimed to explore the molecular mechanisms underlying MU development.

Methods: Six MU and matched para-cancerous scar tissues were collected and subjected to transcriptome sequencing. Differentially expressed genes (DEGs) were screened, followed by function enrichment analysis. Gene set enrichment analysis (GSEA) was performed to investigate changes in pathways. Expression patterns of hub genes were confirmed by quantitative RT–PCR.

Results: A total of 1109 DEGs were identified between MU and para-cancerous scar tissues, including 471 upregulated and 638 downregulated genes. These DEGs mainly participated in extracellular matrix organization-related pathway, cytokine−cytokine receptor interaction, cytokine activity, and IL-17 signalling pathway. GSEA revealed 48 significantly altered pathways. Four common pathways were identified in different enrichment analyses, including PPAR signalling pathway. Among the genes in the PPAR pathway, 18 DEGs were found. Upstream transcription factors RXRG and PPARG were downregulated, while downstream target genes FABP6, OLR1, and MMP1 were upregulated, and FABP3, LPL, and ADIPOQ were downregulated. These results were confirmed by quantitative RT–PCR.

Conclusions: PPARγ-mediated PPAR signalling pathway is important during the development of MU. This study promotes future research on the molecular mechanisms of MU and provides novel insights for the prevention and treatment of this disease.

Highlights:

1. A total of 1109 DEGs were identified in MU tissue.

2. These DEGs were mainly involved in extracellular matrix organization and inflammation-associated pathways.

3. Various enrichment analyses identified four common pathways, including the PPAR signalling pathway.

KEYWORDS:

• Marjolin ulcer
• Scar squamous cell carcinomas
• PPAR signalling pathway

Introduction

Marjolin ulcers (MU) are a well-defined but rare cutaneous malignancy that usually occurs in traumatized skin, longstanding scars, and chronic wounds (Esther et al., 1999), and its malignant development generally originates in burn wounds. Many previous studies have reported that other types of non-healing wounds are also implicated in the malignant development of MU, such as traumatic wounds, osteomyelitis, and fistulas (Pavlovic et al., 2011). Well-differentiated squamous cell carcinoma (SCC) is the most common histopathologic type and tends to degenerate malignantly after approximately 25 years on average (Grenier et al., 1993; Baldursson et al., 1995). The slow transformation of MU results in inadequate treatment, which leads to poor prognosis (Pavlovic et al., 2011). Radical surgical excision remains the best choice of treatment (Xiang et al., 2019). However, surgical treatment is not suitable for some patients because of the considerable surgical difficulties and high risks (Xiang et al., 2019). Hence, the management and control of MU development is pivotal. In other words, there is an urgent need to investigate the molecular mechanisms of malignant degeneration and provide timely and effective therapeutic strategies.

The etiopathogenesis of MU has been investigated for over 100 years, but remains indistinct (Shah & Crane, 2022). Several etiological factors have been implicated in the malignant development of MU, such as immune deficiency in chronic scar tissues, which may be vulnerable and prone to malignant transformation (BOSTWICK et al., 1976; Kerr-Valentic et al., 2009). Long-standing chronic irritation and repeated stimulation of the wound may affect cell proliferation and accelerate spontaneous mutations. Chemokines, especially toxins released by surrounding necrotic tissues, may promote cell mutations (Treves, 1930). Mutations in genes involved in cell division and apoptosis are responsible for the increase in the incidence rate of carcinomas. Previous studies have reported that mutations in these genes are detected in patients with MU (Harland et al., 1997; Lee et al., 2000). The microenvironment of chronic wounds, which is modulated by the immune system, is also involved in the malignant development of MU (Bazaliński et al., 2017). In addition, the activity of matrix metalloproteinases and collagen and the function of the epithelium have been reported to be related to MU transformation (Bazaliński et al., 2017). Taken together, these studies indicate that the etiopathogenesis of MU is connected to multiple factors, including environmental conditions, immune systems, and genetics; however, the crucial element remains unclear.

RNA sequencing (RNA-seq) is a time-saving and high-throughput approach for transcriptome analysis that has promoted the precise investigation of eukaryotic transcriptomes (Wang et al., 2009). RNA-seq-based gene expression analysis has provided novel insights into the underlying molecular mechanisms of multiple diseases, such as osteoarthritis (Shi et al., 2019), pediatric acute myeloid leukemia (Shiba et al., 2019), coronary artery disease (Wang et al., 2021) and diabetes (Gong et al., 2017). In this study, we explored the transcriptome profiles of MU cancer and adjacent para-cancerous scar tissues. Genes that were differentially expressed in MU and their functions were investigated. In particular, the key pathways were screened, and the genes in these pathways were investigated. This study could serve as a reference for future research on the degeneration mechanism of MU and provide novel insights into the innovative prevention and therapy of MU.

Materials and methods

Sample collection

From May 2019 to August 2020, ten MU patients who received surgical treatments in the burn department of Beijing Jishuitan Hospital were recruited. Patients were included in this study when met the following two criteria: 1) patients were pathologically diagnosed; 2) patients with squamous cell carcinoma. Patients were excluded when met any of the following criteria: 1) patients with types of cancers other than squamous cell carcinoma; 2) patients who dropped out; 3) patients who failed to sign a consent form. Finally, six patients were included in this study. The details of the included patients were listed in Table S1. A total of 12 tissue samples consisting of six MU cancer tissues and six matched para-cancerous scar tissues were collected from the six MU patients. Informed consent was obtained from all patients. This study was conducted in accordance with the Declaration of Helsinki and was authorized by the Ethics Committee of Beijing Jishuitan Hospital (Approval No. 202203-69). All tissue samples were immediately frozen in liquid nitrogen and stored at −80°c after excision.

RNA extraction and library construction

Total RNA was isolated from all tissue samples using TRIzol® Reagent (Invitrogen), as described in the manufacturer's manual. RNA quantity and quality were assessed using a NanoDrop 2000 (Thermo Fisher Scientific) and 2100 Bioanalyzer (Agilent Technologies), respectively. Agarose gel electrophoresis was used to confirm RNA integrity. All RNA samples were stored at −80°C. The TruSeq® RNA Library Prep Kit v2 (Illumina) was used to construct the RNA-seq libraries, and subsequent paired-end sequencing of these libraries was conducted on the HiSeq/MiSeq Illumina sequencing platform.

Raw data filtering and transcriptome assembly

To obtain high-quality clean reads, SeqPrep software (Sturm et al., 2016) and Sickle software (Piel et al., 2017) were employed to remove the low-quality raw reads using the following four strategies: 1) reads with adapters and adapter-polluted reads were removed; 2) after trimming the low quality bases (Q < 30) from the 3 ‘end of the reads, there were still bases with quality score Q less than 10 in the reads, and the reads were excluded; 3) reads with a frequency of unrecognized bases of more than 10% were filtered out; and 4) reads with a length of less than 50 base pairs after trimming were removed. The remaining reads were considered as clean reads, and the quality of these reads was evaluated using fastx_toolkit_0.0.14 software (version: 0.0.14). The clean reads were aligned with the genome of Homo sapiens (GRCh38.p13), which was regarded as the reference genome, using HISAT software (version 2.1.0) (Kim et al., 2015) under default parameters. StringTie software (version 1.3.3b) (Pertea et al., 2015) was used to further reconstruct the transcriptome of the RNA-seq reads. Protein databases, including NCBI_NR, Swiss-Prot, Pfam, COG, GO, and KEGG, were used to annotate the assembled transcriptome.

Identification of DEGs

To obtain the expression profiles, the read counts of individual genes were calculated using RSEM software (version: 1.3.1) (Li & Dewey, 2011). The transcripts per kilobase million (TPM) were calculated based on the read counts, which represent the expression level of each gene. The DEGs between the cancer tissue samples and the para-cancerous scar tissue samples were identified on the basis of the TPM via the edgeR R package (version: 3.28.0) (Robinson et al., 2010). The absolute value of the logarithmic fold change (|logFC|) >2 and a p-value <0.05 indicated that the corresponding gene was differentially expressed between the cancer tissue samples and the para-cancerous scar tissue samples.

Functional enrichment analysis

KEGG (Kanehisa & Goto, 2000) and GO (2019) are two powerful and pragmatic databases that uncover the biological functions of the genes of interest. The DEGs between the cancer tissue samples and para-cancerous scar tissue samples were submitted to the KEGG and GO databases for enrichment analysis using the clusterProfiler tool (version: 3.16.0) (Yu et al., 2012). Enriched gene counts ≥2 and p <0.05 were considered to be significantly enriched pathways.

Gene set enrichment analysis (GSEA)

Based on the reference KEGG pathway gene sets (c2.cp.kegg.v7.2.symbols.gmt) from the MSigDB database, GSEA was performed for all genes (ranked by LogFC descending) using the clusterProfiler tool, followed by Benjamini-Hochberg corrections. KEGG pathways with an adjusted p <0.05 was selected.

Construction of the PPI network

To explore the protein–protein interactions (PPI) network in the intricate biological system, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (version: 11.0) (Keshava Prasad et al., 2009) was employed to uncover potential interactions among the DEGs. The threshold of the PPI confidence score was set to 0.4, which is a commonly used value. Based on the data provided by the STRING database, the Cytoscape software (version 3.6.1) (Shannon et al., 2003) was used to construct the PPI network.

Verification of the transcriptome analysis

Total RNA was isolated from all tissue samples using TRIzol® Reagent (Invitrogen), as described in the manufacturer's manual. First-strand cDNA was synthesized using the PrimeScript™ RT Master Mix (Perfect Real Time) (Takara Biomedical Technology), as described in the manual. Gene-specific primers are listed in Table S2. PowerUp SYBR Green PCR Master Mix (Thermo) was used to conduct quantitative PCR as described in the manual. Quantitative PCR was carried out on an Applied Biosystems 7900HT Fast Real-Time PCR system (Applied Biosystems, USA) with 10-µL volumes comprised of 5 µL PowerUp SYBR Green PCR Master Mix, 0.2 µL forward primer, 0.2 µL reverse primer, 1 µL cDNA, and 3.6 µL RNase free water. The assay was conducted in triplicate for each sample. GAPDH was used as a reference gene.

Statistical analysis

R version 4.0.5 was used for all bioinformatics analyses. Quasi-likelihood (QL) F-test provided in edgeR package was used for differential expression analysis. DEGs was selected using |log2 FC|>2 and P <0.05. The ‘ggplot2’package and ‘pheatmap’ package were used to generate volcano plot and heatmap, respectively. Benjamini-Hochberg method was used for multiple tests correction in GSEA. Data in quantitate RT–PCR was presented as mean ± standard deviation, and were statistically analyzed using GraphPad Prism 8.0. Expression of genes between two groups were compared using t-test. P < 0.05 was used to show statistical significance.

Results

Mapping results of the RNA-seq data

A total of 606,820,478 paired-end raw reads were obtained for the 12 tissue samples, ranging from 48,738,424–55,553,030 for each tissue sample (Table 1). Among these reads, 308,404,500 raw reads were generated from the six cancer tissue samples, and 298,415,978 raw reads belonged to the six para-cancerous scar tissue samples. On average, 50.6 million raw reads were obtained from each tissue sample. Only 0.7% of the raw reads were ruled out by filtration of the raw data, and the Q30 of each sample was more than 93%. These results suggest that our data were of high quality. Based on the alignment using HISAT software, 95.38% of the clean reads were mapped to the reference genome of Homo sapiens (GRCh38.p13), ranging from 94.16–96.29% of the reads for each sample. Additionally, 91.32% of the clean reads were uniquely mapped to the reference genome (Table 1).

Table 1. Descriptive statistics of RNA-seq data.

Sample	Raw reads	Raw bases	Clean reads	Clean bases	Error rate (%)	Q20^a (%)	Q30^b (%)	GC content (%)	Multiple mapped	Uniquely mapped
AB_125914	48738424	7359502024	48451418	7220055275	0.0249	98.15	94.2	52.31	2303006(4.75%)	43086344(88.93%)
AB_702122	50637594	7646276694	50283434	7489620076	0.0252	97.99	93.85	52.26	1908443(3.8%)	46325788(92.13%)
AB_718697	47730986	7207378886	47411482	7072482876	0.0251	98.03	93.94	52	1464905(3.09%)	44059999(92.93%)
AB_720865	52415652	7914763452	52073008	7753408209	0.0249	98.13	94.15	52.11	2880630(5.53%)	45098071(86.61%)
AB_722938	55553030	8388507530	55215120	8235631436	0.0251	98.04	93.94	51.82	1896353(3.43%)	51229213(92.78%)
AB_723158	53328814	8052650914	52996782	7892072112	0.0251	98.04	93.92	51.1	1909449(3.6%)	49040255(92.53%)
ZC_125914	51318286	7749061186	50932278	7602951165	0.0251	98.04	94	52.59	2265978(4.45%)	46385964(91.07%)
ZC_702122	50757772	7664423572	50330206	7482279273	0.0251	98.02	93.94	51.91	2150774(4.27%)	46311272(92.01%)
ZC_718697	52285000	7895035000	51925416	7677151577	0.0264	97.57	92.7	51.62	2094793(4.03%)	46799812(90.13%)
ZC_720865	48981174	7396157274	48572384	7219646457	0.0254	97.92	93.64	51.4	1852065(3.81%)	44912198(92.46%)
ZC_722938	47692856	7201621256	47267918	7037744544	0.0252	98	93.89	52.54	1626878(3.44%)	43774484(92.61%)
ZC_723158	47380890	7154514390	46926272	6973181076	0.0256	97.85	93.53	51.95	2092785(4.46%)	43023225(91.68%)

Note: ^a: Q20, Phred Quality Score 20, probability of incorrect base call 1 in 100;

^b: Q30, Phred Quality Score 30, probability of incorrect base call 1 in 1000.

Samples with name begin with AB are cancer tissue samples, while samples with name begin with ZC are para-cancerous scar tissue samples.

Identification of DEGs

A total of 1109 DEGs were identified between the cancer tissue samples and the para-cancerous scar tissue samples, consisting of 471 upregulated genes and 638 downregulated genes (Figure 1A, Table S3). The heat map of hierarchical clustering analysis visually showed the detailed expression levels of the DEGs in different samples (Figure 1B).

Figure 1. Identification of the Differentially Expressed Genes. A) Volcano plots of differentially expressed genes (DEGs) between the cancer tissues and the para-cancerous scar tissues. The red and green dots indicate the significantly upregulated and downregulated genes, respectively. B) Hierarchical clustering analysis of the DEGs. AB represents the MU tumor tissue samples, while ZC represents the matched para-cancerous scar tissues.

Functional enrichment for DEGs

In order to investigate the biological functions of these DEGs, enrichment analysis was performed to explore their involved GO annotation terms, and 444 biological process terms, 99 cellular component terms and 145 molecular function terms were obtained (Table S4). The top 10 enriched biological process terms contained muscle system process, muscle contraction, striated muscle cell development and so on (Figure 2A), and the top 10 enriched cellular component terms included contractile fiber, myofibril, sarcomere, sarcoplasm and collagen−containing extracellular matrix and so on (Figure 2B). The top 10 enriched molecular function terms are shown in Figure 2C, including signaling receptor activator activity, receptor ligand/ cytokine/ chemokine/ growth factor activity, etc (Figure 2C). In addition, 42 KEGG pathways were significantly enriched (Table S4). The top 10 KEGG pathways are shown in Figure 2D, including cytokine−cytokine receptor interaction, cardiac muscle contraction, the interleukin (IL)−17 signalling pathway, and the PPAR signalling pathway.

Figure 2. Functional enrichment analysis of the DEGs. A) The top 10 enriched biological processes terms; B) The top 10 enriched cellular component terms; C) The top 10 enriched molecular function terms; D) The top 10 enriched KEGG pathways. The x-axis indicates the number of DEGs, and the y-axis represents the enriched pathways. The adjusted p-value was represented by the color key on the right.

GSEA

GSEA revealed that 48 KEGG pathways showed significant differences between MU and para-cancerous scar tissues (Table S5), including 18 negatively correlated pathways (NES <0) and 30 positively correlated pathways (NES >0). Among them, the cell cycle, cytokine-cytokine receptor interaction, and primary immunodeficiency pathways were the top three positively correlated pathways, whereas cardiac muscle contraction, dilated cardiomyopathy, and hypertrophic cardiomyopathy were the top three negatively correlated pathways (Table 2).

Table 2. The significant enriched top 5 positively and negatively correlated KEGG pathways in GSEA.

Description	setSize	enrichmentScore	NES	pvalue	p.adjust
KEGG_CELL_CYCLE	124	0.632190771	2.3079156	0.0024691	0.015144
KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	262	0.534936818	2.1392087	0.0029155	0.0163512
KEGG_PRIMARY_IMMUNODEFICIENCY	35	0.704029322	2.076579	0.0020921	0.015144
KEGG_GRAFT_VERSUS_HOST_DISEASE	37	0.697821921	2.0708643	0.0021186	0.015144
KEGG_P53_SIGNALING_PATHWAY	67	0.617181805	2.0595748	0.002193	0.015144
KEGG_CARDIAC_MUSCLE_CONTRACTION	77	−0.703274411	−2.267061	0.0017921	0.015144
KEGG_DILATED_CARDIOMYOPATHY	90	−0.622646456	−2.056833	0.0017857	0.015144
KEGG_HYPERTROPHIC_CARDIOMYOPATHY_HCM	83	−0.616946677	−2.000995	0.0018083	0.015144
KEGG_CALCIUM_SIGNALING_PATHWAY	177	−0.550730098	−1.97924	0.0016502	0.015144
KEGG_DRUG_METABOLISM_CYTOCHROME_P450	62	−0.633979472	−1.958604	0.0018416	0.015144

Identification of key pathways and DEGs

Based on the enrichment results of GO, KEGG, and GSEA analyses, Pathview Web (Luo et al., 2017) was employed to explore the common significant enrichment pathway. Four pathways, including the hsa03320 PPAR signalling pathway, hsa04260 cardiac muscle contraction, hsa04020 calcium signalling pathway, and hsa05414 dilated cardiomyopathy, were identified (Figure 3A).

Figure 3. Identification of key pathways and DEGs. A) GSEA enrichment score curves of the four common pathways. The top panel shows the enrichment scores of each pathway. The middle vertical bars display the position of each hit. The bottom panel indicates the ranking metric scores. B) Protein-protein interactions between the important pathways and genes. The different pathways are represented by different big light color ellipses and the hub genes are represented by the small dark color ellipses. The grey lines indicated the associations between the hub genes.

PPI analysis was performed for all genes in these four pathways. As shown in Figure 3B, the genes in the calcium signalling, cardiac muscle contraction, and dilated cardiomyopathy pathways showed close interactions with each other, and ATP2A1, CACNA1S, and TNNC1 were common between these pathways. No genes were common to either the PPAR signalling pathway or the other three pathways. Increasing evidence indicates that the PPAR signalling pathway could be a potential molecular target for the prevention or treatment of SCC (Ren & Konger, 2019). Peroxisome proliferator activated receptor gamma (PPARG), a crucial gene of PPAR pathway, has been linked to mediate the invasion, differentiation and tumor growth of various SCC (Masuda et al., 2005; Nijsten et al., 2005; Konger et al., 2019). Moreover, change of PPARG expression has been found to drive multiple molecules to involve the chemosensitivity of hypopharyngeal SCC (Lian et al., 2020; Lian et al., 2021). Therefore, we focused on the PPAR pathway in this study.

Among the genes in the PPAR signalling pathway, 18 were differentially expressed in MU (Figure 4). Fatty acid-binding protein 6 (FABP6), oxidized low-density lipoprotein receptor 1 (OLR1), and matrix metallopeptidase 1 (MMP1) were upregulated in MU. The other genes were downregulated in MU, including FABP3, lipoprotein lipase (LPL), retinoid X receptor gamma (RXRG), PPARG, and adiponectin, C1Q, and collagen domain containing (ADIPOQ). Of these, FABP3, LPL, ADIPOQ, OLR, and MMP1 are downstream target genes, while RXRG and PPARG might be upstream transcription factors in the PPAR signalling pathway.

Figure 4. Identification of key pathways and DEGs. A) The diagram of the PPAR signalling pathway. B) The heat map of the DEGs in the PPAR signalling pathway. AB represents the MU tumor tissue samples, while ZC represents the matched para-cancerous scar tissues.

Validation of hub genes in the PPAR signalling pathway

We performed quantitative RT–PCR to confirm the changes in the expression of several genes in the PPAR signalling pathway, including OLR1, FABP3, PPARG, MMP1, RXRG, ADIPOQ, and LPL. As shown in Figure 5, LPL, RXRG, PPARG, ADIPOQ, and FABP3 were significantly downregulated in MU cancer tissues, while the expression level of OLR1 and MMP1 were significantly upregulated in cancer tissues. The expression patterns of these genes were consistent with the transcriptome data.

Figure 5. qRT-PCR analysis of the hub genes in the PPAR signalling pathway. GAPDH was used as the reference gene. ** indicates that the gene expression was significantly different between the MU and the para-cancerous scar tissues (p <0.01).

Discussion

In the current study, we performed systematic transcriptome analysis based on RNA-seq data of MU cancer tissues and para-cancerous scar tissues. A total of 1109 DEGs, including 471 upregulated and 638 downregulated genes, were identified. Further functional enrichment analysis revealed that these DEGs mainly participated in the extracellular matrix organization-related pathway, cytokine−cytokine receptor interaction, cytokine activity, and the IL-17 signalling pathway. Chronic inflammation is considered a common etiological factor for MU (Pekarek et al., 2011). Severe inflammation causes DNA damage and mutagenic assault, which are responsible for tumor transformation (Coussens & Werb, 2002). Similar with our findings, previous studies also revealed a large-scale of expression changes of genes in the development of MU (Sinha et al., 2017; Liu et al., 2018). Sinha et al. identified 665 DEGs between MU tumor cells and epithelial cells and 1673 DEGs between SCC cells and epithelial cells (Sinha et al., 2017). Further functional enrichment analysis of these DEGs suggested that mainly three pathways were affected, including the extracellular matrix organization/degradation-related pathways, DNA damage-related pathways, and cytokine signalling pathways (Sinha et al., 2017). Liu et al. revealed that there were 9867 genes and 7130 long noncoding RNAs differentially expressed among normal skin, scar, and MU tissues, and they were mainly implicated in p53 signal pathway, cell cycle and immune response (Liu et al., 2018). These results suggest that there was a large-scale of changes at molecular level during the scar carcinogenesis and development of MU, and the dysregulation at the molecular level might contribute to scar carcinogenesis and development of MU by regulating various important biological functions and pathways.

Based on the integration of the enrichment results, four KEGG pathways, including the PPAR signalling pathway, were obtained. The peroxisome proliferator-activated receptor (PPAR) is an important family of nuclear transcription factors (TFs), including PPARα, PPARβ/δ, and PPARγ (Berger and & Moller, 2002; Michalik et al., 2006). A variety of genes are regulated by PPARs, including sugar-, lipid-, and energy metabolism-associated genes (Fucci et al., 2012). Four PPARG transcriptional splice variants under the control of the corresponding promoters were identified: PPARG1–4 (Fucci et al., 2012). Interestingly, the proteins encoded by PPARG1, PPARG3, and PPARG4 were identical (PPARγ1), whereas the PPARγ2 protein was encoded by PPARG2 (Sahu et al., 2012). To bind to DNA targets, PPARγ forms a heterodimer with the retinoid X receptor (RXR) (Ren & Konger, 2019). PPARγ is considered essential for regulating the expression of genes associated with cell proliferation, differentiation, and inflammatory responses in human and mouse skin (Man et al., 2006). To bind to DNA targets, PPARγ forms a heterodimer with the retinoid X receptor (RXR) (Ren & Konger, 2019). Later studies have demonstrated that PPARγ agonists play a protective role in the development of oral cavity and oesophageal SCC (Sawayama et al., 2014; Tseng, 2014). Our transcriptional analysis and quantitative RT–PCR results suggested that PPARG and its potential chaperone RXRG were extremely downregulated in MU tissue, suggesting that dysregulation of the PPAR signalling pathway is involved in the progression of MU.

During epidermal differentiation, lipids are synthesized in keratinocytes and squeezed into the extracellular region to form lipid-enriched extracellular layers, which represent a part of the skin barrier (Jensen & Proksch, 2009). Skin diseases can occur when the skin barrier is damaged. Among the genes in the PPAR signalling pathway, 18 genes were differentially expressed in MU, including the downstream target genes FABP3, LPL, and OLR1. FABP3, a member of the FABP family, has been reported to participate in the uptake, intracellular metabolism, and transport of long-chain fatty acids (Kleine et al., 1992). In the PPAR signalling pathway, FABP3 is the target of PPAR that is involved in the lipid metabolism (Tsukahara et al., 2014). LPL encodes lipoprotein lipase, a crucial lipid-modulating enzyme that participates in lipid metabolism and transport (Xu et al., 2018). OLR encodes a low-density lipoprotein receptor that regulates lipid metabolism (Khaidakov et al., 2011). Therefore, we speculated that dysregulation of the PPAR signalling pathway plays a role in the progression of MU by regulating lipid metabolism-related genes that damage the skin barrier.

This study has some limitations. 1) the low number of samples was the major limitation of this study, with only six samples in each group. We advocate that analysis based on a high number of samples may lead to more convincing conclusions, in spite of this, many studies were performed using a small sample size (Chen et al., 2020; Xu et al., 2021). 2) Only MU cancer tissues and matched para-cancerous scar tissues were analyzed, and normal skin should also be analyzed together to obtain more comprehensive and reliable results. 3) This study proposed that the PPAR signalling pathway might be important in the development and progression of MU. However, the role of the PPAR signalling pathway has not been effectively investigated and further study is required.

In conclusion, systematic transcriptome analysis of cancer tissues and para-cancerous scar tissues of MU was performed. Changes in the expression patterns of large-scale genes were observed in MU. These genes are involved in the extracellular matrix organization and inflammation-associated pathways. This study found that the PPARγ-mediated PPAR signalling pathway is a key pathway in the development and progression of MU.

Authors’ contributions

Cheng Wang and Lin Cheng conceived and designed the study. Ying Zhang and Huijun Zhang conducted the data analysis and interpretation. Xiaozhuo Zhao participated in the study design and performed the statistical analysis. Yuming Shen and Cheng Wang conceived the study, participated in its design and coordination, and helped draft and revise the manuscript for important intellectual content. All authors have read and approved the final manuscript.

Data availability statement

The raw sequencing data used in this study were deposited in the National Center for Biotechnology Information (NCBI) Sequence ReadArchive (SRA) database with accession number PRJNA821549 (http://www.ncbi.nlm.nih.gov/bioproject/821549). The data for RT-qPCR has been uploaded to 4TU. ResearchData repositories with doi:10.4121/19375049 (https://figshare.com/s/fa15a0a8ef60ee34f0e8).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Beijing Jishuitan Hospital Elite Young Scholar Program (Program no. XKGG202114), and by Beijing JST Research Funding (Program No. ZR202208).