# Putative Genes and Pathways Involved in the Acne Treatment of Isotretinoin via Microarray Data Analyses

Authors: Biao Chen; Peishan Li; Jun Li; Jinping Chen
Journal: BioMed Research International
 (2020)

---

## Abstract

Acne is the eighth most common disease worldwide. Disease burden of acne such as anxiety, reduced self-esteem, and facial scarring lowers the life quality of acne patients. Isotretinoin is the most potent treatment for moderate-severe acne. However, the adverse events of isotretinoin especially teratogenicity limit its use. This study aims at investigating the therapeutical mechanisms of isotretinoin using bioinformatics analysis. Differentially expressed genes (DEGs) were filtered from microarray datasets GSE10432, GSE10433, and GSE11792. Functional and pathway enrichment analyses of DEGs were performed. Protein–protein interaction (PPI) network and module analyses were also conducted based on DEGs. Using isotretinoin for 1 week,LCN2, PTGES, and GDF15 were upregulated and might mediate sebocytes apoptosis and thus decreased sebum production; CCL2 originated from activated TNF signaling pathway and S100A7 could be related with “acne-flare”. While treating with isotretinoin for 8 weeks, key genes were downregulated, including HMGCS1, HMGCR, FDFT1, MVD, IDI1, and FDPS, which may be associated with decreased sebum synthesis; HMGCS1, HMGCR, and FDFT1 also probably associated with apoptosis of sebocytes. There were only two common genes including ACSBG1 and BCAT2 which worked in both 1 week and 8 weeks and could associate with decreased sebum synthesis and apoptosis of sebocytes, respectively. These results indicate potential therapeutics and side effect mechanisms of isotretinoin in the acne treatment and provide a research direction to further investigate the therapeutic mechanism of isotretinoin and thus develop retinoid-like compounds with similar curative effect and without teratogenicity.


---

## Body

## 1. Introduction

Acne is the eighth most common disease worldwide with a morbidity of 94% according to the Global Burden of Disease Project [1]. Disease burden of acne includes anger, depression, anxiety, reduced self-esteem, facial scarring, and impaired social interaction, which usually negatively impact the quality of life of acne patients [2]. Treatments of acne include topical medication such as antibiotics, benzoyl peroxide, and retinoids; and systemic medication such as isotretinoin, tetracyclines, and spironolactone [3]; as well as other miscellaneous therapies including chemical peels, intralesional injection of triamcinolone acetonide, and laser and light therapies [4]. Among them, isotretinoin is the first line and the most potent treatment for moderate-severe acne [4, 5]. The adverse events of isotretinoin include elevated liver enzymes [6] and serum lipids, xerosis, cheilitis, eye dryness, irritation and conjunctival injection, fatigue, depressed mood, arthralgia, and teratogenicity [7]. Among them, teratogenicity is the main and disastrous side effect of isotretinoin. Considering the therapeutic effect and adverse events, there is an emerging need to explore the mechanism of acne treatment with isotretinoin, and thus develop retinoid-like compounds with similar curative effect and without teratogenicity.Microarray is a high-throughput technology which exhibits the overall gene expression data of recruited specimens from different diseases. NCBI-Gene Expression Omnibus database (NCBI-GEO) (https://www.ncbi.nlm.nih.gov/geo) is a free public website and archives for functional genomics datasets, in which microarray data are usually deposited and available. Nelson et al. performed microarray analyses of skin biopsies from acne patients treated with isotretinoin for 1 week and cultured human sebaceous gland cells (SEB-1 sebocytes) incubated with isotretinoin for 72 hours and found that lipocalin 2 is the most highly upregulated gene induced by isotretinoin, which mediates the apoptosis of sebaceous glands of acne patients [8]. Nelson et al. also conducted microarray analyses of skin biopsies from acne patients treated with isotretinoin for 8 weeks and found that isotretinoin significantly reduces the size of the sebaceous gland at 8 weeks with a trend observed at 1 week, and only 3 common DEGs were found in samples of acne patients treated with isotretinoin at 8 weeks and 1 week, indicating that isotretinoin induces temporal changes in gene expression of acne patients’ skin [9]. Bioinformatical analyses of multiple microarray datasets derived from different studies contribute to identify the hub genes and pathways which participate in the acne treatment of isotretinoin, thus provide potential therapeutic targets of acne.In the present study, three microarray data including GSE10432, GSE10433, and GSE11792 were obtained from NCBI-GEO. Common DEGs between GSE10432 and GSE10433, DEGs in GSE11792 were filtered via the online tool GEO2R. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs were conducted using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/). PPI networks of DEGs were constructed via the free public online website STRING (https://string-db.org/, version 11) and free software Cytoscape; in which the hub genes and modules were further assessed to demonstrate the functions of the DEGs. This study aimed at obtaining deep insights into the therapeutic mechanism of acne treatment using isotretinoin.


## 2. Materials and Methods

### 2.1. Microarray Data and Identification of DEGs

GSE10432, GSE10433, and GSE11792 were downloaded from NCBI-GEO to identify the DEGs between SEB-1 treated with vehicle control and isotretinoin or between acne patients before and after treatment with isotretinoin. The microarray profile GSE10432 was based on the GPL8300 platform (Affymetrix Human Genome U95 Version 2 Array, Palo Alto, CA, USA) and comprised 3 biological replicates of SEB-1 treated with isotretinoin or vehicle control (Submission date: Feb 07, 2008) [8]. GSE10433 was based on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array, Palo Alto, CA, USA) and included 6 acne patients at baseline and treated with isotretinoin for 1 week (Submission date: Feb 07, 2008) [8]. GSE11792 was based on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array, Palo Alto, CA, USA) and consisted of 8 acne patients at baseline and treated with isotretinoin for 8 weeks (Submission date: Jun 16, 2008) [9].The DEGs in GSE10432, GSE10433, and GSE11792 were analyzed using the online tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). The adjusted P value and ∣logFC∣ were calculated; and P<0.05 and ∣logFC∣≥0.263034405833794, 0.263034405833794, and 0.584962501 were used as a cut-point for the statistically significant DEGs in GSE10432, GSE10433, and GSE11792, respectively. The Benjamini and Hochberg false discovery rate method was used as a correction factor for the adjusted P value in GEO2R.


