Peripheral Blood Gene Expression Profile of Infants with Atopic Dermatitis

To enhance the understanding of molecular mechanisms and mine previously unidentified biomarkers of pediatric atopic dermatitis, PBMC gene expression profiles were generated by RNA sequencing in infants with atopic dermatitis and age-matched controls. A total of 178 significantly differentially expressed genes (DEGs) (115 upregulations and 63 downregulations) were seen, compared with those in healthy controls. The DEGs identified included IL1β, TNF, TREM1, IL18R1, and IL18RAP. DEGs were validated by real-time RT- qPCR in a larger number of samples from PBMCs of infants with atopic dermatitis aged <12 months. Using the DAVID (Database for Annotation, Visualization and Integrated Discovery) database, functional and pathway enrichment analyses of DEGs were performed. Gene ontology enrichment analysis showed that DEGs were associated with immune responses, inflammatory responses, regulation of immune responses, and platelet activation. Pathway analysis indicated that DEGs were enriched in cytokine‒cytokine receptor interaction, immunoregulatory interactions between lymphoid and nonlymphoid cells, hematopoietic cell lineage, phosphoinositide 3-kinase‒protein kinase B signaling pathway, NK cell‒mediated cytotoxicity, and platelet activation. Furthermore, the protein‒protein interaction network was predicted using the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database and visualized with Cytoscape software. Finally, on the basis of the protein‒protein interaction network, 18 hub genes were selected, and two significant modules were obtained. In conclusion, this study sheds light on the molecular mechanisms of pediatric atopic dermatitis and may provide diagnostic biomarkers and therapeutic targets.


INTRODUCTION
Atopic dermatitis (AD) is the most common chronic inflammatory skin disease in early childhood. It affects children with a prevalence of up to 20% and adults with prevalence rates of 7-10% (Weidinger et al., 2018). AD is a complex multifactorial disease, thought to result from an interplay between environmental factors, an impaired skin barrier, and immune dysfunctions; however, the overall pathologic mechanisms are still not fully understood (Langan et al., 2020). Although much has been learned about the molecular basis of AD, most investigations have focused on adult AD with years of disease activity that is remarkably different from those of early-onset AD in children. Few studies have profiled skin tissue in infants with AD (Brunner et al., 2019Cole et al., 2014;Esaki et al., 2016). Over the last decade, RNA sequencing (RNA-seq)-based transcriptome profiles have been implemented in identifying transcripts and pathways in many diseases; however, limited studies using this method were performed on skin tissues in AD, particular in children with AD (Brunner et al., 2019;Cole et al., 2014). Only one study conducted has compared transcriptome profiles of both blood and skin tissue in children with AD at various ages up to age 5 years (Brunner et al., 2019).
Given that AD is an early childhood disease that generates a systemic immunological response and that 60% of all cases of AD begin during the first year of life (Bieber, 2008), we aimed to discover signature biomarkers of AD in infants that might help to identify new diagnostic biomarkers and molecular targets for treatment modalities in pediatric AD.
For this purpose, we performed an integrative study comprising RNA-seq transcriptome profile of peripheral blood cells obtained from infants with AD or healthy infants, quantitative RT-PCR, and systems biology analysis.

