Bioinformatics Integrative Analysis of Circadian Rhythms Effects on Atopic Dermatitis and Dendritic Cells

Introduction

Atopic dermatitis (AD) is a chronic inflammatory skin disease that affects approximately 20% of the global population.¹ It manifests as dry skin and severe nocturnal pruritus,² affecting sleep continuity and quality.³ AD also deteriorates the patients’ quality of life and is currently posing a heavy social burden.⁴ Its molecular pathogenetic mechanism, which remains poorly understood, involves complex interactions among genetic, environmental, and immunological factors. The main factors causing disease development include skin barrier dysfunction, immune dysregulation, and microbiology.^5,6 AD is associated with circadian mechanisms, and aberrations in the circadian rhythm homeostasis may play a role in its pathogenesis.⁷

The circadian rhythm is a 24-h endogenous timing oscillating system that regulates various biological functions, including cell proliferation and metabolism and immune cell functions.^8,9 Circadian rhythm disturbances are usually associated with an increased risk for cancer progression and poor outcomes, indicating a tumor-suppressive function of the circadian rhythm homeostasis. In addition, the circadian rhythm plays a role in the tumor immune microenvironment.^10–12 Although circadian rhythm has been widely investigated in the field of cancer in recent decades, its effect on AD has rarely been explored. Circadian rhythm comprises a central circadian rhythm, located in the suprachiasmatic nucleus of the hypothalamus, and peripheral circadian rhythm, which is synchronized with the central rhythm but can also modulate locally on its own.^3,13 Some physiological functions of the cutaneous organs, such as hydration, transepidermal water loss (TEWL), sebum levels, capillary blood flow, temperature, and surface pH, also exhibit circadian rhythm patterning,¹⁴ which are also the physiological basis of some pruritus skin diseases aggravated at night like AD, urticaria, psoriasis, and eczema. Summarizing the current evidence, compared to healthy people, skin barrier impairment (poorer hydration, lower sebum levels, and higher TEWL), higher nocturnal cutaneous blood flow, immune cell and cytokines dysfunction and melatonin, cortisol secretion changes exacerbate pruritus and disrupt sleep in AD patients, which cause the vicious itch/cycle and are closely related to circadian rhythm disturbances.^3,7 However, the specific molecular mechanisms underlying circadian rhythm involvement in AD pathogenesis remain unclear.

This study aimed to clarify the relationship between circadian rhythm and AD and explore potential biomarkers in the etiology of AD by comprehensively analyzing bulk and single-cell RNA-sequencing data from the GSE121212^15,16 (see Table S1), GSE120721^17,18 (see Table S2) and GSE153760^19,20 (see Table S3) datasets.

Materials and Methods

Public Datasets

Bulk and single-cell RNA-sequencing data from AD patients were downloaded from the Gene Expression Omnibus (GEO), including GSE121212, GSE120721, and GSE153760 datasets.

Differential Expression Analysis

Differentially expressed genes (DEGs) were calculated using the R package “edgeR”. “EdgeR” provides a procedure to normalize the data, filter under-expressed genes, and implement a series of statistical methods including empirical Bayes estimation, exact tests, generalized linear models, and quasi-likelihood tests to calculate DEGs. The selection criteria for DEGs were |log2FC| > 1 and p < 0.05.

Profiling of Immune Cells

Immune cells were computed using the computational algorithm CIBERSORT,²¹ which allows the precise calculation of immune cell proportions using bulk RNA-sequencing data.²² The efficacy of CIBERSORT was confirmed by wet experiments.

Circadian Rhythm Gene Expression Levels

The circadian rhythm gene expression levels were estimated using single-sample gene set enrichment analysis (ssGSEA).²³ ssGSEA is a non-parametric, unsupervised approach for estimating the relative levels of a gene list. The specific parameters were set as follows: kcdf = “Gaussian”, tau = 0.25, and abs.ranking = TRUE.

Characterizations of Biological Function and Signaling Pathway

Biological function and signaling pathway were analyzed using the R package “clusterProfiler” (version: 3.18.1),²⁴ which provides a series of approaches to investigate the biological function and signaling pathway in patients. Functional annotation included Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. We also performed a gene set enrichment analysis (GSEA) using “clusterProfiler”. GSEA can identify signaling pathways that are missed in the GO and KEGG analyses. The important parameters included the Benjamini-Hochberg procedure for p-value adjustment, a cutoff p-value of 0.05, a cutoff q-value of 0.2, an nperm of 1000, a minGSSize of 10, and a maxGSSize of 500.

