Transcriptomic Profiling of Peripheral Edge of Lesions to Elucidate the Pathogenesis of Psoriasis Vulgaris

Elucidating transcriptome in the peripheral edge of the lesional (PE) skin could provide a better understanding of the molecules or signalings that intensify inflammation in the PE skin. Full-thickness biopsies of PE skin and uninvolved (UN) skin were obtained from psoriasis patients for RNA-seq. Several potential differentially expressed genes (DEGs) in the PE skin compared to those in the UN skin were identified. These DEGs enhanced functions such as angiogenesis, growth of epithelial tissue, chemotaxis and homing of cells, growth of connective tissues, and degranulation of myeloid cells beneath the PE skin. Moreover, the canonical pathways of IL-17A, IL-6, and IL-22 signaling were enriched by the DEGs. Finally, we proposed that inflammation in the PE skin might be driven by the IL-36/TLR9 axis or IL-6/Th17 axis and potentiated by IL-36α, IL-36γ, IL-17C, IL-8, S100A7, S100A8, S100A9, S100A15, SERPINB4, and hBD-2. Along with IL-36α, IL-17C, and IκBζ, ROCK2 could be an equally important factor in the pathogenesis of psoriasis, which may involve self-sustaining circuits between innate and adaptive immune responses via regulation of IL-36α and IL-36γ expression. Our finding provides new insight into signaling pathways in PE skin, which could lead to the discovery of new psoriasis targets.


Introduction
Psoriasis is a chronic skin disease influenced by the interplay of genetic predisposition, environmental triggers, and immune response. The classic feature of chronic plaque or psoriasis vulgaris is a well-demarcated salmon-pink plaque covered with silvery scales [1]. In addition to psoriasis vulgaris, inverse, guttate, pustular, and erythrodermic forms have been reported as the different phenotypes of psoriasis, which may be associated with comorbidities such as arthritis, metabolic syndrome, and cardiovascular diseases [1][2][3]. Various genes involved in antigen presentation (HLA-Cw6), skin barrier formation, innate immunity, and adaptive immunity have been identified as psoriasis-susceptible genes that increase the risk of developing the disease [4,5]. Individuals who carry these susceptible genes will be more likely to develop psoriasis under the influence of environmental factors, which initiate an inflammatory response leading to psoriasis phenotypes [1,[5][6][7]].

DEGs in the PE Skin
A total of 1202 DEGs were identified. Of these, 653 (54%) were upregulated, and 549 (46%) were downregulated in PE skin compared to the uninvolved (UN) skin ( Figure 1). A heatmap of the differentially regulated transcripts between the PE and UN skin is shown in Supplementary Figure S1. The top 50 upregulated and downregulated DEGs according to the log2FC value are listed in Table 1 Based on QIAGEN Knowledge Base as of October 2021, the blue molecules have been previously mentioned in psoriasis, either their mRNAs or their proteins. * Molecule types are classified into complex, cytokine, enzyme, fusion gene/product, G-protein coupled receptor, group, growth factor, ion channel, kinase, ligand-dependent nuclear receptor, mature microRNA, peptidase, phosphatase, transcription regulator, translation regulator, transmembrane receptor and transporter. None of these were classified as other. FC, fold change.  The 1202 DEGs were analyzed for downstream effects in terms of associated diseases; among them, 135 genes have previously been reported to be associated with psoriasis based on the Ingenuity Knowledge Base, 86 of which showed upregulated in the PE skin. Thus, we hypothesized that the proteins encoded by these 86 DEGs might be involved in inflammation in the PE skin to contribute to psoriasis (Figure 1 and Supplementary Table S3). Furthermore, the 1202 DEGs were analyzed for their functions. These DEGs enhanced functions are to enrich angiogenesis (EDNRB, LAMA2, S100A8, S100A9, SPINK5, TNFRSF1A, and VAV3), growth of epithelial tissue (ADAM17, KRT16, IKBKB, IL36G, ITGB1, RICTOR, and XDH), chemotaxis and homing of cells (A2M, C10orf99, DEFB4A/DEFB4B, GFRA1, PLEC, S100A7, and S100A7A), growth of connective tissue (ADAR, ELN, IL36A, KLK6, PLAT, STAT1, and STAT3) and degranulation of myeloid cells (CD59, CXCL8, DSC1, DSG1, LCN2, and SERPINB3). The list of all DEGs in each function is provided in Supplementary  Table S4. In addition, the upregulated-and downregulated DEGs were analyzed separately as pathway and function enrichments using Metascape. Notably, cornified envelope formation, VEGFA-VEGFR2 signaling, interferon signaling, and neutrophil degranulation were enriched by upregulated DEGs (Figure 2A), whereas epithelial differentiation was enriched by downregulated DEGs ( Figure 2B) such as KRT10, KRT2, LCE1C, LCE1D, LCE2B, and LORICRIN. The list of all DEGs in each DEG-enriched pathway and function is provided in Supplementary Tables S5 and S6. We reviewed and selected the potential DEGs based on their biofunctions and psoriasis associations; they are shown in Table 2 and Figure 3. Moreover, we compared and mapped PE vs. UN skin DEGs from previous studies [11][12][13][14][22][23][24]29,30]; the details are provided in Supplementary Table S7.      These DEGs were selected due to their biofunction and psoriasis association. Based on QIAGEN Knowledge Base as of October 2021, the bold molecules have not been previously mentioned in psoriasis. FC, fold change.

