Skip to content

Advertisement

  • Research article
  • Open Access

Transcriptional landscape of psoriasis identifies the involvement of IL36 and IL36RN

  • 1, 2,
  • 3, 4Email author,
  • 3, 4,
  • 3, 4,
  • 1, 2 and
  • 1, 2
BMC Genomics201516:322

https://doi.org/10.1186/s12864-015-1508-2

  • Received: 28 September 2014
  • Accepted: 30 March 2015
  • Published:

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.

Keywords

  • Psoriasis
  • Transcriptome
  • Gene expression profiling
  • High-throughput nucleotide sequencing
  • Functional genomics

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

Patient

Age, years

Sex

AoO (duration), years

PASI

Nail involvement

PsA

Scalp involvement

Flexural ininvolvement

Treatment

P847

19

M

18 (0.5)

12.5

No

No

Yes

No

No

P652

25

M

13 (12)

18.4

Yes

No

Yes

Yes

No

P844

27

M

18 (9)

8.8

No

No

Yes

No

No

P848

29

M

26 (3)

3.7

No

No

Yes

No

No

P840

49

M

15 (34)

6.0

No

No

Yes

Yes

Calcipotriol + betamethasone

P845

52

M

22 (30)

14.0

Yes

No

Yes

Yes

No

P851

60

M

56 (4)

10.8

No

No

Yes

No

No

P853

28

F

28 (0.5)

23.3

No

No

Yes

No

No

P849

37

F

30 (7)

4.7

No

No

No

No

Topical steroid

P856

54

F

14 (40)

7.3

No

Yes

Yes

Yes

Topical steroid

P843

57

F

53 (4)

12.6

No

No

No

No

Topical steroid

P846

58

F

57 (1)

15.2

No

No

Yes

No

Topical 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
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)

Symbol

logFC

logCPM

PValue

FDR

Entrez Gene Name

OAS2

4.43

7.16

9.26E-70

2.14E-65

2'-5'-oligoadenylate synthetase 2

PLA2G4E

2.55

7.60

2.74E-68

3.16E-64

Phospholipase A2, group IVE

KLK9

5.20

4.01

2.20E-67

1.69E-63

Kallikrein-related peptidase 9

S100A12

9.00

4.41

1.25E-65

7.20E-62

S100 calcium binding protein A12

LCE3E

6.99

4.45

6.02E-65

2.78E-61

Late cornified envelope 3E

TGM1

3.42

6.51

4.94E-61

1.90E-57

Transglutaminase 1

OAS3

2.99

7.25

1.01E-59

3.33E-56

2'-5'-oligoadenylate synthetase 3

PARP9

2.20

7.11

2.80E-59

8.08E-56

Poly (ADP-ribose) polymerase, 9

CRABP2

2.49

5.65

1.99E-57

5.10E-54

Cellular retinoic acid binding protein 2

PLA2G4D

4.33

6.10

3.03E-56

7.01E-53

Phospholipase A2, group IVD

IL1F9

5.64

6.61

2.26E-55

4.74E-52

Interleukin 36, gamma

ALOX12B

3.01

7.30

7.12E-55

1.37E-51

Arachidonate 12-lipoxygenase

SAMD9

3.01

6.23

8.59E-55

1.53E-51

Sterile alpha motif domain

IL1F5

3.59

7.91

1.14E-53

1.87E-50

Interleukin 36 receptor antagonist

C10orf99

6.22

5.44

3.20E-52

4.93E-49

Chrom 10 open reading frame 99

DEFB4A

10.43

6.59

6.30E-52

9.09E-49

Defensin. beta 4A

AKR1B10

6.22

4.39

4.02E-51

5.46E-48

Aldo-keto reductase family 1

PAPL

4.01

4.60

1.56E-50

2.01E-47

Iron/zinc purple acid protein

GLTP

2.01

7.42

5.00E-50

6.08E-47

Glycolipid transfer protein

KDM6B

−1.36

11.51

1.48E-49

1.71E-46

