Identification of potential hub genes associated with atopic dermatitis-like recombinant human epidermal model using integrated transcriptomic and proteomic analysis

Atopic dermatitis (AD) is a severe inflammatory skin disorder, characterized by elevated levels of proinflammatory cytokines that fuel a vicious cycle of inflammation. While inflammatory recombinant human epidermal (RHE) models relevant to AD have been established, comprehensive understanding remains limited. To illuminate changes and identify potential hub genes involved in AD-related inflammation, RHE models, stimulated by an inflammatory cocktail including polyinosinic-polycytidylic acid, tumor necrosis factor-alpha (TNF-α ), interleukin 4 (IL-4), and interleukin 13 (IL-13), were constructed and examined using tandem mass tags-proteomic coupled with RNA-seq transcriptomic analyses. Principal component analysis (PCA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional enrichment were employed for the analysis of related genes and proteins. Protein – protein interaction networks helped identify hub genes, which were further confirmed by qPCR and western blot. We observed high expression of thymic stromal lymphopoietin in the inflammatory RHE. Our study identified 2369 differentially expressed genes and 880 differentially expressed proteins in the cocktail-induced group versus the normal control group. A total of 248 overlapping symbols were enriched in various biological processes and signaling pathways, including cornification envelope, cell – cell junction, calcium ion binding, extracellular matrix receptor, terpenoid backbone biosynthesis, and peroxisome proliferator-activated receptors signaling pathway, among others. Among the 248 overlapping symbols, CytoHubba identified 10 hub molecules, namely, signal transducer and activator of transcription 3 (STAT3), integrin subunit beta 1 (ITGB1), filaggrin (FLG), involucrin (IVL), DEAD (Asp-Glu-Ala-Asp) box polypeptide 58 (DDX58), small proline-rich protein 1B (SPRR1B), interferon induced with helicase C domain 1 (IFIH1), desmoglein 1 (DSG1), collagen type XVII alpha 1 chain (COL17A1), and integrin subunit alpha 6 (ITGA6), based on the degree. These integrated results offer valuable insights into the molecular mechanisms of AD and present potential tools for screening cosmetic formulations intended for the treatment of AD.


Introduction
Atopic dermatitis (AD), also known as eczema or atopic eczema, is a common inflammatory disease that affects up to 20% of children and 3% of adults [1,2]. Unfortunately, the worldwide incidence of AD has increased remarkably and alarmingly, although it is stable in high-income countries. The disease has a heterogeneous clinical presentation among individuals and is characterized by persistent pruritus and recurrent eczematous lesions [3]. Common clinical signs observed in patients with AD include generalized skin dryness, early onset of the disease, and a family history of atopic diseases, such as allergic rhinitis and asthma [4]. The pathophysiology of AD is characterized by a complicated interplay of multiple factors, including a dysfunctional epidermal barrier, abnormalities in the skin microbiome, immunologic dysregulation, inflammation, neuroimmune interactions, and genetic risk factors [5]. Mutations in the filaggrin gene (FLG) leading to alterations in the skin barrier protein FLG and profilaggrin represent the best known genetic risk factor associated with abnormalities in lipid chain length [6], water loss, and finally inflammation [7,8]. AD is classified as a type 2 immune-mediated disease [9]. Cutaneous inflammation plays a key role in the pathogenesis of AD, and overexpression of inflammatory factors is also one of the main features of AD. Lesional AD skin exhibits aberrant expression patterns of proinflammatory factors associated with keratinocytes and T cells, particularly high expression of interleukin 4 (IL-4), interleukin 13 (IL-13), interleukin 4 receptor (IL-4R), interleukin 33 (IL- 33), and thymic stromal lym-phopoietin (TSLP) [10]. TSLP, a member of the IL-2 cytokine family, is expressed in epidermal cells when the skin barrier is compromised by scratching, antigens (pollen and mites), and cytokines, such as tumor necrosis factor-alpha (TNF-α) [11]. TSLP triggers the activation of dendritic cells (from the human thymic stroma) that upregulate costimulatory molecules, such as CD86, which in turn promotes the proliferation and differentiation of CD4+ T cells and enhances the immune response and the production of IL-4 and IL-13 [12]. In type 2 inflammation, IL-4 and IL-3 play major roles in opening cellular tight junctions [13], reducing structural epidermal proteins such as FLG [14], and defects in endogenous protease inhibitors [15].
Many studies have used tissue engineering and mice to develop in vitro models with the characteristics of AD to investigate disease mechanisms and find potential therapeutic targets [16]. AD mouse models are classified into three groups: 1) inbred models, 2) genetically engineered models, 3) models induced by exogenous agents [17]. Due to The Organization for Economic Cooperation and Development's promotion of alternative methods to animal testing in cosmetic research, two-dimensional (2D) and three-dimensional (3D) models have become important tools. The simplest model is a 2D cell system such as 2,4-dinitrochlorobenzene (DNCB)-induced AD in HaCaT cells [18]. Recombinant 3D skin constructs can mimic the real skin dermal structures and barrier function. In the context of AD, some recombinant human epidermal (RHE) models have been specifically designed to mimic FLG deficiency [19,20], while others use TNF-α or polyinosinic-polycytidylic acid (poly I:C) alone or in combination with Th2 cytokines to induce AD-like features in RHE models [21,22]. Previous studies on AD-related inflammation models have been limited at the transcriptome and morphology levels, and there is no integrated transcriptomic and proteomic study to gain a comprehensive understanding of the AD-like inflammatory model. In view of the above considerations, we prepared the construction of AD-like inflammatory RHE models by adding an inflammatory cocktail (poly I:C, TNF-α, IL-4, and IL-13). Subsequently, RNA-seq and tandem mass tags (TMTs) proteomic technology were performed to examine gene expression and protein profiles of the RHE models under inflammatory and normal conditions (NC). Integrated transcriptomic and proteomic analysis was performed to investigate hub genes. A better understanding of AD-like RHEs will provide deep thinking for further search for mechanization of AD and can be used as a relevant tool for screening cosmetic formulations to alleviate AD.

