Transcriptional landscape of psoriasis identifies the involvement of IL36 and IL36RN.
Journal: 2016/January - BMC Genomics
ISSN: 1471-2164
Abstract:
BACKGROUND
In present study we performed whole transcriptome analysis in plaque psoriasis patients and compared lesional skin with non-lesional skin and with the skin from healthy controls. We sequenced total RNA from 12 lesional (LP), 12 non-lesional (NLP) and from 12 normal (C) skin biopsies.
RESULTS
Compared with previous gene expression profiling studies we had three groups under analysis - LP, NLP and C. Using NLP samples allows to see the transcriptome of visually normal skin from psoriasis patient. In LP skin S100A12, S100A7A, LCE3E, DEFB4A, IL19 were found up regulated. In addition to already these well-described genes, we also found several other genes related to psoriasis. Namely, KLK9, OAS2, OAS3, PLA2G, IL36G, IL36RN were found to be significantly and consistently related to the psoriatic lesions and this finding is supported also by previous studies. The genes up-regulated in the LP samples were related to the innate immunity, IL17 and IL10 networks. In NLP samples innate immunity and IL17 network were activated, but activation of IL10 network was not evident. The transcriptional changes characteristic in the NLP samples can be considered as a molecular signature of "dormant psoriasis".
CONCLUSIONS
Taken together, our study described the transcriptome profile characteristic for LP and NLP psoriatic skin. RNA profile of the NLP skin is in between the lesional and healthy skin, with its own specific pattern. We found that both LP and NLP have up-regulated IL17 network, whereas LP skin has up regulated IL10 related cytokines (IL19, IL20, IL24). Moreover, IL36G and IL36RN were identified as strong regulators of skin pathology in both LP and NLP skin samples, with stronger influence in LP samples.
Relations:
Content
Citations
(15)
References
(42)
Diseases
(1)
Chemicals
(2)
Genes
(2)
Organisms
(1)
Processes
(3)
Anatomy
(1)
Affiliates
(2)
Similar articles
Articles by the same authors
Discussion board
BMC Genomics. Dec/31/2014; 16(1)
Published online Apr/18/2015

Transcriptional landscape of psoriasis identifies the involvement of IL36 and IL36RN

Abstract

Background

In present study we performed whole transcriptome analysis in plaque psoriasis patients and compared lesional skin with non-lesional skin and with the skin from healthy controls. We sequenced total RNA from 12 lesional (LP), 12 non-lesional (NLP) and from 12 normal (C) skin biopsies.

Results

Compared with previous gene expression profiling studies we had three groups under analysis - LP, NLP and C. Using NLP samples allows to see the transcriptome of visually normal skin from psoriasis patient. In LP skin S100A12, S100A7A, LCE3E, DEFB4A, IL19 were found up regulated. In addition to already these well-described genes, we also found several other genes related to psoriasis. Namely, KLK9, OAS2, OAS3, PLA2G, IL36G, IL36RN were found to be significantly and consistently related to the psoriatic lesions and this finding is supported also by previous studies. The genes up-regulated in the LP samples were related to the innate immunity, IL17 and IL10 networks. In NLP samples innate immunity and IL17 network were activated, but activation of IL10 network was not evident. The transcriptional changes characteristic in the NLP samples can be considered as a molecular signature of “dormant psoriasis”.

Conclusions

Taken together, our study described the transcriptome profile characteristic for LP and NLP psoriatic skin. RNA profile of the NLP skin is in between the lesional and healthy skin, with its own specific pattern. We found that both LP and NLP have up-regulated IL17 network, whereas LP skin has up regulated IL10 related cytokines (IL19, IL20, IL24). Moreover, IL36G and IL36RN were identified as strong regulators of skin pathology in both LP and NLP skin samples, with stronger influence in LP samples.

Electronic supplementary material

The online version of this article (doi:10.1186/s12864-015-1508-2) contains supplementary material, which is available to authorized users.

Background

Psoriasis is one of the most prevalent chronic inflammatory disease affecting skin and joints [1]. The disease affects 2–3% of the population worldwide and it can have variable clinical course and severity [2]. The commonest form of psoriasis is plaque psoriasis and chronic, symmetrical, silvery-scaled plaques characterize this form [1,3]. Although the cause of psoriasis is unknown, it is complex disease with multifactorial pathogenesis where several genes interact and eventually induce active disease [4-6].

Whole-genome analysis of RNA provides useful tool to study complex disorders and to find common pathways related to the phenotype. Studies of the psoriatic transcriptome with microarrays and RNA-seq technologies have revealed very large number of differentially expressed genes (DEGs) in lesional skin [7-12]. With microarrays more than 1000 genes have been found differentially regulated between psoriatic and normal control skin [8,10]. These studies have resulted in the new genes and candidates involved in the pathogenesis of psoriasis. A more recent study has found that 80% of significantly elevated genes in psoriasis lesions are related to the keratinocyte activity and infiltration by T-cells and macrophages [9]. In addition to genechips, recently evolved RNA-seq provides more comprehensive overview about the transcriptional landscape. RNA-seq is able to detect new transcripts and splicing forms that are undetectable with other tools. There have been only two RNA-seq based studies on psoriasis so far [12,13]. Both of these studies described several new transcripts involved in the psoriasis. In one study three pairs of lesional and normal skin samples from psoriatic individuals were analysed [13]. In the most recent study based on biopsies of 92 psoriatic patients and 82 normal individuals 3,577 DEGs between lesional and normal skin were described [12]. In present study, we performed RNA-seq transcriptome analysis of 12 paired lesional (LP) and non-lesional skin (NLP) samples with 12 samples from healthy controls (C). Differential gene expression analysis was combined with functional network annotation. As a result we found gene expression pattern characteristic for the psoriatic lesions and described the molecular signature of the non-lesional skin of psoriatic patient.

Methods

Patients and controls