### 2.2. GO and Pathway Enrichment Analyses of DEGs in Acne Treated with Isotretinoin

GO analyses, consisting of molecular function (MF), cellular component (CC), and biological process (BP), were used to elucidate gene functions. DAVID is a free online website available for analyzing the biological meaning behind a large list of genes. In our study, GO and KEGG pathway enrichment analyses of significant DEGs, including the common DEGs between GSE10432 and GSE10433 and the DEGs in GSE11792, were conducted by DAVID, with a threshold ofP value < 0.05 and enrichment gene count>2.


### 2.3. Establishment of PPI Network and Modular Analysis

The common DEGs between GSE10432 and GSE10433 and DEGs in GSE11792 were analyzed using the online tool STRING, with 0.400 (medium confidence) and 0.700 (high confidence) as the minimum required interaction score, respectively. The PPI network of the above DEGs was established using Cytoscape, respectively. The node degree of the PPI network was calculated by Network Analyzer in Cytoscape. The hub genes of these PPI networks were filtered by CytoHubba, which is one of the Apps in Cytoscape. Modular analysis of the PPI network of GSE11792 was conducted using MCODE, with aDegreeCutoff=2, NodeScoreCutoff=0.2, K−Core=2, and MaxDepth=100. Finally, KEGG pathway enrichment analyses of these candidate genes in every cluster of PPI networks of GSE11792 were conducted by DAVID.


## 2.1. Microarray Data and Identification of DEGs

GSE10432, GSE10433, and GSE11792 were downloaded from NCBI-GEO to identify the DEGs between SEB-1 treated with vehicle control and isotretinoin or between acne patients before and after treatment with isotretinoin. The microarray profile GSE10432 was based on the GPL8300 platform (Affymetrix Human Genome U95 Version 2 Array, Palo Alto, CA, USA) and comprised 3 biological replicates of SEB-1 treated with isotretinoin or vehicle control (Submission date: Feb 07, 2008) [8]. GSE10433 was based on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array, Palo Alto, CA, USA) and included 6 acne patients at baseline and treated with isotretinoin for 1 week (Submission date: Feb 07, 2008) [8]. GSE11792 was based on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array, Palo Alto, CA, USA) and consisted of 8 acne patients at baseline and treated with isotretinoin for 8 weeks (Submission date: Jun 16, 2008) [9].The DEGs in GSE10432, GSE10433, and GSE11792 were analyzed using the online tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). The adjusted P value and ∣logFC∣ were calculated; and P<0.05 and ∣logFC∣≥0.263034405833794, 0.263034405833794, and 0.584962501 were used as a cut-point for the statistically significant DEGs in GSE10432, GSE10433, and GSE11792, respectively. The Benjamini and Hochberg false discovery rate method was used as a correction factor for the adjusted P value in GEO2R.


## 2.2. GO and Pathway Enrichment Analyses of DEGs in Acne Treated with Isotretinoin

GO analyses, consisting of molecular function (MF), cellular component (CC), and biological process (BP), were used to elucidate gene functions. DAVID is a free online website available for analyzing the biological meaning behind a large list of genes. In our study, GO and KEGG pathway enrichment analyses of significant DEGs, including the common DEGs between GSE10432 and GSE10433 and the DEGs in GSE11792, were conducted by DAVID, with a threshold ofP value < 0.05 and enrichment gene count>2.


## 2.3. Establishment of PPI Network and Modular Analysis

The common DEGs between GSE10432 and GSE10433 and DEGs in GSE11792 were analyzed using the online tool STRING, with 0.400 (medium confidence) and 0.700 (high confidence) as the minimum required interaction score, respectively. The PPI network of the above DEGs was established using Cytoscape, respectively. The node degree of the PPI network was calculated by Network Analyzer in Cytoscape. The hub genes of these PPI networks were filtered by CytoHubba, which is one of the Apps in Cytoscape. Modular analysis of the PPI network of GSE11792 was conducted using MCODE, with aDegreeCutoff=2, NodeScoreCutoff=0.2, K−Core=2, and MaxDepth=100. Finally, KEGG pathway enrichment analyses of these candidate genes in every cluster of PPI networks of GSE11792 were conducted by DAVID.


