Multi-omics integration reveals a core network involved in host defence and hyperkeratinization in psoriasis

Abstract

Objectives: The precise pathogenesis of psoriasis remains incompletely explored. We aimed to better understand the underlying mechanisms of psoriasis, using a systems biology approach based on transcriptomics and microbiome profiling.

Methods: We collected the skin tissue biopsies and swabs in both lesional and non-lesional skin of 13 patients with psoriasis, 15 patients with psoriatic arthritis and healthy skin from 12 patients with ankylosing spondylitis. To study the similarities and differences in the molecular profiles between these three conditions, and the associations between the host defence and microbiota composition, we performed high-throughput RNA-sequencing to quantify the gene expression profile in tissues. The metagenomic composition of 16S on local skin sites was quantified by clustering amplicon sequences and counted into operational taxonomic units. We further analysed associations between the transcriptome and microbiome profiling.

Results: We found that lesional and non-lesional samples were remarkably different in terms of their transcriptome profiles. The functional annotation of differentially expressed genes showed a major enrichment in neutrophil activation. By using co-expression gene networks, we identified a gene module that was associated with local psoriasis severity at the site of biopsy. From this module, we found a 'core' set of genes that was functionally involved in neutrophil activation, epidermal cell differentiation and response to bacteria. Skin microbiome analysis revealed that the abundances of Enhydrobacter, Micrococcus and Leptotrichia were significantly correlated with the genes in core network.

Conclusions: We identified a core gene network that associated with local disease severity and microbiome composition, involved in the inflammation and hyperkeratinization in psoriatic skin.

Keywords: gene network; hyperkeratinization; multi-omics; psoriasis.

FIGURE 1

Differential expression profile in skin transcriptome: (A) principal component analysis (PCA) of all samples; (B) volcano plot of differentially expressed genes (DEGs) in psoriatic lesion compared to non‐lesion. The red dots are the genes with both FDR smaller than .05 and log2FoldChange greater than 1. The green dots are the genes with both FDR smaller than .05 and log2FoldChange smaller than −1; (C) expression profile of psoriasis markers DEGs. Gene expression values were normalized with median of ratios method; (D) top 10 gene ontology (GO) annotations for the DEGs; (E) top 10 Reactome annotations for the DEGs; (F) gene set enrichment analysis (GSEA) of the gene expression profile for the enrichment of neutrophil counts and function; (G) heat map of the gene expression profile of neutrophil activation; (H) expression profile of neutrophil markers; (I) histological section of psoriatic lesion (left: 10×; right 40×). Haematoxylin and eosin (H&E) staining demonstrates an accumulation of neutrophils in the skin.

FIGURE 2

Weighted gene co‐expression network analysis (WGCNA) construction for gene expression profile and module–trait relationship in lesional samples: (A) hierarchical clustering tree (dendrogram) and gene co‐expression module definition in lesional samples. A total of eight modules were identified; (B) gene ontology (GO) annotation for the modules; (C) table of module–trait (clinical parameter) correlations and p values. Each cell reports the correlation (and p value) resulting from correlating module eigengenes (rows) to traits (columns). The table is colour‐coded by correlation according to the colour legend; (D) the significant relationship between modules and phenotypic variables Psoriasis Area Severity Index (PASI) score and PSI score. The results are expressed as the means ± SD. GS, gene significance

FIGURE 3

Core network exploration in green module: (A) core network in green module: The clusters with the highest adjacency in heat map (two triangles) were extracted and visualized with cytoscape. Genes for neutrophil activation (light blue), epidermal cell differentiation (green) and response to bacteria (yellow); (B) protein‐protein interactions network for core network genes; (C) gene regulatory network of CRABP2 and its target genes only showed the genes overlap with core network. Gene regulatory analysis was based on random forest algorithm.

FIGURE 4

Validation of core network signature with public datasets: (A) normalized enrichment score (NES) of core network across six independent bulk RNA‐Seq datasets. *p < .05, **p < 0.01, ****p < 0.0001; (B) volcano plots of six independent bulk RNA‐Seq datasets. The red dots are the genes with both FDR smaller than .05 and log2FoldChange greater than 0. The green dots are the genes with both FDR smaller than .05 and log2FoldChange smaller than 0. The genes in core network were labelled with their names when they were significantly upregulated; (C) upset plot visualized intersections of genes in core network as differentially expressed genes (DEGs) across six independent bulk RNA‐Seq datasets

FIGURE 5

Core network validation in single‐cell level: (A) Heat map of correlation between the genes in core network in three subsets of keratinocyte (KC) in a single‐cell dataset of psoriatic skin. * Adjusted p < .05; (B) gene expression of CRABP2 in different cell‐types across different conditions. FB, fibroblast; Inf., inflammatory; LC, Langerhans cell; LE, lymphatic endothelium; Mac, macrophage; Mig., migratory; MoDC, monocyte‐derived dendritic cell; Tc, cytotoxic T cell; Th, T helper cell; VE, vascular endothelium

FIGURE 6

Microbiota profile and host–commensal interaction: (A) principal coordinates analysis (PCoA) of skin microbiota based on operational taxonomic units (OTUs) in genus level; (B) the Shannon diversity of skin microbiota in different cohorts; (C) the association between the Shannon diversity and psoriasis duration (left) or age (right); (D) heat map of the significant associations between skin microbiota abundance and clinical parameters for all samples (upper) and for only lesional samples (lower). Colour for association significance (p < .05). Sign ‘+’ and ‘−’ for coefficient; (E) the association between the abundance of Enhydrobacter, Micrococcus, Leptotrichia and core network genes. Size of ribbon size encodes R values. All ribbons shown are for correlations with statistical significance. *p < .05, **p < .01