Cell cultures
Normal human epidermal keratinocytes (NHEKs) were obtained from Guangdong BioCell Biotechnology Co., Ltd. NHEK cells were cultured in a monolayer with epilife growth medium supplemented with human keratinocyte growth supplement (HKGS) (Thermofisher Scientific, Waltham, MA, USA) until they were 80% confluent. Second passage NHEKs were used to generate RHE models.
Establishing 3D reconstructed human epidermis models The inserts precoated with collagen I (Corning, NY, USA) were placed in the lowest position in the respective wells. NHEK cells were seeded into the inserts at a density of 3.5 × 105 cells/0.5 cm 2 . The NHEKs were seeded with 0.5 mL epilife growth medium supplemented with HKGS, 10 ng/mL KGF, 140 μM calcium chloride (CaCl 2 ), and 50 μg/mL ascorbic acid (Sigma-Aldrich, Burlington, MA, USA) in the lower and upper compartments for two days. The culture was raised to the air-liquid interface and then the appropriate volume of medium was added with an additional 1.5 mM CaCl 2 .
After 16 days in culture, an inflammatory cocktail containing poly I:C 10 μg/mL, TNF-α 10 ng/mL, IL-4 50 ng/mL, and IL-13 50 ng/mL (R&D, Minneapolis, MN, USA) was added to the medium and incubated for 72 h. The 3D RHE models were divided into cocktail-induced positive control RHEs (PC group) and normal control RHEs (NC group) depending on the treatment. These RHE models were used for the following experiments.

