Multifaceted analysis of cross-tissue transcriptomes reveals phenotype-endotype associations in atopic dermatitis

Abstract

Atopic dermatitis (AD) is a skin disease that is heterogeneous both in terms of clinical manifestations and molecular profiles. It is increasingly recognized that AD is a systemic rather than a local disease and should be assessed in the context of whole-body pathophysiology. Here we show, via integrated RNA-sequencing of skin tissue and peripheral blood mononuclear cell (PBMC) samples along with clinical data from 115 AD patients and 14 matched healthy controls, that specific clinical presentations associate with matching differential molecular signatures. We establish a regression model based on transcriptome modules identified in weighted gene co-expression network analysis to extract molecular features associated with detailed clinical phenotypes of AD. The two main, qualitatively differential skin manifestations of AD, erythema and papulation are distinguished by differential immunological signatures. We further apply the regression model to a longitudinal dataset of 30 AD patients for personalized monitoring, highlighting patient heterogeneity in disease trajectories. The longitudinal features of blood tests and PBMC transcriptome modules identify three patient clusters which are aligned with clinical severity and reflect treatment history. Our approach thus serves as a framework for effective clinical investigation to gain a holistic view on the pathophysiology of complex human diseases.

Fig. 1. Summary of study design.

a This study consists of two parts, i) a cross-sectional part (n; Atopic dermatitis: AD = 115, healthy = 14) and ii) a longitudinal part (n; AD = 30) to elucidate endotypes that are associated with phenotypes in AD. b We focused of two classes of disease phenotypes highlighted by clinical data; i) skin manifestation and ii) longitudinal disease course along with medication history, that were examined in association with endotypes in cross-sectional and longitudinal analysis, respectively. EMR electronic medical records, EASI Eczema Area and Severity Index.

Fig. 2. Compositional analysis of clinical scores highlighted two distinct skin manifestations in AD.

a Separation pattern by multidimensional scaling (MDS) on individual components of EASI across AD patients. Components that are correlated with each other (Pearson r > 0.40) were connected with gray lines. Two major clusters were identified in the aspect of key signs of eczema, among which erythema and induration/papulation are two primary skin manifestations that bear the distinction. Clinical pictures (b) and immunohistochemistry of skin tissue for CD4 (c, target protein was stained in red) in two representative patients who have a score composition that are skewed to either of erythema (upper) or induration/papulation (lower). Upper: a 51-year-old male patient who has erythema-skewed EASI composition (total = 19.6, erythema = 5.2, papulation = 3.4). Lower: a 50-year-old male patient who has papulation-skewed EASI composition (total = 21.0, erythema = 3.0, papulation = 8.4). One slide per patient was stained for one marker protein in histological analysis. Assays for other markers in the same patient samples are shown in Supplementary Figs. 6 and 7.

Fig. 3. General transcriptional characteristics of skin and PBMC in AD.

a Volcano plot with significantly differentially expressed genes (|log2 fold change (log2FC)| ≧ 2 and false discovery rate (FDR) < 0.01 for skin and |log2FC| ≧1 and FDR < 0.05 for PBMC) highlighted in red (up-regulated in AD) and blue (down-regulated in AD) compared to healthy controls (n; AD = 115, healthy = 14). b Gene ontology (GO) terms enriched in differentially expressed genes in AD (FDR < 0.1 and Enrichment score > 0.5). Enrichment score was obtained based on the size of a given gene set in GO terms (see “Methods”). Source data are provided as a Source Data file.

Fig. 4. Inference in ligand-receptor coupling suggests augmented skin-PBMC crosstalk in AD patients.

a Connection map of cytokine–receptor coupling across skin and PBMC in a representative healthy control (left) and AD patient (right). Genes that code cytokines and receptors are aligned along the perimeter of the circles. From the outer layer to the center is the tissue expressing the genes (either skin or PBMC), inferred cell specificity, classification of cytokine or receptor, and the connections between cytokines and its matching receptors. The connections were indicated in different colors according to the classification of direction, i.e. in which tissue the cytokines are produced and on which tissue they act. VEC: vascular endothelia cell, vSMC: vascular smooth muscle cell. b Number of active connections between cytokines and receptors. Connections were enumerated according to 4 classes defined by a sender organ and a receiver organ. Boxplots show median and first and third quartiles, whiskers extending to the highest and lowest values no further than 1.5*interquartile range. Brunner-Munzel rank test, two-sided, **p < 0.01, NS: not significant. c Pearson correlation between number of active connections and clinical index, two-sided. N; Atopic dermatitis: AD = 115, healthy = 14 (biologically independent samples). Source data are provided as a Source Data file.

Fig. 5. Identification and characterization of transcriptional modules from skin/PBMC RNA-seq data.

a Cluster dendrograms of transcripts produced by implementation of WGCNA. Color indicates separation of transcriptional module. Cell type expression and GO enrichment in skin tissue (b) and PBMC (d) analyzed by referring public database. Visualization of gene-gene networks in PC1 top 30 genes from each transcriptome module in skin (c) and PBMC (e). Genes that have eigengene-based connectivity > 0.65 were connected with lines. sModu skin transcriptome module, pModu PBMC transcriptome module, vSMC vascular smooth muscle cell, IRS/seba inner root sheeth/sebaceous gland, VEC vascular endothelia cell.

Fig. 6. Regression analysis revealed differential patterns of modular involvement in erythema and papulation skin manifestation in AD.

a Regression models for the prediction of clinical phenotypes. Adjustment was made for R² in training set with the number of prediction variables. b Predicted dysregulated networks of blood tests and skin/PBMC transcriptome modules contributing to distinct phenotypes. Node size and node frame color represent size and the sign of coefficients for each variable predicted by elastic net regression. sXX skin transcriptome module XX, pXX PBMC transcriptome module XX, ALT alanine transaminase, BUN blood urea nitrogen.

Fig. 7. Prediction performance of regression models on longitudinal dataset.

a Performance of elastic net regression models to predict general disease severity (log2(EASI.total+1)). Models were trained with cross-sectional patient dataset and tested on longitudinal dataset. Adjustment was made for R² in training set with the number of prediction variables. b Trajectories of observed and predicted disease severity (log2(EASI.total+1)) in two representative patients both with successful prediction outcome (left, r = 0.81, p = 2.4E−3) and with unsuccessful prediction outcome (right, r = −0.44, p = 0.20) assessed by two-sided Pearson correlation. Source data are provided as a Source Data file.

Fig. 8. Time series features of disease severity, clinical lab and PBMC transcriptome in each patient in association with history of systemic therapy.

a Schematic of extraction of time series features in 30 AD patients. b Silhouette width plot for identifying the optimal number of patient clusters based on time series features. c–f PCA on 30 AD patients using time series features of blood tests and PBMC transcriptome modules. Color indicates patient clusters defined by k-means (c), the intensity of time series feature (upper; mean, lower; MAC) of 5 variables normalized among patients (d), time series features of clinical severity (e), and history of internal medication (f). g Dynamics of EASI (total), pModu07, pModu09, lymphocyte, neutrophil and WBC as well as period of internal medication in representative patients. MAC mean absolute change, pModu PBMC transcriptome module, WBC white blood cell. Source data are provided as a Source Data file.