The Ethics Review Committee on Human Research of the University of Tartu approved the protocols and informed consent forms used in this study. All the participants signed a written informed consent. The patients and control subjects in the study were unrelated Caucasians living in Estonia. Unrelated patients with plaque psoriasis from the Dermatology Clinic of Tartu University Hospital were included in the study. Sex- and age-matched (+/− 10 years) control subjects were recruited among patients with melanocytic nevi at the dermatologic outpatient clinic. Healthy volunteers were free from inflammatory dermatoses and without a positive family history of psoriasis. The main characteristics of the psoriasis patients are presented in Table 1.

Table 1
Characteristics of psoriasis patients in present study
PatientAge, yearsSexAoO (duration), yearsPASINail involvementPsAScalp involvementFlexural ininvolvementTreatment
P84719M18 (0.5)12.5NoNoYesNoNo
P65225M13 (12)18.4YesNoYesYesNo
P84427M18 (9)8.8NoNoYesNoNo
P84829M26 (3)3.7NoNoYesNoNo
P84049M15 (34)6.0NoNoYesYesCalcipotriol + betamethasone
P84552M22 (30)14.0YesNoYesYesNo
P85160M56 (4)10.8NoNoYesNoNo
P85328F28 (0.5)23.3NoNoYesNoNo
P84937F30 (7)4.7NoNoNoNoTopical steroid
P85654F14 (40)7.3NoYesYesYesTopical steroid
P84357F53 (4)12.6NoNoNoNoTopical steroid
P84658F57 (1)15.2NoNoYesNoTopical steroid
AoO = Age of Onset, M = male, F = female, PASI = Psoriasis Area and Severity Index, PsA = psoriatic arthritis.

Two 4 mm punch biopsy specimens were taken from patients with psoriasis, one from the lesional (LP sample) and another from distant non-sun-exposed nonlesional (NLP sample) skin. One 4 mm punch biopsy specimen was taken from non-sun-exposed skin of each control subject (C sample). All biopsy specimens were instantly frozen in liquid nitrogen and stored at -80C° until RNA extraction.

RNA sequencing

Precellys24 tissue homogenizer (Bertin Technologies, France) with Cryolys adaptor was used to homogenise biopsy specimens. Total RNA was extracted with RNeasy Fibrous Tissue Mini Kit (Qiagen, Valencia CA, USA) according to the manufacturer’s protocol. Isolated RNA was dissolved in RNase free water and stored in -80C°. The quality of total RNA was evaluated with Agilent 2100 Bioanalyzer and the RNA 6000 Nano Kit (Agilent Technologies Inc., CA, USA). The average RNA integrity number (RIN) of samples was ≥7.

50 ng of total RNA was amplified by applying Ovation RNA-Seq System V2 (NuGen, Emeryville, CA, USA) after which the resulting cDNAs were pooled in equal amount and the pool was used to prepare the DNA fragment library with SOLiD System chemistry (Life Technologies Corp., Carlsbad, CA, USA). Sequencing was performed using SOLiD 5500W platform and DNA sequencing chemistry (Life Technologies Corp., Carlsbad, CA, USA). Raw reads (75 bp) were color-space mapped to the human genome hg19 reference using Maxmapper algorithm implemented in the Lifescope software (Life Technologies, Ltd). Mapping to multiple locations was permitted. The quality threshold was set to 10, giving the mapping confidence was more than 90. Reads with score less than 10 were filtered out. Average mapping quality was 30. Analysis of the RNA content and gene-based annotation was done with whole-transcriptome workflow. Raw sequencing data with appropriate experimental information is available in the NCBI Gene Expression Omnibus (GEO) repository under the accession number GSE66511 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66511).

Statistical analysis

Patients and controls were tested for potential bias caused by covariates.

Non-normalized raw counts were used for the EdgeR package to perform differential gene expression analysis after quality control of samples. EdgeR performs model-based scale normalization, estimates dispersions and applies negative binomial model. EdgeR is very flexible tool for RNAseq data analysis to find differentially expressed genes [14,15]. It implements negative binomial model fitting followed by testing procedures for determining differential expression.

As our sample contains paired-samples (lesional and non-lesional skins are form the same persons) we used two different approaches. For the paired samples we applied general linear modelling where the subjects were added to the contrast matrix. GLM likelihood ratio test was applied for LP-NLP comparison. Different approach was applied for comparisons between LP-C and NLP-C. In this approach we used group-wise comparisons where negative binomial fitting was followed by exact test. False discovery rate (FDR) adjustment was used for multiple testing correction [16]. FDR threshold 0.1 for statistical significance was applied. Genes with larger differential expression were defined with logFC threshold 0.5 (i.e. 50% change between experimental conditions).

Functional annotation of differentially expressed genes

The functional analysis of a gene network was used to identify the biological functions that are most significantly related to the molecules in the network. To define the functional networks of differentially expressed genes, data were analyzed by using the Ingenuity Pathway Analysis (IPA, Ingenuity Systems, www.ingenuity.com) that calculates a significance score (network score) for each network. This score indicates whether the likelihood that the assembly of a set of focus genes in a network could be explained by random chance alone (e.g., score of 2 indicates that there is a chance of 1 in 100 that the focus genes are together in a network due to random chance). A data set containing the gene identifiers and their corresponding fold change (log2) values were uploaded into the IPA software. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base to identify molecules whose expression was significantly differentially regulated (focus genes or Networks Eligible molecules). These focus genes were overlaid onto a global molecular network developed from information contained in the Ingenuity Knowledge Base. Networks of these focus genes were then algorithmically generated based on their connectivity. A network represents the molecular relationships between genes or gene products, which are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). All edges are supported by at least one reference from the literature, or from canonical information stored in the Ingenuity Pathways Knowledge Base. Fold change difference threshold for the functional analysis was set on 2 (log2 transformed).

Quantitative real-time PCR

