Abstract
Biologic drugs targeting the IL-23/IL-17 axis have transformed outcomes for people living with psoriasis, allowing most patients to achieve skin clearance. Despite the potential to reveal fundamental regulators of skin homeostasis, the early mechanisms of action of these drugs remain poorly understood. Here, we performed longitudinal single-cell RNA sequencing in individuals with psoriasis who had been successfully treated with IL-23 inhibitor therapy. By profiling skin at baseline, day 3 and day 14 of treatment, we demonstrated that IL-23 blockade causes marked gene expression shifts, with fibroblast and myeloid populations displaying the most extensive changes at day 3. We identified a transient WNT5A+/IL24+ fibroblast state, which was only detectable in lesional skin. The study of ligand-receptor interactions revealed that pro-inflammatory signals stemming from these WNT5A+/IL24+ fibroblasts upregulated multiple genes in differentiated keratinocytes. Importantly, the abundance of WNT5A+/IL24+ fibroblasts was significantly reduced after treatment. This observation was validated in-silico, by deconvolution of multiple transcriptomic datasets, and experimentally, by RNA in situ hybridization of additional skin samples. These findings demonstrate that the evolution of inflammatory fibroblast states is a key feature of resolving psoriasis skin.
INTRODUCTION
Skin is an archetypal barrier tissue which protects from external insults, while orchestrating innate and adaptive responses to microbial agents (1). Our understanding of these functions has been transformed by advances in single-cell transcriptomics, which have uncovered key differences between the cellular ecosystems of healthy and inflamed skin (2–5). Profiling of skin lesions, however, provides only a snapshot of end-stage disease, with limited insight into the sequence of events leading to the onset and resolution of pathology. In this context, a longitudinal analysis of resolving disease can illuminate early mechanisms of drug action, identifying the fundamental events that promote a return to homeostasis. While this potential is being realised in oncology (e.g., by studying evolving cell states during immunotherapy) (6, 7), inflammatory conditions remain relatively underexplored (8, 9).
Psoriasis is a chronic, immune-mediated skin disorder affecting up to 2% of individuals of European descent. It presents with red, disfiguring skin plaques and is associated with multimorbidity, poor quality of life and increased cardiovascular mortality (10). Excessive activation of the IL-23/IL-17 axis is a key disease driver. IL-23 produced by dendritic cells and macrophages promotes the terminal differentiation of T17 cells, leading to the release of IL-17. This cytokine, which acts in synergy with TNF, causes keratinocytes to proliferate and produce chemokines that attract further T17 cells to sites of inflammation. The ensuing feedback loop sustains chronic inflammation and epidermal hyperplasia, two hallmarks of psoriasis (11).
The IL-23/IL-17 immune axis has been successfully targeted by powerful biologics inhibiting IL-17 or IL-23. IL-23 antagonists have proven particularly efficacious, with rates of skin clearance (remission) exceeding 60% after 1 year (12). However, the early events leading to drug-induced remission are poorly understood. Skin biopsies have been profiled during trials of IL-17 and IL-23 antagonists (13, 14), but these studies have been limited by the poor resolution of bulk RNA-sequencing and the use of samples obtained weeks after the onset of clinical benefit.
Here, we sought to elucidate the early events underlying the effects of IL-23 antagonists, with a view to uncovering key processes regulating skin homeostasis. We performed longitudinal single-cell RNA-sequencing (scRNA-seq) of psoriasis skin, in individuals receiving the IL-23 inhibitor risankizumab. We generated a high-resolution atlas of early disease resolution, identifying a psoriasis specific WNT5A+/IL24+ fibroblast state. This population, which shows a marked upregulation of IL-17 signalling, starts to decline in abundance 3 days after treatment initiation. Importantly, an early reduction in WNT5A+/IL24+ fibroblasts was also detectable in patients receiving other systemic and topical therapeutics, suggesting that the normalisation of this cell state is critical to the resolution of psoriatic inflammation.
RESULTS
Single-cell RNA sequencing identifies key cell populations responding to IL-23 blockade
To investigate the early effects of IL-23 inhibition, we obtained serial, full-thickness skin biopsies from five individuals with psoriasis who had achieved skin clearance (remission) by 16 weeks of risankizumab treatment (Supplementary Table 1, Fig. 1, A and B). We sampled lesional skin at baseline (day 0, when we also obtained non-lesional biopsies) and at two early timepoints (day 3 and day 14 of treatment) before full resolution of inflammation was clinically apparent. Following scRNA-seq, we generated expression profiles from 164,553 viable cells (Supplementary Table 2). We analysed this dataset with Seurat (15), identifying cell clusters corresponding to the main cell populations found in the epidermis (keratinocytes and melanocytes), dermis (fibroblasts, vascular and lymphatic endothelial cells, pericytes, epithelial cells) and in the psoriasis immune infiltrate (mast cells, myeloid cells and T cells) (Fig. 1, C and D; Supplementary Fig. 1).
(A) Experimental design of the study, showing timepoints for the sampling of non-lesional (day 0) and lesional (day 0, day 3, day 14) skin biopsies. (B) Left: Percentage reduction of disease severity (measured by the psoriasis area and severity index, PASI) after 3 and 14 days of risankizumab treatment. Right: representative skin images showing clinical improvement during risankizumab treatment. (C) Unifold Manifold Approximation and Projection (UMAP) of 164,553 single cells (lesional and non-lesional skin), visualised as clusters (left) or skin cell types (right). (D) Dot plot showing the expression of marker genes used for the annotation of cell identities. (E) UMAP visualization of lesional skin showing the number of differentially expressed genes (DEG) observed in each cluster, after 3 (left) and 14 (right) days of treatment. (F) Heatmap depicting the number of ligand-receptor interactions between cell types. The cell types expressing the ligand are labelled on the y-axis, whereas those expressing the receptor are on the x-axis. KC, keratinocytes; FB, fibroblasts; PE, pericytes; LEC, lymphatic endothelial cells; VEC, vascular endothelial cells; EpC, epithelial cells; MEL, melanocytes; MAST, mast cells; MYEL, myeloid cells; TC, T cells.
Differential expression analysis (day 3 vs baseline) showed that the majority of day 3 changes were observed in myeloid cells (25 differential expressed genes, DEG) and fibroblasts (25 DEG), identifying these populations as early responders to IL-23 inhibition (Fig. 1E). Ligand-receptor analyses undertaken with CellChat (16) demonstrated a widespread reduction in the number of cell-cell interaction after 3 days of treatment (Fig. 1F). In keeping with the results of the differential expression analysis, the number of interactions originating from fibroblasts and myeloid cells was markedly reduced.
By day 14, myeloid cells and keratinocytes were displaying the strongest gene expression changes (121 and 115 DEG, respectively) (Fig. 1E). The ligand-receptor analyses were broadly consistent with these observations, showing that the number of interactions driven by myeloid cells was further reduced after 14 days of treatment (Fig. 1F).
Thus, myeloid cells and fibroblasts are the cell types showing the earliest response to IL-23 inhibition, with keratinocyte changes becoming apparent shortly after. Importantly, these cellular shifts are observed before clinical remission becomes apparent.
IL-23 blockade affects cell proliferation and IL-17 signalling in keratinocytes
To investigate the effects of IL-23 inhibition in further detail, we sought to define population subsets for each cluster. In the epidermis, melanocytes could not be further subdivided, but keratinocytes could be classified in five distinct subclusters (Fig. 2A and Supplementary Fig. 2A). These corresponded to basal (KRT15+), follicular (KRT17+), proliferating (MKI67+), spinous (KRT1+/SPRR2E-) and supraspinous (KRT1+/SPRR2E+) populations (Fig. 2B). In keeping with the early clinical improvement of skin lesions (Fig. 1B), the abundance of proliferating keratinocytes was significantly reduced by day 14 (Fig. 2C). The largest number of DEG was observed in supraspinous keratinocytes (n=228 at day 14 vs day 0) and spinous keratinocytes (n=193 at day 14 vs day 0) (Fig. 2D). Pathway enrichment analysis of the DEG identified a significant over-representation of IL-17 related genes (“Role of IL-17A in psoriasis”; FDR <10-4 in supraspinous keratinocytes, FDR<10-2 in spinous keratinocytes), confirming the suppressive effect of IL-23 inhibition on T17 activation. Likewise, an upstream regulator analysis demonstrated a significant enrichment for genes that are induced by IL-17 and its synergistic partner TNF (z< -2.50 for both in supraspinous keratinocytes, FDR<10-5). Conversely, the impact of risankizumab on IFN-γ (another cytokine induced downstream of IL-23) was limited, as the activation score (z= -1.97) did not reach statistical significance (z< -2.0).
(A) UMAP of 68,695 keratinocytes forming five distinct subclusters. (B) Dot plot showing the expression of marker genes used for the annotation of cell identities. (C) Abundance of proliferating cells (expressed as keratinocyte percentage) at different time points after risankizumab treatment. *P<0.05, **P<0.01 (repeated measures one-way ANOVA with Dunnett’s post-test). (D) UMAP visualization of keratinocytes from lesional skin, showing the number of differentially expressed genes (DEG) observed in each cluster, after 3 (left) and 14 (right) days of treatment. KC, keratinocytes; L, lesional skin; NL, non-lesional skin.
Thus, the main effect of IL-23 blockade in the epidermis is a marked downregulation of IL-17 signalling in spinous and supraspinous keratinocytes.
IL-23 inhibition downregulates cytotoxic T cell activity and increases dendritic cells (DC)-3 signalling
Our initial clustering identified three immune cell types: mast cells, T cells and dendritic cells (DC). While mast cells could be subdivided into three subclusters (TPSAB1-/IL1RL1+, TPSAB1+/IL1RL1+, TPSAB1+/IL1RL1-) (Supplementary Fig. 3, A and B), the frequencies of these populations did not change with treatment (Supplementary Fig. 3C). The number of DEG observed in each cluster was also low (<15 at day 14), suggesting that IL-23 inhibition has minimal effects on mast-cell activity.
T cells formed six distinct populations, including T helper (CD40LG+), Treg (CTLA4+/FOXP3+), proliferating Treg (MKI67+/CTLA4+/FOXP3+), terminally differentiated effector memory (TEMRA) T cells (CD8+/CCL5+), cytotoxic T17 cells (CD8+/IL17A+) and activated PLCG2+ cells. A CD3-/XCL1+ cluster was also identified, corresponding to natural killer/innate lymphoid cells (NK/ILC) (Fig. 3, A and B, and Supplementary Fig. 3, D and E). Differential expression analysis identified TEMRA T cells as the population showing the strongest response to IL-23 blockade, with 28 DEG detected at day 3 and 30 at day 14 (Fig. 3C). These included several genes related to effector function (e.g., GZMA and GZMB at day 3; NKG7, KLRB1 and IL7R at day 14), suggesting a marked reduction in cytotoxic activity (Fig. 3D). The analysis of the myeloid compartment identified 10 LYZ+ subclusters. These broadly recapitulated the main mononuclear phagocyte populations described by Reynolds et al in psoriatic skin (5). In fact, we observed cDC1 (CLEC9A+) and proliferating cDC1 (CLEC9A+/MKI67+) cells, two cDC2 populations (CLEC10A+/CD1C+), two DC3 populations (IL1B+/IL23A+), mregDC (LAMP3+/BIRC3+), Langerhans cells (CLEC4A+/CD207+) and two macrophage populations (CD68+/CD163+). A LYZ-subcluster was also detected and identified as a plasmacytoid DC (pDC) population, based on IRF7 expression (Fig. 3, E and F, and Supplementary Fig. 3, F and G).
(A) UMAP of 6,263 lymphocytes forming seven distinct clusters. (B) Dot plot showing the expression of marker genes used for the annotation of T cell identities. (C) UMAP visualization of lymphocytes from lesional skin, showing the number of differentially expressed genes (DEG) observed in each cluster, after 3 (left) and 14 (right) days of treatment. (D) Downregulation of cytotoxicity markers in TEMRA T cells, following treatment with risankizumab. (E) UMAP of 18,544 myeloid cells, forming 11 distinct clusters. (F) Dot plot showing the expression of marker genes used for the annotation of myeloid cell identities. (G) UMAP visualization of myeloid cells from lesional skin, showing the number of differentially expressed genes (DEG) observed in each cluster, after 3 (left) and 14 (right) days of treatment.
Day 14 vs day 0 differential expression analysis identified the two DC3 subsets as the cells showing the strongest response to risankizumab (180 DEG detected in DC3-2 and 155 in DC3-1) (Fig. 3G). DEG included IL23A, CXCL3 and IL1B, all of which were upregulated at day 14. Notably, the expression of TNFAIP3 (a negative regulator of IL-17 and NF-κB signalling) was also increased (Supplementary Fig. 3H), suggesting that the elevation in cytokine gene expression is a transient phenomenon, likely caused by the modulation of negative feedback loops regulating IL-17/IL-23 signalling (17).
IL-23 blockade downregulates mediators of leukocyte diapedesis in pericytes and endothelial cells
While fibroblasts (described below) were the most abundant cell type in the dermis (Supplementary Fig 1B), we also observed epithelial cells, lymphatic endothelial cells, vascular endothelial cells and pericytes. Neither epithelial cells nor lymphatic endothelial cells could be further subdivided, while vascular endothelial cells formed five subclusters (Fig. 4A, Supplementary Fig. 4, A and B). These included arterioles (HEY1+/S100A4+), capillaries (PLVAP+/RGCC+; ACKR1+/FABP4+), post-capillary venules (ACKR1+/SELE+) and venules (ICAM1+/VWF+) (Fig. 4B). Although the number of DEG observed at day 14 was <40 in all populations, there was evidence for a significant downregulation of IL-17 signalling in lymphatic endothelial cells and ACKR1+/FABP4+ capillaries FDR<10-7 in both).
(A) UMAP of 16,420 vascular endothelial cells forming five clusters. (B) Dot plot showing the expression of marker genes used for the annotation of vascular endothelial cell identities. (C) Risankizumab-induced downregulation of proinflammatory molecules and collagen genes, in PLVAP+/RGCC+ capillaries. (D) UMAP of 13,507 pericytes forming two clusters. (E) Dot plot showing the expression of marker genes used for the annotation of pericyte cell identities. (F) UMAP visualization of vascular endothelial cells from lesional skin, showing the number of differentially expressed genes (DEG) observed in each cluster, after 3 (left) and 14 (right) days of treatment. (G) Genes promoting leukocyte extravasation are downregulated in pericyte 1, following risankizumab treatment.
Downregulation of genes promoting leukocyte chemotaxis (S100A7, S100A8, S100A9) and extravasation (COL4A1, COL4A2) was also observed in PLVAP+/RGCC+ capillaries (Fig. 4C).
Pericytes formed two clusters (pericyte 1 and 2), distinguishable based on PDGFRB and MYL9 expression (Fig. 4, D and E, and Supplementary Fig. 4, C and D). At day 14, pericyte 1 cells showed the largest number of DEG (n=63, vs day 0) (Fig. 4F), including several genes implicated in leukocyte chemotaxis (CXCL8, S100A7, S100A8, S100A9) and diapedesis (VIM, ACTA2, SPARC) (Fig. 4G).
Thus, IL-23 blockade causes an early downregulation of the intermediate filaments, extracellular matrix proteins and chemokines that facilitate the migration of leukocytes towards sites of skin inflammation.
IL-23 inhibition reduces the abundance of a psoriasis specific, IL24+ fibroblast population
Fibroblasts formed 11 distinct subclusters (Fig. 5A and Supplementary Fig. 5, A and B). These included CCN5+/PI16+, COL11A1+/DPEP1+ and ANGPTL7+ populations, as well as three COMP+ subclusters, three APOE+ populations and a WNT5A+/IL24+ subcluster (Fig. 5B). While CCN5+/PI16+ cells showed a strong response to treatment at day 3 (143 DEG), the number of differentially expressed genes at day 14 was <60 in all subclusters (Supplementary Fig. 5C). Conversely, IL-23 inhibitor treatment had a marked impact on cell abundance, causing a significant decrease in the frequency of WNT5A+/IL24+ cells. This phenomenon mirrors the absence of WNT5A+/IL24+ cells from non-lesional skin (Fig. 5C), identifying the subcluster as a psoriasis-specific population. Accordingly, an analysis of the 237 genes that were differentially expressed in WNT5A+/IL24+ cells compared to all other fibroblast clusters demonstrated a marked enrichment for the IL-17 signalling pathway (FDR=10-5).
(A) UMAP of 31,765 fibroblasts forming eleven clusters. (B) Dot plot showing the expression of marker genes used for the annotation of cell identities. (C) Plot showing the decline in WNT5A+/IL24+ cell abundance, following treatment with risankizumab. Every line represents a patient. *P<0.05; **P<0.01 (repeated measures ANOVA, with Dunnett’s post-test). (D) Inferred interaction between IL-24 produced by WNT5A+/IL24+ fibroblasts and its receptor (encoded by IL20RA and IL20RB) on spinous keratinocytes. (E) Plot showing the SCENIC regulons detected in WNT5A+/IL24+ fibroblasts. The three regulons with the highest specificity scores are indicated by labels, with their size (number of regulated genes) in brackets. (F) Dot plot showing the expression of the top regulons in the various fibroblast clusters. (G) Boxplots showing fibroblast clusters ordered by pseudotime. L, lesional skin; NL, non-lesional skin; Fib, fibroblast; KC, keratinocyte.
In keeping with the specificity of WNT5A+/IL24+ cells, IL-24 signalling was only detectable in day 0 lesional skin, where ligand-receptor analysis uncovered an interaction between the cytokine (originating from WNT5A+/IL24+ fibroblasts) and its receptor on spinous keratinocytes (Fig. 5D). The effects of risankizumab on IL-24 signalling were independently validated with NicheNet (18), which showed that reduced IL24 levels at day 14 correlated with the downregulation of IL-24 target genes in spinous keratinocytes (Supplementary Fig. 5D). Of note, the DEG included several transcription factors regulating keratinocyte proliferation and differentiation (JUN, FOS, ID1, STAT1), as well as inflammatory mediators (S100A8, S100A9) and structural proteins (DST, KRT15). These observations suggest that IL-24 produced by WNT5A+/IL24+ fibroblasts promotes the activation of spinous keratinocytes and that IL-23 inhibition reverses this process.
To further characterise WNT5A+/IL24+ fibroblasts, we carried out an analysis of gene regulatory networks with SCENIC (19). This demonstrated an enrichment for FOXQ1, TBX3 and PRDM1 regulons in WNT5A+/IL24+ cells (Fig. 5E). Further investigation of fibroblasts trajectories using Monocle 3 (20) identified the three COMP+ subsets as the populations mapping closest to WNT5A+/IL24+ fibroblasts (Fig. 5G). This is in line with the progressively decreasing expression of WNT5A (Fig. 5B) and FOXQ1-regulated genes in COMP+3, COMP+2 and COMP+1 cells (Fig. 5F) and their absence from all other clusters. Given the role of FOXQ1 in activating Wnt signalling (21), this suggests that the drug-induced reduction in WNT5A+/IL24+ cell numbers reflects a progressive downregulation of Wnt-related genes, leading to the re-establishment of a physiological COMP+ phenotype. Accordingly, we found that IL-23 inhibitor treatment caused a significant increase in the abundance of COMP+ fibroblasts (Supplementary Fig. 5E).
Taken together, these observations identify WNT5A+/IL24+ fibroblasts as a psoriasis-specific, IL-17 dependent cell state, that is rapidly normalised following IL-23 blockade.
A reduction in WNT5A+/IL24+ cell numbers is an early event in the response to psoriasis drugs
To validate the pathogenic involvement of WNT5A+/IL24+ fibroblasts, we sought to detect these cells in other psoriasis datasets generated in the first two weeks of treatment. We retrieved three transcriptomic datasets from public repositories and used a deconvolution algorithm (CIBERSORTx (22)) to estimate the frequency of WNT5A+/IL24+ fibroblasts in patient skin. We then compared the inferred cell abundance in different sample groups.
The first dataset we examined was generated in individuals with psoriasis who were receiving guselkumab, a biologic therapy that blocks the same IL-23 subunit targeted by risankizumab (23). The analysis of these samples confirmed that WNT5A+/IL24+ fibroblasts are only detectable in lesional psoriasis skin, where their abundance begins to decrease after one week of treatment (Fig. 6A). The second dataset was obtained in a placebo-controlled clinical trial of the IL-17 inhibitor ixekizumab (24). Again, our analysis showed that the decline of WNT5A+/IL24+ fibroblasts occurred early (2 weeks) after treatment. Importantly, this reduction in cell numbers was not observed among the individuals receiving placebo (Fig. 6B). Finally, we examined RNA-seq data generated in individuals receiving a topical glucocorticoid (halometasone monohydrate 0.05% cream) for the treatment of psoriasis (25). We observed that the abundance of WNT5A+/IL24+ fibroblasts in lesional skin was significantly reduced after two weeks of effective treatment. Of note, this trend was reversed upon treatment withdrawal and clinical relapse (Fig. 6C).
(A-C) Deconvolution of transcriptome datasets shows that the frequency of WNT5A+/IL24+ fibroblasts is significantly reduced by treatment with guselkumab (A), ixekizumab (B) or topical glucocorticoid (halomethasone monohydrate 0.05%), with the latter effect going into reverse upon drug withdrawal and clinical relapse (C). (D) The percentage of WNT5A+ fibroblasts (PDGFRA+ cells) detected by RNA in-situ hybridisation decreases following treatment with risankizumab. (E) Representative confocal microscopy image showing the results of RNA in-situ hybridization (10x magnification, scale bar 100 µm). Arrows indicate cells triple-positive for PDGFRA (green), COMP (red) and WNT5A (magenta). (F) Diagram showing the early effects of therapeutic IL-23 inhibition (schematic created using BioRender.com). Left: in lesional skin, IL-23 drives the differentiation of T17 cells. These produce IL-17, leading to the activation of keratinocytes and the evolution of COMP+ fibroblasts into a WNT5A+/IL24+ state. The ensuing IL-24 production further reinforces keratinocyte activation. Right: following IL-23 blockade, T17 cell differentiation and IL-17 production are inhibited, while WNT5A+/IL24+ fibroblasts revert to their COMP+ state. As a result of diminished IL-17 and IL-24 production, keratinocyte homeostasis is re-established. Ixe, ixekizumab; L, lesional skin; NL, non-lesional skin; ns, not significant; *P<0.05; **P<0.01; ***P<0.0001 (repeated measures ANOVA, with Dunnett’s post-test (A, C, D) or unpaired t-test (B)).
To complement these in-silico observations and add spatial context to our findings, we carried out RNA in-situ hybridization in the skin of three newly ascertained patients, who had been successfully treated with risankizumab. We sampled lesional psoriasis skin at baseline (day 0, when non-lesional biopsies were also obtained) and at day 14 of treatment. The analysis of these samples confirmed that our population of interest was present in the upper dermis of lesional psoriasis skin, where its abundance was markedly reduced by IL-23 blockade (Fig. 6, D and E).
Taken together, these observations identify the reduction of WNT5A+/IL24+ fibroblasts as an early event mediating the resolution of skin inflammation in psoriasis, following systemic or topical treatment (Fig. 6F).
DISCUSSION
We present a comprehensive single-cell atlas of resolving psoriasis, which captures the dynamic cell states mediating a return to skin homeostasis. We specifically investigated the effects of risankizumab, as this highly effective drug targets the central IL-23/T17 pathogenic axis (12).
We readily detected a marked downregulation of IL-17 signalling in keratinocytes, fibroblasts and vascular endothelial cells sampled after 14 days of risankizumab treatment. At this time point, the impact of IL-23 blockade on IFN-γ activity was modest, with inhibition scores only trending towards statistical significance in keratinocytes. These results complement those of a recent study where the systemic effects of IL-23 inhibition were examined by phenotyping individuals with IL23R deficiency. Interestingly, the analysis of this rare condition uncovered prominent defects of IFN-γ responses and only a limited impairment of IL-17 activity (26). Thus, the consequences of IL-23 inhibition may be organ specific. This warrants the study of IL-23 antagonists in other tissues targeted by these drugs (e.g., the intestinal mucosa (27) and the synovium (28)).
A key result from our study is the identification of a WNT5A+/IL24+ inflammatory fibroblast state, which is expanded in lesional psoriasis skin and rapidly shrinks after treatment initiation. Importantly, our interrogation of external datasets indicates that the decline in WNT5A+/IL24+ cell numbers is a drug-agnostic phenomenon, which is also observed after IL-17 inhibitor or topical corticosteroid therapy. Inflammatory fibroblasts are typified by their ability to produce and respond to cytokines. As single-cell transcriptomics have illuminated the functional diversity of stromal cells, these inflammatory subsets have been implicated in the pathogenesis of several immune mediated conditions (3, 29, 30). For example, an integrated analysis of four immune-mediated disorders (rheumatoid arthritis, inflammatory bowel disease, interstitial lung disease, and Sjogren’s syndrome) identified two shared pro-inflammatory states (CXCL10+/CCL19+ and SPARC+/COL3A1+ fibroblasts) that were also detectable in atopic dermatitis (31). While we observed both cell states in our samples (corresponding to APOE+/CCL19+ and COL11A1+ fibroblasts), neither showed increased abundance in lesional skin or declined in frequency following IL-23 inhibitor treatment (Supplementary Fig. 5B). A WNT5A+/POSTN+ fibroblast subset has recently been identified in the inflammatory skin disease prurigo nodularis (32). In affected individuals, these cells were shifted towards a cancer-associated phenotype and showed a prominent activation of type I interferon responses. Neither of these features were observed in WNT5A+/IL24+ cells, which therefore represent a distinct, disease-specific population.
Both WNT5A and IL24 are overexpressed in psoriasis skin lesions (33, 34). Furthermore, the upregulation of IL-24 has specifically been detected in human dermal fibroblasts exposed to IL-17 cytokines (34), which is in keeping with our observation of enhanced IL-17 signalling in WNT5A+/IL24+ cells. Notably, prior research has demonstrated that both Wnt5a and IL-24 promote keratinocyte proliferation and cytokine release (33, 34). Accordingly, our ligand-receptor analyses showed that WNT5A+/IL24+ fibroblasts communicate with spinous keratinocytes in lesional skin. We also found that IL-23 inhibition abolishes this interaction, causing the downregulation of genes mediating keratinocyte proliferation and activation. Importantly, our in-situ hybridization experiments mapped WNT5A/IL24+ fibroblasts to the upper dermis, providing further support for their crosstalk with keratinocytes.
Our work also sheds new light on the relationships between fibroblast populations, suggesting that the WNT5A+/IL24+ state originates from COMP+ cells. As dermal fibroblast plasticity has been implicated in wound healing (35), fibrosis (36) and tumour responses (37), these results add to emerging evidence that the evolution of dynamic fibroblast states is also important for the onset and resolution of chronic inflammation.
Our study has some limitations. We were unable to obtain viable neutrophils and only captured a relatively small number of T cells (∼6,000 in total), which prevented us from resolving rarer subsets (e.g., MAIT cells) and cell states. Similar issues have been reported in other scRNA-seq studies (38) and likely reflect cell damage caused by the skin dissociation protocol. Given our stringent inclusion criteria, our cohort only included individuals of European descent without comorbidities. Thus, further investigations will be required to establish whether our findings can be generalised to other demographics.
Since the aim of this study was to investigate the early effects of IL-23 inhibition, our patient sample collection deliberately targeted the first two weeks of treatment. Nonetheless, our analysis of publicly available transcriptome data indicates that the abundance of WNT5A+/IL24+ cells rapidly increases upon treatment withdrawal and clinical relapse. If this finding is confirmed in future cohorts, the frequency of WNT5A+/IL24+ cells could be monitored in the increasing number of patients receiving biologics, with the potential to inform personalized, drug tapering strategies in psoriasis remission.
METHODS
Participants
Written informed consent was obtained from all participants, in line with approval from the London - Westminster Research Ethics Committee (REC ref 11/H0802/7). All participants were biologic-naïve, adult males of European descent, with a dermatologist confirmed diagnosis of severe psoriasis (PASI > 10) (Supplementary Table 1). To minimise confounders, we also applied these additional inclusion criteria: no comorbidities, no systemic immune-modifying therapy for at least 4 weeks prior to risankizumab treatment, no topical steroids to the sampled area for at least 5 days prior to sampling. Six to eight mm full thickness punch biopsies were obtained from participants at baseline (day 0; lesional and non-lesional skin) and after 3 and 14 days of treatment (lesional skin only). All biopsies were taken from site-matched areas.
Whole skin dissociation, scRNA-seq library preparation and scRNA-sequencing
Skin biopsies were processed immediately after collection and incubated in 5U/ml dispase (Stem Cell Technologies) at 4°C overnight, to separate the epidermis from the dermis. The dermis was minced and then processed using the Whole Skin Dissociation Kit (Miltenyi Biotec) and a gentleMACS™ Tissue Dissociator (Miltenyi Biotec). The epidermis was minced and digested with a 50:50 trypsin-EDTA: Versene solution (Gibco) at 37°C for 15 minutes. Dissociated cells were then filtered through a 100μm cell strainer. Following incubation with Fc block (Biolegend), TotalSeq™ Hashtag (Biolegend) antibodies were added to delineate lesional dermis, lesional epidermis, non-lesional dermis and non-lesional epidermis cell suspensions. Following sorting of viable cells on a FACSAria™ flow cytometer (BD), dermal and epidermal suspensions were mixed at a 1:1 ratio. Libraries were prepared with the Chromium Single Cell 3’ kit (10X Genomics) and sequenced using a NextSeq 2000 instrument (Illumina) to generate 150 base pair paired end reads.
scRNA-seq data processing
FASTQ files were aligned against the human GRCh38 reference genome using Cell Ranger v4.0.0 (10X Genomics). Cells expressing <300 or >4000 unique genes, or with greater than 20% of genes of mitochondrial origin, were removed. The scDblFinder package (v3.1.6) (39) was used to examine each sample and cells identified as doublets were also excluded. To avoid batch effects, data generated in different NextSeq runs was integrated using the reciprocal principal component analysis workflow provided by the Seurat package (v4.1.0) (15). The data were normalized and then scaled by regressing out the following variables: number of unique molecular identifiers per cell, percentage of mitochondrial genes, day of single-cell dissociation, patient identifier, and difference between G2M and S phase score.
Cell clustering and cell type annotation
The batch-corrected coordinate space was used for linear dimensional reduction with Seurat. A K-Nearest Neighbour graph was constructed with the FindNeighbours function and unsupervised cell clustering was implemented using the FindClusters function at a resolution of 0.4). Subclustering was performed using the same approach, except the resolution was optimised separately for each population, using Clustree (40).
Genes that were differentially expressed (fold change >1.5; Bonferroni adjusted p-value <0.05) in a cluster compared to all the others were identified using the FindAllMarkers function. Cell identities were then annotated by cross referencing these cluster-specific genes with published signature genes for the relevant cell type. The same approach was applied to subclusters, except small subclusters expressing multiple lineage markers were considered as doublets and removed.
Differential expression analysis
Genes that were differentially expressed between day 0 and day 3 or day 0 and day 14 were identified with the FindMarkers function. Genes detected in <5% of cells from either sample group were removed. A negative binomial regression was applied, including patient id as a latent variable.
Identification of enriched pathways, upstream regulators and gene regulatory networks
The Ingenuity Pathway Analysis package (Qiagen) was used to identify canonical pathways and upstream regulators that were enriched (FDR<0.05) among the DEG identified with the FindAllMarkers and FindMarkers functions. DEG networks centred around IL-17, IL-17R and IL-22 were visualised with the igraph (v1.2.6) R Package (41).
The gene regulatory networks that were active in lesional fibroblast subclusters at day 0 were identified using the workflow implemented by SCENIC (v1.3.0) (19). Briefly, co-expression modules between transcription factor and target genes (regulons) were derived using GRNBoost (a scalable alternative to GENIE3), then RcisTarget was used to refine the regulons by inferring direct targets of the transcription factors. Finally, regulon activity scores were calculated for each cell, using AUCell.
Cell-cell communication analysis
Receptor-ligand interaction analysis was performed using CellChat (v1.4.0) (16). The main analysis was based on major cell type annotations. A separate run was performed for each sample group and the number of significant interactions was calculated for each cell type pair. For the analysis of IL-24 signalling, keratinocytes, fibroblasts, T cells, and myeloid cells were divided into subclusters. Ligand-receptor interactions originating from WNT5A+/IL24+ fibroblasts were then queried at day 0 and day 14. The interaction between IL-24 and its receptor was also confirmed by using NicheNet (v1.1.1) (18) to analyse day 0 and day 14 lesional skin data. WNT5A+/IL24+ cells were set as source and the various keratinocyte subclusters were in turn set as recipient cells. Results were generated using the ‘nichenet_seruatobj_aggregate’ function.
Pseudotime analysis
Seurat fibroblast data from day 0 lesional skin were converted and imported into Monocle 3 (v1.3.1) (20) using the cell_data_set function from SeuratWrappers. Cells were separated into clusters by Leiden community detection (cluster_cells function). A cell trajectory was then fitted with the learn_graph function. Root nodes were selected within the WNT5A+/IL24A+ population and the order_cells function was applied.
Deconvolution of transcriptome datasets
The scRNA-seq dataset produced in this study was down-sampled to 5000 cells and used to generate a reference signature matrix for CIBERSORTx (22). The transcriptomic datasets generated by Sofen et al., (23) Krueger et al., (24) and Cai et al., (25) were retrieved from the Gene Expression Omnibus repository (https://www.ncbi.nlm.nih.gov/gds/), using their respective identifiers (GSE51440, GSE31652 and GSE114729). The reference matrix was then used to infer the relative abundance of WNT5A+/IL24+ fibroblasts.
RNA fluorescence in situ hybridization
16μm sections were cut from skin biopsies frozen in optimal cutting temperature compound (VWR International). RNA fluorescence in situ hybridization of Hs-PDGFRA, Hs-COMP-C2 and Hs-WNT5A-C4 probes (ACD Biologics) was performed using RNAscope® Multiplex Fluorescent Reagent Kit v2 (ACD Biologics) and 4-plex Ancillary kit (ACD Biologics). The manufacturer’s instructions were followed with minor modifications. Fixation in 4% paraformaldehyde/PBS (Santa Cruz Biotechnology) was extended to 60 minutes at 4°C, and TSA VIVID™ Fluorophore 520, 570 and 650 (ACD Biologics) were diluted at 1:1000. At least two representative images were captured for each section, using an Eclipse Ti-2 inverted microscope (Nikon).
Statistical analysis
All statistical tests were implemented in GraphPad Prism (version 9.3.0 for Windows). For the analysis of scRNA-seq and RNA in situ hybridisation results, the relative abundance of cell types in the four sample groups (non-lesional, lesional day 0, lesional day 3, lesional day 14) was compared using repeated measures analysis of variance (ANOVA), with a Dunnett’s post-test. For the analysis of cell fractions derived by deconvolution, the relative abundance of WNT5A+/IL24+ fibroblasts was compared across treatment groups, using a repeated measures ANOVA with a Dunnet’s post-test or a paired t-test, as appropriate. Details of all statistical tests are reported in the relevant figure legends.
Data availability
The scRNA-seq data associated with this manuscript have been deposited in the NCBI Gene Expression Omnibus (series id: GSE228421). The rest of the data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials.
Author contributions
F.C. and S.K.M conceived the study, designed the experiments and obtained the necessary funding. Further funding was provided by S.V. S.K.M supervised the patient recruitment, which was also facilitated by C.H.S and J.N.B. Y.K. prepared the single-cell libraries under P.D.’s supervision. L.F. carried out the remaining experimental work with input from S.K.M., C.G., J.G. and X.D.H. L.F. and D.Mc. carried out the computational analyses under F.C.’s supervision. F.C. and S.K.M. drafted the manuscript, which was then reviewed by all co-authors.
Competing interests
SV is a Boehringer-Ingelheim employee. JNB has attended advisory boards and/or spoken at sponsored symposia and/or received research funding from: AbbVie, Almirall, Amgen, Boehringer-Ingelheim, Bristol Myers Squibb, Celgene, Janssen, Leo, Lilly, Novartis, Samsung, Sun Pharma. CHS reports departmental research funding as investigator in EU-IMI consortia involving multiple industry partners (see biomap-imi.eu and hippocrates-imi.eu for details). FC has received grant support and consultancy fees from Boehringer Ingelheim. SKM reports departmental income from Abbvie, Almirall, Eli Lilly, Leo, Novartis, Sanofi, UCB, outside the submitted work.
Acknowledgments
We are grateful to all the patients who made this study possible. We also wish to acknowledge the support of King’s College London Advanced Cytometry Platform (Isabel Correa, Richard Ellis, Leanne Farnan, Anna Rose) and Genomics Core (Ulrich Kadolsky, Shichina Kannambath, Michelle Kleeman, Athul Menon, Rosamond Nuamah, Heli Vaikkinen), and Paola Di Meglio for her input on the manuscript.
This research was supported by the National Institute for Health and Care Research (NIHR) Biomedical Research Centre based at Guy’s and St Thomas’ NHS Foundation Trust and King’s College London (guysbrc-2012-1). The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care. This work was funded by the Psoriasis Association (PhD studentship ST2/21 RE18411 supporting LF and Grant BSTOP50/5 to CHS) and the Wellcome Trust (grant 096540/Z/11/Z). DMc was supported by the Medical Research Council (grant MR/R015643/1) and King’s College London as a member of the MRC Doctoral Training Partnership in Biomedical Sciences. SKM is funded by a NIHR Advanced Fellowship (NIHR302258). CHS is supported by a NIHR Senior Investigator Award.
References
Subject Area
Reviews and Context
0
Comment
0
TRIP Peer Reviews
1
Community Reviews
0
Automated Services
0
Blogs/Media
Author Videos