Identifying the Potential Therapeutic Targets for Atopic Dermatitis Through the Immune Infiltration Analysis and Construction of a ceRNA Network.

PURPOSE: This study was meant to analyze immune infiltration and construct a ceRNA network to explore the new therapeutic targets for atopic dermatitis (AD) through bioinformatics way. PATIENTS AND METHODS: We downloaded the AD patients' RNA expression profile datasets (GSE63741, GSE124700) from the Gene Expression Omnibus (GEO) database, which were analyzed through the GEO2R. We explored the hub genes by the enrichment analysis and the protein-protein interaction (PPI) analysis. Moreover, we estimated immune cell types and their proportions by ImmucellAI. GSE121212 dataset validation was performed to verify the robustness of the hub genes. Then, a ceRNA network was constructed by the miRWalk, miRNet, miRDB, DIANA, TargetScan, and starbase database. Finally, gene expression analysis was performed by using RT-qPCR.
RESULTS: In total, we detected 22 differentially expressed genes (DEGs), which contained 8 downregulated genes and 14 upregulated genes. There were 5 hub genes confirmed as key genes through PPI network analysis and the ROC curves. KEGG pathway analysis revealed that they were significantly enriched in the IL-17 signaling pathway and GO analysis showed mainly in the immune cell chemotaxis. The immune infiltration profiles were different between normal controls and AD, and each of the key genes (S100A7, S100A8, S100A9, and LCE3D) was significantly correlated with the main infiltration cell of AD. A lncRNA-miRNA-mRNA ceRNA network containing the key genes was constructed, and NEAT1 and XIST, the core of ceRNA network, were significantly overexpressing verified by RT-qPCR in AD patients.
CONCLUSION: Altogether, the key genes and their ceRNA network provided a novel perspective to the immunomodulation of AD, which may be potential and new therapeutic targets for AD.

Introduction

Atopic dermatitis (AD) is a chronic, recurrent, inflammatory skin disease, which is one of the core challenges for dermatological research in the world.1 With increasing the risk of allergic rhinitis, asthma, and other immune-mediated inflammatory diseases, AD is considered as a systemic disease.2 Over the past 30 years, the prevalence of AD has been increasing worldwide, which has reached 10%~20% in developed countries, and children’s are higher than adults’.3,4 The persistence and recurrence of AD not only has a great impact on patients’ life and mental health but also brings a certain economic burden to both their families and society.5 Although current kinds of research reveal a combination of genetic, environmental, and immunological factors are considered to take part in the pathogenesis of AD, its pathogenesis has been still unclear.6 Therefore, it is important and urgently required to identify the causes, molecular biomarkers and molecular mechanisms for early diagnosis and personalized therapy.

Although the pathogenesis of AD has not been fully understood, more and more studies suggest that immune dysregulation is a potentially key contributor to the pathogenesis of this disease.7,8 It has been found that AD is the non-specific inflammation of the skin caused by the interaction of various immune cells (especially T cells, eosinophils, mast cells, and a small amount of monocytes/macrophages) and their chemokines, cytokines, and inflammation-related factors, such as IL-4, IL-5, IL-6, IFN-γ, CCL22 and so on.1,8 In this setting, there are a few immunotherapy medicine strategies such as Dupilumab, the first systemic biological medicine ratified to AD by Food and Drug Administration (FDA), could block IL-4 related to IL-13 signaling; Crisaborole, a new no steroidal drug, blocks phosphodiesterase 4 (PDE4) and inhibits proinflammatory cytokines, which was approved for AD by FDA in 2016; Nemolizumab, a humanized anti-IL-31 receptor antibody, could reduce the pruritus and the extension of moderate-to-severe AD. These strategies seem to be promising options for AD.9 However, reducing the relapse of AD remains a worldwide challenge.10 Therefore, fuller and further understanding of the immune response mechanism of AD has important theoretical significance for the clinical diagnosis as well as the research of immunotherapy.