## 3. Results

### 3.1. The Volcano Plots of Genes Expression and the Common DEGs of GSE10432, GSE10433, and GSE11792

Three microarray datasets including GSE10432, GSE10433, and GSE11792 were downloaded from GEO. Basing on a threshold ofP value < 0.05 and ∣logfoldchangeFC∣≥0.263034405833794, namely, ∣FC∣≥1.2, 638 upregulated and 775 downregulated genes in GSE10432, and 176 upregulated and 142 downregulated genes in GSE10433 were obtained. There were 62 upregulated and 249 downregulated genes in GSE11792 with a cut off value of P value < 0.05 and ∣logFC∣≥0.584962501, namely, ∣FC∣≥1.5. The volcano plots of gene expression of GSE10432 (Figure 1(a)), GSE10433 (Figure 1(b)), and GSE11792 (Figure 1(c)) were generated using a free online website imageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/Volcanoplot.html). There were 27 upregulated (Figure 1(d)) and 9 downregulated (Figure 1(e)) genes in GSE10432 and GSE10433 concurrently (Table S1). 62 upregulated genes and 249 downregulated genes were observed in GSE11792 (Table S2). There were consistently 0 upregulated (Figure 1(f)) and 2 downregulated (Figure 1(g)) genes, including ACSBG1 and BCAT2, among GSE10432 (SEB-1 sebocytes incubated with isotretinoin for 72 hours), GSE10433 (acne patients treated with isotretinoin for 1 week), and GSE11792 (acne patients treated with isotretinoin for 8 weeks).Figure 1
Volcano plots of genes expression and the common DEGs in GSE10432, GSE10433, and GSE11792. Volcano plots of genes expression in microarray data GSE10432 (a), GSE10433 (b), and GSE11792 (c). In volcano plots, green represents downregulated genes, red represents upregulated genes, and gray means no significant DEGs. UP: upregulated DEGs; DW: downregulated DEGs; NoDiff: nondifferentially expressed genes. (d) The common upregulated DEGs of GSE10432 and GSE10433. (e) The common downregulated DEGs of GSE10432 and GSE10433. (f) The common upregulated DEGs among GSE10432, GSE10433, and GSE11792. (g) The common downregulated DEGs among GSE10432, GSE10433, and GSE11792. Different colors indicate different datasets. The cross part represents the common DEGs.
(a)(b)(c)(d)(e)(f)(g)


### 3.2. GO, KEGG Pathway Enrichment Analyses and PPI Network of the Common DEGs in GSE10432 and GSE10433

GSE10432 and GSE10433 represented that acne and SEB-1 sebocytes treated with isotretinoin for 1 week and 72 hours, respectively. GO and KEGG pathway enrichment analyses of the consistently upregulated or downregulated DEGs in GSE10432 and GSE10433 were clustered via DAVID. Figures2 and 3 showed the results of GO and KEGG pathway enrichment analyses. As far as the upregulated DEGs of GSE10432 and GSE10433 were concerned, the upregulated DEGs mainly participated in cell adhesion and inflammatory response in term of BP (Figure 2(a)); located in the cytoplasm and extracellular exosome in term of CC (Figure 2(b)); involved in protein binding in term of MF (Figure 2(c)); enriched in TNF signaling pathway in term of KEGG (Figure 2(d)); all details can be seen in Table S3. In the matter of the downregulated DEGs of GSE10432 and GSE10433, the downregulated DEGs mainly involved in long-chain fatty-acyl-CoA biosynthetic process and definitive hemopoiesis in term of BP (Figure 3(a)); located in the external side of the plasma membrane and focal adhesion in term of CC (Figure 3(b)); enriched in fatty acid metabolism in term of KEGG (Figure 3(c)); all details can be seen in Table S4. The PPI network of the common DEGs in GSE10432 and GSE10433 was showed in Figure 3(d), in which the hub genes included CCL2, LCN2, S100A7, PTGES, and GDF15 (Table 1), and CCL2 participated in TNF signaling pathway (Table S3).Figure 2
GO and KEGG pathway enrichment analyses of the upregulated DEGs in GSE10432 and GSE10433. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)Figure 3
GO and KEGG pathway enrichment analyses of the downregulated DEGs, and the PPI network in GSE10432 and GSE10433. (a), (b), and (c) represent BP, CC, and KEGG, respectively. (d) is on behalf of the PPI network of GSE10432 and GSE10433; red represents the upregulated genes, and green represents the downregulated genes. Count: number of DEGs.
(a)(b)(c)(d)Table 1
The hub genes in the PPI network of GSE10432 and GSE10433.
IDDegreeIDDegreeIDDegreeIDDegreeCCL28GDF153PTP4A31CRABP21LCN25CX3CL12ITGA11RARRES11S100A73ALOX52IVL1PLAT1PTGES3S100A92ANXA51CTNND11Note: Degree: the number of interactions between two genes in the PPI network. Abbreviations: PPI: protein–protein interaction.


### 3.3. GO and KEGG Pathway Enrichment Analyses of the DEGs in GSE11792