250 ng of total RNA was used with High Capacity RNA-to-cDNA Kit (Life Technologies Co, USA) for cDNA synthesis. QuantiTect Reverse Transcription Kit (Qiagen) was applied for cDNA synthesis in case of g1 or s1 assays. Both cDNA synthesis kits were used according to the manufacturer’s protocol. cDNA was used as a template for TaqMan qRT-PCR in 7900 Fast QRT-PCR System (Life Technologies Co). Two primers and labelled probe were used to detect the mRNA expression level of the reference gene hypoxanthine phosphoribosyl-transferase-1 (HPRT1; primer sequences available upon request). Expression levels of IL36G, IL36RN, LCE3D, IFI6, TGM1, SPRR2B, OAS2, CRABP2, S100A12 and PARP9 were detected with TaqMan Gene Expression Assays (Life Technologies Co) Hs00219742_m1, Hs01104220_g1, Hs00754375_s1, Hs00242571_m1, Hs01070310_m1, Hs01595682_s1, Hs00942643_m1, Hs00275636_m1, Hs00942835_g1 and Hs00967084_m1, respectively.

The relative gene expression levels from qRT-PCR were calculated relative to the reference gene HPRT-1 using 2-ΔCT method. The normality of the data was tested with Shapiro test and based on its results unpaired t-test or Mann–Whitney test was used. We also performed correlation analysis to find the strength of expressional co-regulation of selected genes.

Results

The effect of covariates

Present study compared skin RNA expressional profiles in psoriasis patients with healthy controls (C). In addition, we compared the lesional (LP) and non-lesional (NLP) skin biopsies from psoriasis patients in order to find activity related genes. In order to exclude most common confounding variables, we tested for the differences in age, body weight, height, BMI and smoking status. There was no statistical difference between any of the studied covariate. In case of age (38 versus 41) the p-value was 0.61, for height (173 versus 174 cm) p-value was 0.77, for body weight (79 kg versus 87 kg) p-value was 0.28, for BMI (26.2 versus 29.0) p-value was 0.33. Also, the smoking status did not differ between the study groups (p-value 0.4).

RNA sequencing results

RNA sequencing experiment gave high quality reads with good similarity between different samples (Figure 1a). Multidimensional scaling analysis indicated good separation of study groups based on the gene expression fold-change difference (Figure 1a). In Figure 1a the clear separation of LP and C samples is evident, whereas NLP samples are between LP and C. This is expected, as in the LP skin active inflammation is evident. NLP is a skin without inflammation, but nevertheless, non-lesional skin of psoriasis patients has molecular differences compared to the normal skin. Indeed, differential expression analysis found significant differences between all study groups.

Figure 1

General illustration of study groups. a. Multidimensional scaling plot of the original data indicates clear separation of experimental groups – controls (C), non-lesional psoriatic skin (NLP) and lesional psoriatic skin (LP). This plot verifies good sampling and confirms reliability of collected data. b. Venn diagram for different comparisons between study groups.

We had three different groups in our model: lesional psoriatic skin (LP, n = 12) paired with the non-lesional psoriatic skin (NLP, n = 12) and normal skin from healthy controls (C, n = 12). We applied generalized linear model (GLM) for paired sample analysis. Comparison between other groups was done with the exact test (NLP versus C and LP versus C).

In the lesional psoriatic skin compared to the skin of healthy controls (LP-C), 7,932 genes were differentially expressed with the confidence level of FDR < 0.1 (for 30 genes Table 2, for all genes Additional file 1: Table S1). Applying of LogFC threshold verified that 5,853 genes are at least 50% different between LP-C (FDR < 0.1). 7,539 genes were differentially expressed (False Discovery Rate, FDR < 0.1) after pair-wise comparison of lesional psoriatic skin to the non-lesional skin of psoriasis patients (LP-NLP) (Additional file 2: Table S2). After logFC filtering 4,626 genes remained differentially expressed at least by 50% between LP-NLP. Comparison between non-lesional psoriatic skin and the skin of healthy controls (NLP-C) identified 2,595 genes to be differentially expressed with FDR < 0.1 (Additional file 3: Table S3). Additional logFC threshold indicated that 970 genes were 50% differentially expressed between NLP and C. The number of differentially expressed genes between different groups is illustrated in the Venn diagram (Figure 1b).

Table 2
Differentially expressed genes in psoriatic lesions compared to healthy controls (LP-C)
SymbollogFClogCPMPValueFDREntrez Gene Name
OAS24.437.169.26E-702.14E-652'-5'-oligoadenylate synthetase 2
PLA2G4E2.557.602.74E-683.16E-64Phospholipase A2, group IVE
KLK95.204.012.20E-671.69E-63Kallikrein-related peptidase 9
S100A129.004.411.25E-657.20E-62S100 calcium binding protein A12
LCE3E6.994.456.02E-652.78E-61Late cornified envelope 3E
TGM13.426.514.94E-611.90E-57Transglutaminase 1
OAS32.997.251.01E-593.33E-562'-5'-oligoadenylate synthetase 3
PARP92.207.112.80E-598.08E-56Poly (ADP-ribose) polymerase, 9
CRABP22.495.651.99E-575.10E-54Cellular retinoic acid binding protein 2
PLA2G4D4.336.103.03E-567.01E-53Phospholipase A2, group IVD
IL1F95.646.612.26E-554.74E-52Interleukin 36, gamma
ALOX12B3.017.307.12E-551.37E-51Arachidonate 12-lipoxygenase
SAMD93.016.238.59E-551.53E-51Sterile alpha motif domain
IL1F53.597.911.14E-531.87E-50Interleukin 36 receptor antagonist
C10orf996.225.443.20E-524.93E-49Chrom 10 open reading frame 99
DEFB4A10.436.596.30E-529.09E-49Defensin. beta 4A
AKR1B106.224.394.02E-515.46E-48Aldo-keto reductase family 1
PAPL4.014.601.56E-502.01E-47Iron/zinc purple acid protein
GLTP2.017.425.00E-506.08E-47Glycolipid transfer protein
KDM6B−1.3611.511.48E-491.71E-46Lysine (K)-specific demethylase 6B
KYNU4.494.722.54E-482.79E-45Kynureninase
RRM22.685.473.29E-483.46E-45Ribonucleotide reductase M2
ZC3H12A3.064.882.30E-472.31E-44Zinc finger CCCH-type cont 12A
SDR9C72.705.006.01E-475.79E-44Short chain dehydrogenase/reductase
HPSE3.474.171.31E-461.21E-43Heparanase
APOL62.236.521.91E-461.69E-43Apolipoprotein L, 6
NIPAL41.399.136.90E-465.90E-43NIPA-like domain containing 4
DMD−1.817.469.69E-467.99E-43Dystrophin
S100A7A10.0210.302.84E-452.26E-42S100 calcium binding protein A7A