Statistics Analysis

All statistical analyses were performed using R software (version 4.0.1). According to the homogeneity of variance and distribution of data, the independent-samples, t-test or Wilcoxon signed rank test was selected. Pearson’s correlation coefficient was used to compute the correlation between two continuous variables. A p-value < 0.05 was considered statistically significant.

Results

Association Between Circadian Rhythms and AD

Circadian rhythm disruption has been widely reported to involve in comprehensive biological processes, including immune cell function, cell division and cycle, and metabolism.^8,9 Thus, we wondered whether the circadian rhythm was associated with AD. To identify the relationship between circadian rhythms and AD, we compared circadian rhythm gene expression levels between lesional skin (AD group) and site-matched non-lesional skin (control group) samples from 21 AD patients in the GSE121212 and found that the circadian rhythm gene expression levels were significantly lower in lesional samples than in non-lesional samples (p < 0.05; Figure 1A) (see Tables S4 and S5). Specifically, the circadian rhythm gene expression levels refer to the overall expression levels of the panel of circadian rhythm-related genes, and was quantified using the ssGSEA algorithm. Furthermore, we compared circadian rhythm gene expression levels between lesional skin samples from patients with AD and matched normal skin from healthy volunteers, and observed that the circadian rhythm gene expression levels were indeed decreased in lesional samples than in normal skin samples (p < 0.05; Figure 1B).

Figure 1 Association between circadian rhythm and atopic dermatitis. (A) The circadian rhythm gene expression levels are significantly higher in the non-lesional samples than in the lesional samples. (B) The circadian rhythm gene expression levels are higher in the normal samples than in the lesional samples. (C) The volcano plot shows 828 differentially expressed genes between the control and atopic dermatitis groups, including 516 upregulated and 312 downregulated genes (|log2FC| > 1, p < 0.05). (D) The heatmap analysis shows different expression modes between the control and atopic dermatitis groups. (E) The principal component analysis shows different expression modes between the control and atopic dermatitis groups. *p < 0.05.

To identify potential key genes involved in AD development, we conducted a differential expression analysis between the control and AD group from the GSE121212 dataset and acquired 828 DEGs, including 516 upregulated and 312 downregulated genes (|log2FC| > 1, p < 0.05; Figure 1C). Notably, we observed that three circadian genes were significantly expressed in AD. Among them, Aryl hydrocarbon receptor nuclear translocator like 2 (ARNTL2) and nocturnin (NOCT) were upregulated, whereas retinoid-related orphan receptor C (RORC) was downregulated in lesional skin of AD. To test whether a distinct expression mode existed between the control and AD group, we carried out a heatmap analysis and principle component analysis (PCA). Both of them showed distinct expression patterns in the control and AD groups (Figure 1D and E).

Characterization of Differentially Enriched Signaling Pathways in AD

Since previous results showed a distinct expression pattern in AD, we wondered to profile the signaling pathways enriched in the development of AD by performing GO and KEGG pathway enrichment analysis for 516 upregulated genes and 312 downregulated genes between the control and AD groups. The top 10 enriched biological functions were cell chemotaxis, T cell activation, cytokine-mediated signaling pathways, positive regulation of cell-cell adhesion, and leukocyte migration (Figure 2A). Consistently, the top 10 upregulated signaling pathways comprised cytokine-cytokine receptor interaction; IL-17, chemokine, and NF-κB signaling pathways, and Th1 and Th2 cell differentiations (Figure 2B). GSEA demonstrated that allograft rejection, graft-versus-host disease, IL-17 signaling pathway, and primary immunodeficiency were enriched (Figure 2C and D), suggesting that AD is an autoimmune condition.

Figure 2 Characterization of differentially enriched signaling pathways in atopic dermatitis. (A) Top 10 overexpressed and downregulated biological functions in atopic dermatitis samples. (B) Top 10 overexpressed and downregulated signaling pathways in atopic dermatitis samples. (C and D) Top 10 overexpressed signaling pathways generated by the gene set enrichment analysis.

Profiling of Immune Microenvironment of AD

Since the results of characterization of differentially enriched signaling pathways indicated that immunity disorders occurred in AD, we analyzed changes in the abundance of immune cells between the control and AD samples using RNA-sequencing data from the GSE121212 dataset (see Table S6). We estimated the expression values of 22 subsets of immune cells, deleted the subsets of cells with low expression values in both groups, and retained 7 subsets of cells with high expression values that may have potential biological significance for further analysis. The results showed that only myeloid dendritic cells were significantly upregulated in the AD samples compared with in the control samples (p < 0.001, Figure 3A). The heatmap also showed different expression patterns in dendritic cells between the control and AD groups (Figure 3B). The correlation analysis showed that dendritic cells were positively correlated with monocytes and NK cells (Figure 3C). Collectively, our results revealed that myeloid dendritic cells were abnormally upregulated in AD.