Canonical Pathways
We analyzed the known canonical pathways enriched by the DEGs of PE skin and the activity prediction of these pathways based on the QIAGEN Knowledge Base. To include both significant and highly predicted activity pathways, we determined the threshold of minimum significance (p-value of overlap > 0.05 or −log (p-value) > 1.3) and absolute z-score (|z-score| ≥ 2). Only seven canonical pathways are shown in Figure 4A,B. It is noteworthy that the role of IL-17A in psoriasis, eukaryotic initiation factor 2 (eIF2) signaling, and the coronavirus pathogenesis pathway had the top three highest ratios (mapped DEGs divided by total molecules in a given pathway = 0.36, 0.30, and 0.25, respectively). Moreover, these three pathways are involved in the stress response and inflammation, especially the coronavirus pathogenesis pathway, the hallmark of which is the cytokine storm [31]. The ratios of all seven pathways are listed in Supplementary Table S8. storm [31]. The ratios of all seven pathways are listed in Supplementary Table S8.

Upstream Regulators
We further identified the upstream regulators of the DEGs in the PE skin. We used the upstream analysis function according to a statistically significant overlap between the dataset genes and the genes that are regulated by a potential regulator (p-value overlap < 0.05) and |z-score| ≥ 2. The lists of potential regulators with their predicted activity are provided in Supplementary Table S9. Among the potential regulators, seven molecules were identified as DEGs, including CCND1, DDX21, IL36A, KDM5A, NUPR1, RICTOR, To explore the role of cytokine signaling in the PE skin in further detail, we analyzed the canonical pathway of cytokines, with a −log (p-value) > 1.3 indicating significance. A total of 11 cytokine signaling pathways were among the DEG-enriched pathways ( Figure 4C,D), including the IL-8, IL-17A, IL-6, C-C motif chemokine receptor 3(CCR3), IL-15, granulocytemacrophage colony-stimulating factor (GM-CSF), IL-22, and formyl-methionyl-leucylphenylalanine (fMLP) signaling pathways.

Upstream Regulators
We further identified the upstream regulators of the DEGs in the PE skin. We used the upstream analysis function according to a statistically significant overlap between the dataset genes and the genes that are regulated by a potential regulator (p-value overlap < 0.05) and |z-score| ≥ 2. The lists of potential regulators with their predicted activity are provided in Supplementary Table S9. Among the potential regulators, seven molecules were identified as DEGs, including CCND1, DDX21, IL36A, KDM5A, NUPR1, RICTOR, and STAT1. Some relationships between the potential regulators and their target DEGs are shown in Figure 5 and Supplementary Figure S2.