GSE11792 represented that acne treated with isotretinoin for 8 weeks. GO and KEGG pathway enrichment analyses of the DEGs in GSE11792 were also clustered via DAVID. The upregulated DEGs of GSE11792 mainly involved in cell adhesion, extracellular matrix organization in term of BP (Figure4(a)); located in the extracellular region and extracellular exosome in term of CC (Figure 4(b)); participated in integrin binding, serine-type endopeptidase activity, and calcium ion binding in term of MF (Figure 4(c)); enriched in focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway in term of KEGG (Figure 4(d)); all details can be seen in Table S5. The downregulated DEGs of GSE11792 mainly involved in the oxidation-reduction process, cholesterol biosynthetic process, and lipid metabolic process in term of BP (Figure 5(a)); located in integral component of the membrane and extracellular exosome in term of CC (Figure 5(b)); participated in protein homodimerization activity, identical protein binding, and catalytic activity in term of MF (Figure 5(c)); enriched in metabolic pathways, biosynthesis of antibiotics, and peroxisome in term of KEGG (Figure 5(d)); all details can be seen in Table S6.Figure 4
GO and KEGG pathway enrichment analyses of the upregulated DEGs in GSE11792. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)Figure 5
GO and KEGG pathway enrichment analyses of the downregulated DEGs in GSE11792. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)


### 3.4. PPI Network of the DEGs of GSE11792 and Identification of Hub Genes

The PPI network of the DEGs of GSE11792 was generated via STRING and Cytoscape. There were 303 nodes, including 34 upregulated and 269 downregulated DEGs, and 478 edges in the PPI network (Figure6(a)). The term “Degree” denotes the number of interactions between two genes and was used to screen the hub genes with a threshold of Degree≥10. A total of 30 hub genes were identified (Table 2), and the top 10 hub genes included HMGCS1, HMGCR, SLC27A2, FDFT1, MVD, IDI1, IDH1, FDPS, ACAT2, and ACAA1, in which all of them belonged to the downregulated DEGs of GSE11792.Figure 6
The PPI network of DEGs and module analysis in GSE11792. (a) The PPI network of GSE11792 contains 303 nodes and 478 edges; five parts of the PPI network encompassed by five circles represent five modules filtered by app MCODE in Cytoscape. (b–f) represent modules 1, 2, 3, 4, and 7, respectively, which were extracted from the PPI network. Red represents the upregulated genes, and green represents the downregulated genes.
(a)(b)(c)(d)(e)(f)Table 2
The hub genes in the PPI network of GSE11792.
IDDegreeIDDegreeIDDegreeIDDegreeHMGCS121ACAT215MSMO113TM7SF210HMGCR18ACAA115GNPAT13PC10SLC27A217FN114ACOX213G6PD10FDFT117DHCR714ACLY13FBN110MVD15SC5D13ACAA212DHCR2410IDI115PEX1313HADH11CRAT10IDH115NSDHL13DBT11FDPS15MVK13ACADM11Note: Degree: the number of interactions between two genes in the PPI network. Abbreviations: PPI: protein–protein interaction.


### 3.5. Module Analysis of the PPI Network of GSE11792

Five modules were created from the PPI Network of GSE11792 with a cut-off point ofMCODEscore≥4 and nodes≥6, including module 1 (Score=11.273, nodes=12) (Figure 6(b)), module 2 (Score=9, nodes=9) (Figure 6(c)), module 3 (Score=7, nodes=7) (Figure 6(d)), module 4 (Score=5, nodes=15) (Figure 6(e)), and module 7 (Score=4.8, nodes=6) (Figure 6(f)). Most of the hub genes with Degree≥10 belonged to the downregulated DEGs of GSE11792 and mainly enriched in modules 1, 2, and 3, which were mainly discussed for their KEGG pathway involved. As shown in Table 3, module 1 mainly participated in metabolic pathways and biosynthesis of antibiotics. Module 2 was involved in peroxisome and primary bile acid biosynthesis. Module 3 was mainly enriched in valine, leucine and isoleucine degradation, biosynthesis of antibiotics, metabolic pathways, and fatty acid degradation.Table 3
KEGG pathway enrichment analysis of genes of modules 1, 2, and 3 in the PPI network of GSE11792.
TermDescriptionCountP valueGenesModule 1hsa01100Metabolic pathways125.21E-09SC5D, MSMO1, MVD, HMGCR, DHCR7, FDPS, HMGCS1, MVK, ACAT2, IDI1, NSDHL, FDFT1hsa01130Biosynthesis of antibiotics116.70E-15SC5D, MSMO1, MVD, HMGCR, FDPS, HMGCS1, MVK, ACAT2, IDI1, NSDHL, FDFT1hsa00900Terpenoid backbone biosynthesis72.32E-13MVD, HMGCR, FDPS, HMGCS1, MVK, ACAT2, IDI1hsa00100Steroid biosynthesis51.69E-08SC5D, MSMO1, DHCR7, NSDHL, FDFT1hsa00072Synthesis and degradation of ketone bodies20.015886HMGCS1, ACAT2hsa00650Butanoate metabolism20.042368HMGCS1, ACAT2Module 2hsa04146Peroxisome93.18E-16ACOX2, PECR, AMACR, HAO2, GNPAT, IDH1, PEX13, CRAT, SLC27A2hsa00120Primary bile acid biosynthesis20.01961ACOX2, AMACRModule 3hsa00280Valine, leucine and isoleucine degradation67.14E-11ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa01130Biosynthesis of antibiotics61.55E-07ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa01100Metabolic pathways68.88E-04ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa00071Fatty acid degradation44.18E-06ACAA2, ACADM, HADH, ACAA1hsa01212Fatty acid metabolism46.29E-06ACAA2, ACADM, HADH, ACAA1hsa00640Propanoate metabolism32.37E-04ACADM, ECHDC1, PCCBhsa00062Fatty acid elongation20.021616ACAA2, HADHNotes: Count: the number of DEGs. Abbreviations: DEGs: differentially expressed genes; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein–protein interaction.