Enzyme-linked immunosorbent assay (ELISA)
The lower compartment medium was collected and stored at −20°C. TSLP was quantified using the human TSLP kit (Thermofisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions.
RNA extraction and sequencing RNA extraction and sequencing were performed by TIANGEN Biotech (Beijing) Co, Ltd. RHEs were cut from the inserts with a surgical knife and soaked in trizol (Thermofisher Scientific, Waltham, MA, USA) for RNA extraction. The concentration and purity of the RNA were determined using the Nanodrop 2000 spectrophotometer (Thermofisher Scientific, Waltham, MA, USA), and the integrity of the RNA was analyzed using the Agilent 2100 bioanalyzer.
RNA containing polyA in total RNA was enriched using the TIANSeq mRNA Capture Kit (TIANGEN, Beijing, China). TIANSeq Fast RNA Library was used to construct the transcriptome sequencing libraries sequenced on Illumina Novaseq 6000 platform, obtaining 150 bp paired-end reads. HISAT2 [23] was used to map the clean reads obtained from the raw files in fast format to the reference genome. The fragments per kilobase of the exon model per million mapped fragments (FPKM) were calculated using the featureCounts software [24] for subsequent analysis. The detailed versions of the software used in the study are listed in Table S1. RNA-seq data are available through NCBI with the identifier SRR24941377-82.
Protein preparation and tandem mass tag quantification RHE samples were quickly frozen in liquid nitrogen, placed in PASP lysis buffer (100 mM NH4HCO3, 8 M urea, pH 8), and ultrasonicated on ice for 5 min. The lysate was centrifuged at 12000 g for 15 min at 4°C and supernatant was added 10 mM DTT for 1 h followed by alkylation with IAM for 1 h in the dark. Samples were mixed with 4 times the volume of cooled acetone and incubated overnight at −20°C, then centrifuged at 12000 g for 15 min at 4°C and the precipitate was collected, which was dissolved in dissolution buffer (8M urea, 100 mM TEAB, pH 8.5). The amount of sample solution was measured by the Bradford method and the quality of protein was determined by 12% SDS-PAGE.
The dissolved protein samples were collected and the volume was made up to 100 μL with dissolution buffer. The sample was mixed with trypsin and 100 mM TEAB buffer and digested overnight at 37°C. Labeling samples were desalted and lyophilized after being treated with acetonitrile-dissolved TMT labeling reagent. A C18 column (Waters BEH C18, 4.6 × 250 mm, 5 μm) on a Rigol L3000 HPLC system was used to separate the fractions. Finally, fractions were analyzed using the EASY-nLC 1200 UHPLC system coupled to the Q Exactive HF-X mass spectrometer operating in data-dependent acquisition (DDA) mode. Proteome Discoverer 2.4 was used to identify and quantify the peptides and proteins. Mass spectrometry proteomics data were deposited with the ProteomeXchange Consortium via the partner repository PRIDE [25] with dataset identifier PXD043085.

Identification of differentially expression genes (DEGs) and proteins (DEPs)
The DESeq2 package was used to identify DEGs with criteria |log2Foldchange|>2 and adjusted P value < 0.05. DEPs were analyzed with a score of logFoldchange > 1.2 or logFoldchange < 0.83 and false discovery rate < 0.05. All DEGs and DEPs were performed using normal RHEs as reference.

Bioinformatic analysis
Principal component analysis (PCA) was performed using the Factoextra R packages with FPKM values for RNA-Seq and mass spectrometry intensities for TMT. The clusterProfier R package was used for gene enrichment analysis of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways based on the org.Hs.eg.db database. Metascape (https:// metascape.org/) was used for the analysis of protein GO and KEGG. Heatmap was created using complexHeatmap and Circlize packages, while other plots were created using ggplot2. The detailed versions of the packages used in the study are listed in Table S2.
Protein-protein interaction (PPI) network and hub gene analysis A PPI network of selected symbols was plotted using the Search Tool for the Retrieval of Interacting Genes (STRING) (http:// string-db.org) to assess all protein interaction relationships. The PPI network was visualized using Cytoscape software (version 3.9.2). The Cytohubba Apps plugin was used to analyze the hub genes that play key roles in the PPI network by calculating the hubba top 10 degree nodes. The PPI network graph was drawn and optimized using Cytoscope.
Quantitative real-time PCR (qRT-PCR) Relative mRNA expression of specific genes was determined by qRT-PCR. Total RNA was isolated using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Isolated RNA was reverse transcribed into cDNA using the FastKing cDNA RT kit (TIANGEN, Beijing, China), and aliquots were stored at −20°C. Quantitative PCR was performed using commercially available reagents and the LC 96 qPCR system with specific primers (Table S3). The results obtained were normalized to the level of the glyceraldehyde-3-phosphate dehydrogenase (GADPH) gene, and the relative expression of the gene was calculated by 2 −ΔΔCt .

Statistical analysis
Differences between the means of variables containing ELISA, qPCR, and western blot measured in triplicates of normal and cocktail-induced RHE models were tested with a student's t-test using GraphPad Prism 9.0. P < 0.05 was considered statistically significant in each experiment.

Results
Compromised reconstructed epidermis and ELISA assay To mimic the inflammatory milieu observed in AD, we selected four factors (poly I:C, TNF-α, IL-4, and IL-13) as the inflammatory cocktail according to previous literature [21]. Our previous results of hematoxylin-eosin (HE) stain indicated that the living layers of the cocktail-induced PC group were looser and showed spongiosis compared with normal RHEs, NC group [26]. To confirm the relationship between damaged RHEs and the inflammatory content observed in AD patients, we performed ELISA analysis to quantify the secretion of TSLP. As shown in Figure 1, the inflammatory cocktail significantly increased the level of TSLP in the PC group compared with the NC group (P < 0.01).

Sample distribution and correlation
The RNA sequencing depth of this study was approximately 6G with three replicates per group. The total quantification of peptides in TMT proteomic analysis was 46,530. The violin plot and box plot showed the distribution of gene and protein expression using FPKM values and relative quantitative peptide spectrum matches (PSM) (Figures S1 and S2). The correlation coefficients of FPKM values and quantitative PSM values in each sample were calculated. The similarities between different groups were visualized using the pheatmap and circlize packages. As shown in Figures S1 and S2, the RHEs exhibited excellent biological duplication under rooster tail-induced (PC group) and NC group, indicating the high reliability of the experimental design.
PCA analysis of recombinant human epidermal (RHE) samples PCA was used to analyze the differences between the six RHE samples. Based on genes and proteins, the distribution plots of PCA for the different samples were shown in Figure 2A and 2B, respectively. The variance contributions of the first principal components were 88.26% and 65.39% at the gene and protein levels, respectively. The RHEs were grouped in two separate clusters in both graphs, which clearly separated the NC group from the PC group and also showed a clear grouping of biological triplicates. This result suggests that the cocktail-induced RHEs have significant differences in terms of gene and protein expression.

Analysis of transcriptomic sequencing expression
Based on the threshold criteria of |log2foldchange|>1 and padj < 0.05, we identified 2963 DEGs from the screening of 13,482 genes in this study, with 1178 upregulated genes and 1785 downregulated genes ( Figure 3A). Some representative genes, such as FLG, IVL, loricrin (LOR), AQP3, and keratin 10 (KRT10) were labeled. The heatmap diagram of DEGs showed optimal clustering and well-controlled biological replications among the PC group and NC group ( Figure 3B). As shown in Figure 3C, the 30 terms of GO were selected to create a GOcircle plot of the biological process (BP), molecular function (MF), and cell component (CC). At BP, DEGs were predominantly associated with extracellular matrix (ECM) organization, extracellular structure organization, skin development, and epidermis development. At MF, DEGs were predominantly associated with growth factor binding, growth factor receptor binding, and receptor ligand activity. At CC, DEGs were predominantly associated with collagen-containing ECM, basement membrane, and cornified envelope. For KEGG assessment of DEGs, the clusterProfile package was used, and the top 20 KEGG pathways were selected based on their importance and relevance to the study ( Figure 3D). The phosphatidylinositol 3-kinase protein kinase B (PI3K-Akt) pathway, cytokine-cytokine receptor interaction, focal adhesion, Janus kinases signal transducer and activator of transcription proteins (JAK-STAT), peroxisome proliferator-activated receptors (PPAR) pathway, and ECM-receptor interaction were enriched at the gene level.

Differential analysis of protein expression
After scanning the protein database, we identified a total of 4673 proteins, including 643 upregulated proteins and 237 downregulated proteins, using fold-change > 1.2 or foldchange < 0.83 and false discovery rate < 0.05 as screening standards. The volcano and heatmap plot showed the different distribution of DEPs in Figure 4A and 4B. Some representative proteins, such as FLG, IVL, AQP3, and S100 calcium-binding protein A9, had been labeled. The top 20 terms of GO showed that DEPs were enriched in mitochondrial protein-containing complex, amide metabolic process, membrane organization, cell-substrate compound, fatty acid metabolic process, and skin development, as shown in Figure 4C. The dot plot of KEGG showed that DEPs were mainly enriched in oxidative phosphorylation, ribosome, nucleotide metabolism, T-cell receptor signaling pathway, ECM-receptor interaction, PPAR signaling pathway, and tricarboxylic acid (TCA) cycle ( Figure 4D).

Integrated transcriptomic and proteomic analysis
We combined the results of transcriptomic and proteomic data to gain a deeper understanding of the AD in RHE models. Compared with the normal RHEs, a total of 13,482 genes and 4673 proteins were identified in the PC groups by the differential expression analysis of transcriptome and proteome, of which 4295 factors were found to be regulated at both the RNA and protein levels ( Figure 5A). Of the 4295 factors, 248 factors were consistent at both the RNA and protein levels, including 118 upregulated and 130 downregulated. At the protein level, there were 520 upregulated and 103 downregulated, and at the RNA level, there were 1059 upregulated and 1650 downregulated ( Figure 5B). As shown in Figure 5C, the Pearson correlation coefficient of the log2 fold-change values at the protein and gene levels was 0.61, indicating a high correlation trend at the mRNA level and the corresponding protein level. There were multiple differences (log2 values) at different levels, showed in nine quadrant diagrams ( Figure 5C). Red dots indicated the same trend in upregulation and downregulation of proteins/genes. The presence of purple dots indicated differential changes at the protein level without corresponding changes at the gene level. The green dots indicated differential changes at the gene level without corresponding changes at the protein level. The gray dots represent situations in which no significant difference was observed at both the protein and gene levels. In the results of the integrated transcriptomic and proteomic analysis results, the top 20 terms, including BP, CC, and MF, were shown in Figure 5D. The genes and proteins with a similar trend were enriched in the cornified envelope, cell-cell junction, cell-substrate junction, defense response   regulation, calcium ion binding, organic hydroxy compound metabolic process, and lipid biosynthesis process. After KEGG pathway enrichment analysis and annotation with the Metascope database, the top eight pathways were identified from the results of proteomic and transcriptomic analysis ( Figure 5E). ECM-receptor interaction, terpenoid backbone biosynthesis, fluid shear stress and atherosclerosis, glycerophospholipid metabolism, glycerolipid metabolism, PPAR signaling pathway, Qiao   and amino acid metabolism were enriched with genes showing the same expression trends in the proteome and transcriptome.

PPI network and hub gene identification
The PPI network with overlapping symbols containing 118 common upregulated genes and 130 common downregulated genes was generated from the STRING database. The hub genes from the PPI network were analyzed using the degree algorithm in the CytoHubba app plugin, and the top 10 hub genes were ranked based on the degree scores (Table 1), including signal STAT3, integrin subunit beta 1 (ITGB1), FLG, IVL, DEAD (Asp-Glu-Ala-Asp) box polypeptide 58 (DDX58), small proline-rich protein 1B (SPRR1B), interferon induced with helicase C domain 1 (IFIH1), desmoglein 1 (DSG1), collagen type XVII alpha 1 chain (COL17A1), and integrin subunit alpha 6 (ITGA6). As shown in Figure 6, the top ten hub genes have been highlighted in red according to their degree ranking (Table 1). Other genes are distinguished based on their log2foldchange values, with circles representing upregulation and diamonds representing downregulation of genes and proteins.

Verification of gene and protein expression
To further validate the accuracy of the omics-related data and the above results, we then used qPCR and western blot methods to detect the expression of relevant genes and proteins under inflammatory cocktail conditions and under NC. The  expression of 13 genes was validated by qPCR. As shown in Figure 7A, the expression of differentiation-associated genes IVL, LOR, FLG, keratin 1 (KRT1), KRT10, SPRR1B, and DSG1 was significantly downregulated (P < 0.01). The expression of toll-like receptor 2 (TLR2), TLR3, STAT3, and IFIH1 genes associated with inflammation and immunity was significantly upregulated ( Figure 7B). In addition, other epidermal-related genes AQP3, periostin (POSTN), and DDX58 were also significantly upregulated (P < 0.01). The observed changes in these genes were consistent with the transcriptomic results, which are consistent with previous reports in AD patients. Compared with the NC group, protein expression of FLG (P < 0.05) and IVL (P < 0.01) was decreased in the positive control group of RHE models. Compared with the NC group, the protein expression of AQP3 (P < 0.05), STAT3, and p-STAT3 (P < 0.01) was significantly increased ( Figure 7C and 7D).

Discussion
Replication in vitro of the relative gene expression changes associated with the pathogenesis of AD is a major challenge, mainly because of the complicated interplay of the immune system and environmental factors in individuals. However, the reproducibility and convenience of in vitro models are irreplaceable while maintaining a high degree of representativeness. It is critical to recognize the limitations of RHE models in accurately reproducing the complex immunological response.
First, in the current experiment, we generated AD-like RHE models by adding the inflammatory cocktail (IL-3, IL-4, TNF-α, and poly I:C) to the culture medium for 72 h. The inflammatory factors and their respective concentrations were obtained from the literature [21], and the stimulation duration was based on our previous study [26]. The non-lesional in keratinocytes caused spongiosis in acute eczema [27].
In the transcriptomic results, the mRNA expression of HAS3 was also significantly upregulated, which is consistent with the literature. In the lesional skin of AD, dysregulated expression of a wide range of genes primarily related to keratinocyte activity and T-cell infiltration was observed. In particular, genes associated with Th2 response (IL-4, IL-10, and IL-13) and Th22 response (IL-22) showed prominent alterations in their expression levels [28,29]. TSLP was strongly expressed by keratinocytes in AD patients with skin lesions. The combination of Th2-and TNF-α-induced TSLP production in full-thickness human skin explants, whereas TNF-α alone had no significant effect on TSLP levels [30]. Poly I:C as an agonist of TLR3 ligand mimicking viral double-stranded RNA (dsRNA) could increase TSLP concentration, which could directly activate dendritic cells to secrete Th2-recruiting chemokines [31]. The ELISA result of TSLP showed that the combined inflammatory cocktail effectively stimulated secretion. Secondly, it should be noted that RHEs, being derived solely from keratinocytes differentiation, lack immune system and the ability to distinguish the skin lesional and non-lesional states. The result of volcano showed that downregulation of FLG, IVL, ITGB1, ITGA3, KRT10, KRT1, LOR, DSG1, and DSG2, and upregulation of carbonic anhydrase 2 (CA2), C-C motif chemokine ligand 26 (CCL26), TLR2, suppressor of cytokine signaling 3 (SOCS3), tenascin-C (TNC), and AQP3. These genes have been associated with skin barrier function and signal transduction [32]. GO functional classification annotation showed that the DEGs were enriched in ECM, structural organization, skin development, growth factor, and receptor binding ( Figure 3C). Interestingly, Christian Cole observed GO enrichment in the areas of extracellular space, receptor binding, and defense response. He compared 26 AD cases with ten control cases [7]. Cytokine-cytokine receptor interaction signaling molecules and interaction and signatures associated with diseases such as measles. Downregulated DEGs were enriched in KEGG pathway signaling pathways, including ECM-receptor interaction, focal adhesion cellular community, and PPAR signaling pathway ( Figure S3) [33,34].
Third, the correlation between mRNA and corresponding protein expression levels depends on various regulatory factors and metabolic processes [35]. The volcano proteomics result also showed that the presence of the inflammatory cocktail caused a delay in the synthesis of barrier proteins, including FLG, IVL, and LOR, which was also observed in other AD-like models and in individuals with AD [36,37]. Suprabasin (SBSN), a novel precursor of the epidermal differentiation envelope precursor, showed low expression of SBSN in the stratum corneum of AD patients. Masahiro discovered that there was no effect of IL4/IL-13 on SBSN expression in RHE, but our proteomic data showed downregulation of SBSN, although there was no difference at the mRNA level. This discrepancy could be attributed to the influence of TNF-α and poly I:C [38]. The GO and KEGG of DEPs were performed using the Metascape database. It was found that the downregulated DEPs were enriched in skin development, lipid biosynthetic biosynthesis process, cell-cell adhesion, and desmosome desmosomes ( Figure S4), which are considered important features of AD [39]. In contrast to the transcriptomes, the upregulated DEPs were enriched in the GO terms mitochondrial protein-containing complex, mitochondrial matrix, and peptide biosynthetic process.
Leman et al. found that mitochondrial activity is upregulated in non-lesional AD (NLAD). NLAD keratinocytes exhibited increased mitochondrial oxidation of long-chain fatty acids, and metabolomic analysis revealed upregulation of TCA cycle turnover [40]. Interestingly, the upregulated pathways identified in KEGG pathway analysis, such as oxidative phosphorylation, ribosome, citrate cycle (TCA cycle), and amino acid degradation, are consistent with Geraldine, further supporting our findings. PPARs were a group of nuclear hormone receptors that can be divided into three subclasses: PPARα, PPARβ/δ, and PPARγ [41]. Studies had demonstrated the beneficial effects of PPARα agonists on the maintenance of barrier homeostasis including the upregulation of molecules associated with epidermal differentiation and lipid synthesis [42].
Finally, the integrated analysis of DEPs and DEGs identified 248 overlapping genes that were enriched in GO terms, including cornified envelope, cell-substrate junction [37], regulation of defense response, and calcium ion binding [43]. The KEGG pathway contains ECM-receptor interaction [44], Terpenoid backbone biosynthesis, and lipid metabolism. The elongase of very long chainfatty acid 1 (ELOVL1) was involved in the biosynthesis of long-chain free fatty acids (FFAs). Interestingly, there was a discrepancy between the mRNA and protein expression of ELOVL1, as the mRNA showed no significant changes while the protein was downregulated [22]. In AD, AQP3 is upregulated to maintain barrier stability when FLG expression is low (Figure 7). The top 10 hub gene include STAT3, ITGB1, FLG, IVL, DDX58, SPRR1B, IFIH1, DSG1, COL17A1, and ITGA6. IL-13/ IL-4 signaling is now considered to be the essential core of the pathogenesis of AD [16,32]. The cytokines IL-13 and IL-4 bind to the IL-4Rα/ IL-13Rα1 receptor, triggering downstream JAK1/JAK2 and STAT6/STAT3. Activation of the STAT6/STAT3 axis leads to downregulation of FLG expression, disruption of the skin barrier function, and an increase in TSLP production [45]. In the epidermal keratinocytes of individuals with AD, IL-4 binds to its receptors, IL-4α, and triggers activation of STAT3. This activation subsequently stimulates the skin inflammation and contributes to the progression of AD [46]. In this study, both STAT3 and p-STAT3 showed the significant upregulation induced with triggered by the inflammatory cocktail, which in turn led to downregulation of FLG and IVL, as well as upregulation of AQP3 (Figure 7). Surprisingly, IL-13 did not significantly affect the expression levels of IL-4R and IL-13Rα1 in keratinocytes, but it induced the upregulation of IL-13Rα2 in cells exposed to IL-13 [47]. Under physiological conditions in AD, the homeostasis of skin barrier function depends on the precise expression of barrier-related proteins, intercellular lipids, and desmosomes within the granular and keratinized layers of the epidermis [48]. FLG and IVL play key roles in the skin barrier. SPRR proteins have been identified as crucial components of the cornified envelope, serving as cross-bridging that connect LOR [49]. In our experiment, the addition of poly I:C resulted in the upregulation of DDX58 and interferon induced with IFIH1 [50], both of which are dsRNA sensors and play a critical role in the immunological function of the epidermal barrier. Corneodesmosomes play a central role in maintaining the structural integrity of the stratum corneum, while desmosomal cadherins function as intercellular adhesive junctions that connect keratinocytes through their association with intermediate filaments. Research suggests that IL-4 and IL-13 downregulate DSG1 expression in both the cultured human keratinocyte (HaCaT) cell line and the skin of people with AD, which may contribute to the impaired skin barrier observed in individuals with AD [51]. In patients with AD, ITGA6 shows marked expression on endothelial cells in both lesional and non-lesional skin [52]. It is interesting to observe that inflammatory differentiated keratinocytes show a remarkable decrease in the expression of both undifferentiated markers, tumor protein P63 (TP63) and ITGA6, and differentiated markers, including KRT1 and KRT10, at the transcriptional level [53]. COL17A1 is a protein involved in the formation of hemidesmosomes, which are structures that anchor the epidermis. Understanding the role of COL17A1 in AD may have implications for the development of targeted therapies aimed at restoring the integrity of the epidermal barrier [54]. The comprehensive integration of IL-4, IL-13, TNF-α, and ploy I:C into the RHE model enabled a multidimensional simulation of the immunological aspects of AD, which includes TLR2-, TLR3and TLR4-mediated inflammation. The integrated results provided a deeper understanding of the RHE inflammation model and offered valuable insights for further research and practical applications.

Conclusion
This study provided comprehensive transcriptome and protein profiling in RHE models that mimic the inflammation associated with AD. Integrated transcriptomic and proteomic analysis further revealed alterations in signaling pathways and potential hub genes. These results provide valuable insights into the molecular mechanisms of AD as well as potential biomarkers for screening cosmetic formulations for the treatment of AD.

Conflicts of interest: Authors declare no conflicts of interest.
Funding: Authors received no specific funding for this work.