Regulator Effect
Next, we established a causal hypothesis between upstream regulators and their downstream effects using the regulator effect function to correlate the significant upstream regulators with their downstream functions. The useful predicted networks were prioritized and ranked according to their consistency scores. The two highest consistency score networks are shown in Figure 6. Notably, IL-36α, IL-17C, nuclear factor-kappa B inhibitor zeta (IκBζ), and Rho-associated kinase 2 (ROCK2) emerged as potentially important regulators in the network. and STAT1. Some relationships between the potential regulators and their target DEGs are shown in Figure 5 and Supplementary Figure S2.

Regulator Effect
Next, we established a causal hypothesis between upstream regulators and their downstream effects using the regulator effect function to correlate the significant upstream regulators with their downstream functions. The useful predicted networks were prioritized and ranked according to their consistency scores. The two highest consistency score networks are shown in Figure 6. Notably, IL-36α, IL-17C, nuclear factor-kappa B inhibitor zeta (IκBζ), and Rho-associated kinase 2 (ROCK2) emerged as potentially important regulators in the network.

Validation of Selected DEGs by Multiplex Real-Time PCR
Finally, to validate the expression of DEGs such as SERPINB4, IL36A, PLSCR1, CXCL8, DEFB4A/B, DEFB103A/B, S100A7, S100A7A, S100A8, and S100A9, we performed quantitative multiplex Real-Time PCR. Furthermore, we compared the mRNA expression of upstream regulators IL-17C, IκBζ, and ROCK2 in the PE skin to those in the UN skin. The results are presented in Figure 7.

Discussion
To the best of our knowledge, this is the first study to establish the transcriptomic profile at the edge of psoriatic lesions. We revealed that the increased functions in the PE skin were involved in angiogenesis, growth of epithelial tissue, chemotaxis and homing of cells, growth of connective tissue, and the degranulation of myeloid cells. The upregulated DEGs were enriched in functions such as cornified envelope formation, VEGFA- The fold change of significantly altered genes. All graphs were designated by geometric mean and standard deviation. The significance was determined using *; * p < 0.05, ** p < 0.01 and *** p < 0.001. n ≥ 10 for each gene. CXCL8, C-X-C motif chemokine ligand 8; DEFB103A/B, defensin beta 103A/B; DEGs, differentially expressed genes; IL17C, interleukin 17C; NFKBIZ, nuclear factor kappa B inhibitor zeta; PE skin, peripheral edge of lesional skin; PLSCR1, phospholipid scramblase 1; ROCK2, rho-associated kinase 2; S100A7, S100 calcium binding protein A7; SERPINB4, serpin family B member 4; UN skin, uninvolved skin.