In recent years, microarray technology has been used in scientific research, especially in life sciences. Bioinformatics data mining of microarray data such as non-coding RNA and mRNA expression is a valuable tool for exploring the key genes involved in disease pathogenesis and provides useful insights and a basis for further new study. Moreover, Wang et al have suggested that the competing endogenous RNA (ceRNA) network played a part of pivotal regulatory in the development of various skin diseases including AD.11 And the ceRNA network is also involved in various pathophysiological processes such as immunomodulation, apoptosis and so on.12,13 In this research, through using bioinformatic methods, we analyzed the data (GSE63741, GSE124700, etc.) generated by microarray technology from the Gene Expression Omnibus (GEO), to identify the hub genes and construct a ceRNA network, and further investigate their potential role in AD immunomodulation.

Materials and Methods

Microarray Data

The AD patients’ RNA expressed datasets (GSE63741, GSE124700) were downloaded from the GEO database ().14 GSE63741 included an RNA expression profile of 150 samples, containing 30 control samples and 30 AD samples, while GSE124700 included RNA expression datasets for 2 control samples and 2 AD samples (6 subjects). All control samples were assembled from healthy volunteers’ skin biopsies, and AD samples were harvested from lesions of AD skin biopsies (Table 1).

Table 1

Details of the GEO AD Data

Dataset	Platform	Number of Samples (AD/Control, Subjects)
GSE63741	GPL19471 PIQOR (TM) Skin 2.0 Microarray, human, antisense (591)	150(30/30,60)
GSE124700	GPL10558 Illumina HumanHT-12 V4.0 expression beadchip	4(2/2,6)

Abbreviation: AD, atopic dermatitis.

Recognization of Differentially Expressed Genes (DEGs)

The online tool, GEO2R (),15 was employed in preprocessing data and screening for mRNAs that were differentially expressed between control tissue samples and AD. p<0.05 and |log FC| >1 were selected as the difference.

Function and Pathway Enrichment Analyses

GSEA was performed to analyze the biological function of the DEGs. The idiographic flow acted under the recommended protocol of the WEB-based GEne SeT AnaLysis Toolkit database (),16 and set FDR q-value < 0.05 as the cut-off criteria. GO (included molecular function (MF), cellular component (CC), and biological process (BP)), and KEGG pathway analysis were performed through Metascape tool ().17 P< 0.01 was recognized as the cut-off criterion.

Protein–Protein Interaction (PPI) Network and Module Analysis

The STRING database ()18 was carried out to explore a PPI network, and the Cytoscape software19 to build the PPI network and analyze the reciprocities of the DEGs. The hub genes were identified by the cytoHubba plug-in. The modules of the PPI network were developed by the Molecular Complex Detection (MCODE) plug-in. The ClueGO plugin was used for enrichment analysis of the hub genes.

Data Verification

The GSE121212 downloaded from the GEO database was utilized for verifying the hub genes’ reliability, which contained 27 non-lesion samples and 27 AD lesion samples. In addition, GraphPad Prism 8.0.2 (Graph-Pad Software Inc., San Diego, CA, USA) was used to construct the ROC curves and estimate the areas under the curve (AUCs).

Immune Infiltration and Association Analysis

The immune infiltration was calculated based on the gene expression data using the web tool ImmucellAI (),20 of which 24 immune cells are composed of 6 immune cell types and 18 T-cell subtypes. The difference of immune cell infiltration between the tissue of AD and control was analyzed through the principal component analysis (PCA). The proportion of immune cells of the two groups and the correlation between the hub genes and immune cells in AD were analyzed in GraphPad Prism 8.0.2.

The ceRNA Network Analysis

To explore the lncRNA-miRNA-miRNA ceRNA network, DIANA Tools (),21 the miRNet (),22 TargetScan (),23 miRWalk ()24 and miRDB ()25 online databases were performed to identify possible miRNA targeting hub genes, and those identified in ≥3 object databases were considered as true. The starbase (),26 miRNet, lncbase v.3(),27 lncbase v.2 ()28 online databases were utilized for predicting potential lncRNA targeting the miRNA, and those identified in ≥4 object databases were considered as true. The web-based tools, the Cytoscape software, Wei Sheng Xin () and Draw Venn Diagram () were used for data visualization.

The RT-qPCR Analysis