## 3.1. The Volcano Plots of Genes Expression and the Common DEGs of GSE10432, GSE10433, and GSE11792

Three microarray datasets including GSE10432, GSE10433, and GSE11792 were downloaded from GEO. Basing on a threshold ofP value < 0.05 and ∣logfoldchangeFC∣≥0.263034405833794, namely, ∣FC∣≥1.2, 638 upregulated and 775 downregulated genes in GSE10432, and 176 upregulated and 142 downregulated genes in GSE10433 were obtained. There were 62 upregulated and 249 downregulated genes in GSE11792 with a cut off value of P value < 0.05 and ∣logFC∣≥0.584962501, namely, ∣FC∣≥1.5. The volcano plots of gene expression of GSE10432 (Figure 1(a)), GSE10433 (Figure 1(b)), and GSE11792 (Figure 1(c)) were generated using a free online website imageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/Volcanoplot.html). There were 27 upregulated (Figure 1(d)) and 9 downregulated (Figure 1(e)) genes in GSE10432 and GSE10433 concurrently (Table S1). 62 upregulated genes and 249 downregulated genes were observed in GSE11792 (Table S2). There were consistently 0 upregulated (Figure 1(f)) and 2 downregulated (Figure 1(g)) genes, including ACSBG1 and BCAT2, among GSE10432 (SEB-1 sebocytes incubated with isotretinoin for 72 hours), GSE10433 (acne patients treated with isotretinoin for 1 week), and GSE11792 (acne patients treated with isotretinoin for 8 weeks).Figure 1
Volcano plots of genes expression and the common DEGs in GSE10432, GSE10433, and GSE11792. Volcano plots of genes expression in microarray data GSE10432 (a), GSE10433 (b), and GSE11792 (c). In volcano plots, green represents downregulated genes, red represents upregulated genes, and gray means no significant DEGs. UP: upregulated DEGs; DW: downregulated DEGs; NoDiff: nondifferentially expressed genes. (d) The common upregulated DEGs of GSE10432 and GSE10433. (e) The common downregulated DEGs of GSE10432 and GSE10433. (f) The common upregulated DEGs among GSE10432, GSE10433, and GSE11792. (g) The common downregulated DEGs among GSE10432, GSE10433, and GSE11792. Different colors indicate different datasets. The cross part represents the common DEGs.
(a)(b)(c)(d)(e)(f)(g)


## 3.2. GO, KEGG Pathway Enrichment Analyses and PPI Network of the Common DEGs in GSE10432 and GSE10433

GSE10432 and GSE10433 represented that acne and SEB-1 sebocytes treated with isotretinoin for 1 week and 72 hours, respectively. GO and KEGG pathway enrichment analyses of the consistently upregulated or downregulated DEGs in GSE10432 and GSE10433 were clustered via DAVID. Figures2 and 3 showed the results of GO and KEGG pathway enrichment analyses. As far as the upregulated DEGs of GSE10432 and GSE10433 were concerned, the upregulated DEGs mainly participated in cell adhesion and inflammatory response in term of BP (Figure 2(a)); located in the cytoplasm and extracellular exosome in term of CC (Figure 2(b)); involved in protein binding in term of MF (Figure 2(c)); enriched in TNF signaling pathway in term of KEGG (Figure 2(d)); all details can be seen in Table S3. In the matter of the downregulated DEGs of GSE10432 and GSE10433, the downregulated DEGs mainly involved in long-chain fatty-acyl-CoA biosynthetic process and definitive hemopoiesis in term of BP (Figure 3(a)); located in the external side of the plasma membrane and focal adhesion in term of CC (Figure 3(b)); enriched in fatty acid metabolism in term of KEGG (Figure 3(c)); all details can be seen in Table S4. The PPI network of the common DEGs in GSE10432 and GSE10433 was showed in Figure 3(d), in which the hub genes included CCL2, LCN2, S100A7, PTGES, and GDF15 (Table 1), and CCL2 participated in TNF signaling pathway (Table S3).Figure 2
GO and KEGG pathway enrichment analyses of the upregulated DEGs in GSE10432 and GSE10433. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)Figure 3
GO and KEGG pathway enrichment analyses of the downregulated DEGs, and the PPI network in GSE10432 and GSE10433. (a), (b), and (c) represent BP, CC, and KEGG, respectively. (d) is on behalf of the PPI network of GSE10432 and GSE10433; red represents the upregulated genes, and green represents the downregulated genes. Count: number of DEGs.
(a)(b)(c)(d)Table 1
The hub genes in the PPI network of GSE10432 and GSE10433.
IDDegreeIDDegreeIDDegreeIDDegreeCCL28GDF153PTP4A31CRABP21LCN25CX3CL12ITGA11RARRES11S100A73ALOX52IVL1PLAT1PTGES3S100A92ANXA51CTNND11Note: Degree: the number of interactions between two genes in the PPI network. Abbreviations: PPI: protein–protein interaction.


