2D Visualization of the Psoriasis Transcriptome Fails to Support the Existence of Dual-Secreting IL-17A/IL-22 Th17 T Cells

The present paradigm of psoriasis pathogenesis revolves around the IL-23/IL-17A axis. Dual-secreting Th17 T cells presumably are the predominant sources of the psoriasis phenotype-driving cytokines, IL-17A and IL-22. We thus conducted a meta-analysis of independently acquired RNA-seq psoriasis datasets to explore the relationship between the expression of IL17A and IL22. This analysis failed to support the existence of dual secreting IL-17A/IL-22 Th17 cells as a major source of these cytokines. However, variable relationships amongst the expression of psoriasis susceptibility genes and of IL17A, IL22, and IL23A were identified. Additionally, to shed light on gene expression relationships in psoriasis, we applied a machine learning nonlinear dimensionality reduction strategy (t-SNE) to display the entire psoriasis transcriptome as a 2-dimensonal image. This analysis revealed a variety of gene clusters, relevant to psoriasis pathophysiology but failed to support a relationship between IL17A and IL22. These results support existing theories on alternative sources of IL-17A and IL-22 in psoriasis such as a Th22 cells and non-T cell populations.


INTRODUCTION
Psoriasis is a chronic inflammatory skin condition with nail and systemic manifestations that affects ∼3% of the general United States population. It is commonly associated with psoriatic arthritis and is likely linked to other comorbidities, such as cardiovascular disease and metabolic syndrome (1)(2)(3)(4).
Of the many clinical variants, plaque psoriasis (psoriasis vulgaris) is the most common, accounting for ∼80-90% of cases (1,5). It is also the most well-characterized histologically and genetically. Plaque psoriasis was initially proposed to be driven by hyperproliferative keratinocytes.
However, in 1890, neutrophil involvement was suggested after histologic evaluation revealed early neutrophil accumulation within the dermis and epidermis (i.e., microabscesses of Munro and pustules of Kogoj, respectively) (6).
Despite the clear existence of neutrophils in lesional skin, the role of the adaptive immune system in psoriasis pathophysiology became the main focus of the field after the T cell-targeting agent, cyclosporine, was shown to be an effective treatment (7)(8)(9). Thus, psoriasis researchers became very quickly focused on characterizing CD4 + and CD8 + T cell responses in normal and diseased human skin (10)(11)(12). Subsequently, experiments performed in animal models were also developed that supported the T cell-centric view of psoriasis. For example, it was demonstrated that a psoriasis-like phenotype could be induced following adoptive transfer of dysregulated CD4 + T cells (13). With this knowledge came the development of the next generation of T cell-targeting therapeutics (alefacept, efalizumab) (14)(15)(16)(17), which further corroborated the essential role of T cells in psoriasis pathophysiology. At the time T cells became the focus of psoriasis, adaptive immune responses were typically divided into two types, T helper type 1 (Th1) and T helper type 2 (Th2) responses. In psoriasis, the absence of Th2-defining cytokines [interleukin (IL)-4, IL-5, and IL-10] (18) and the increased presence of Th1 cytokines (interferon gamma (IFN-γ), tumor necrosis factor (TNF) and  prompted researchers to classify psoriasis as a Th1-mediated disease (18). Soon thereafter, however, it became increasingly apparent that IL-17-secreting T cells (Th17 cells) played a major role in disease pathogenesis, not only in psoriasis, but also across a wide spectrum of animal models of autoimmunity (19)(20)(21)(22).
Psoriasis is now thought to be a predominantly Th17-driven disease (23,24) that is maintained by the key Th-17-supporting cytokine, IL-23 (25,26). The dominant role of the IL-23/IL-17A axis in psoriasis is also evident by the overwhelming clinical success of newly developed IL-23/IL-17A axis-targeting biologics, which could induce near complete resolution of psoriasis, even in the most severely affected individuals (27)(28)(29). IL-22 is also a highly investigated cytokine involved in psoriasis pathophysiology. It is thought to be the primary promoter of keratinocyte hyperproliferation (30,31). The predominant view is that this cytokine is secreted by IL-17A/IL-22 dual-secreting Th17 cells (32).
However, the observed pathogenicity of IL-17A/IL-22 dualsecreting Th17 cells has never been formally demonstrated in vivo. In fact, the vast majority of evidence in support of these cells have come from animal studies and in vitro analysis of human T cells cultured under extreme polarizing conditions (32)(33)(34)(35). Even when studied directly ex vivo, the dual secretion is usually seen only after non-physiologic T cell stimulation (36,37). Since naturally processed autoimmune epitopes are difficult to identify (38), it is challenging to study cytokine secretion using more physiologic stimuli.
Thus, we sought evidence for the existence of dual secreting IL-17A/IL-22 Th17 cells within the psoriasis transcriptome. Weighted gene co-expression networks analysis (WGCNA) (39) have previously been used to analyze gene-gene correlations within RNA-Seq datasets. While this strategy has certain advantages, it is not ideally suited to explore gene relationships across multiple RNA-Seq datasets. Herein, we conduct metaanalyses of RNA-seq datasets to directly evaluate the current hypothesis that dual-secreting IL-17A/IL-22 Th17 cells are the dominant effector population in psoriasis. We also used this strategy to correlate the expression of IL17A, IL22, and IL23A with genes linked to psoriasis susceptibility identified through genome-wide association studies (GWAS). Finally, to explore the gene expression profile of IL17A, IL22, and IL23A in relation to other genes expressed in psoriatic plaques, we utilized a machine learning nonlinear dimensionality reduction strategy to visualize the entire psoriasis transcriptome as a 2-dimensional (2D) image. This allowed us to clearly visualize the relationship between IL17A, IL22, and IL23A and all other genes that are expressed in psoriatic skin.

Correlations
Correlation analyses of gene expressions were performed on read counts of each identified gene normalized with DESeq2 package (43). Values were subsequently log transformed and winsorized when necessary. Spearman's correlation coefficients were calculated (r s ) using the cor.test function in R (44). P values of the correlations were estimated by algorithm AS 89.

2D Visualization of the Psoriasis Transcriptome
We computed the gene pairwise distance using a formula, 1-r 2 , where r represents Pearson's correlation. A visual representation of the gene co-expression network was created using a dimensionality reduction technique, t-Distributed Stochastic Neighbor Embedding (t-SNE), calculated with Rtsne package (45).

Gene Selection
A Pubmed search was performed to identify genes linked to psoriasis through GWAS.

Meta-Analysis
Meta-analysis was completed using the R package "metafor" (46). A weighted random-effects model was used to estimate a summary effect size. To estimate between-study variance, a restricted maximum-likelihood estimator was applied. A weighted estimation with inverse-variance weights was used to fit the model.

RESULTS
The Expression of IL17A and IL22 Do Not Strongly Correlate With One Another in Psoriatic Plaques We hypothesized that if a significant amount of IL-17A and IL-22 is produced by IL-17A/IL-22 dual secreting Th17 cells in psoriasis, then the gene expression of these two cytokines should correlate with one another. In theory, their expression would be directly linked to the number of dual-secreting Th17 cells in a psoriasis plaque. Their gene expression should also correlate with IL23A, which activates and maintains Th17 cells.
To test this hypothesis, gene expression of IL17A vs. IL22 was graphed and the Spearman's correlation coefficient (r s ) was calculated ( Figure 1A). These correlative studies demonstrated that the expression of IL22 does not strongly correlate with the expression of IL17A (r s = 0.04, p = 0.67). To obtain a weighted average across all four independently acquired psoriasis datasets, a meta-analysis was performed and the resulting Forest plot ( Figure 1B) demonstrated again that IL17A and IL22 do not strongly correlate with one another [r s = 0.18, with a confidence interval that crosses 0 (−0.05, 0.41)] (Supplemental Figure 1).
Similarly, IL22 gene expression did not correlate well with IL23A (Figures 2A,B). In contrast, the expression of IL17A did correlate very well with IL23A ( Figure 3A), a result that was consistent amongst a majority of datasets. The weighted average of this correlation across all psoriasis datasets was highly significant, [r s = 0.31 (0.12, 0.51); p = 0.0014], with no FIGURE 4 | 2D gene coexpression network illustrates gene expression relationships of the psoriasis plaque. Points closer together represent genes that have higher correlation coefficients, while genes at further distances generally do not correlate well. Gene clusters appear to correspond to neutrophil-recruiting cytokines, neutrophil effector molecules, T cell, B-cell, and psoriasis-associated genes in lesional psoriatic skin. FIGURE 5 | IL22 is located at a relatively great distance away from IL17A and IL23A in the psoriasis transcriptome, as illustrated in the 2D gene expression network map of psoriasis. These results suggest that IL22 is not highly linked to either IL17A or IL23A in psoriasis pathophysiology. evidence (p = 0.33) of any substantial residual heterogeneity (i.e., there was no remaining variability in effect sizes that was unexplained) ( Figure 3B). Our data confirms the wellcharacterized dependency of IL-17A on IL-23A. However, IL-22 was not found to have a similar dependency on IL-23A, casting doubt on the theory that IL-17A and IL-22 are secreted mainly by the same dual-secreting cell.

2D Visualization of the Psoriasis Transcriptome Reveals T Cell, B Cell, Inflammatory Cytokines, Neutrophil-Recruiting, and Neutrophil Gene Clusters
To determine how genes expressed in psoriatic plaques are related to one another, correlation coefficients were calculated for all pairwise comparisons. The distances between each gene pair was calculated as described in the methods. The resulting distance matrix was then used to construct a 2D image using t-SNE.
In the 2D plot (Figure 4), genes that highly correlate with one another tend to be located in the same region, known as a cluster. Genes that do not cluster near each other do not correlate well.  With respect to other cytokines and chemokines involved in psoriasis pathophysiology, IL22 is located peripherally at a relatively great distance away on the 2D plot of the psoriasis transcriptome ( Figure 5). This supports our results from the meta-analyses and suggests that IL22 does not correlate well with IL23A. Interestingly, IL22 does not cluster well with any of the most commonly implicated cytokines in psoriasis pathophysiology. In contrast, IL17A clusters together with IL23A and the other cytokines thought to be involved in psoriasis pathophysiology.

IL22 Correlates With Keratins
Several studies have demonstrated that IL-22 stimulates keratinocytes. There is a variety of evidence, including data obtained from in vitro studies with skin-like organoid cultures, that support IL-22 as the main cytokine responsible for epidermal hyperplasia, a hallmark of psoriasis (47,48).
Because IL22 failed to strongly correlate with IL17A and IL23A (Figures 1A, 2A), the relationships between IL22 and keratin genes were explored across the four independently acquired RNA-Seq psoriasis datasets. Again, Spearman's correlation coefficients (r s ) were calculated for each keratin gene's relationship with IL17A, IL22, and IL23A. These correlative studies demonstrate that the expression of IL22 did indeed strongly correlate with the expression of the different keratin genes (Figure 6), especially KRT6C (keratin 6C) (r s = 0.32, p = 0.0011). To obtain a weighted average across all four independent psoriasis datasets, a meta-analysis was performed and the resulting Forest plots (Figure 7) confirm the close relationship between IL22 expression and keratin gene expression [KRT6C: r s = 0.34, with a confidence interval that did not cross 0 (0.18-0.50)]. The weighted average of this correlation across all psoriasis datasets was highly significant (p = 0.000025), with no evidence (p = 0.56) of any substantial residual heterogeneity (i.e., there was no remaining variability in effect sizes that was unexplained). Additional genes that were found to positively correlate with IL22 expression are listed in Supplemental Figure 2.
In contrast to IL22's relationship with keratin gene expression, IL17A did not correlate well with the keratins (Figure 6), a finding that was confirmed by a meta-analysis across all four RNA-Seq datasets.

IL23A Correlates With Other Genes Besides IL17A
IL-23A is known for its ability to support Th17 T cells but it likely has a variety of functions independent of this role. To investigate this, IL23A's ability to independently correlate with other immune-relevant genes was explored. Figure 8 reveals that IL23A correlates with several genes unrelated to IL17A, a finding confirmed by meta-analyses across all psoriasis RNA-Seq datasets. Included in the analysis were genes identified by GWAS to be linked to psoriasis. Of these genes, CARM1, KRT14, KRT37, TNFAIP3, UBE2L3 are elevated in psoriasis plaques compared to control healthy skin ( Table 1). Thus, IL23A appears to be linked to other genes putatively involved in the pathophysiology of psoriasis that are unrelated to IL17A.

IL17A, IL22, and IL23A Expression Correlates With Psoriasis Susceptibility Genes
A variety of genes have been linked to psoriasis susceptibility through GWAS (49)(50)(51)(52)(53)(54)(55). Table 2 demonstrates that many of these genes are differentially regulated in the setting of psoriasis. We thus sought to determine how the expression of genes located at psoriasis susceptibility loci correlated with the expression of IL17A, IL22, and IL23A, genes known to be linked to the pathophysiology of psoriasis. For this analysis, the expression of IL17A, IL22, and IL23A was plotted against the expression of each of the genes identified through GWAS studies. Spearman's correlation coefficients (r s ) were then calculated, which demonstrated a variety of significant correlations ( Table 2) between GWAS-identified genes and IL17A, IL22, and IL23A. Correlation values between atopic dermatitis GWAS-identified genes and IL17A, IL22 and IL23A expression in psoriasis samples were also obtained for comparison ( Table 3). To obtain a weighted average across all four independent psoriasis data sets, meta-analysis was performed. The resulting Forest plots are depicted in Figure 9, which confirm the close relationship between IL17A, IL22, and IL23A and the different genes linked to psoriasis susceptibility. These results support a direct or indirect link between IL17A, IL22, and IL23A and these genes. Of note, the genes that significantly correlated with IL17A, IL22, and IL23A varied for each cytokine. These results will hopefully help investigators better understand the pathophysiology of psoriasis.

DISCUSSION
Investigators have employed numerous genetic strategies to characterize the immune response in the setting of psoriasis. Microarray and RNA-Seq have provided insight into the psoriatic transcriptome, identifying thousands of differentially expressed genes (40). However, differential expression alone does not necessarily mean that the gene is involved in psoriasis pathogenesis. For example, a gene that is normally downregulated in psoriatic T cells may actually appear falsely upregulated in psoriasis simply because there are more T cells in a psoriatic plaque. With the rising popularity of single cell sequencing, investigators are now focused on re-characterizing the psoriasis transcriptome at a greater cellular resolution, not previously obtained with whole tissue transcriptomics.
However, single cell sequencing is also not without its drawbacks. Purifying immune cell populations from skin biopsy specimens can alter their transcriptome, especially for cells isolated by positive selection or flow cytometry. Furthermore, immune cells within the skin will undoubtedly have different purification yields. T cells in particular are especially difficult to analyze because once purified, they require additional non-physiologic in vitro stimulation with lectins or anti-CD3/anti-CD28 antibodies to identify their cytokine secretion profiles. How closely the garnered information from these studies will relate to in vivo cellular function remains unclear. Although each technique will yield important discoveries, none can perfectly decipher the in vivo pathogenic immune response.
With these limitations in mind, we have focused on developing new methods to characterize immune responses from whole tissue RNA-Seq (12,62). We view this strategy to be an important complement to work currently being conducted by other investigators. The main advantage being that the data is not subject to experimentally-induced changes in gene expression. Its main disadvantage, however, is that it cannot discriminate between direct or indirect correlations between genes of interest.
In our current study, we utilize a machine-learning 2D visualization strategy, t-SNE, to characterize IL17A, IL22, and IL23A gene expression in the context of the entire psoriatic transcriptome. The 2D map of the psoriatic transcriptome revealed distinct gene clusters corresponding to common immune cell types (e.g., B cells, T cells, neutrophils).
Our data did not support the existence of a dual-secreting IL-17A/IL-22 Th17 cell as the major source of these cytokines in psoriasis. In fact, in the 2D model, these genes are located far from one another. As such, IL22 correlated with several genes that did not appear to have a relationship with IL17A. In addition, a set of genes identified to be involved in psoriasis pathophysiology (CARD14, CXCL5, CXCR2, DDX58, IFIH, PTPN22, and TNFRSF9) correlated with IL17A and IL23A, but did not correlate with IL22.
Though, IL-22 is commonly considered a hallmark Th17 cytokine (63), our results are in line with studies demonstrating  the existence of uniquely secreting IL-17 and IL-22 T cells or the existence of other cytokine-secreting phenotypes (48,(64)(65)(66)(67)(68)(69), although these other studies usually relied upon non-physiologic ex vivo T cell stimulation. Another possibility is that other cell types, such as γδ T cells or mast cells, contribute to the IL-22 production in psoriasis (48,65,70). Even neutrophils have been implicated as major producers of IL-22 and IL-17A (71) and recent animal models have re-explored their role as effector cells in psoriasis pathophysiology (22,72,73). Indeed, there are numerous studies supporting a key function of these cells (71,(74)(75)(76)(77)(78). Single cell sequencing may provide information to verify the relationship between IL17A and IL22 expression.
Although it is possible that dual-secreting IL-17A/IL-22 Th17 cells exist, our results suggest that they are not a major source of IL-22.
Although we did not find evidence for a strong link between IL22 with IL17A or IL23A, our results do support a strong correlation between the expression of IL22 and the keratin genes, such as KRT6A, KRT6B, KRT6C, KRT16, KRT17, a finding in accord with IL-22's ability to induce epidermal hyperproliferation (48). IL-22 is clearly a major cytokine involved in psoriasis pathophysiology. In animal models, it has been demonstrated to simulate psoriasis-like epidermal changes (47,79) and elevated levels of IL-22 positively correlate with disease severity in humans, as measured by Psoriasis Area Severity Index (PASI) scores (80)(81)(82)(83).

CONCLUSION
Although dual-secreting T cells may exist, our results demonstrate that it is unlikely that the classical Th17 cells (IL-17A/ IL-22 dual-secreting T cells) play a universal role in  psoriasis pathophysiology. RNA-Seq analysis revealed that the expression of these cytokines seems to be largely unrelated to one another in the psoriasis transcriptome. However, the expression of IL17A did correlate with IL23A but, interestingly, unique relationships between IL23A and genes unrelated to IL17A were also established, supporting a broad function of IL-23. Taken together, these results do not support the current dogma that IL-17A/IL-22 dual-secreting Th17 T cells are the major driver of psoriasis pathophysiology. In addition, our results support unique functions of IL-23 that are unrelated to its known role in supporting Th17 responses. Finally, we demonstrate that the expression of genes linked to psoriasis susceptibility also correlate with expression of either IL17A, IL22, or IL23. This supports the aforementioned cytokines' involvement in multiple avenues of psoriasis susceptibility.
2D mapping of inflammatory transcriptomes is an exciting innovative modality that may help us visualize relationships of all genes expressed in a disease process. When applied to gene expression relationships in psoriatic lesional skin, distinct clusters of cell lineage genes could be identified, supporting the presence of a complex crosstalk among separate cell lines in disease development. In the near future, single cell transcriptome analysis will provide additional insight into psoriasis pathogenesis. Identifying the cells responsible for the psoriasis phenotype will bring us one step closer to developing a cure for psoriasis.

DATA AVAILABILITY
Publicly available datasets were analyzed in this study. This data can be found here: "http://www.ncbi.nlm.nih.gov/Traces/sra".

AUTHOR CONTRIBUTIONS
EM, AAM, AIM, GL, MS, SR, and SH contributed to the conception and design of the study; LT and JG organized the database; AAM performed the statistical analysis; SL, JW, and CA performed data mining, SL and EM wrote the first draft of the manuscript; SL, IA, and EM wrote sections of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu. 2019.00589/full#supplementary-material Supplemental Figure 1 | Funnel plot representation demonstrating RNA-seq datasets analyzed in the meta-analysis of IL17A and IL22. All data points, representing individual data sets, fall within the 95% confidence interval. In our meta-analysis, the p value for residual heterogeneity did not reach significance (p = 0.34), indicating that all datasets are within the variation that is expected for this particular meta-analysis.
Supplemental Figure 2 | Additional genes that positively correlate with IL22 expression.