To confirm the core lncRNAs of ceRNA network, serums from 3 patients with OA and 3 healthy control were harvested for RT-qPCR validation. The experiment was executed following the authorization of the ethics committee of The First Affiliated Hospital of Hunan Normal University (Hunan Province People’s Hospital) (2021 Scientific Research Ethics Review NO: 05). All patients signed the informed consent. The TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) was used to extract the total RNA from the serum. The cDNA was created by reverse transcription with RNA samples from total RNA, and RT-qPCR was performed by using the Revert Aid First Strand cDNA Synthesis Kit (Kangweishiji, Beijing, China). β-actin was set as the internal reference. All primers were engineered and compounded by Sheng Gong Company (Shanghai, China) (Table 2). Relative lncRNA expression was computed using the 2−ΔΔCt method.29 The experiment was done independently 3 times.

Table 2

The Primer Used in RT-qPCR

Gene	Primer Sequence (5ʹ – 3ʹ)
H-actin	F:ACCCTGAAGTACCCCATCGAG
	R:AGCACAGCCTGGATAGCAAC
H-NEAT1	F:CCTGCCTTCTTGTGCGTTT
	R:TAGCACAACACAATGACACCCT
H-XIST	F: ACTGCCACCCATATATAAGCTA
	R:AGTAATCACCATTCAGTAAGCCA

Statistical Analysis

Two groups’ data were analyzed through the unpaired Student’s t-test. P<0.05 indicated to be significantly different. GraphPad Prism 8.0.2 (Graph-Pad Software Inc., San Diego, CA, USA) was recognized as the statistical software.

Results

Identification of DEGs in AD

As shown in Figure 1, the AD RNA expression profile datasets (GSE63741, GSE124700) were normalized. The GSE63741 dataset included 66 different genes, and the GSE124700 dataset contained 434 differential genes, whose volcano plot are shown in Figure 2A and B. Utilizing integrated bioinformatics analysis, altogether consistent 22 DEGs (containing 8 consistently downregulated genes and 14 consistently upregulated genes) were identified (Figure 2C and D, and Table 3). Figure 3 shows the heat map of DEGs.

Figure 1

Data after normalization. (A) Standardization of selected samples of GSE63741. (B) Standardization of GSE124700. Green bars represent AD, blue bars represent control.

Figure 2

Identification of DEGs in AD. (A) The differentially expressed genes of GSE63741. (B) The differentially expressed genes of GSE124700. (C) 8 DEGs were consistently downregulated in the two datasets. (D) 14 DEGs were consistently upregulated in the two datasets. The green point represents downregulated genes, the red point represents upregulated. |log FC|>1 and p-value <0.05 were set as the difference.

Table 3

The DEGs of AD

Regulation	Common DEGs of GSE63741 and GSE124700
Upregulated(n=14)	CCL18, CCL22, AKR1B10, CCL19, S100A7, LCE3D, IFI27
	CHI3L2, AQP3, LY6D, S100A8, RAB31, MMP12, S100A9
Downregulated(n=8)	FLG2, CST6, SERPINA12, HBB, KRT15, CCND1, LCE1B, HBA1

Abbreviations: AD, atopic dermatitis; DEGs, differentially expressed genes.

Figure 3

The heat map of the DEGs. The color from green to red indicates the expression of genes from low to high.

Function and Pathway Enrichment Analyses of the DEGs

GSEA was used for enrichment analysis of the DEGs to explore their biological roles. The results by GSEA show that in the biological processes categories, most of DEGs were enriched in response to stimulus, biological regulation and developmental process etc; In the cellular component categories, the DEGs were mostly involved in membrane and extracellular space; In the molecular function categories, the major DEGs were involved in enzyme regulator activity and protein binding (Figure 4A). The results shown in Figure 4B indicated that the significantly enriched pathway was the IL-17 signaling pathway. Utilizing Metascape online database to further validate our results, the top ten GO and KEGG items of DEGs enriched showed that they were obviously involved in epidermis development, granulocyte chemotaxis, lymphocyte chemotaxis, and monocyte chemotaxis in BPs; collagen-containing extracellular matrix and extracellular matrix in CCs; cytokine activity, cytokine receptor binding, and chemokine receptor binding in MFs; IL-17 signaling pathway, Cytokine-cytokine receptor interaction, and Chemokine signaling pathway in KEGG pathway (Table 4 and Figure 5).

Figure 4