Discussion
To the best of our knowledge, this is the first study to establish the transcriptomic profile at the edge of psoriatic lesions. We revealed that the increased functions in the PE skin were involved in angiogenesis, growth of epithelial tissue, chemotaxis and homing of cells, growth of connective tissue, and the degranulation of myeloid cells. The upregulated DEGs were enriched in functions such as cornified envelope formation, VEGFA-VEGFR2 signaling, interferon signaling, and neutrophil degranulation. Moreover, we showed that the genes involved in epithelial differentiation were downregulated in the PE skin compared with the UN skin, which might limit the normal differentiation capacity of KC and contribute to plaque formation. We discussed some potential DEGs that were related to AMPs, angiogenesis, autoantigens, and chemotaxis that might play a pivotal role in the pathogenesis of PE skin.
AMPs activate the innate immune response to create an inflammatory environment in the skin through antimicrobial, chemotactic, angiogenic, and proinflammatory properties [32,33]. We identified various AMP-coding genes that were upregulated in the PE skin. We confirmed that DEFB4A/B mRNA was significantly upregulated in PE skin and its fold change was the highest among the AMP-coding genes. It was previously reported to be expressed in psoriatic plaques [34,35] and more highly in lesional skin than in UN skin [36]. Its protein, human β-defensin 2 (hBD-2), promotes pDCs to uptake self-DNA and secrete IFN-α in a TLR-9-dependent manner [37]. Moreover, we confirmed that S100A7, S100A7A, S100A8, and S100A9 mRNA were significantly upregulated in PE skin compared to those in UN skin. Interestingly, among the psoriasis-associated S100 protein-coding genes, S100A7A (also known as S100A15) was the most upregulated gene in the PE skin. Its protein koebnerisin is expressed in psoriatic lesions along with psoriasin (encoded by S100A7), S100A8, and S100A9 [38,39] by induction of psoriatic-derived Th17/Th22 cytokines (IL-17, IL-22, and TNF-α) [40][41][42]. Of particular interest, apart from their antimicrobial properties of them [33,43], both koebnerisin and psoriasin have been proposed to be responsible for the inflammation priming of KCs [41], as well as S100A8 and S100A9 also function as proangiogenic factors [39]. However, further study is needed to determine their role in the PE skin, in particular the very-high-fold change of S100A7A.
We also identified the upregulation of SERPINB3 and SERPINB4 in the PE skin and confirmed that SERPINB4 was significantly upregulated in the PE skin. Their proteins, serpin family B members 3 and 4, were reported to be expressed in psoriatic lesions [44] by induction of IL-17 and IL-22 [45]. Pso p27 antigen, a proteolyzed product of SERPINB4 and SERPINB3 [46,47], was detected in the PE skin, which is considered to be equivalent to its detection on the lesion, but was not detected on the UN skin [48]. The pso p27 forms an immune complex with anti-pso p27 to activate the complement cascade at the psoriatic scale [49]. However, the roles of SERPINB3 and SERPINB4 as antigens in psoriasis are controversial and require further investigation, especially SERPINB4, which was highly upregulated in the PE skin.
We further used the canonical pathway, regulator analysis, and regulator effect functions of QIAGEN's IPA tool to generate the network among the molecules. The proposed pathogenic network is schematically displayed in Figure 8. We found that IL36A and IL36G expression were upregulated in the PE skin and then confirmed that IL36A expression was significantly upregulated in the PE skin. Their proteins IL-36α and IL-36γ are expressed in psoriatic plaque [50][51][52], and their mRNA levels were reported to be increased in KCs after injury [53,54]. Moreover, IL-36α and IL-36γ alone or synergistically with IL-17A could amplify their own mRNA and proteins in KCs [55]. We identified and confirmed the significant upregulation of PLSCR1 (encoding a TLR-9 transporter) in the PE skin. A previous RNA-seq analysis of whole blood from psoriasis patients implicated the pDC-mediated release of IFN-α after induction by IL-36 in a TLR-9-dependent manner, evidenced by the upregulation of PLSCR1 in pDCs [56]. The pDC-derived IFN-α promotes mDC maturation and function, which amplifies and sustains T cells for the psoriatic inflammatory response [7,8]. Therefore, we suggest that the IL-36/TLR-9 axis, which upregulates IFN-α production, might also drive psoriatic inflammation in the PE skin (Figure 8). . IL-36α and IL-36γ activate myeloid dendritic cells (mDCs) to mature and secrete IL-6, which drives the differentiation of T helper (h) cells (IL-6/Th17 axis). ROCK2 (rhoassociated kinase 2) regulates CD4 + T cells to secrete IL-17. Note that inflammation in the PE skin is potentiated by IL-36α, IL-36γ, IL-17C, IL-8, S100A7, S100A8, S100A9, S100A15, SERPINB4, and hBD-2 and that a self-sustaining circuit between the innate and adaptive immune response in the PE skin is mediated through IL-36α and IL-36γ.
IL-36α and IL-36γ also activate mDCs to mature and secrete IL-6 [57]. Dermal DCderived IL-6 then drives Th17 and Th22 cell differentiation [7,58]. Furthermore, IL-6 is also secreted by KCs after IFN-γ treatment via the suppression of miR-149, which amplifies skin inflammation [59]. Consistently, we found that the DEGs in the PE skin were enriched in the IL-6 canonical pathway, suggesting that the IL-6/Th17 axis might be another axis driving psoriatic inflammation in the PE skin ( Figure 8). However, despite the lack of IL-6, KCs expressed various additional proinflammatory cytokines after the transient suppression of psoriasiform skin manifestation, delineating the ineffective treatment of psoriasis plaque by IL-6 blockage [60]. Our results showed that the DEGs were enriched in the IL-17A and IL-22 canonical pathways and that IL-22 was an activated potential regulator, thus suggesting the crucial roles of IL-17A and IL-22 in the PE skin. IL-17A and IL-22 potentiate the psoriatic milieu by inducing AMPs, recruiting other immune cells, and increasing KC proliferation [7,8]. In addition, a synergistic effect of the Th17 cytokines IL-17A and IL-22 induces the expression of IL-36α and IL-36γ mRNA in KCs [50,61]. IL-22, in particular, promotes KC proliferation and suppresses KC terminal differentiation [42]. Previous studies comparing LS skin to healthy skin suggested that IL-22 functions better when the levels of soluble scavenging receptor IL-22 binding protein (encoded by IL22RA2) are decreased. However, we could not identify the alteration of IL22RA2 as DEG of PE vs. UN skin. Furthermore, serum IL-22 level was higher in psoriasis patients compared with healthy individuals and was also positively correlated with disease severity [42,62].
IL-36α and IL-36γ also activate mDCs to mature and secrete IL-6 [57]. Dermal DCderived IL-6 then drives Th17 and Th22 cell differentiation [7,58]. Furthermore, IL-6 is also secreted by KCs after IFN-γ treatment via the suppression of miR-149, which amplifies skin inflammation [59]. Consistently, we found that the DEGs in the PE skin were enriched in the IL-6 canonical pathway, suggesting that the IL-6/Th17 axis might be another axis driving psoriatic inflammation in the PE skin ( Figure 8). However, despite the lack of IL-6, KCs expressed various additional proinflammatory cytokines after the transient suppression of psoriasiform skin manifestation, delineating the ineffective treatment of psoriasis plaque by IL-6 blockage [60]. Our results showed that the DEGs were enriched in the IL-17A and IL-22 canonical pathways and that IL-22 was an activated potential regulator, thus suggesting the crucial roles of IL-17A and IL-22 in the PE skin. IL-17A and IL-22 potentiate the psoriatic milieu by inducing AMPs, recruiting other immune cells, and increasing KC proliferation [7,8]. In addition, a synergistic effect of the Th17 cytokines IL-17A and IL-22 induces the expression of IL-36α and IL-36γ mRNA in KCs [50,61]. IL-22, in particular, promotes KC proliferation and suppresses KC terminal differentiation [42]. Previous studies comparing LS skin to healthy skin suggested that IL-22 functions better when the levels of soluble scavenging receptor IL-22 binding protein (encoded by IL22RA2) are decreased. However, we could not identify the alteration of IL22RA2 as DEG of PE vs. UN skin. Furthermore, serum IL-22 level was higher in psoriasis patients compared with healthy individuals and was also positively correlated with disease severity [42,62].
Considering IL-36, IL-36α was identified as an activated potential regulator in the PE skin, which induces the expression of IL36G, IL17C, S100A9, and NFKBIZ in KCs [63]. Notably, the mouse skin showed enhanced chemokine expression, leukocyte infiltration, and acanthosis after the injection of IL-36α [57]. Moreover, imiquimod induced much milder psoriasiform dermatitis in Il36a −/− mice compared with that induced in wild-type mice [64].
IL-17C was identified as an activated potential regulator. We established that its mRNA expression was not significantly upregulated in the PE skin because similar relative expression levels were also detected in the UN skin. These results suggest that IL-17C plays a role during the early stages of psoriasis inflammation or inflammatory priming for plaque formation. A previous report demonstrated that IL-17C transgenic mice developed psoriasiform skin. Moreover, IL-17C and its mRNA significantly increase in UN and LS skin compared to healthy skin [65]. IL-17C is mainly secreted from KCs and not leukocytes [66]. It has also been detected in psoriatic plaques at a higher concentration than IL-17A [65,67]. IL-17C stimulates the release of TNF-α from monocytes [68] and also simulates neutrophil migration, especially in synergy with TNF-α [67]. IL-17C could amplify itself in KCs, especially in synergy with TNF-α [65,67]. In addition, IL-36α and IL-36γ alone or synergistically with IL-17A increased IL-17C levels from KCs [55]. Moreover, IL-17C induced the expression of S100A15, S100A7, S100A8, S100A9, DEFB4A, CXCL8, IL36G, NFKBIZ, and TNIP3 in KCs [66,67]. In particular, IL-17C synergistically with TNF-α induced the expression of S100A7, S100A8, S100A9, LCN2 (coding lipocalin-2), CXCL8, IL36G, and DEFB4 in KCs [65].
We identified significant upregulation of CXCL8 in the PE skin. It encodes IL-8, which attracts T lymphocytes to the psoriasis plaque [73], promotes angiogenesis [74,75], and promotes epidermal proliferation [76]. Interestingly, IL-8 is secreted by accumulated neutrophils (Munro's microabscess) and is the chemokine for other neutrophils to aggregate into the plaque in an autocrine manner [77]. Moreover, IL-8 is secreted from KCs after scratching [54] or the induction of IL-22 or induction of IL-36α and IL-36γ, especially in synergy with IL-17A [55,57,61,78]. Particularly, KC-derived IL-8 is capable of neutrophil chemotaxis, but this capacity was reduced when miR-146a, an IL-17 counteractant, which is expressed in psoriasis plaque, was overexpressed [79].
We identified that ROCK2, as an activated potential regulator, could be equally important as IL-36α, IL-17C, and IκBζ in the network. ROCK2 is located in the nucleus and cytoplasm of the human T-cell [80]. It plays an important role in keratinocyte differentiation [81] as well as in angiogenesis [82]. Moreover, ROCK2 was found to regulate CD4 + T cells to secrete IL-17 [83]. These processes are important for the formation of psoriasis plaques [7]. In 2017, the activity of ROCK2 was studied in moderately severe plaque psoriasis patients using selective oral ROCK2 inhibitors; 71% of the patients achieved a reduction of 50% in the psoriasis area and severity index (PASI 50) at 12 weeks, peripheral blood IL-17 and IL-23 significantly decreased, and ROCK2 staining decreased with decreasing epidermal thickness [84]. As the predicted regulator in the PE skin, ROCK2 was associated with the upregulation of nine DEGs (CXCL8, DSC1, DSC2, DSG1, DSG3, SPINK5, PI3, S100A7, and S100A8) ( Figure 5). We further investigated the ROCK2 mRNA expression in PE skin and compared it to those in the UN skin. Our results showed that ROCK2 mRNA was significantly downregulated. However, activation of the ROCK2 kinase protein depends on post-translational modification steps. Previous reports suggest that activation of ROCK2 kinase needs phosphorylation [85,86]. Interestingly, its structure is capable of self-phosphorylation (auto-phosphorylation) [87]. Thus, we suggest that the downregulation of ROCK2 mRNA in PE skin might be the result of autoregulation overfunction. In addition, ROCK2 kinase has been proposed to cooperate with phosphorylated signal transducer and activator of transcription 3 (STAT3) to regulate Th17 function [80]. Due to the complexity of psoriasis, further investigation of ROCK2 mRNA expression in the PE skin compared to the mature lesion or healthy skin could provide new information toward a better understanding of the mechanism of oral ROCK2 inhibitor action and the role of ROCK2 in the pathogenesis of psoriasis.
Finally, based on our results and these previous findings, we propose that IL-36α, IL-36γ, IκBζ, IL-17C, IL-8, S100A7, S100A8, S100A9, S100A15, SERPINB4, hBD-2, and ROCK2 potentiate inflammation in the PE skin and that IL-36α and IL-36γ may run a self-sustaining circuit between the innate and adaptive immune response, creating a psoriatic milieu in the PE skin ( Figure 8).
The limitation of this study is the small sample size. Ideally, more samples of a larger size need to be evaluated to determine the PE skin transcriptomic profile.
In conclusion, this study identified a set of upregulated DEGs, including PLSCR1, IL36A, CXCL8, S100A7, S100A8, S100A9, S100A15, DEFB4A/B, and SERPINB4, and the molecules that may regulate these DEGs including IL-17C, IL-36α, IκBζ, and ROCK2. These DEGs and upstream regulators may potentiate the inflammation in PE skin as the proposed network might be different from that seen in the LS skin. This study provides new insight into the molecular mechanisms occurring in the PE skin and could lead to the identification of specific molecules that could serve as therapeutic targets in psoriasis vulgaris. Further, it would be helpful to compare these molecules with those observed in the center of LS skin, normal-appearing skin in patients with psoriasis, and healthy skin, to further confirm changes in the levels of other DEGs and regulators of interest at both the mRNA and protein levels, and to further examine therapeutic modalities that might target these molecules by comparing the before and after treatment transcript levels in the PE skin. These insights could provide a better understanding of the pathogenesis of psoriasis and lead to the discovery of novel therapeutic strategies.

Patients
We enrolled voluntary patients with plaque-type psoriasis (aged ≥18 years) with moderate-to-severe psoriasis [Psoriasis Area and Severity Index (PASI) score ≥ 10] diagnosed by either classical clinical features or a confirmed pathological examination. Patients who had been treated with topical agents (steroids, vitamin D analogs) within 2 weeks or had undergone phototherapy and systemic medication (prednisolone, methotrexate, biologics) within 4 weeks of the date of tissue sampling were excluded. Based on previous RNA-seq profiling studies in human psoriasis skin, three samples were deemed to be adequate for analysis [11,88]; therefore, three eligible patients were included in this study. The PASI scores of the three patients are shown in Supplementary Table S10. All three patients were informed of all study protocols and signed the consent forms. The protocol was designed in accordance with the principles of the Declaration of Helsinki and was approved by the ethics committee of Thammasat University, Thailand.

Tissue Sampling
Full-thickness skin biopsies were collected from the three patients using a 6 mm punch biopsy under local anesthesia (2% lidocaine with adrenaline). The PE skin was taken at the edge of the plaque (Figure 9), and UN skin was collected from phenotypically normal skin a minimum of 10 cm away from the plaque from the same patient. The biopsies were preserved in RNAlater (ThermoFisher, AMBION, Austin, TX, USA) at −80 • C until gene expression profiling.

High-Throughput Next-Generation Sequencing (NGS)
Total RNA was extracted from each skin biopsy of the volunteers (PE skin, n = 3, UN skin, n = 2) using TRIzol Reagent (Ambion, Austin, TX, USA). NGS was performed by Macrogen (Seoul, Korea). Briefly, the purity and structural integrity of RNA were analyzed using an Agilent 2100 Bioanalyzer (RNA integrity number > 7). Sequence libraries were constructed using the SMARTer Universal Low Input RNA Kit and TruSeq RNA Sample Prep Kit v2 (pair-end). NGS was performed using the NovaSeq 6000 S4 Reagent Kit on the NovaSeq 6000 System (Illumina Inc., San Diego, CA, USA) with 100-bp pairedend reads. The raw data generated in this study have been deposited in the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo, accessed on 8 September 2021) and are accessible through the GEO Series accession number GSE183732.