Lysine (K)-specific demethylase 6B

KYNU

4.49

4.72

2.54E-48

2.79E-45

Kynureninase

RRM2

2.68

5.47

3.29E-48

3.46E-45

Ribonucleotide reductase M2

ZC3H12A

3.06

4.88

2.30E-47

2.31E-44

Zinc finger CCCH-type cont 12A

SDR9C7

2.70

5.00

6.01E-47

5.79E-44

Short chain dehydrogenase/reductase

HPSE

3.47

4.17

1.31E-46

1.21E-43

Heparanase

APOL6

2.23

6.52

1.91E-46

1.69E-43

Apolipoprotein L, 6

NIPAL4

1.39

9.13

6.90E-46

5.90E-43

NIPA-like domain containing 4

DMD

−1.81

7.46

9.69E-46

7.99E-43

Dystrophin

S100A7A

10.02

10.30

2.84E-45

2.26E-42

S100 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
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
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 Diapedesis

19.60

0.31

Agranulocyte Adhesion and Diapedesis

12.20

0.25

Role of IL-17A in Psoriasis

10.10

0.79

Atherosclerosis Signaling

9.75

0.25

Pathogenesis of Multiple Sclerosis

9.73

0.90

Differential Regulation of Cytokine Production in Intestinal Epithelial Cells by IL-17A and IL-17F

9.54

0.61

Differential Regulation of Cytokine Production in Macrophages and T Helper Cells by IL-17A and IL-17F

8.91

0.67

Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza

8.38

0.39

LXR/RXR Activation

6.29

0.20

Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses

6.13

0.22

Role of Cytokines in Mediating Communication between Immune Cells

5.89

0.31

Altered T Cell and B Cell Signaling in Rheumatoid Arthritis

5.85

0.22

T Helper Cell Differentiation

5.66

0.26

Hepatic Fibrosis / Hepatic Stellate Cell Activation

5.03

0.18

Communication between Innate and Adaptive Immune Cells

4.86

0.19

IL-10 Signaling

3.95

0.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 Psoriasis

7.64

0.50

γ-linolenate Biosynthesis II (Animals)

4.05

0.21

Granulocyte Adhesion and Diapedesis

3.80

0.08

Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza

3.62

0.15

Atherosclerosis Signaling

3.10

0.08

Altered T Cell and B Cell Signaling in Rheumatoid Arthritis

3.09

0.09

Agranulocyte Adhesion and Diapedesis

3.02

0.07

Role of Cytokines in Mediating Communication between Immune Cells

3.02

0.13

Differential Regulation of Cytokine Production in Macrophages and T Helper Cells by IL-17A and IL-17F

2.82

0.22

LPS/IL-1 Mediated Inhibition of RXR Function

2.72

0.06

Stearate Biosynthesis I (Animals)

2.51

0.10

Primary Immunodeficiency Signaling

2.44

0.09

Differential Regulation of Cytokine Production in Intestinal Epithelial Cells by IL-17A and IL-17F

2.41

0.17

IL-17 Signaling

2.34

0.09

Communication between Innate and Adaptive Immune Cells

2.33

0.07

Uracil Degradation II (Reductive)

2.32

0.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).

Abbreviations

LP: 

Lesional psoriatic skin

NLP: 

non-lesional psoriatic skin

C: 

normal skin of healthy controls

Declarations

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).

Authors’ Affiliations

(1)
Department of Dermatology, University of Tartu, Tartu, Estonia
(2)
Department of Dermatology, Tartu University Hospital, Tartu, Estonia
(3)
Department of Pathophysiology, Centre of Translational Medicine, University of Tartu, 19 Ravila Street, 50411 Tartu, Estonia
(4)
Department of Reproductive Biology, Estonian University of Life Sciences, Tartu, Estonia