LogFC is fold changes differences in log2 scale and it describes how many times gene expression differs between groups. Positive values indicate up-regulation in psoriasis. LogCPM (log2 counts-per-million) is average gene expression signal in all samples. FDR is genome-wide corrected P-value.

Figure 2
Heatmap of the 50 genes with largest fold change differences between NLP and C samples. Violet bar is for control samples, red bar is for non-lesional samples and green bar is for lesional skin samples.

To illustrate the observed differences and their relations to the skin conditions, we performed cluster analysis of the entire sample. Heatmap (Additional file 4: Figure S2) of top 200 differing genes between LP and C samples with the lowest FDR values were generated. Moreover, we generated heatmaps from smaller sample subset based on their gene expression fold-changes (Figure 2, Additional file 5: Figure S1). Again, our samples were distinguishable based on the gene expression data what confirms the quality of RNAseq data. An interesting pattern was found in Figure 2. This is comparison of NLP-C and therefore should illustrate the genes that are different in non-lesional and control sample. NLP sample has genes that have similar expressional pattern to the lesional sample (central block of genes in Figure 2) and about the same amount of genes that have similar expressional pattern to the control sample (two blocks in the upper and lower part of Figure 2). The genes with expressional pattern similar to the lesional sample, represent the genes that form the molecular signature of psoriasis in non-lesional skin. These genes can be described as “dormant pathology” or “background inflammation”, indicating certain standby situation in non-lesional skin. Interestingly, among these genes are IL1F6 and IL6. IL6 is well known to be involved in the pathogenesis of psoriasis [17,18]. IL1F6 (IL36A) is a gene that was recently found to be part of the signalling system what is active in psoriasis [19,20].

Figure 3
Quantitative real-time PCR analysis confirmed transcriptional differences found with RNA-seq. In case of TGM1, LCE3D and SPRR2B there is no difference between C and NLP. *** - p < 0.001. ** - p < 0.01.

Functional analysis of gene expression networks

Ingenuity pathway analysis software was used for functional analysis. We analysed the genes differentially expressed between all three groups and performed three separate analyses. This enabled us to describe the functional networks that characterize non-lesional skin and lesional skin. Comparison between lesional skin (LP) and healthy controls (C) revealed activation of the “Granulocyte adhesion and diapedesis” (Table 3) as the top enriched canonical pathway (p = 2.81E-20). In addition, “Agranylocyte adhesion and diapedesis” (p = 6.08E-13) and “Role of IL-17A in Psoriasis” (p = 8.52E-11) were two other highly enriched canonical pathways in LP compared to C sample. In addition to the IL-17 signalling, canonical pathway “IL-10 signalling” was also significantly enriched in the lesional samples compared to control samples. Activation of IL-10 signalling means that IL19, IL20, IL24 and IL26 were significantly up regulated, whereas IL-10 was not. Activation of IL10 canonical signalling was not found in NLP samples.

Table 3
Activated canonical pathways in LP-C comparisons
Ingenuity canonical pathways-log(p-value)Ratio
Granulocyte Adhesion and Diapedesis19.600.31
Agranulocyte Adhesion and Diapedesis12.200.25
Role of IL-17A in Psoriasis10.100.79
Atherosclerosis Signaling9.750.25
Pathogenesis of Multiple Sclerosis9.730.90
Differential Regulation of Cytokine Production in Intestinal Epithelial Cells by IL-17A and IL-17F9.540.61
Differential Regulation of Cytokine Production in Macrophages and T Helper Cells by IL-17A and IL-17F8.910.67
Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza8.380.39
LXR/RXR Activation6.290.20
Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses6.130.22
Role of Cytokines in Mediating Communication between Immune Cells5.890.31
Altered T Cell and B Cell Signaling in Rheumatoid Arthritis5.850.22
T Helper Cell Differentiation5.660.26
Hepatic Fibrosis / Hepatic Stellate Cell Activation5.030.18
Communication between Innate and Adaptive Immune Cells4.860.19
IL-10 Signaling3.950.21
Differentially expressed genes were used as a “signature” to find what biological function is changed in the skin of psoriasis patients.