Sequencing Data Analysis
Basic data analysis was conducted by Macrogen. Briefly, overall read qualities, total bases, total reads, guanine-cytosine content (%GC content), and basic statistics were calculated. Trimmed reads were mapped to the reference genome using HISAT2, a spliceaware aligner. The transcripts were assembled using StringTie with aligned reads. Expression profiles were repeated as reading counts and normalized based on the transcript length and depth of coverage. The fragments per kilobase of transcript per million mapped reads (FPKM) value was used as a normalization value. Furthermore, we used the read counts to statistically analyze the expression profile with edgeR Bioconductor statistical library version 3.13 [89,90] on R Studio [91]. The DEGs between the PE skin and UN skin were determined according to p < 0.01 (0.026 ≤ FDR ≤ 0.194) and fold change (FC) ratio (|log2FC|) ≥ 1.5.

DEGs Analysis
All DEGs with expression profiles were uploaded into QIAGEN's Ingenuity Pathway Analysis (IPA) software (QIAGEN Redwood City, CA, USA; www.qiagen.com/ingenuity, accessed on 1 October 2021) using the "core analysis" of DEGs. Disease and biofunction, canonical pathways, upstream regulators, and regulator effects of the DEGs were analyzed based on the Ingenuity Knowledge Base. Pathway and function enrichment analysis