References

  1. Nestle FO, Kaplan DH, Barker J. Psoriasis. N Engl J Med. 2009;361(5):496–509.View ArticlePubMedGoogle Scholar
  2. Lebwohl M. Psoriasis. Lancet. 2003;361(9364):1197–204.View ArticlePubMedGoogle Scholar
  3. Lowes MA, Suarez-Farinas M, Krueger JG. Immunology of psoriasis. Annu Rev Immunol. 2014;32:227–55.View ArticlePubMed CentralPubMedGoogle Scholar
  4. Liu Y, Helms C, Liao W, Zaba LC, Duan S, Gardner J, et al. A genome-wide association study of psoriasis and psoriatic arthritis identifies new disease loci. PLoS Genet. 2008;4(3):e1000041.View ArticlePubMed CentralPubMedGoogle Scholar
  5. Kõks S, Kingo K, Rätsep R, Karelson M, Silm H, Vasar E. Combined haplotype analysis of the interleukin-19 and −20 genes: relationship to plaque-type psoriasis. Genes Immun. 2004;5(8):662–7.View ArticlePubMedGoogle Scholar
  6. Kõks S, Kingo K, Vabrit K, Rätsep R, Karelson M, Silm H, et al. Possible relations between the polymorphisms of the cytokines IL-19, IL-20 and IL-24 and plaque-type psoriasis. Genes Immun. 2005;6(5):407–15.View ArticlePubMedGoogle Scholar
  7. Bowcock AM, Shannon W, Du F, Duncan J, Cao K, Aftergut K, et al. Insights into psoriasis and other inflammatory diseases from large-scale gene expression studies. Hum Mol Genet. 2001;10(17):1793–805.View ArticlePubMedGoogle Scholar
  8. Gudjonsson JE, Ding J, Johnston A, Tejasvi T, Guzman AM, Nair RP, et al. Assessment of the psoriatic transcriptome in a large sample: additional regulated genes and comparisons with in vitro models. J Invest Dermatol. 2010;130(7):1829–40.View ArticlePubMed CentralPubMedGoogle Scholar
  9. Swindell WR, Johnston A, Voorhees JJ, Elder JT, Gudjonsson JE. Dissecting the psoriasis transcriptome: inflammatory- and cytokine-driven gene expression in lesions from 163 patients. BMC Genomics. 2013;14:527.View ArticlePubMed CentralPubMedGoogle Scholar
  10. Gudjonsson JE, Ding J, Li X, Nair RP, Tejasvi T, Qin ZS, et al. Global gene expression analysis reveals evidence for decreased lipid biosynthesis and increased innate immunity in uninvolved psoriatic skin. J Invest Dermatol. 2009;129(12):2795–804.View ArticlePubMed CentralPubMedGoogle Scholar
  11. Suarez-Farinas M, Lowes MA, Zaba LC, Krueger JG. Evaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (GSEA). PLoS One. 2010;5(4):e10247.View ArticlePubMed CentralPubMedGoogle Scholar
  12. Li B, Tsoi LC, Swindell WR, Gudjonsson JE, Tejasvi T, Johnston A, et al. Transcriptome Analysis of Psoriasis in a Large Case–control Sample: RNA-Seq Provides Insights into Disease Mechanisms. J Invest Dermatol. 2014;134(7):1828–38.View ArticlePubMed CentralPubMedGoogle Scholar
  13. Jabbari A, Suarez-Farinas M, Dewell S, Krueger JG. Transcriptional profiling of psoriasis using RNA-seq reveals previously unidentified differentially expressed genes. J Invest Dermatol. 2012;132(1):246–9.View ArticlePubMed CentralPubMedGoogle Scholar
  14. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40(10):4288–97.View ArticlePubMed CentralPubMedGoogle Scholar
  15. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.View ArticlePubMed CentralPubMedGoogle Scholar
  16. Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A. 2003;100(16):9440–5.View ArticlePubMed CentralPubMedGoogle Scholar
  17. Lo YH, Torii K, Saito C, Furuhashi T, Maeda A, Morita A. Serum IL-22 correlates with psoriatic severity and serum IL-6 correlates with susceptibility to phototherapy. J Dermatol Sci. 2010;58(3):225–7.View ArticlePubMedGoogle Scholar
  18. Bonifati C, Solmone M, Trento E, Pietravalle M, Fazio M, Ameglio F. Serum interleukin-6 levels as an early marker of therapeutic response to UVB radiation and topical steroids in psoriatic patients. Int J Clin Lab Res. 1994;24(2):122–3.View ArticlePubMedGoogle Scholar
  19. Johnston A, Xing X, Guzman AM, Riblett M, Loyd CM, Ward NL, et al. IL-1F5, −F6, −F8, and -F9: a novel IL-1 family signaling system that is active in psoriasis and promotes keratinocyte antimicrobial peptide expression. J Immunol. 2011;186(4):2613–22.View ArticlePubMed CentralPubMedGoogle Scholar
  20. Russell CB, Rand H, Bigler J, Kerkof K, Timour M, Bautista E, et al. Gene expression profiles normalized in psoriatic skin by treatment with brodalumab, a human anti-IL-17 receptor monoclonal antibody. J Immunol. 2014;192(8):3828–36.View ArticlePubMedGoogle Scholar
  21. Ainali C, Valeyev N, Perera G, Williams A, Gudjonsson JE, Ouzounis CA, et al. Transcriptome classification reveals molecular subtypes in psoriasis. BMC Genomics. 2012;13:472.View ArticlePubMed CentralPubMedGoogle Scholar
  22. Zhou X, Krueger JG, Kao MC, Lee E, Du F, Menter A, et al. Novel mechanisms of T-cell and dendritic cell activation revealed by profiling of psoriasis on the 63,100-element oligonucleotide array. Physiol Genomics. 2003;13(1):69–78.View ArticlePubMedGoogle Scholar
  23. Reischl J, Schwenke S, Beekman JM, Mrowietz U, Sturzebecher S, Heubach JF. Increased expression of Wnt5a in psoriatic plaques. J Invest Dermatol. 2007;127(1):163–9.View ArticlePubMedGoogle Scholar
  24. Yao Y, Richman L, Morehouse C, de los Reyes M, Higgs BW, Boutrin A, et al. Type I interferon: potential therapeutic target for psoriasis? PLoS One. 2008;3(7), e2737.View ArticlePubMed CentralPubMedGoogle Scholar
  25. Zaba LC, Suarez-Farinas M, Fuentes-Duculan J, Nograles KE, Guttman-Yassky E, Cardinale I, et al. Effective treatment of psoriasis with etanercept is linked to suppression of IL-17 signaling, not immediate response TNF genes. J Allergy Clin Immunol. 2009;124(5):1022–1010 e1021-1395.Google Scholar
  26. Suarez-Farinas M, Li K, Fuentes-Duculan J, Hayden K, Brodmerkel C, Krueger JG. Expanding the psoriasis disease profile: interrogation of the skin and serum of patients with moderate-to-severe psoriasis. J Invest Dermatol. 2012;132(11):2552–64.View ArticlePubMed CentralPubMedGoogle Scholar
  27. Tian S, Krueger JG, Li K, Jabbari A, Brodmerkel C, Lowes MA, et al. Meta-analysis derived (MAD) transcriptome of psoriasis defines the "core" pathogenesis of disease. PLoS One. 2012;7(9):e44274.View ArticlePubMed CentralPubMedGoogle Scholar
  28. Joyce CE, Zhou X, Xia J, Ryan C, Thrash B, Menter A, et al. Deep sequencing of small RNAs from human skin reveals major alterations in the psoriasis miRNAome. Hum Mol Genet. 2011;20(20):4025–40.View ArticlePubMed CentralPubMedGoogle Scholar
  29. Swindell WR, Xing X, Voorhees JJ, Elder JT, Johnston A, Gudjonsson JE. Integrative RNA-seq and microarray data analysis reveals GC content and gene length biases in the psoriasis transcriptome. Physiol Genomics. 2014;46(15):533–46.View ArticlePubMedGoogle Scholar
  30. Kulka M, Calvo MS, Ngo DT, Wales SQ, Goswami BB. Activation of the 2-5OAS/RNase L pathway in CVB1 or HAV/18f infected FRhK-4 cells does not require induction of OAS1 or OAS2 expression. Virology. 2009;388(1):169–84.View ArticlePubMedGoogle Scholar
  31. Grammatikos AP, Kyttaris VC, Kis-Toth K, Fitzgerald LM, Devlin A, Finnell MD, et al. A T cell gene expression panel for the diagnosis and monitoring of disease activity in patients with systemic lupus erythematosus. Clin Immunol. 2014;150(2):192–200.View ArticlePubMed CentralPubMedGoogle Scholar
  32. Wolk K, Witte K, Witte E, Raftery M, Kokolakis G, Philipp S, et al. IL-29 is produced by T(H)17 cells and mediates the cutaneous antiviral competence in psoriasis. Sci Transl Med. 2013;5(204):204ra129.View ArticlePubMedGoogle Scholar
  33. Swindell WR, Johnston A, Xing X, Voorhees JJ, Elder JT, Gudjonsson JE. Modulation of epidermal transcription circuits in psoriasis: new links between inflammation and hyperproliferation. PLoS One. 2013;8(11):e79253.View ArticlePubMed CentralPubMedGoogle Scholar
  34. Sun GY, Chuang DY, Zong Y, Jiang J, Lee JC, Gu Z, et al. Role of Cytosolic Phospholipase A in Oxidative and Inflammatory Signaling Pathways in Different Cell Types in the Central Nervous System. Molecular neurobiology. 2014;50(1):6–14.View ArticlePubMedGoogle Scholar
  35. Kumar S, McDonnell PC, Lehr R, Tierney L, Tzimas MN, Griswold DE, et al. Identification and initial characterization of four novel members of the interleukin-1 family. J Biol Chem. 2000;275(14):10308–14.View ArticlePubMedGoogle Scholar
  36. Bachmann M, Scheiermann P, Hardle L, Pfeilschifter J, Muhl H. IL-36gamma/IL-1F9, an innate T-bet target in myeloid cells. J Biol Chem. 2012;287(50):41684–96.View ArticlePubMed CentralPubMedGoogle Scholar
  37. Dinarello CA. Interleukin-1 in the pathogenesis and treatment of inflammatory diseases. Blood. 2011;117(14):3720–32.View ArticlePubMed CentralPubMedGoogle Scholar
  38. Vigne S, Palmer G, Lamacchia C, Martin P, Talabot-Ayer D, Rodriguez E, et al. IL-36R ligands are potent regulators of dendritic and T cells. Blood. 2011;118(22):5813–23.View ArticlePubMedGoogle Scholar
  39. DaErme AM, Wilsmann-Theis D, Wagenpfeil J, Holzel M, Ferring-Schmitt S, Sternberg S, et al. IL-36gamma (IL-1F9) is a Biomarker for Psoriasis Skin Lesions. J Invest Dermatol. 2014;135(4):1025–32.View ArticleGoogle Scholar
  40. Swindell WR, Stuart PE, Sarkar MK, Voorhees JJ, Elder JT, Johnston A, et al. Cellular dissection of psoriasis for transcriptome analyses and the post-GWAS era. BMC Med Genet. 2014;7:27.Google Scholar
  41. Garlanda C, Dinarello CA, Mantovani A. The interleukin-1 family: back to the future. Immunity. 2013;39(6):1003–18.View ArticlePubMed CentralPubMedGoogle Scholar
  42. Martin DA, Towne JE, Kricorian G, Klekotka P, Gudjonsson JE, Krueger JG, et al. The emerging role of IL-17 in the pathogenesis of psoriasis: preclinical and clinical findings. J Invest Dermatol. 2013;133(1):17–26.View ArticlePubMed CentralPubMedGoogle Scholar

Copyright

© Keermann et al.; licensee BioMed Central. 2015

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Advertisement