Comparison between LP and NLP samples (LP-NLP) indicated the enrichment of the “Granulocyte adhesion and diapedesis” as the top canonical pathway (p = 1.13E-13, Additional file 6: Table S4), similarly to the LP-C analysis. LP samples also exhibited significant activation of the canonical pathway “Role of cytokines in mediating communication between immune cells” (p = 2.04E-11). Comparison between NLP and C samples found up-regulation of canonical pathway “Role of IL-17A in Psoriasis” (p = 2.31E-08) (Additional file 7: Table S5). However, the enrichment was not that large as in case of lesional skin and the number of differentially expressed genes was reduced on NLP-C comparison (Additional file 3: Table S3). Therefore, differences were more quantitative than qualitative (Additional file 8: Figure S4 and Additional file 9: Figure S3. However, there were some gene sets forming quite distinctive pattern or footprint characteristic only for the LP sample (Additional file 9: Figure S3). For instance, the genes LCE3E, SPRR2B, DEFB4A, S100A7A, S100A8, SERPINB4 and S100A12 were specific only for the lesional and not found in non-lesional sample (Figure 3). Activation of “IL-10 signalling” was also found only in the LP samples. On the other hand, IL1F6, IL6, GPR151, WIF1 and TNFRSF10C were activated in both lesional and non-lesional samples (Figure 3). Therefore, there are genes characteristic for the lesional skin and there are genes what are activated also in the skin without psoriatic lesion. The functional impact of non-lesional transcriptional signature needs further studies.

Confirmative quantitative real-time PCR

In order to verify the RNA seq findings we performed quantitative real-time PCR analysis for 10 genes. We analysed the expression of IL1F5 (IL36RN) and IL1F9 (IL36G) genes and found that their expression profile correlated very well with the inflammation in the skin. The results from qRT-PCR were similar with RNA-seq data and the differences in gene expression levels followed completely the skin status – C, NLP or LP. Only in case of TGM1, LCE3D and SPRR2B genes there wasn’t statistically significant difference between NLP and C, but still the difference between LP-C and LP-NLP was significant (Figures 3). These results support that LCE3D and SPRR2B are genes more characteristic for the psoriatic lesion, than for the non-lesional skin. Correlation analysis indicated very high and significant correlation between the analysed genes (Additional file 10: Figure S6).

Discussion

In present study we performed whole transcriptome analysis of lesional and non-lesional psoriasis skin samples and compared them with the healthy control samples. There are several previous studies where transcriptome changes in psoriasis patients have been analysed [8-11,21-27]. In most of these studies hybridization based genechips were used. However, recently several additional papers reported the results of RNA-seq based transcriptome analysis [9,12,13,28,29]. One of these studies was focused on the differences in the small RNAs [28]. The other studies analysed whole genome transcriptome [9,12,13,29]. One major difference between RNA-seq and genechip based transcriptome studies is the substantially larger number of differentially expressed genes and the appearance of the genes previously not discovered in similar studies. This difference can be explained by the lower detection limit and wider dynamic range that is characteristic for RNA-seq. However, the majority of changes described in the genechip based studies and in published RNA-seq studies are quite similar. Most of the findings described in our paper coincide with previously published large scale RNA-seq study [12]. However, we applied slightly different approach and we also found some interesting targets not described in previous studies. We used total RNA (not only polyA RNA), we applied only gene-level analysis and functional annotation of the networks related to the differential transcriptome signatures. Also, the statistical approaches in previous studies have been slightly different. However, despite these technical differences, the results of our study and previous studies are generally similar.

Our study revealed significant and substantial differences between all three groups in their transcriptome profile. Most remarkably, oligoadenylate synthetases (OAS2 and OAS3) and phospholipases from A2 (PLA2G4E and PLA2G4D) group were found up regulated in lesional psoriatic skin. OAS2 and OAS3 are related to the innate immunity and antiviral response [30]. Recent studies have shown the association of OAS genes with psoriasis and systemic lupus erythematosus [31,32]. OAS genes have thoroughly been discussed as part of the IFN signature in several studies focusing on psoriasis [24,33]. Phospholipases A2 are enzymes related to the metabolism of fatty acids in membrane phospholipids and related to the inflammatory signalling pathways [34]. Lipid metabolism and antimicrobial defence pathways are shown activated during psoriasis [10]. Moreover, OAS and PLA genes have also been found in other psoriasis transcriptome studies [12,29]. Therefore, these genes are clearly linked with the psoriasis and their general role in the chronic inflammatory conditions is shown in previous papers.

Another interesting set of genes we found activated during psoriasis are interleukin-36 (IL36) cytokines that belong to interleukin-1 family (IL1F). Namely, with very high confidence (FDR below 1E-50) we found IL1F9 (IL36G) and IL1F5 (IL36RN) up regulated in lesional skin samples. Quantitative real-time PCR verified our findings and even confirmed IL36G and IL36RN expression to be significantly associated with the skin condition. Expression was statistically different even between C and NLP skin samples, in LP samples the expression was even more elevated. IL36G is a recently discovered novel member of IL-1 family of cytokines [35]. IL36G has been shown to be involved in innate immunity and all IL36 cytokines (except IL36RN) have proinflammatory activity [36-38]. We can consider the activation of IL36G as the most significant finding and propose the family of these genes to be important biomarkers for psoriasis. Indeed, one recent study identified IL36G as the most outstanding biomarker for psoriasis [39]. Several other studies confirm the involvement of IL36 cytokines in the psoriasis [19,36,40]. More precisely, IL36G was recently identified as one of the 13 hallmark psoriasis genes universally or near-universally up-regulated in psoriasis lesions [40]. IL-1 family members play a key role in the function of innate and adaptive immunity and are new promising targets for immunpathologies [41]. In our study we only found IL36 cytokines (IL36G, IL36RN and IL36A) to be related to psoriasis. IL37 (IL1F7) was up-regulated only in the non-lesional psoriatic skin, which makes its role in psoriasis very intriguing. In another previous study using RNA-seq for psoriasis samples the IL37 was found to be down-regulated in lesional skin [12]. Therefore, our results confirm previous studies, that IL36 genes have significant role in regulation of psoriasis and psoriasis activity.

Pathway analysis found clear activation of innate immune response and inflammation. Pathway “Role of IL-17A in Psoriasis” was enriched in both LP-C and NLP-C comparisons (Table 3 and Table 4). The finding that the IL-17 pathway is activated during psoriasis is not new and suits with findings from previous studies [42]. Moreover, the IL36 cytokine discussed in previous section is also linked to the IL-17 network. Granulocyte and agranulocyte adhesion canonical pathways illustrate the activation of inflammatory pathways and inflammatory conditions. Additional file 8: Figure S4 illustrates gene expression profiles compared between LP and C samples. It is visible that the majority of genes are up-regulated and they are involved in inflammation. Comparison of Additional file 8: Figures S4, Additional file 9: Figure S3 and Additional file 11: Figure S5 illustrates how the number of activated genes is decreasing as inflammatory conditions getting milder (from lesional to non-lesional and healthy skin). Additional file 9: Figure S3 has remarkably less activated genes from the “Psoriasis” genetic network. On the other hand, the Additional file 10: Figure S6 illustrates genes that are related to the background inflammation in the non-lesional psoriatic skin. Most well known psoriasis related genes activated in non-lesional skin are IL6, IL22, IL36A, IL36G, IL19, IL20, S100A7 and S100A12. One previous study has also described similar “pre-psoriatic” gene expression signature [10]. In accordance with this study, we found that the gene expression differences between LP and C samples were more prominent than differences between NLP and C samples. We can conclude that non-lesional skin in psoriasis patients, while considered normal, still contains molecular signature characteristic for psoriasis.

Table 4
Activated canonical pathways in NLP-C comparisons
Ingenuity canonical pathways-log(p-value)Ratio
Role of IL-17A in Psoriasis7.640.50
γ-linolenate Biosynthesis II (Animals)4.050.21
Granulocyte Adhesion and Diapedesis3.800.08
Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza3.620.15
Atherosclerosis Signaling3.100.08
Altered T Cell and B Cell Signaling in Rheumatoid Arthritis3.090.09
Agranulocyte Adhesion and Diapedesis3.020.07
Role of Cytokines in Mediating Communication between Immune Cells3.020.13
Differential Regulation of Cytokine Production in Macrophages and T Helper Cells by IL-17A and IL-17F2.820.22
LPS/IL-1 Mediated Inhibition of RXR Function2.720.06
Stearate Biosynthesis I (Animals)2.510.10
Primary Immunodeficiency Signaling2.440.09
Differential Regulation of Cytokine Production in Intestinal Epithelial Cells by IL-17A and IL-17F2.410.17
IL-17 Signaling2.340.09
Communication between Innate and Adaptive Immune Cells2.330.07
Uracil Degradation II (Reductive)2.320.18
Differentially expressed genes were used as a “signature” to find what biological function is changed in the skin of psoriasis patients.

Difference between NLP and C samples may be caused by the confounding covariates that should be considered. After analysis of the distribution of covariates in study groups we found no statistical difference in age (p = 0.61), body weight (p = 0.30), height (p = 0.77), BMI (p = 0.32) and smoking status (p = 0.41) between controls and patients. Based on these data we can assume, that the influence of covariates is minimal and the transcriptional differences are caused by the psoriasis.

Taken together, the most of the genes we found have also been identified in earlier studies. Here we described results of complex analysis of the transcriptional signatures characteristic for lesional and non-lesional psoriatic skin and normal skin. There are several published results on the transcriptional changes of the psoriasis and in large extend the results are coinciding. However, differences exist in technical details (total RNA versus polyA RNA), statistical models and study designs (lesional versus non-lesional or lesional-non-lesional-normal). Our study design allowed us to use wider model (comparison between three conditions) and to gain more detailed information about the psoriasis. We found significant activation of IL36 cytokines, what has not been described in previous studies. Moreover, we found that IL17 related cytokines are also activated in non-lesional psoriatic biopsies indicating more general immune pathology in the skin. Therefore, even the skin looks noormal, immune system inside the skin is still altered.

Conclusions

We found the genes of IL36 cytokine family to be involved in the pathogenesis of psoriasis. These genes have received very little attention in previous research and therefore, further exploration of their regulation can give new insights to the psoriasis.

Availability of supporting data

The data set supporting the results of this article is available in the NCBI Gene Expression Omnibus (NCBI-GEO) repository. The accession number for dataset is GSE66511 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66511).