## 3.3. GO and KEGG Pathway Enrichment Analyses of the DEGs in GSE11792

GSE11792 represented that acne treated with isotretinoin for 8 weeks. GO and KEGG pathway enrichment analyses of the DEGs in GSE11792 were also clustered via DAVID. The upregulated DEGs of GSE11792 mainly involved in cell adhesion, extracellular matrix organization in term of BP (Figure4(a)); located in the extracellular region and extracellular exosome in term of CC (Figure 4(b)); participated in integrin binding, serine-type endopeptidase activity, and calcium ion binding in term of MF (Figure 4(c)); enriched in focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway in term of KEGG (Figure 4(d)); all details can be seen in Table S5. The downregulated DEGs of GSE11792 mainly involved in the oxidation-reduction process, cholesterol biosynthetic process, and lipid metabolic process in term of BP (Figure 5(a)); located in integral component of the membrane and extracellular exosome in term of CC (Figure 5(b)); participated in protein homodimerization activity, identical protein binding, and catalytic activity in term of MF (Figure 5(c)); enriched in metabolic pathways, biosynthesis of antibiotics, and peroxisome in term of KEGG (Figure 5(d)); all details can be seen in Table S6.Figure 4
GO and KEGG pathway enrichment analyses of the upregulated DEGs in GSE11792. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)Figure 5
GO and KEGG pathway enrichment analyses of the downregulated DEGs in GSE11792. (a), (b), (c), and (d) represent BP, CC, MF, and KEGG, respectively. Count: number of DEGs.
(a)(b)(c)(d)


## 3.4. PPI Network of the DEGs of GSE11792 and Identification of Hub Genes

The PPI network of the DEGs of GSE11792 was generated via STRING and Cytoscape. There were 303 nodes, including 34 upregulated and 269 downregulated DEGs, and 478 edges in the PPI network (Figure6(a)). The term “Degree” denotes the number of interactions between two genes and was used to screen the hub genes with a threshold of Degree≥10. A total of 30 hub genes were identified (Table 2), and the top 10 hub genes included HMGCS1, HMGCR, SLC27A2, FDFT1, MVD, IDI1, IDH1, FDPS, ACAT2, and ACAA1, in which all of them belonged to the downregulated DEGs of GSE11792.Figure 6
The PPI network of DEGs and module analysis in GSE11792. (a) The PPI network of GSE11792 contains 303 nodes and 478 edges; five parts of the PPI network encompassed by five circles represent five modules filtered by app MCODE in Cytoscape. (b–f) represent modules 1, 2, 3, 4, and 7, respectively, which were extracted from the PPI network. Red represents the upregulated genes, and green represents the downregulated genes.
(a)(b)(c)(d)(e)(f)Table 2
The hub genes in the PPI network of GSE11792.
IDDegreeIDDegreeIDDegreeIDDegreeHMGCS121ACAT215MSMO113TM7SF210HMGCR18ACAA115GNPAT13PC10SLC27A217FN114ACOX213G6PD10FDFT117DHCR714ACLY13FBN110MVD15SC5D13ACAA212DHCR2410IDI115PEX1313HADH11CRAT10IDH115NSDHL13DBT11FDPS15MVK13ACADM11Note: Degree: the number of interactions between two genes in the PPI network. Abbreviations: PPI: protein–protein interaction.


## 3.5. Module Analysis of the PPI Network of GSE11792