The GSEA analysis of the DEGs. (A) Biological processes, cellular components and molecular functions. (B) The signaling pathways.

Table 4

The Top 10 Function and Pathway Enrichment Items of the DEGs

Category	Description	P value	Genes
GO BP	GO:0008544 epidermis development	5.45551E-08	AQP3,CST6,KRT15,S100A7,LCE3D,LCE1B,FLG2
GO BP	GO:0071621 granulocyte chemotaxis	5.44051E-10	S100A7,S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0097530 granulocyte migration	1.54858E-09	S100A7,S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0097529 myeloid leukocyte migration	1.64556E-08	S100A7,S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0030595 leukocyte chemotaxis	2.08055E-08	S100A7,S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0060326 cell chemotaxis	1.17284E-07	S100A7,S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0030593 neutrophil chemotaxis	1.53171E-08	S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:1990266 neutrophil migration	3.74245E-08	S100A8,S100A9,CCL18,CCL19,CCL22
GO BP	GO:0048247 lymphocyte chemotaxis	1.87602E-07	S100A7,CCL18,CCL19,CCL22
GO BP	GO:0002548 monocyte chemotaxis	2.25311E-07	S100A7,CCL18,CCL19,CCL22
GO CC	GO:0060205 cytoplasmic vesicle lumen	5.6214E-06	HBA1,HBB,S100A7,S100A8,S100A9
GO CC	GO:0031983 vesicle lumen	5.7024E-06	HBA1,HBB,S100A7,S100A8,S100A9
GO CC	GO:0031012 extracellular matrix	0.000923351	MMP12,S100A7,S100A8,S100A9
GO CC	GO:0034774 secretory granule lumen	0.001974261	S100A7,S100A8,S100A9
GO CC	GO:0062023 collagen-containing extracellular matrix	0.004403432	S100A7,S100A8,S100A9
GO CC	GO:0030139 endocytic vesicle	0.001821123	HBA1,HBB,RAB31
GO MF	GO:0005509 calcium ion binding	0.000192896	MMP12,S100A7,S100A8,S100A9,FLG2
GO MF	GO:0043177 organic acid binding	3.1967E-05	HBA1,HBB,S100A8,S100A9
GO MF	GO:0050786 RAGE receptor binding	6.92211E-08	S100A7,S100A8,S100A9
GO MF	GO:0048020 CCR chemokine receptor binding	6.67899E-06	CCL18,CCL19,CCL22
GO MF	GO:0008009 chemokine activity	7.58116E-06	CCL18,CCL19,CCL22
GO MF	GO:0042379 chemokine receptor binding	2.22851E-05	CCL18,CCL19,CCL22
GO MF	GO:0048306 calcium-dependent protein binding	3.99042E-05	S100A7,S100A8,S100A9
GO MF	GO:0016209 antioxidant activity	4.13255E-05	HBA1,HBB,S100A9
GO MF	GO:0005125 cytokine activity	0.000799355	CCL18,CCL19,CCL22
GO MF	GO:0005126 cytokine receptor binding	0.001205692	CCL18,CCL19,CCL22
KEGG pathway	hsa04657 IL-17 signaling pathway	5.56373E-05	S100A7,S100A8,S100A9
KEGG pathway	hsa04062 Chemokine signaling pathway	0.000515188	CCL18,CCL19,CCL22
KEGG pathway	hsa04060 Cytokine-cytokine receptor interaction	0.002080696	CCL18,CCL19,CCL22

Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological processes; CC, cellular components; MF, molecular functions.

Figure 5

The top 10 items of GO and KEGG enrichment analysis of DEGs. (A) GO Biological processes. (B) GO cellular components. (C) GO molecular functions. (D) KEGG signaling pathways. **p<0.01; ***p<0.001.

PPI Network and Module Analyses

Based on the STRING database, the PPI network of the DEGs was built through Cytoscape software, including 15 edges and 14 nodes. Furthermore, among the 14 nodes, 9 nodes were upregulated, while 5 were downregulated (Figure 6A). As shown in Figure 6B, the most important module containing 6 nodes and 9 edges was selected. Using two algorithms (degree and MCC), the cytoHubba was used to detect hub genes. The 5 highest-scoring genes were recognized as hub genes of the PPI network: HBB, LCE3D, S100A7, S100A8, and S100A9 (Figure 6C and Table 5).