Additional files

Additional file 1: Table S1.

The list of 500 most significant differentially expressed genes based on comparison of LP to the C. LogFC is fold changes differences in log2 scale and it describes how many times gene expression differs between groups. Positive values indicate up-regulation in psoriasis. LogCPM (log2 counts-per-million) is average gene expression signal in all samples. FDR is genome-wide corrected P-value.

Additional file 2: Table S2.

The list of 500 most significant differentially expressed genes based on comparison of LP to the NLP skin. LogFC is fold changes differences in log2 scale and it describes how many times gene expression differs between groups. Positive values indicate up-regulation in psoriasis. LogCPM (log2 counts-per-million) is average gene expression signal in all samples. FDR is genome-wide corrected P-value.

Additional file 3: Table S3.

The list of 500 most significant differentially expressed genes based on comparison of NLP skin to the C skin. LogFC is fold changes differences in log2 scale and it describes how many times gene expression differs between groups. Positive values indicate up-regulation in psoriasis. LogCPM (log2 counts-per-million) is average gene expression signal in all samples. FDR is genome-wide corrected P-value.

Additional file 4: Figure S2.
Heatmap of the 50 genes with largest fold change differences after pairwise comparison of LP and NLP samples. Violet bar is for control samples, red bar is for non-lesional samples and green bar is for lesional skin samples.
Additional file 5: Figure S1.
Heatmap of the 50 genes with largest fold change differences between LP and C samples. Violet bar is for control samples, red bar is for non-lesional samples and green bar is for lesional skin samples.
Additional file 6: Table S4.
Activated canonical pathways in LP sample compared to NLP sample.
Additional file 7: Table S5.
Activated canonical pathways in NLP-C comparions.
Additional file 8: Figure S4.
Activated gene expression network related to psoriasis based on the differential gene expression profile between LP and NLP. The gene expression network was constructed based on their biological role. The figure illustrates expressional changes (red genes are higher, green genes are lower expressed) of particular genes and their cellular location (extracellular, membrane, cytoplasm or nucleus) in relation with other genes in particular network.
Additional file 9: Figure S3.
Activated psoriasis-related gene expression network based on the differential gene expression profile between LP and C. The gene expression network was constructed based on their biological role. The figure illustrates expressional changes (red genes are higher, green genes are lower expressed) of particular genes and their cellular location (extracellular, membrane, cytoplasm or nucleus) in relation with other genes in particular network.
Additional file 10: Figure S6.
Correlation analysis of the genes analysed with quantitative real-time PCR.
Additional file 11: Figure S5.
Activated psoriasis gene expression network based on the differential gene expression profile between NLP and C. The gene expression network was constructed based on their biological role. The figure illustrates expressional changes (red genes are higher, green genes are lower expressed) of particular genes and their cellular location (extracellular, membrane, cytoplasm or nucleus) in relation with other genes in particular network.