Five modules were created from the PPI Network of GSE11792 with a cut-off point ofMCODEscore≥4 and nodes≥6, including module 1 (Score=11.273, nodes=12) (Figure 6(b)), module 2 (Score=9, nodes=9) (Figure 6(c)), module 3 (Score=7, nodes=7) (Figure 6(d)), module 4 (Score=5, nodes=15) (Figure 6(e)), and module 7 (Score=4.8, nodes=6) (Figure 6(f)). Most of the hub genes with Degree≥10 belonged to the downregulated DEGs of GSE11792 and mainly enriched in modules 1, 2, and 3, which were mainly discussed for their KEGG pathway involved. As shown in Table 3, module 1 mainly participated in metabolic pathways and biosynthesis of antibiotics. Module 2 was involved in peroxisome and primary bile acid biosynthesis. Module 3 was mainly enriched in valine, leucine and isoleucine degradation, biosynthesis of antibiotics, metabolic pathways, and fatty acid degradation.Table 3
KEGG pathway enrichment analysis of genes of modules 1, 2, and 3 in the PPI network of GSE11792.
TermDescriptionCountP valueGenesModule 1hsa01100Metabolic pathways125.21E-09SC5D, MSMO1, MVD, HMGCR, DHCR7, FDPS, HMGCS1, MVK, ACAT2, IDI1, NSDHL, FDFT1hsa01130Biosynthesis of antibiotics116.70E-15SC5D, MSMO1, MVD, HMGCR, FDPS, HMGCS1, MVK, ACAT2, IDI1, NSDHL, FDFT1hsa00900Terpenoid backbone biosynthesis72.32E-13MVD, HMGCR, FDPS, HMGCS1, MVK, ACAT2, IDI1hsa00100Steroid biosynthesis51.69E-08SC5D, MSMO1, DHCR7, NSDHL, FDFT1hsa00072Synthesis and degradation of ketone bodies20.015886HMGCS1, ACAT2hsa00650Butanoate metabolism20.042368HMGCS1, ACAT2Module 2hsa04146Peroxisome93.18E-16ACOX2, PECR, AMACR, HAO2, GNPAT, IDH1, PEX13, CRAT, SLC27A2hsa00120Primary bile acid biosynthesis20.01961ACOX2, AMACRModule 3hsa00280Valine, leucine and isoleucine degradation67.14E-11ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa01130Biosynthesis of antibiotics61.55E-07ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa01100Metabolic pathways68.88E-04ACAA2, DBT, ACADM, HADH, PCCB, ACAA1hsa00071Fatty acid degradation44.18E-06ACAA2, ACADM, HADH, ACAA1hsa01212Fatty acid metabolism46.29E-06ACAA2, ACADM, HADH, ACAA1hsa00640Propanoate metabolism32.37E-04ACADM, ECHDC1, PCCBhsa00062Fatty acid elongation20.021616ACAA2, HADHNotes: Count: the number of DEGs. Abbreviations: DEGs: differentially expressed genes; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein–protein interaction.


## 4. Discussion