Analysis of gene expression by RNA-seq
A total of 100 infants with moderate or severe AD in the first year of life and 20 age-matched healthy control infants were initially recruited (McAleer et al., 2019). PBMCs were isolated from 42 patients and 19 controls. RNA samples were extracted, and only 27 samples from patients with AD and 17 controls passed quality control and were used in this study. The use of samples is presented in a schematic flow chart (Figure 1). We performed RNA-seq profiles on PBMCs from randomly selected infants with AD (n ¼ 8) and controls (n ¼ 5) using the Illumina platform. Differential expression analysis was conducted to identify differentially expressed genes (DEGs) between AD and controls on the basis of the following criteria: false discovery rate <0.05 and fold change !1.5. We identified a total of 178 significantly DEGs with 115 upregulations and 63 downregulations in AD PBMCs when compared with those in control PBMCs. Among highly upregulated genes, we identified IL1b, previously shown to be upregulated in the serum of adult patients with AD (Thijs et al., 2018); TNF, a proinflammatory cytokine whose role in the pathogenesis of AD is well-known (Jacobi et al., 2005;Sumimoto et al., 1992); and early growth response genes EGR2 and EGR3, known to have a crucial role in the regulation of the immune system (Li et al., 2012). Other upregulated genes included TREM1, previously shown to be elevated in lesional skin and serum in patients with AD (Suarez-Farinas et al., 2015), and CXCL5, an inflammatory chemokine found to be at elevated levels in the blood of patients with AD (Brunner et al., 2017). Among downregulated genes, we identified IL18R1 and IL18RAP found to be associated with AD (Hirota et al., 2012). We summarized the 10 DEGs randomly selected from our study and compared their expression with published data (Table 1). A complete list of DEGs in blood cells is shown in Table 2.
Validation of RNA-seq data by RT-qPCR To confirm the results of RNA-seq, real-time RT-qPCR was performed to detect the mRNA expression of five randomly selected DEGs in PBMCs from controls (n ¼ 17) and infants with AD (n ¼ 27). As shown in Figure 2a-e, the mRNA levels of all five genes-IL18RAP, IL1b, TNF, TREM1, and EGR3had significant differences between the AD and control groups in accordance with RNA-seq results.

Effect of infant's age on gene expression in the blood of infants with AD
We wondered whether the infant's age has an effect on gene expression in the blood of patients with AD in the first year of life. Previous studies have shown that differences in skin microbiome depend on infant's age in healthy infants (Capone et al., 2011) and AD infants (Nakamura et al., 2020). To elucidate the effect of age, we stratified patients with AD and healthy controls accordingly (0-6 months and 7-12 months) and performed differential expression analysis on RNA-seq data. Table 3 summarizes the DEGs between four groups: patients with AD aged >6 months (n ¼ 5) versus agematched healthy controls (n ¼ 2), patients with AD aged <6 months (n ¼ 3) versus age-matched healthy controls (n ¼ 3), patients with AD aged >6 months versus patients with AD aged <6 months, and healthy controls aged >6 months versus healthy controls aged <6 months. Interestingly, four DEGs that have been identified between AD infants and healthy controls and validated by RT-qPCR in this study (IL1B, TNF, TREM1, and EGR3) were differentially expressed in patients with AD aged <6 months compared with those in age-matched healthy controls, suggesting its unique differential expression in the first 6 months of life in patients with AD. IL18RAP has been shown to be differentially expressed in both age groups between patients with AD and healthy controls, suggesting no effect of age stratification on this DEG in patients with AD (Table 4). Data were further validated by . Four DEGs (IL1B, TNF, TREM1, and IL18RAP) showed significant differential expression affected by age in accordance with RNA-seq results. Expression of EGR3 showed a trend to be affected by age; however, it was not significant as shown by RNA-seq analysis. Altogether, these data indicate that identified DEGs in patients with AD in the first year of life could be affected by age; however, more samples are required to approve this effect.

Gene ontology and pathway analysis of DEGs
To gain insight into the gene ontology (GO) categories of DEGs between the AD group and control group, all DEGs were uploaded to the DAVID (Database for Annotation, Visualization and Integrated Discovery) database. GO analysis contained three categories: biological process (BP), cellular component (CC), and molecular function (MF). In total, 86 significant GO terms were enriched for DEGs identified in the AD group, of which 57 were within the BP category, 21 were within the CC category, and 8 were within the MF category. The enriched BP categories included immune response, inflammatory response, regulation of immune response, leukocyte migration, positive regulation of NF-kB transcription factor activity, cell adhesion, cell surface receptor signaling pathway, and many others. In the category CC, the DEGs were significantly enriched in the plasma membrane, extracellular region, extracellular exosome, cell surface, and cytoplasm. Furthermore, in the category MF, DEGs were mainly enriched in antigen binding, receptor activity, cytokine activity, enzyme binding, and protein binding. The selected pathways significantly enriched in the AD group included hematopoietic cell lineage, NK cell-mediated cytotoxicity, phosphoinositide 3-kinase-protein kinase B signaling pathway, cytokine-cytokine receptor interaction, extracellular matrix-receptor interaction, and immunoregulatory interactions between a lymphoid and a nonlymphoid cell (Figure 4a and b). Further details of the results of the GO enrichment and pathway analyses are provided in Table 5.
Recruitment of patients with AD (n = 100) and healthy controls (n = 20) under 12 months of age PBMCs isolation from patients with AD (n = 42) and healthy controls (n = 19) RNA sequencing profile of patients with AD (n = 8) and healthy controls (n = 5) mRNA validation for patients with AD (n = 27) and healthy controls (n = 17) RNA extraction from patients with AD (n = 36) and healthy controls (n = 19) QC passed RNA from patients with AD (n = 27) and healthy controls (n = 17) Figure 1. The overall framework of study design. AD, atopic dermatitis; QC, quality control.