High-Throughput Next-Generation Sequencing (NGS)
Total RNA was extracted from each skin biopsy of the volunteers (PE skin, n = 3, UN skin, n = 2) using TRIzol Reagent (Ambion, Austin, TX, USA). NGS was performed by Macrogen (Seoul, Korea). Briefly, the purity and structural integrity of RNA were analyzed using an Agilent 2100 Bioanalyzer (RNA integrity number > 7). Sequence libraries were constructed using the SMARTer Universal Low Input RNA Kit and TruSeq RNA Sample Prep Kit v2 (pair-end). NGS was performed using the NovaSeq 6000 S4 Reagent Kit on the NovaSeq 6000 System (Illumina Inc., San Diego, CA, USA) with 100-bp paired-end reads. The raw data generated in this study have been deposited in the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo, accessed on 8 September 2021) and are accessible through the GEO Series accession number GSE183732.

Sequencing Data Analysis
Basic data analysis was conducted by Macrogen. Briefly, overall read qualities, total bases, total reads, guanine-cytosine content (%GC content), and basic statistics were calculated. Trimmed reads were mapped to the reference genome using HISAT2, a splice-aware aligner. The transcripts were assembled using StringTie with aligned reads. Expression profiles were repeated as reading counts and normalized based on the transcript length and depth of coverage. The fragments per kilobase of transcript per million mapped reads (FPKM) value was used as a normalization value. Furthermore, we used the read counts to statistically analyze the expression profile with edgeR Bioconductor statistical library version 3.13 [89,90] on R Studio [91]. The DEGs between the PE skin and UN skin were determined according to p < 0.01 (0.026 ≤ FDR ≤ 0.194) and fold change (FC) ratio (|log2FC|) ≥ 1.5.