There are mainly four key factors involved in the pathogenesis of acne, including follicular colonization by Propionibacterium acnes (P. acnes), increased sebum production, infundibular hyperkeratinization of the pilosebaceous unit, and inflammation [10]. Isotretinoin is the most potent treatment for moderate-severe acne, which targets these four mechanisms of acne pathogenesis [11]. Because of the adverse effects, especially teratogenicity of isotretinoin, multiple studies have been conducted to explore the therapeutical mechanism of isotretinoin and thus develop retinoid-like alternatives without teratogenicity in their side effects’ profile. Dermcidin is an antimicrobial peptide and manifests as anti-inflammatory properties in acne vulgaris; isotretinoin can ameliorate acne vulgaris via enhancing Dermcidin level [12]. Isotretinoin increases transcription factor p53 expression, which induces expression of sestrin 1, sestrin 2, FoxO1, FoxO3, p21, and TRAIL, and suppresses levels of androgen receptor and IGF-1 receptor [13]. Studies of Nelson et al. performed microarray analyses of skin biopsies from acne patients or SEB-1 sebocytes treated with isotretinoin, explored the biological effect of lipocalin 2, which is the most highly upregulated gene [8]; however, there were other significantly upregulated or downregulated genes which might also have key functions in the acne treatment of isotretinoin. Nelson et al. also conducted GO and pathway annotations of GSE10433 and GSE11792, respectively [9], without GSE10432; also PPI networks of GSE10432, GSE10433, and GSE11792 were not constructed.Using integrated bioinformatics analyses, there were 27 upregulated and 9 downregulated common DEGs observed in datasets of GSE10432 and GSE10433. 62 upregulated and 249 downregulated genes were found in GSE11792. The hub genes of the PPI network of the common DEGs in GSE10432 and GSE10433 includedCCL2, LCN2, S100A7, PTGES, and GDF15; in which all of them were upregulated and CCL2 mainly participated in TNF signaling pathway. The hub genes of PPI network of DEGs in GSE11792 included HMGCS1, HMGCR, SLC27A2, FDFT1, MVD, IDI1, IDH1, FDPS, ACAT2, and ACAA1; in which all of them were downregulated. There were mainly three modules in the PPI network of GSE11792. Module 1 mainly participated in metabolic pathways. Module 2 was mainly involved in peroxisome. Module 3 was mainly enriched in metabolic pathways and fatty acid degradation. These three modules support the effects of isotretinoin on sebum decreased production. There were only 2 DEGs including ACSBG1 and BCAT2, which are downregulated among GSE10432, GSE10433, and GSE11792.As mentioned before, the hub genes of the PPI network of the common DEGs in GSE10432 (SEB-1 sebocytes incubated with isotretinoin for 72 hours) and GSE10433 (acne patients treated with isotretinoin for 1 week) were upregulated.CCL2, namely, C-C Motif Chemokine Ligand 2, is a chemotactic factor which attracts monocytes and basophils but not for neutrophils or eosinophils and involves in immunoregulatory and inflammatory processes. A study of Schuster et al. showed that S1pr4 activation increases CCL2 production which promotes macrophage infiltration in a murine psoriasis model [14]. TNF signaling pathway involves in a wide range of biological processes including cell survival, apoptosis, immunity, and inflammation. A previous study showed that the aryl hydrocarbon receptor promotes secretion of inflammatory factors TNF-α and IL-8 in human SZ95 sebocytes, which indicates a possible mechanism in TLR2-mediated acne [15]. Compared to normal skin, the expression of TNF in acne lesions is increased; TNF also induces the formation of lipid droplets in sebocytes [16, 17]. The overexpression of CCL2 originated from activated TNF signaling pathway and S100A7 with chemotactic activity probably associates with some patients’ side effects of “acne-flare” in early stage of isotretinoin treatment [9]. Suppression of the TNF signaling pathway and S100A7 might be an alternative option for inhibiting “acne-flare” of acne patients while using isotretinoin. Previous studies have showed the proapoptotic effect of isotretinoin on sebocytes via mediating associated gene expression. Isotretinoin increases the levels of nuclear FoxO1 and FoxO3 proteins and participates in isotretinoin-induced proapoptotic signaling in sebocytes [18]. TRAIL is increased and involved in the apoptosis of human sebaceous gland cells mediated by isotretinoin [19]. Upregulated LCN expression promotes sebaceous gland cell apoptosis induced by isotretinoin [8]. In the present study, PTGES is a glutathione-dependent prostaglandin E synthase. Interleukin 1 beta or TP53 mediate PTGES expression, which probably participates in TP53 induced apoptosis. The study of Seo et al. showed that the expression of PTGES associates with the prognosis of colorectal cancer [20]. GDF15 belongs to TGF-beta superfamily of proteins. A previous study showed that GDF15 induces cytotoxicity and apoptosis in A549 cells [21]. Our study showed that the proapoptosis effect of isotretinoin on sebocytes might depend on mediating PTGES and GDF15 expression except for LCN2.As has been noted, the hub genes of the PPI network of DEGs in GSE11792 were downregulated, includingHMGCS1, HMGCR, FDFT1, MVD, IDI1, and FDPS. All of these genes serve as critical enzymes in squalene synthesis in the biosynthetic pathway of cholesterol. Downregulated expressions of these genes were in line with the known effect of isotretinoin on decreasing sebum production. In addition, HMGCS1 and HMGCR have been reported to promote cell proliferation and induce apoptosis. HMGCS1 is upregulated and promotes cell proliferation in colon cancer tissues [22]. Overexpression of HMGCR is observed in human lung adenocarcinoma, and knockdown of HMGCR suppresses growth and promotes apoptosis of malignant cells [23]. Study of Ashida et al. showed that the knockdown of endogenous HMGCS1 or HMGCR in prostate cancer (PC) cells significantly decreases PC cell viability [24]. FDFT1 is a membrane-associated enzyme with the ability of stimulating cell growth which associates with disease development and progression in some tumors [25]. This study demonstrated that isotretinoin might also induce apoptosis of sebocytes via mediating HMGCS1, HMGCR, and FDFT1 expression besides LCN2.There were only two genes in common among datasets of GSE10432, GSE10433, and GSE11792, indicating that the effects of isotretinoin on gene expression are time-dependent in most cases. These two genes wereACSBG1 and BCAT2 which were downregulated. ACSBG1 involves in S1P metabolism, which functions as a bioactive lipid molecule [26]. BCAT2 plays an important role in the production of the branched-chain amino acids; the study showed that the depletion of BCAT2 enzyme impairs myoblast survival, indicating a suppressing apoptosis property of BCAT2 [27]. Isotretinoin induced downregulated expression of ACSBG1 and BCAT2 might be related to decreased sebum production and apoptosis of sebocytes, respectively.In conclusion, this study exhibits that the upregulated expression ofLCN2, PTGES, and GDF15 might mediate sebocytes apoptosis and thus decreased sebum production; while CCL2 originated from activated TNF signaling pathway and S100A7 could be related with “acne-flare” using isotretinoin for 1 week. Downregulated genes, including HMGCS1, HMGCR, FDFT1, MVD, IDI1, and FDPS may be associated with decreased sebum synthesis; HMGCS1, HMGCR, and FDFT1 also probably associated with apoptosis of sebocytes treated with isotretinoin for 8 weeks. ACSBG1 and BCAT2 work in both 1 week and 8 weeks and could associate with decreased sebum synthesis and apoptosis of sebocytes, respectively. These results indicate potential therapeutic mechanisms in the treatment of acne using isotretinoin; provide a possible target for avoiding side effect of “acne-flare”. Microarray data from acne treated with isotretinoin for 20 weeks will be helpful to investigate the mechanism of long-term remission of acne induced by isotretinoin.


## 5. Conclusions

Based on microarray datasets GSE10432, GSE10433, and GSE11792, upregulated expression ofLCN2, PTGES, and GDF15 might mediate sebocytes apoptosis; CCL2 and S100A7 could be related with “acne-flare” using isotretinoin for 1 week. Downregulated genes, including HMGCS1, HMGCR, FDFT1, MVD, IDI1, and FDPS may be associated with decreased sebum synthesis; HMGCS1, HMGCR, and FDFT1 also probably associate with the apoptosis of sebocytes treated with isotretinoin for 8 weeks. ACSBG1 and BCAT2 work in both 1 week and 8 weeks and could associate with decreased sebum synthesis and apoptosis of sebocytes, respectively.



---

*Source: 5842795-2020-06-29.xml*