Figure 6

The PPI network analysis and identification of hub genes. (A) PPI network of the DEGs. The sizes of the edges and nodes represent the degree, which is the bigger size representing the higher degree. The green node represents downregulated gene and the red represents upregulated. (B) The first module of the PPI network. The green node represents downregulated gene and the red represents upregulated. (C) 5 hub genes were identified through the cytoHubba plug-in. (D) The biological process of hub genes through ClueGO plug-in. (E) Signaling pathways of hub genes through ClueGO plug-in.

Table 5

The Top 5 Hub Genes

Genes	Description	Degree	MCC	Expression Change
HBB	Hemoglobin subunit beta	4	4	Downregulated
LCE3D	Late cornified envelope 3D	3	3	Upregulated
S100A7	S100 calcium binding protein A7	4	4	Upregulated
S100A9	S100 calcium binding protein A9	3	4	Upregulated
S100A8	S100 calcium binding protein A8	3	4	Upregulated

Abbreviation: MCC, maximal clique centrality.

Enrichment Analysis of the Hub Genes

Due to the fact that the most closely-knit genes in a network are keys to regulatory, the ClueGO plugin was used to further investigate the functions of the hub genes. As Figure 6D and E exhibited, they are still mostly associated with immunomodulation including antimicrobial humoral response, regulation of chemotaxis, leukocyte chemotaxis, granulocyte chemotaxis, and IL-17 signaling pathway (Figure 6D and E).

Data Validation of the Hub Genes

The validation data GSE121212 was used to ensure the reliability of the hub genes. The gene-level of the hub genes in AD and control tissue exhibited in Figure 7A indicated that, compared with the control, HBB is significantly lower expression while S100A7, S100A8, S100A9 and LCE3D are significantly overexpression in AD, which is the same as the above results. GraphPad Prism 8.0.2 was used to perform ROC curve analysis, and the results exhibited the hub genes correlated with AD were identified as possible biomarkers for AD diagnosis, which suggested that they are the key genes of AD (Figure 7B and Table 6).

Figure 7

Data verification of hub genes through GSE121212. (A) Gene expression of each hub gene. (B) The ROC curve of each hub gene. **p<0.01; ****p<0.0001.

Table 6

Results of AUCs for Hub Genes

Hub Genes	AUC	95% Confidence Interval	P value
LCE3D	0.9342	0.8710 to 0.9974	<0.0001
S100A8	0.9328	0.8717 to 0.9938	<0.0001
S100A9	0.9246	0.8586 to 0.9905	<0.0001
S100A7	0.9122	0.8383 to 0.9861	<0.0001
HBB	0.7956	0.6720 to 0.9193	0.0002

Abbreviation: AUC, area under the curve.

Immune Infiltration and Association Analyses

To reveal the landscape of immune infiltration in AD, the ImmucellAI algorithm was used to investigate the difference between normal tissues and AD in immune infiltration. The proportions of immune cells from the tissues of AD and normal control displayed group-bias clustering by PCA (Figure 8A). In comparison to control tissue, AD tissue owned a higher proportion of CD4 T cell, Th2 cell, DC cell, Tfh cell, Tr1 cell, iTreg cell, nTreg cell, CD8 T cell, and cytotoxic cell, but the proportions of NKT cell, Th17 cell, CD8 naive cell, and neutrophil cell were relatively lower (Figure 8B). Subsequently, we studied the correlation between hub genes expression and the abundance of the immune cell based on the Pearson’s correlation coefficient. The results shown in Figure 9 indicated that each of S100A7, S100A8, S100A9, and LCE3D was significantly positively correlated with the expression of Th2 cell, DC cell, Tfh cell, Tr1 cell, nTreg cell, iTreg cell, and CD8+T cell, but negatively correlated with the expression of CD8 naive cell, NKT cell, and Th17 cell. HBB was just significantly correlated with neutrophil cell, Tfh cell, and Th17 cell.

Figure 8

The difference of immune infiltration between AD and controls. (A) PCA performed on all samples. (B) The proportion of immune cells in AD and controls. *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001; blue represents control and red represents AD.