Construction of protein-protein interaction network and module analysis
To explore interactions among the DEG genes, STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) analysis was applied, and the most important modules were then screened and visualized using Cytoscape software. A protein-protein interaction (PPI) network containing 82 connected nodes (proteins) and 194 interaction edges (interactions of proteins), where the average degree of connectivity (i.e., the average number of neighbors) was 4.732, is presented in Figure 5a. The hub nodes with the greatest number of neighbors (!8), such as TNF, IL-1b, Von Willebrand factor, and ITGB3, were identified (labeled in red in Figure 5a) and analyzed by GO enrichment and pathway analyses ( Table 6). The Kyoto Encyclopedia of Genes and Genomes pathway analysis revealed that the hub genes were involved in the cytokine-cytokine receptor interaction, hematopoietic cell lineage, extracellular matrix-receptor interaction, platelet activation, cell division, and other pathways. In addition, two significant modules with 10 nodes were obtained from the PPI network of DEGs using Molecular Complex Detection (Figure 5b and c). Enrichment analysis suggested that the genes in the first significant module ( Figure 5b) were mainly associated with functional terms in the category BP, including cell division, cell proliferation, and mitotic nuclear division. In the category CC, the genes in this module were significantly enriched in cytosol and nucleus, and in the category MF, the genes were mainly enriched in protein and adenosine triphosphate binding. The genes in the second module ( Figure 5c) were significantly enriched in inflammatory response, chemokine-mediated signaling, platelet degranulation and activation, immune response, and signal transduction in the category BP. In the category CC, the genes were significantly enriched in the extracellular region, extracellular space, and platelet alpha granule lumen, and in the category MF, the genes were mainly enriched in chemokine activity and CXCR chemokine receptor binding. Furthermore, results from the Kyoto Encyclopedia of Genes and Genomes analysis showed that the genes in this significant module were associated with chemokine signaling pathway and cytokine-cytokine receptor interaction (Table 7).