DEGs Analysis
All DEGs with expression profiles were uploaded into QIAGEN's Ingenuity Pathway Analysis (IPA) software (QIAGEN Redwood City, CA, USA; www.qiagen.com/ingenuity, accessed on 1 October 2021) using the "core analysis" of DEGs. Disease and biofunction, canonical pathways, upstream regulators, and regulator effects of the DEGs were analyzed based on the Ingenuity Knowledge Base. Pathway and function enrichment analysis was additionally analyzed using Metascape Online Tool (https://metascape.org, accessed on 1 October 2021) [92].

Quantitative Multiplex Real-Time PCR (Multiplex qPCR)
For validation experiments, UN and PE skin samples were biopsied from moderate to severe psoriasis vulgaris patients. RNA was extracted and reverse transcribed to complementary DNA using the ImProm-II TM Reverse Transcription Kit (Promega Corp., Madison, WI, USA). TagMan ® Gene Expression Assay predesigned qPCR primers of the selected DEGs (IL36A, CXCL8, PLSCR1, SERPINB4, DEFB4A/B, DEFB103A/B, S100A7, S100A7A, S100A8, and S100A9), upstream regulators (IL17C, NFKBIZ, and ROCK2), and of housekeeping gene (GAPDH) for qPCR were purchased from Gene Plus Co., Ltd. (ThermoFisher Scientific, West Sacramento, CA, USA). Their assay IDs are listed in Supplementary Table S11. The multiplex qPCR was performed using the QuantiStudio TM 6 Flex Real-Time PCR system (Applied Biosystem). A statistical difference in the mean of ∆Ct (targeted gene vs GAPDH) between the UN and PE skin groups was calculated by a paired t-test. Graphs created by GraphPad Prism 9.3.1 were shown for relative mRNA expression using the 2 −∆Ct method and for fold change using the 2 −∆∆Ct method [93].

Data Availability Statement:
The raw data generated in this study have been deposited in the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo (accessed on 1 March 2022)) and are accessible through the GEO Series accession number GSE183732.