Figure 9

The correlation between hub genes and immune cells. (A) HBB. (B) LCE3D. (C) S100A7. (D) S100A8. (E) S100A9. P< 0.05 was considered as statistical significance.

The lncRNA-miRNA-mRNA ceRNA Network of AD

LncRNA can act as a competitive endogenous RNA to play an opposite role with miRNA in regulating mRNA because of containing the miRNA binding sites.30 Five online databases were carried out to predict potential miRNA targeting hub genes, and the Venn diagram to get the intersection between each of the databases, and the results were shown in Figure 10. Being identified in ≥3 object databases was set as the cut-off criteria, and the miRNAs are shown in Table 7. Subsequently, we used four online databases to predict potential lncRNA targeted by the miRNA and the Venn diagram to get the intersection between each of the databases, and the result as shown in Figure 11. Being identified in ≥4 object databases were set as the cutoff criteria and there are a total of 13 lncRNAs identified. Finally, the ceRNA network containing S100A9/hsa-miR-588/NEAT1, S100A9/hsa-miR-588/NEAT1, LCE3D/hsa-miR-1224-5p/NEAT1, S100A8/hsa-miR-98-5p/XIST, S100A9/hsa-miR-766-5p/XIST, S100A9/hsa-miR-588/XIST and hence a total of 18 axes, was constructed (Figure 12 and Table 8).

Figure 10

Venn diagram showed the miRNAs targeting each hub gene which were predicted in the online database. (A) HBB. (B) LCE3D. (C) S100A7. (D) S100A9. (E) S100A8. Predicted in ≥3 object databases were considered as true.

Table 7

The miRNAs Targeting Hub Genes

Genes	The miRNAs
HBB	hsa-miR-6851-3p hsa-miR-3922-5p hsa-miR-361-3p hsa-miR-1910-3p
	hsa-miR-6511a-5p hsa-miR-3914 hsa-miR-33a-3p hsa-miR-5009-5p
	hsa-miR-6822-5p hsa-miR-4307 hsa-miR-4423-3p hsa-miR-6750-5p
	hsa-miR-6886-3p hsa-miR-8058
LCE3D	hsa-miR-193b-5p hsa-miR-1224-5p hsa-miR-3170 hsa-miR-4700-5p
	hsa-miR-4320 hsa-miR-6763-5p hsa-miR-4667-5p hsa-miR-3690
	hsa-miR-4471 hsa-miR-513a-5p hsa-miR-4269 hsa-miR-4314
	hsa-miR-6750-5p hsa-miR-6855-5p hsa-miR-4691-5p hsa-miR-6792-3p
	hsa-miR-5004-5p
S100A7	hsa-miR-30c-2-3p hsa-miR-1343-3p hsa-miR-6788-5p hsa-miR-6852-5p
S100A8	hsa-miR-98-5p
S100A9	hsa-miR-4701-5p hsa-miR-132-5p hsa-miR-588 hsa-miR-1204
	hsa-miR-3677-5p hsa-miR-940 hsa-miR-6878-3p
	hsa-miR-766-5p hsa-miR-92a-1-5p hsa-miR-660-3p

Figure 11

Venn diagram showed the lncRNAs which were predicted to target the identified miRNA in online database. (A) hsa-miR-361-3p. (B) hsa-miR-513a-5p. (C) hsa-miR-1224-5p. (D) hsa-miR-98-5p. (E) hsa-miR-132-5p. (F) hsa-miR-766-5p. (G) hsa-miR-588. Predicted in ≥4 object databases were considered as true.

Figure 12

The lncRNA–miRNA–mRNA ceRNA network of AD. (A) The ceRNA network was constructed through Cytoscape software. The triangle represents the hub gene, and the circle represents the miRNA and the rhombus represents the lncRNA. The lncRNA whose degree was highest was painted red. (B) The alluvial diagram of the ceRNA network.

Table 8

The ceRNA Network of AD