DISCUSSION
AD is a complex disease associated with immunological and epidermal barrier dysfunctions. Most of our knowledge in the field of AD is based on studies performed in adult patients with AD, although remarkable differences between pediatric and adult AD have been shown recently. Therefore, it is of great importance to identify the molecular basis of pediatric AD and elucidate the biomarkers that could help to identify young patients at risk at an earlier stage of life and to explore new therapies in pediatric AD. Given that more than half of all cases of AD begin during the first year of life, we aimed to   discover signature biomarkers of AD in infants. Considering that skin biopsies are very difficult to obtain at such a young age, that AD generates a systemic immunological response, and that blood is a noninvasive source of biological tissue, we analyzed blood profiles of pediatric patients with AD in the first year of life. Using RNA-seq transcriptome profile of peripheral blood cells obtained from infants with AD or healthy infants, we identified 178 genes differentially expressed in pediatric patients with AD: 115 were upregulated, and 63 were downregulated. To further investigate the functions of the DEGs, GO functional annotation and pathway enrichment analysis were used on the basis of the DAVID database. The GO analysis showed that DEGs were associated with immune responses, inflammatory responses, regulation of immune responses, and platelet activation, which are all known to be AD related. The results of the pathway analysis indicated that the DEGs were enriched in immunoregulatory interactions between a lymphoid and a nonlymphoid cell, hematopoietic cell lineage, phosphoinositide 3-kinase-protein kinase B signaling pathway, cytokine-cytokine receptor interaction, NK cell-mediated cytotoxicity, and platelet activation.
Randomly selected DEGs were further validated in a larger number of samples collected from patients with AD in the first year of life using RT-qPCR. Among highly upregulated genes, we identified IL1b, previously shown to be upregulated in the serum of adult AD (Thijs et al., 2018) and stratum corneum of our pediatric AD collection as reported previously (McAleer et al., 2019). It has been shown to be involved in AD development (Bernard et al., 2017). IL-1b is a potent proinflammatory cytokine that can mediate inflammatory responses by supporting T-cell survival, upregulation of the IL-2 receptor on lymphocytes, enhancing antibody production of B cells, and promoting B-cell proliferation and T helper 17 cell differentiation (Lamkanfi et al., 2011). IL-1b activity is regulated at multiple levels, one of which is controlled by inflammasomes (Schroder and Tschopp, 2010). Recent findings suggest that inflammasome-dependent IL-1b activation plays a role in a variety of disorders, including AD. Of note, among upregulated DEGs, we identified NLRP3, one of the important inflammasome proteins.
Another interesting upregulated DEG in our pediatric patients was TREM1. Recently, it has been reported to be highly expressed in lesional skin and serum of adult AD (Suarez-   Farinas et al., 2015). It has also been reported to be expressed in psoriasis and has been suggested to be a therapeutic target to modify the effects of inflammatory myeloid dendritic cells in psoriasis (Hyder et al., 2013). TREM1 (CD354) is a cell surface receptor that is expressed on various types of cells: monocytes, neutrophils NK cells, dendritic cells, and B and T cells and has been implicated in innate and adaptive immune responses. Activation of TREM1 was shown to result in the production of a variety of inflammatory cytokines, including TNF, IL6, MCP1, and IL-1b, and amplification of toll-like receptor-initiated inflammation (Roe et al., 2014). Of interest, TNF was highly expressed in the blood cells of our pediatric AD collection. It has been shown to be involved in inflammatory processes in AD (Jacobi et al., 2005;Sumimoto et al., 1992). Furthermore, TNF together with the T helper 2 cytokines induced AD-like features in a skin model (Danso et al., 2014). In addition, TNF together with TNF-like weak inducer of apoptosis induced keratinocyte apoptosis in AD skin (Zimmermann et al., 2011). Another interesting group of genes found to be upregulated in the blood of pediatric patients with AD included early growth response genes (EGR1, EGR2, and EGR3), a family of zinc-finger transcription factors. EGR1, an important player in the regulation of cell growth, differentiation, cell survival, and immune responses, has been reported to be upregulated in psoriatic skin lesions (Jeong et al., 2014). EGR2/3 is known to play a crucial role in the regulation of the immune system. They control inflammation, regulate B-and T-cell function in adaptive immune responses, and have been suggested to be involved in preventing the development of autoimmune disease (Morita et al., 2016) and limiting immunopathology during productive adaptive immune responses (Li et al., 2012). Notably, EGR2 is located in a susceptibility locus for AD identified by GWAS in the Japanese population (Hirota et al., 2012).
Among the downregulated genes, we identified IL18R1 and IL18RAP, also previously found to be associated with AD (Hirota et al., 2012). IL18RAP enhances the IL18-binding activity of the IL18 receptor (IL18R1) and plays a role in signaling by IL-18, a pleiotropic immune regulator. IL-18 plays a strong proinflammatory role by inducing IFN-g. Previous studies have implicated IL-18 in the pathogenesis of AD. It has been shown to contribute to the spontaneous development of AD-like skin lesions in a transgenic mouse model (Konishi et al., 2002). It has been reported to be elevated in skin lesions of adults with AD. In our previous study, we analyzed plasma and stratum corneum biomarkers in this collection of patients and showed that IL-18 was observed in very high levels in the stratum corneum of pediatric patients; however, no difference was observed in IL-18 plasma levels (McAleer et al., 2019). Another study showed that PBMCs from patients with AD have a decreased IL-18 expression and capacity to produce IFN-g, which is inversely correlated with serum IgE concentrations (Higashi et al., 2001), suggesting an IL-18 role in the skewing of the immune system in patients with AD. Validation of RNA-sequencing data by RT-qPCR. RT-qPCR analyses for five genes from the top 10 differentially expressed genes identified by highthroughput RNA sequencing: (a) IL18RAP, (b) IL1b, (c) TNF, (d) TREM1, and (e) EGR3 in children with AD (n ¼ 27) and healthy controls (n ¼ 17). Fold change was calculated by 2-DDCT method. The normalized expression data were log 2 transformed and shown as the means AE SD. Significant difference among groups was calculated by unpaired t test with Welch's correction for normal distribution or with Mann-Whitney rank-sum test for non-normal distribution data. *P < 0.05, **P < 0.01, and ***P < 0.001                                      Another gene significantly downregulated in AD infants was GLDC. GLDC, glycine metabolism and the metabolic enzyme glycine decarboxylase, is a key enzyme of the mitochondrial glycine cleavage system (Hiraga and Kikuchi, 1980). GLDC plays important role in many human cancers (Zhang et al., 2012). It has been shown to be differentially expressed in psoriatic skin (Rittie et al., 2016). Interestingly, GLDC has been reported to be differentially expressed in AD-like-reconstructed human epidermis (Evrard et al., 2021), suggesting its involvement in AD development.
A PPI network among the screened DEGs was predicted. The PPI analysis allowed us to determine significant modules and hub genes. In the resulting PPI network, 18 hub genes with the highest degree of connectivity were selected, which included IL1b, Von Willebrand factor gene VWF, PF4, ITGB3, ITGA2B, APP, F5, AURKB, SKA3, MELK, CDC20, PPBP, NCAPG, GTSE1, KIF2C, GP1BA, UBE2C, and TNF. Pathway analysis revealed that the hub genes were involved in the cytokine-cytokine receptor interaction, hematopoietic cell lineage, extracellular matrix-receptor interaction, cell division, platelet activation, and other pathways. In addition, two significant modules were identified. The genes in the first significant module were mainly associated with cell division, cell proliferation, and mitotic nuclear division. The genes in the second module were significantly enriched in inflammatory response, chemokine-mediated signaling, platelet degranulation and activation, immune response, and signal transduction and were associated with chemokine signaling pathway and cytokine-cytokine receptor interaction.  Patients with AD and HCs were stratified by age (0-6 months and 7-12 months), and differential expression analysis on RNA-seq data was performed between the four groups: children with AD aged >6  We wondered whether the hub genes could be linked to AD or other skin inflammatory diseases. The important role of IL-1b and TNF in AD has been shown earlier in this report.
Von Willebrand factor, a key player in hemostasis, has been reported in relation to cutaneous inflammation (Hillgruber et al., 2014). Increased expression of PF4 has been Shown is the effect of infant's age on the expression of selected genes that were differentially expressed in infants with AD and validated by RT-qPCR. Differential expression analysis between patients with AD and HCs was performed using DESeq2 on RNA-seq data after stratifying the dataset by age (0-6 months and 7-12 months). Abbreviations: AD, atopic dermatitis; AIF, apoptosis-inducing factor; DEG, differentially expressed gene; FC, fold change; FDR, false discovery rate; HC, healthy control; LFC, log 2 fold change; OSM, oncostatin M; RNA-seq, RNA sequencing; XAF1, XIAP-associated factor 1. , and (e) EGR3 in children with AD (n ¼ 27) and healthy controls (n ¼ 17) after stratifying the dataset by age (0-6 months and 7-12 months). Fold change was calculated by 2-DDCT method. The normalized expression data were log 2 transformed and shown as the means AE SD. Significant difference among groups was calculated by unpaired t test with Welch's correction for normal distribution or with Mann-Whitney rank-sum test for non-normal distribution data. *P < 0.05, **P < 0.01, ***P < 0.001. AD, atopic dermatitis; HC, healthy control.