Figure 3 Profiling of immune microenvironment of atopic dermatitis. (A) Comparison of the proportion of immune cells between the control and atopic dermatitis samples. (B) The heatmap shows different expression patterns in the dendritic cells between the control and atopic dermatitis groups. (C) The correlation analysis shows that the dendritic cells are positively correlated with the monocyte and natural killer cells. ***p < 0.001. Abbreviation: ns, not significance.

Validation of Three Key Genes in an External Cohort

In the previous step, we revealed that the dendritic cell was aberrantly overexpressed in AD, and we first wondered to identify its relationship with circadian rhythm gene expression levels. As we expected, the results showed that dendritic cells were negatively correlated with circadian rhythm gene expression levels (R = − 0.4, p = 0.0092; Figure 4A). Next, we assessed the association of the three key genes with dendritic cells, and found that ARNTL2 and NOCT were critically positively correlated with dendritic cells (R = 0.83 and R = 0.68, respectively; Figure 4B and C), whereas RORC was negatively associated with dendritic cells (R = − 0.44; Figure 4D).

Figure 4 Validation of three key genes in an external cohort. (A) Dendritic cells are negatively correlated with the circadian rhythm levels (R = − 0.4, p = 0.0092). (B) ARNTL2 expression levels are positively correlated with dendritic cells (R = 0.83, p = 1e-11). (C) NOCT expression levels are positively correlated with the dendritic cells (R = 0.68, p = 6e-07). (D) RORC expression levels are negatively correlated with dendritic cells (R = − 0.44, p = 0.0032). (E and F) ARNTL2 and NOCT expression levels are significantly upregulated in atopic dermatitis group. (G) RORC expression levels are significantly downregulated in atopic dermatitis group. (H and I) ARNTL2 and NOCT expression levels are positively correlated with the severity of atopic dermatitis. (J) RORC expression levels are unchanged in the GSE120721 database. (K–M) The receiver operating characteristic curve shows robust predictive abilities of ARNTL2, NOCT, and RORC genes, with areas under the curve of 0.907, 0.866, and 0.943, respectively. **p < 0.01; ****p < 0.0001. Abbreviations: ns, not significance; AUC, area under the curve; NN, normal tissue; NL, non-lesional atopic dermatitis; LS, lesional atopic dermatitis.

Further, we analyzed the ARNTL2, NOCT, and RORC expression levels in the GSE121212 and GSE120721datasets, respectively. Surprisingly, the expression levels of ARNTL2, NOCT, and RORC were differently expressed in the boxplot (Figure 4E–G). Besides, ARNTL2 and NOCT expression levels elevated as the AD severity increased (Figure 4H–I), although RORC expression levels did not change in the GSE120721dataset (Figure 4J).

We further assessed the predictive abilities of the three key genes for the AD diagnosis using a receiver operating characteristic (ROC) curve. Impressively, three genes all displayed robust predictive abilities in ROC curve, with areas under the curve of 0.907, 0.866, and 0.943, respectively (Figure 4K–M).

Validation of the Relationship Between ARNTL2, NOCT, and RORC and the Dendritic Cells Using Single-Cell RNA-Sequencing Data

To repeatedly validate the association of three key genes with dendritic cells, we investigated their relationships using single-cell RNA-sequencing data of four paired of AD and matched healthy control (HC) samples from the GSE153760 dataset. We clustered 16,435 cells into 10 types of cells, including epithelial cells, monocytes, dendritic cells, ciliated cells, basal cells, NK cells, memory CD4 cells, fibroblast, plasma cells, and goblet cells (Figure 5A). Dendritic cells were increased in patients with AD, although not significantly (p = 0.12; Figure 5B), compared with that in healthy controls.

Figure 5 Validation of the relationship of the three key genes with dendritic cells using single-cell RNA-sequencing data. (A) 16,435 cells from 8 samples are clustered into 10 types of cells. (B) Dendritic cells shows an increased tendency in atopic dermatitis than in healthy controls. (C–F) ARNTL2 and NOCT expression levels are significantly correlated with the abundance of dendritic cells. (G and H) RORC expression levels are not associated with the dendritic cells.