mRNA	miRNA	lncRNA
HBB	hsa-miR-361-3p	SNHG1, GUSBP11, MALAT1
LCE3D	hsa-miR-1224-5p	NEAT1
S100A8	hsa-miR-98-5p	NUTM2A-AS1, TRG-AS1, XIST, TTTY15
S100A9	hsa-miR-766-5p	RMRP, STAG3L5P-PVRIG2P-PILRB, XIST
		NEAT1, LINC00894, KCNQ1OT1, NUTM2A-AS1
	hsa-miR-588	XIST,NEAT1
	hsa-miR-132-5p	TTN-AS1

RT-qPCR Validation of the Core lncRNAs of ceRNA Network

NEAT1 and XIST are identified as the core lncRNAs, which have the most degree of connectivity because they both are the intersection of three axes of the ceRNA network (Figure 12). Due to the fact that the downstream target genes of ceRNA (S100A8, S100A9, and LCE3D) are overexpression, they are predicted to be equally overexpressed based on the mechanism of ceRNA. NEAT1 and XIST were validated in AD patients’ serum using RT-qPCR analysis. The results show that NEAT1 and XIST are obviously overexpressed in AD compared to healthy controls (Figure 13), which is consistent with the bioinformatics prediction.

Figure 13

RT-qPCR validation of the core lncRNAs of the ceRNA network. (A) The relative expression of XIST between AD and healthy control. (B) The relative expression of NEAT1 between AD and healthy control. *p<0.05.

Discussion

AD is a common chronic inflammatory skin disease, which has presented a leading non-fatal health burden attributable worldwide.31 Although immunotherapy of AD got some effectiveness, how to alleviate and prevent recurrence more effectively is a worldwide problem for AD at present owing to the fact that there still are many uncertainties in the mechanism.32 With the fast development of science, bioinformatics has provided a powerful strategy for exploring potential molecular markers and pathophysiological mechanisms.33 Therefore, this study aims to develop and validate the potential key genes and molecular mechanisms of AD through bioinformatics analysis, especially in immunomodulation, to propose a new perspective into the pathogenesis and treatment strategies of AD.

In this study, we analyzed three RNA expression profile datasets (GSE63741 and GSE124700 as the training datasets, and GSE121212 as the validation datasets), and only selected the DEGs that overlapped in the two training datasets, which can prevent from resulting in a high probability of false-positive results that may ensue from most of previous research focused on the single dataset. 22 DEGs, which included 14 upregulated genes and 8 downregulated genes, were identified and subjected to gene enrichment analyses. The results of the databases, including GSEA and Metascape, confirmed each other, which showed that these genes were mainly involved in immune response (including leukocyte chemotaxis, granulocyte chemotaxis, lymphocyte chemotaxis and so on), developmental process (including epidermis development), biological regulation (including cytokine activity, chemokine receptor binding and so on), most of which are commonly accepted components of the pathogenesis and pathological changes of AD.34 The researchers found that lesional atopic dermatitis skin shows the dysregulated expression of a wide range of cytokines, mostly related to keratinocyte activity and T-cell infiltration, especially for Th22-associated (IL-22) and Th2-associated (IL-4, IL-5).35,36 The IL−17 signaling pathway was the main result of pathway enrichment. IL-17 was not only produced by Th17 cells, but also excreted by various innate cells such as γδ-T cells, lymphoid tissue inducer, natural killer T, natural killer, dendritic cells (DC), and macrophages.37 The pro-inflammatory effect of IL-17 was initially confirmed in fibroblasts, in which it could stimulate NF-kappa-B and related cytokines.38,39 Epithelial cells could be stimulated by IL-17 to express IL-6 and IL-8 and granulocyte colony-stimulating factor, inducing neutrophils migration and activation, which might mediate AD-related immune dysregulation.40 The DEGs were also involved in cytokine−cytokine receptor interaction and Chemokine signaling pathway. To some extent, these suggest that the DEGs have the function of participating in the immunomodulation of AD.