Patients
We recruited infants aged <12 months with moderate-to-severe AD who were treatment naive (apart from the use of emollients and 1% hydrocortisone cream or ointment) along with age-matched healthy controls. The study was approved by the Research Ethics Committee of Children's Health Ireland at Crumlin (Dublin, Ireland) and was conducted in compliance with the Helsinki Declaration. Written informed consent was given by parents or legal guardians for all study subjects. The age of onset of AD was recorded. Severity was assessed by the SCORing Atopic Dermatitis index. Clinical and demographic features are summarized in Table 8 ICAM1, IGLV1-44, IGKV5-2, NCR1, IGLV2-11, IGLV2-8,  IGLV3-19, IGKV1-5, IGLV2-23, IGKV4-1, KLRF1, TREM1,  KIR2DL3, TREML1, IGLC2,  DEGs with a significant change between children with AD and healthy control children (cutoff FC !1.5 and FDR <0.05) were used for GO enrichment and pathway analyses using DAVID database. KEGG and Reactome pathway analyses were used to determine the pathways of DEGs between two groups. Abbreviations: AD, atopic dermatitis; Akt, protein kinase B; Ca 2þ , calcium ion; DEG, differentially expressed gene; ECM, extracellular matrix; FC, fold change; FCERI, Fc-epsilon receptor; FDR, false discovery rate; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MMP, matrix metalloproteinase; PI3K, phosphoinositide 3-kinase.