Finally, we sought to validate the relationships between three key genes and dendritic cells. Consistent with previous findings generated from bulk RNA-sequencing data, the ARNTL2 and NOCT expression levels were significantly correlated with the abundance of dendritic cells (Figure 5C–F), whereas RORC expression levels were not (Figure 5G and H).

Discussion

This study aimed to clarify the relationship between circadian rhythms and AD and explore potential biomarkers in the etiology of AD. Subsequently, we revealed the association of circadian rhythms with atopic dermatitis and identified three circadian rhythm-related genes, ARNTL2, NOCT, and RORC, involved in the pathogenesis of AD. Moreover, dendritic cells were significantly elevated in AD and positively correlated with ARNTL2 and NOCT, suggesting their roles in the pathogenesis of AD. These results provide insights into the establishment of circadian rhythm-oriented therapies for AD.

Circadian rhythm disruption has been proposed to be implicated in multiple inflammatory conditions.^25–28 Herein, we revealed and verified its effects on AD and correlation with the immune system. The circadian rhythm regulates the immune system, and 8% of macrophages exert their role in a rhythmic manner, exerting a considerable regulatory effect on the immune system.²⁸ Moreover, circadian rhythm-related genes can regulate the secretion of small molecules, such as IL-6,²⁸ TNF-α,²⁹ and the differentiation of various immune cells, including T-helper cells,³⁰ B-cells,^31,32 macrophages and NK cells. AD is an immune-driven inflammatory disease characterized by Th2 inflammation.³³ In AD patients, circadian rhythm disturbances affect the secretion of cortisol, disrupt melatonin levels and worsen barrier function at night, resulting in an aggravated inflammatory response and the natural rhythm of immune cells and inflammatory cytokines.³⁴ At night, T-cell activation is increased, and the expressions of Toll-like receptors, Th1 cytokines, Th2 cytokines, Th22 cytokines, and chemokines are upregulated, further stimulating the production of inflammatory cytokines,³ such as IL-6 and IL-31, which can exacerbate pruritus and disrupt sleep. These processes in AD patients can aggravate the disease and prolong its course. This study also demonstrated that dendritic cells are implicated in AD. Dendritic cells are major effector cells in AD. Skin barrier dysfunction in AD patients facilitates external environmental substances (microorganisms and allergens) to invade the epidermis, which activates inflammatory epidermal dendritic cells and Th2 inflammatory responses; dendritic cells are involved in this process by presenting allergens.^35,36 Our findings further highlighted its potential effect on the pathogenesis of AD and demonstrated that circadian rhythm disturbances are involved in the disruption of immune functions in the skin.

Another main finding of this study was that three circadian rhythm-related genes ARTNL2, NOCT, and RORC were identified as key genes and predictors of AD diagnosis. ARNTL2 is a member of the PAS superfamily and has been comprehensively reported in human malignancies.^37,38 ARNTL2 overexpression is associated with reduced survival and immune disorders in clear cell renal cell carcinoma,³⁹ and can be a predictor of cancer aggressiveness in lung adenocarcinoma,⁴⁰ colorectal cancer,⁴¹ and colon cancer.⁴² In addition, ARTNL2 may be involved in the functional activity of various immune cells, such as neutrophils, dendritic cells, Th1, Th2, Th17, and mast cells.⁴³ However, its role in AD has not been described. Similarly, NOCT and RORC have not been associated with AD.

In this study, we revealed the roles of ARTNL2, NOCT, and RORC and their association with dendritic cells in AD, indicating their potential roles as therapeutic targets. We hypothesized that modulating the abnormal circadian rhythm and immune response of dendritic cells in AD might be possible by intervening in three key genes, thereby repairing skin barrier function. Currently, we can introduce chronotherapy into the treatment of AD, which is a pattern of drug administration that synchronizes with the circadian rhythm of a condition or symptom of the disease. Properly implemented chronotherapy can provide a stable circadian sleep-wake cycle to correct circadian rhythm disturbances, where we can use such as light, melatonin, and regular work and rest guidance to increase the strength and stability of circadian rhythm genes. The results of this study provide a potential therapeutic target for the future precision medicine treatment of AD, but it still needs to be further explored and realized.

This study had some limitations. The key genes and immune cells were identified using bioinformatic approaches, and further in vitro and in vivo studies are required to investigate the molecular mechanisms underlying their roles in AD.

In conclusion, this study revealed the association between circadian rhythms and AD and several genes and immune cells involved in the progression of AD. The results of this study, might facilitate the development of circadian rhythm-oriented therapies for AD.