We identified the five hub genes by the PPI network and module analysis, which were HBB, S100A7, S100A8, S100A9, and LCE3D. ROC curve analysis verified the reliability of their diagnostic value, which suggests that they are the key genes and could be potential prognostic and diagnostic biomarkers. S100A7, S100A8, and S100A9 are known as antimicrobial peptides (AMPs) in the epidermis.41 They can kill microorganisms, such as bacteria, viruses and so on, and also exhibit immunomodulatory activities, inducing altering cytokine/chemokine expression, cell migration, proliferation, and differentiation, as one of the first defensive reaction of the natural immune system.42,43 Gittler et al identified that increased expression of S100A9, S100A8, and S100A7 genes were associated with AD initiation and activation of Th2 and Th22 cytokines.44 Son et al suggested that S100A7-related signaling molecules might be effective targets for restoring skin barrier function in AD where S100A7 is accumulated excessively.45 In agreement with these studies, we also find that S100A7, S100A8, and S100A9 were significantly overexpression in AD, and the results of functional enrichment analysis showed that they mediated the IL−17 signaling pathway, and developed antimicrobial humoral response, leukocyte chemotaxis, granulocyte chemotaxis and so on, which made it more systematically explained that S100A7, S100A8, and S100A9 played a role in AD. HBB and LCE3D have not been found to be associated with AD before. It is well known that defect of the HBB gene leads to the occurrence of β-thalassemia, and LCE3D is associated with susceptibility to psoriasis.46,47 They may also be important for AD because they could, respectively, interact with S100A7, S100A8, and S100A9 as described above, which may be new perspective to the pathogenesis of AD.

To further clarify immune dysregulation of AD, we analysed the immune infiltration. It was determined that there was some diversity between the tissue of normal controls and that of AD patients in immune cell infiltration through the PCA. The basic feature of AD was Th2-type inflammation, and IL 4 and IL 13 were important cytokines mediating the pathogenesis of AD, which were mainly produced by Th2 cells, basophil, and type 2 innate lymphoid cells.48 Brunner et al indicated that skin barrier dysfunction allowed external environmental substances (instance of microorganisms and allergens) to easily enter the epidermis and initiated Th2 type inflammation, which Langerhans cells and dermal dendritic cells participated in by presenting allergens.49 These were further confirmed in our study. We found that CD4+ T cell, Th2 cell, and DC cell were higher proportion in AD. Besides, there were also Tr1 cell, iTreg cell, and nTreg cell polarization, Tfh cell, and other cells infiltration. Our further research showed that each of S100A7, S100A8, S100A9, and LCE3D was significantly and positively correlated with most of the infiltration cells in AD like Th2 cell, DC cell, and iTreg cell, which indicated that S100A7, S100A8, S100A9, and LCE3D correlated with immune dysregulation of AD and might be its the potential critical immunomodulation sites.

Moreover, to screen out key molecules and axes and unveil a systematic regulatory network in AD, a ceRNA network was constructed, in which the S100A9/hsa-miR-588/NEAT1, S100A9/hsa-miR-588/NEAT1, LCE3D/hsa-miR-1224-5p/NEAT1, S100A8/hsa-miR-98-5p/XIST, S100A9/hsa-miR-766-5p/XIST, S100A9/hsa-miR-588/XIST and hence a total of 18 axes might serve critical roles in AD. The fact that NEAT1 and XIST, the core lncRNAs of ceRNA network, were significantly overexpressing was verified in Clinical AD patient serum by RT-qPCR, which was consistent with the bioinformatics prediction based on the ceRNA mechanism. Both NEAT1 and XIST played a pivotal regulatory part in the immunomodulation of various diseases, however, they have not been found to be associated with AD before, so our findings may bring a new perspective for the immunomodulation pathogenesis of AD.

In addition, there are also aware of some limitations in this study. We determined the level of NEAT1 and XIST using Clinical AD patient’ serum but not tissue, which is not as good as possible but still have Clinical representative. And it is insufficient to conduct a secondary analysis of existing data only through bioinformatics.

Conclusion

Based on the bioinformatic analysis, we identified 5 key genes including HBB, S100A7, S100A8, S100A9, and LCE3D, which are mainly involved in immune response in AD. Especially the four of them, S100A7, S100A8, S100A9, and LCE3D, were significantly correlated with the main infiltration cell of AD in particular, indicating they were the potential new therapeutic targets for AD to ameliorate immune dysregulation. Moreover, we set up a ceRNA network containing the key genes and verified the core lncRNAs of the network by RT-qPCR to make it more reliable, which provided a more detailed molecular mechanism to understand AD development. In the future, we are going to do further experimental studies of these genes to further validate these predicted results of bioinformatics analysis.