Footnotes

Competing interests

All authors declare that they have no competing interests.

Authors’ contributions

MK participated in the study design, cohort collection and participated in manuscript writing. SK participated in RNA seq data analysis and manuscript writing. ER performed RNA sequencing. EP performed statistical analysis and participated in manuscript writing. KA participated in cohort collection and clinical data management. KK participated in clinical data collection, data analysis and interpretation and in the writing of manuscript. All authors read and approved the final manuscript.

Acknowledgments

This work was supported by personal research funding PUT177 (Külli Kingo, Kristi Abram) and institutional research funding IUT (IUT20-46, Sulev Kõks, Ene Reimann) of the Estonian Ministry of Education and Research, by the Centre of Translational Genomics of University of Tartu (SP1GVARENG, Sulev Kõks, Ele Prans) and by the European Regional Development Fund (Centre of Translational Medicine, University of Tartu, Sulev Kõks, Maris Keermann).

Abbreviations

LPLesional psoriatic skinNLPnon-lesional psoriatic skinCnormal skin of healthy controls

References

  • 1. NestleFOKaplanDHBarkerJPsoriasisN Engl J Med20093615496509[PubMed][Google Scholar]
  • 2. LebwohlMPsoriasisLancet200336193641197204[PubMed][Google Scholar]
  • 3. LowesMASuarez-FarinasMKruegerJGImmunology of psoriasisAnnu Rev Immunol20143222755[PubMed][Google Scholar]
  • 4. LiuYHelmsCLiaoWZabaLCDuanSGardnerJA genome-wide association study of psoriasis and psoriatic arthritis identifies new disease lociPLoS Genet200843e1000041[PubMed][Google Scholar]
  • 5. KõksSKingoKRätsepRKarelsonMSilmHVasarECombined haplotype analysis of the interleukin-19 and −20 genes: relationship to plaque-type psoriasisGenes Immun2004586627[PubMed][Google Scholar]
  • 6. KõksSKingoKVabritKRätsepRKarelsonMSilmHPossible relations between the polymorphisms of the cytokines IL-19, IL-20 and IL-24 and plaque-type psoriasisGenes Immun20056540715[PubMed][Google Scholar]
  • 7. BowcockAMShannonWDuFDuncanJCaoKAftergutKInsights into psoriasis and other inflammatory diseases from large-scale gene expression studiesHum Mol Genet200110171793805[PubMed][Google Scholar]
  • 8. GudjonssonJEDingJJohnstonATejasviTGuzmanAMNairRPAssessment of the psoriatic transcriptome in a large sample: additional regulated genes and comparisons with in vitro modelsJ Invest Dermatol20101307182940[PubMed][Google Scholar]
  • 9. SwindellWRJohnstonAVoorheesJJElderJTGudjonssonJEDissecting the psoriasis transcriptome: inflammatory- and cytokine-driven gene expression in lesions from 163 patientsBMC Genomics201314527[PubMed][Google Scholar]
  • 10. GudjonssonJEDingJLiXNairRPTejasviTQinZSGlobal gene expression analysis reveals evidence for decreased lipid biosynthesis and increased innate immunity in uninvolved psoriatic skinJ Invest Dermatol2009129122795804[PubMed][Google Scholar]
  • 11. Suarez-FarinasMLowesMAZabaLCKruegerJGEvaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (GSEA)PLoS One201054e10247[PubMed][Google Scholar]
  • 12. LiBTsoiLCSwindellWRGudjonssonJETejasviTJohnstonATranscriptome Analysis of Psoriasis in a Large Case–control Sample: RNA-Seq Provides Insights into Disease MechanismsJ Invest Dermatol20141347182838[PubMed][Google Scholar]
  • 13. JabbariASuarez-FarinasMDewellSKruegerJGTranscriptional profiling of psoriasis using RNA-seq reveals previously unidentified differentially expressed genesJ Invest Dermatol201213212469[PubMed][Google Scholar]
  • 14. McCarthyDJChenYSmythGKDifferential expression analysis of multifactor RNA-Seq experiments with respect to biological variationNucleic Acids Res20124010428897[PubMed][Google Scholar]
  • 15. RobinsonMDMcCarthyDJSmythGKedgeR: a Bioconductor package for differential expression analysis of digital gene expression dataBioinformatics201026113940[PubMed][Google Scholar]
  • 16. StoreyJDTibshiraniRStatistical significance for genomewide studiesProc Natl Acad Sci U S A20031001694405[PubMed][Google Scholar]
  • 17. LoYHToriiKSaitoCFuruhashiTMaedaAMoritaASerum IL-22 correlates with psoriatic severity and serum IL-6 correlates with susceptibility to phototherapyJ Dermatol Sci20105832257[PubMed][Google Scholar]
  • 18. BonifatiCSolmoneMTrentoEPietravalleMFazioMAmeglioFSerum interleukin-6 levels as an early marker of therapeutic response to UVB radiation and topical steroids in psoriatic patientsInt J Clin Lab Res19942421223[PubMed][Google Scholar]
  • 19. JohnstonAXingXGuzmanAMRiblettMLoydCMWardNLIL-1F5, −F6, −F8, and -F9: a novel IL-1 family signaling system that is active in psoriasis and promotes keratinocyte antimicrobial peptide expressionJ Immunol20111864261322[PubMed][Google Scholar]
  • 20. RussellCBRandHBiglerJKerkofKTimourMBautistaEGene expression profiles normalized in psoriatic skin by treatment with brodalumab, a human anti-IL-17 receptor monoclonal antibodyJ Immunol20141928382836[PubMed][Google Scholar]
  • 21. AinaliCValeyevNPereraGWilliamsAGudjonssonJEOuzounisCATranscriptome classification reveals molecular subtypes in psoriasisBMC Genomics201213472[PubMed][Google Scholar]
  • 22. ZhouXKruegerJGKaoMCLeeEDuFMenterANovel mechanisms of T-cell and dendritic cell activation revealed by profiling of psoriasis on the 63,100-element oligonucleotide arrayPhysiol Genomics20031316978[PubMed][Google Scholar]
  • 23. ReischlJSchwenkeSBeekmanJMMrowietzUSturzebecherSHeubachJFIncreased expression of Wnt5a in psoriatic plaquesJ Invest Dermatol200712711639[PubMed][Google Scholar]
  • 24. YaoYRichmanLMorehouseCde los ReyesMHiggsBWBoutrinAType I interferon: potential therapeutic target for psoriasis?PLoS One200837[PubMed][Google Scholar]
  • 26. Suarez-FarinasMLiKFuentes-DuculanJHaydenKBrodmerkelCKruegerJGExpanding the psoriasis disease profile: interrogation of the skin and serum of patients with moderate-to-severe psoriasisJ Invest Dermatol201213211255264[PubMed][Google Scholar]
  • 27. TianSKruegerJGLiKJabbariABrodmerkelCLowesMAMeta-analysis derived (MAD) transcriptome of psoriasis defines the "core" pathogenesis of diseasePLoS One201279e44274[PubMed][Google Scholar]
  • 28. JoyceCEZhouXXiaJRyanCThrashBMenterADeep sequencing of small RNAs from human skin reveals major alterations in the psoriasis miRNAomeHum Mol Genet20112020402540[PubMed][Google Scholar]
  • 29. SwindellWRXingXVoorheesJJElderJTJohnstonAGudjonssonJEIntegrative RNA-seq and microarray data analysis reveals GC content and gene length biases in the psoriasis transcriptomePhysiol Genomics2014461553346[PubMed][Google Scholar]
  • 30. KulkaMCalvoMSNgoDTWalesSQGoswamiBBActivation of the 2-5OAS/RNase L pathway in CVB1 or HAV/18f infected FRhK-4 cells does not require induction of OAS1 or OAS2 expressionVirology2009388116984[PubMed][Google Scholar]
  • 31. GrammatikosAPKyttarisVCKis-TothKFitzgeraldLMDevlinAFinnellMDA T cell gene expression panel for the diagnosis and monitoring of disease activity in patients with systemic lupus erythematosusClin Immunol20141502192200[PubMed][Google Scholar]
  • 32. WolkKWitteKWitteERafteryMKokolakisGPhilippSIL-29 is produced by T(H)17 cells and mediates the cutaneous antiviral competence in psoriasisSci Transl Med20135204204ra129[PubMed][Google Scholar]
  • 33. SwindellWRJohnstonAXingXVoorheesJJElderJTGudjonssonJEModulation of epidermal transcription circuits in psoriasis: new links between inflammation and hyperproliferationPLoS One2013811e79253[PubMed][Google Scholar]
  • 34. SunGYChuangDYZongYJiangJLeeJCGuZRole of Cytosolic Phospholipase A in Oxidative and Inflammatory Signaling Pathways in Different Cell Types in the Central Nervous SystemMolecular neurobiology2014501614[PubMed][Google Scholar]
  • 35. KumarSMcDonnellPCLehrRTierneyLTzimasMNGriswoldDEIdentification and initial characterization of four novel members of the interleukin-1 familyJ Biol Chem2000275141030814[PubMed][Google Scholar]
  • 36. BachmannMScheiermannPHardleLPfeilschifterJMuhlHIL-36gamma/IL-1F9, an innate T-bet target in myeloid cellsJ Biol Chem2012287504168496[PubMed][Google Scholar]
  • 37. DinarelloCAInterleukin-1 in the pathogenesis and treatment of inflammatory diseasesBlood201111714372032[PubMed][Google Scholar]
  • 38. VigneSPalmerGLamacchiaCMartinPTalabot-AyerDRodriguezEIL-36R ligands are potent regulators of dendritic and T cellsBlood201111822581323[PubMed][Google Scholar]
  • 39. DaErmeAMWilsmann-TheisDWagenpfeilJHolzelMFerring-SchmittSSternbergSIL-36gamma (IL-1F9) is a Biomarker for Psoriasis Skin LesionsJ Invest Dermatol20141354102532[PubMed][Google Scholar]
  • 40. SwindellWRStuartPESarkarMKVoorheesJJElderJTJohnstonACellular dissection of psoriasis for transcriptome analyses and the post-GWAS eraBMC Med Genet2014727[Google Scholar]
  • 41. GarlandaCDinarelloCAMantovaniAThe interleukin-1 family: back to the futureImmunity2013396100318[PubMed][Google Scholar]
  • 42. MartinDATowneJEKricorianGKlekotkaPGudjonssonJEKruegerJGThe emerging role of IL-17 in the pathogenesis of psoriasis: preclinical and clinical findingsJ Invest Dermatol201313311726[PubMed][Google Scholar]
Collaboration tool especially designed for Life Science professionals.Drag-and-drop any entity to your messages.