PBMC preparation and RNA isolation
PBMCs were isolated from whole blood as previously described (Nousbeck et al., 2021) using histopaque double-gradient density centrifugation (Sigma-Aldrich, St. Louis, MO) and cryopreserved for further analysis. Total RNA was isolated from PBMCs according to RNeasy Mini Kit protocol (Qiagen, Hilden, Germany). RNA concentrations, RNA integrity, and quality of RNA were evaluated using Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA) and RNA 6000 Nano Lab Chips on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). RNA samples with optimal RNA integrity number values (!8) were considered to construct libraries for sequencing.
RNA-seq, data processing, and differential expression analysis Library preparation (using Illumina TruSeq stranded mRNAseq library kit) and sequencing were conducted by Edinburgh Genomics, The University of Edinburgh (Edinburgh, United Kingdom). The sequencing of the libraries was performed with Illumina NovaSeq 6000 (100 cycles, 50 base pair paired-end sequencing). Sequencing reads showed excellent quality, with the overall Q30 above 94%. After sequencing, reads were trimmed using Cutadapt (Martin, 2011), and clean paired-end reads were mapped to the human reference genome GRCh38 using STAR software (Dobin et al., 2013). The number of reads for each gene was counted using fea-tureCounts (Liao et al., 2014), and the count matrix was used for differential expression analysis. Differential expression was performed using package DESeq2 in R software (version 3.5.2), considering an expression >20 read counts in at least 25% of the samples, a cutoff of at least 1.5-fold change in expression, and a Benjamini-Hochberg-corrected false discovery rate < 0.05.  groups. Any GO terms and pathways with P < 0.05 were considered significantly enriched.

Construction of PPI network and module analysis
Associations between DEGs were investigated using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Szklarczyk et al., 2017) (STRING10.5, http://string-db.org/), and a confidence score >0.6 was considered to indicate significance. Cytoscape software (version 3.6.1) was then used to visualize the PPI network (Shannon et al., 2003). In the network, nodes represented proteins, and lines (edges) represented the interactions. In addition, the most significant modules were identified with the plug-in Molecular Complex Detection (version 1.5.1) with the following settings: degree cutoff of 2, node score cutoff of 0.2, k-core of 2, and a maximum depth of 100, and they were identified by the following criteria: Molecular Complex Detection score >5 and number of nodes >5. Finally, the hub genes in the PPI network were determined, defined as those with a degree of connectivity !8.