Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.

To increase our understanding of psoriasis, we used high-throughput complementary DNA sequencing (RNA-seq) to assay the transcriptomes of lesional psoriatic and normal skin. We sequenced polyadenylated RNA-derived complementary DNAs from 92 psoriatic and 82 normal punch biopsies, generating an average of ∼38 million single-end 80-bp reads per sample. Comparison of 42 samples examined by both RNA-seq and microarray revealed marked differences in sensitivity, with transcripts identified only by RNA-seq having much lower expression than those also identified by microarray. RNA-seq identified many more differentially expressed transcripts enriched in immune system processes. Weighted gene coexpression network analysis (WGCNA) revealed multiple modules of coordinately expressed epidermal differentiation genes, overlapping significantly with genes regulated by the long noncoding RNA TINCR, its target gene, staufen-1 (STAU1), the p63 target gene ZNF750, and its target KLF4. Other coordinately expressed modules were enriched for lymphoid and/or myeloid signature transcripts and genes induced by IL-17 in keratinocytes. Dermally expressed genes were significantly downregulated in psoriatic biopsies, most likely because of expansion of the epidermal compartment. These results show the power of WGCNA to elucidate gene regulatory circuits in psoriasis, and emphasize the influence of tissue architecture in both differential expression and coexpression analysis.

Introduction

Psoriasis is a chronic inflammatory disease of the skin and joints, affecting ~0.2–2% of the world’s population (Gudjonsson and Elder, 2012; Nestle ). Microarray-based studies of the psoriatic transcriptome have revealed a large number of differentially expressed genes (DEGs) in lesional skin (Bowcock ; Gudjonsson ; Reischl ; Suarez-Farinas ; Yao ; Zaba ; Zhou ), both in case-control comparisons and in lesional and normal skin from the same individual. Moreover, one study of non-lesional psoriatic vs. normal skin identified over 200 DEGs (Gudjonsson ). Recently, RNA-seq has emerged as a powerful alternative to microarrays (Garber ; Liu ; Metzker, 2010; Ozsolak and Milos, 2011; Roy ; Wang ) that provides increased accuracy and precision (Roy ), identifies unmapped transcripts and alternative splicing, and markedly increases the dynamic range detected (Liu ; Wang ). To date there has been only one study using RNA-seq to analyze the coding psoriatic transcriptome (Jabbari ), examining 3 pairs of lesional and normal skin samples from psoriatic individuals, and one other study limited to the microRNA genome (Joyce ).

In this paper, we report findings from a large-scale RNA-seq analysis of skin biopsies from 174 individuals. A subset of these samples has been studied previously using microarrays (Gudjonsson ; Gudjonsson ), allowing for a robust comparison of the two technologies. In addition, the large sample size we consider permits new modes of data analysis, such as gene co-expression network analysis (Barabasi and Oltvai, 2004; De Las Rivas and Fontanillo, 2010; Horvath and Dong, 2008; Langfelder and Horvath, 2007, 2008; Yip and Horvath, 2007; Zhang and Horvath, 2005). Our results demonstrate the value of RNA-seq in capturing the full dynamic range of the psoriatic transcriptome, identify key transcription factors involved in psoriasis, and reveal gene co-expression networks illuminating the processes of keratinocyte differentiation, lipid biosynthesis, and the inflammatory interplay between myeloid cells, T-cells, and keratinocytes. We also show that expansion of the epidermal compartment provides an architectural explanation for the apparent down-regulation of dermally-derived transcripts in psoriatic lesions.

Results

Comparison of RNA-seq and microarray results

Using the Illumina GAII platform, we sequenced libraries constructed from cDNA prepared from polyadenylated RNA extracted from skin biopsies donated by 92 psoriatic patients and 82 normal individuals. From reads averaging 80 base pairs in length, we generated on average 38±10 (mean±SD) million reads per sample, yielding approximately 3 gigabases of cDNA sequence per sample. We aligned the reads to the NCBI build 37 reference genome using TopHat (Garber ), yielding an average mapping rate of 80 ± 6 % (mean ± SD). After performing quality control (see Methods), gene expression levels, represented as reads per kilobase per million mapped reads (RPKM) (Mortazavi ), were obtained for 21,099 genes annotated in the RefSeq database (Pruitt ). We identified 3,577 DEGs between lesional and normal skin (|log2FC| ≥ 1 and p<1×10−6, corresponding to family-wise error rate < 0.025), of which 1,049 were up-regulated and 2,528 were down-regulated in psoriasis (the full list of DEGs is provided in Table S1). We then performed analysis to investigate the potential splicing events based on differential usage of exons (Anders ). We identified 343 exons (from 292 genes) with differential usage between the psoriatic and normal skin conditions. Notably, only 95 of the genes were differentially expressed based on the gene-level analysis (Table S1).

In previous studies (Gudjonsson ; Gudjonsson ) we measured gene expression using Affymetrix HU133 Plus 2.0 microarrays on 122 skin samples (58 lesional and 64 normal). Among these, 42 samples overlap with the present experiment (20 lesional and 22 normal). We first examined the expression levels for the 42 samples assayed on both platforms. In general, RNA-seq and microarray measurements were consistent for intermediate- and high-abundance transcripts, whereas discrepancies become noticeable for low-abundance transcripts in both cases and controls (Figure S1). We then compared each gene’s fold change (FC) estimates derived using microarrays with those obtained by RNA-seq (Figure S2a). For medium and high abundance transcripts, FC estimates were relatively consistent between the two platforms (Figure S2, c and d). In contrast, for low-abundance transcripts (Figure S2b), RNA-seq had a wider range of FC estimates, in agreement with other studies (Bradford ; Marioni ).

We next compared our set of DEGs with our full microarray dataset (Gudjonsson ). Over 80% (ie 794 / 987) of the genes identified in the microarray analysis were also flagged by RNA-seq, but only 22% (794 / 3,577) of the genes identified by RNA-seq were flagged by microarray analysis (Figure S3). Again, DEGs detected only by RNA-seq tended to be low-abundance transcripts (Figure S4, panels a and b), whereas those identified by both platforms (panels c and d) or only using microarrays (panels e and f) tended to be intermediate- and high-abundance transcripts. We observed very similar patterns when we restricted the analysis to only those genes and samples that overlapped between the RNA-seq and microarray experiments (Figure S5). We did not observe any difference in the distributions of gene length, number of exons, and number of transcripts in genes identified as DEGs only by RNA-seq, only by microarrays or by both platforms (Figure S6). Using the same statistical criteria that we have used in this analysis, we reprocessed available microarray data in order to compare the number and percentage of overlap of DEGs between our datasets and other independent datasets (Suarez-Farinas ; Zaba ). These results are presented in Table S2).

Analysis of differentially-expressed genes (DEGs)

We performed pathway analyses to identify biological functions/pathways enriched among the DEGs. The most significantly enriched functions among the up-regulated genes included “inflammatory responses” (p=6×10−31), “cytokine-receptor interactions” (p=6×10−27), “cell division” (p=2×10−15), and “keratinization” (p=4×10−14). Using the cytokine-cytokine receptor interaction pathway as an example, 14 of the 17 genes having the highest differences in expression between psoriatic and normal skin were not identified as DEGs in a recent microarray based study of psoriatic skin (Gudjonsson ) (Table S3). We performed QRT-PCR experiments on 8 out of the 17 genes listed in Table S3, and validated that their expression levels are all significantly differentially expressed (Table S3). The most significantly enriched functions among the down-regulated genes included “calcium ion binding” (p=4×10−21), “homophilic cell adhesion” (p=5×10−20), and “muscle contraction” (p=1×10−17). Subsequent analysis revealed that many of these genes were of dermal origin and systematically under-expressed due to expansion of the epidermal compartment in psoriatic lesions (see below). However, the 16 DEGs showing enrichment for “keratinization” that were detected only by RNA-seq comprised a mixture of up-regulated and down-regulated genes, and only 4 of which were represented on the HU133 Plus 2.0 microarray platform. Table S4 shows the full list of significantly enriched functions and the genes contributing to each.

DEGs attributable to Th17 function (Ouyang ) were prominent in our sample, including IL12B, IL17A, IL17F, IL21, IL22, IL26, IL21R and IL23R. With the exception of IL26 and IL21R, increased expression of these genes has previously been described in lesional psoriasis (Caruso ; Lee ; Teunissen ). While their expression in psoriatic skin was quite low, their expression was often undetectable in normal skin (Table S2). Other notable immune-related genes that are significantly up-regulated included IFN-γ (IFNG, >1000-fold), IFN-γ (IFNE, >1000-fold), nitric oxide synthase (NOS2, 50-fold), IL-6 (IL6, 4.3-fold), and IL-24 (IL24, 4-fold); IL-34 is significantly down-regulated (IL34, 7-fold suppressed).

To examine transcriptional mechanisms underlying the differential expression of genes identified by RNA-seq, we performed a transcription factor analysis using Ingenuity Pathway Analysis software. The results revealed significant enrichment for the targets of 61 transcription factors (Table S5). Of these, 15 were predicted to be in an activated state in psoriatic skin, 12 were predicted to be in an inhibited state, and the remaining 34 yielded no prediction regarding activation. Among those predicted to be activated, the most significant (p < 1×10−10) enrichment involved targets of the pro-inflammatory transcription factors STAT3, NF-κB, C/EBPβ, and STAT1, whereas PPAR-γ and FOXA2 were predicted to be inhibited. PPAR-γ is a major regulator of genes involved in lipid metabolism, water transport, and cornified envelope formation during keratinocyte differentiation (Jiang ). Like PPAR-γ, FOXA2 regulates networks of genes controlling complex metabolic functions (Soccio ) and shares many target genes with PPAR-γ (Table S5).

Gene co-expression network analysis

We grouped transcripts with correlated expression levels into gene co-expression modules using the weighted gene co-expression network analysis (WGCNA) approach, which maximizes the network’s scale-free properties (Horvath and Dong, 2008; Langfelder and Horvath, 2008; Yip and Horvath, 2007; Zhang and Horvath, 2005). WGCNA identified 18 gene modules in normal skin and 27 in lesional psoriatic skin. Figure 1 shows the average RPKM, average Spearman correlations, and the proportions of significantly up- and down-regulated genes in each of the modules. We then tested each module for enrichment of specific biological functions or pathways (Table S6). Cluster annotations and FC values for each gene are provided in Table S7.

Figure 1

Averaged RPKM and proportion of DEGs for each module

Averaged RPKM among the genes (top panel), proportions of up- and down-regulated genes (middle panel), and the average Spearman correlation (bottom panel) for the co-expression gene modules constructed from the normal (normal; left panel) and psoriatic (psoriatic; right panel) skin samples, respectively.

Among the normal skin co-expression clusters, modules N10 and N15 had the highest proportion of up-regulated genes (each over 12%), and both were enriched for genes involved in keratinization (p=2×10−10) or located in the epidermal differentiation complex (EDC) (p=2×10−10)(Table S6). Among modules identified in psoriatic skin, modules P20 and P26 had the highest proportion of up-regulated genes (each over 70%). Module P20 was enriched for genes involved in inflammatory and immune responses (p=1×10−19 and 6×10−18, respectively). Over 80% of module P20 genes were up-regulated in lesional skin, averaging ~9-fold difference in expression. Moreover, the average Spearman correlations of expression levels for genes in module P20 was much higher for lesional skin (r = 0.45) than for normal skin (r = 0.07, Figure 1, Table S6). Module P26 was significantly enriched for genes involved in keratinization (p=3×10−20) and located in the EDC (p=7×10−19). Module P23 was not characterized by strong differential expression, but like modules N10, N15 and P26, it was enriched for genes found in the EDC (p=7×10−25) and/or involved in keratinization (p=9×10−25). Consistent with the functional annotation results, modules N10, N15, P23, and P26 are all significantly enriched with the top signature genes of keratinocytes and epidermis (Figure 2). As shown in Figure 3, these four modules also have the highest proportions of overlap with transcripts identified in previous studies as being expressed in the granular layer of the epidermis (Mattiuzzo ; Toulza )., and with genes regulated by the long non-coding RNA TINCR and its target gene, staufen-1 (STAU1), which have been shown to regulate terminal differentiation of the epidermis (Kretz ), as well as the p63 target gene ZNF750 and its target KLF4 (Sen ) (Figure 3).

Figure 2

Proportion of top signature genes in different cell populations for co-expression gene modules constructed from normal (left panel) and psoriatic (right panel) skin samples

From the top to bottom, respectively, the four panels show the proportions of signature genes for different skin cell and tissue compartments (keratinocytes, epidermis, dermis, and adipose tissue); "innate immunity" lymphocytes (gamma-delta T-cells, natural killer (NK) cells, and NK-T cells); "adaptive immunity" lymphocytes (CD4+ T-cells, CD8+ T-cells, and regulatory T-cells); myeloid-derived leukocytes (monocytes, macrophages, and dendritic cells). Asterisks (*) denote significant (ie p< 1 × 10−4) enrichment for the top 5% of cell signature genes in the corresponding module.

Figure 3

Overlap of gene co-expression modules with stratum granulosum and epidermal differentiation transcription factor target genes

The figures show the proportion of genes in each module overlapping with the genes identified in different independent studies. The top two panels depict overlap of each of the normal- and psoriasis- derived modules with stratum corneum genes identified by Toulza et al. (Toulza ), Mattiuzzo et al. (Mattiuzzo ), or the union of both studies. The bottom two panels depict overlap with targets of the epidermal differentiation regulators STAU1 and TINCR (Kretz ) and with the transcription factors ZNF50 and KLF4 (Sen ). The numbers in the parentheses after the names of the studies indicate the number of genes in the corresponding gene sets used in this analysis.

Normal skin-derived modules N6 and N18 were significantly enriched for dermal tissue signature transcripts (Figure 2) related to extracellular space as well as collagen and glycosaminoglycan binding (Table S6). Both of these modules had high proportions of down-regulated genes (Figure 1 and Table S6). The apparent under-representation of dermis-derived transcripts in psoriasis lesions was further assessed with the aid of an existing dataset derived from laser-capture microdissection (LCM) of lesional vs. nonlesional psoriatic skin (Mitsui ). Using this dataset, we identified two sets of transcripts that (a) were specific for either epidermis or dermis (i.e. expression of genes are significantly elevated as compared to other 23 cell types, see Supplementary Information for details), and (b) did not vary substantially in their expression in lesional vs. nonlesional psoriatic skin based on the LCM experiment. We then calculated the FC values (psoriasis vs. normal) for each of these transcripts in our RNA-seq dataset. As shown in Figure S7, the FC distribution of epidermis-specific transcripts was centered on unity (log2 FC = 0), as expected. However, the log2 fold change value of the dermis-specific transcripts was centered around −1. This significant under-representation of dermal transcripts (Wilcoxon rank-sum test p<2×10−16 when comparing dermis-specific genes’ FCs with those from epidermis-specific genes) can be explained by the dramatic expansion of the epidermis that is characteristic of psoriatic lesions (Figure S7e).

The co-expressed gene modules with the most significantly enriched functional annotations were modules N14 and P18 (Table S6), both of which were enriched with the function “keratin filament” (p=3×10−113 for N14; p=8×10−114 for P18). Indeed, 99 genes were common to both of these modules, and the majority (76/99) of them are hair keratins (Schweizer ) or keratin-associated proteins found mainly in hair follicles (Rogers ). Other hair follicle-specific genes mapped to these clusters include S100A3 (Kizawa ), DSG4 (Bazzi ), PADI3 (Dong ), and trichohyalin (TCHH) (Steinert ). Moreover, both modules showed significant (p < 1×10−4) overlap (32% for N14; 26% for P18) with genes identified as down-regulated in a paired microarray comparison of bald vs. normal human scalp (Garza ). Moreover, the 25 most highly down-regulated genes in that study included 18 of the same hair-specific keratins observed in clusters N14 and P18, including 9 hair keratins, 6 keratin-associated proteins, S100A3, and TCHH, all of which were down-regulated by 4-fold or more in bald scalp.

The sebaceous gland and the arrector pili muscle are anatomically joined to the hair follicle to form the pilosebaceous unit (Fujiwara ; Lever and Schaumburg-Lever, 1990). Consistent with this, module N16 was significantly enriched for transcripts engaged in muscle contraction (Table S6), and module N17 was significantly enriched for signature genes expressed in adipose tissue (fat) but not in dermis (Figure 2). Module N17 was also significantly enriched for PPAR signaling transcripts (p=1×10−11) (Table S6), which would be in keeping with the major role of PPAR-γ in adipocyte differentiation (Yessoufou and Wahli, 2010).

Three psoriatic modules (P13, P19, and P21) were significantly enriched (p<1×10−18) for immune-related functions, had high proportions of up-regulated genes, and exhibited higher Spearman correlations in expression in lesional than in normal skin (Table S6). They also contained the highest proportions of signature genes expressed in lymphocytes and/or myeloid cells (Figure 2). As shown in Figure 2, modules P13 and P19 were significantly enriched for the signatures of several types of lymphocytes as well as dendritic cells, whereas module P21 manifested significant overlap for myeloid but not lymphoid signature genes.

In an effort to determine to what extent specific cytokines present in psoriatic lesions might influence the behavior of lesional keratinocytes, we used the signature gene enrichment approach to query co-expression modules expressed by keratinocytes and/or epidermis for genes distinctively expressed by cultured human keratinocytes treated with various cytokines (Swindell ; Swindell ). Interestingly, three (N15, P20, and P26) out of the five modules manifesting a strong keratinocyte gene expression signature (N10, N15, P20, P23, and P26) were significantly enriched for signature genes induced by IL-17 treatment (Table 1). Module P20 was also enriched for signature genes up-regulated by IL1α. This module was also significantly enriched for genes regulated by the transcription factors STAT1, STAT3, and IRF3 (17 genes, p=8×10−11; 11 genes, p=1.4×10−9; and 14 genes, p = 2.7×10−9, respectively) (Table 1 and data not shown).

Table 1

Results of enrichment tests for the modules containing the highest proportions of signature genes in keratinocytes/epidermis. For the cytokine-induced experiments, the dosage (ng/ml) and treatment time are included in the brackets.

			Enrichment test‡
					Gene signatures in cytokine-inducedexperiments
Condition	ModuleLabel	Hub genes#	Functions	TranscriptionFactor	up	down
Normal	10	TMEM54; PKP3;LYPD3; ST14;LAD1	EDC; keratinization;3' UTR-mediatedtranslational regulation
	15	C1orf68;LY6G6C; LCE1B;NIPAL4;SLC46A2	keratinization		IL17A (200ng|24hr)	IL1a (10ng|24hr);IL36g (5Kng|24hr);IL22 (20ng|4d)
psoriatic	20	SERPINB4;SERPINB3;SLC6A14; TGFA;HEPHL1	inflammatory response;immune response;extracellular space	STAT1; STAT3;IRF3	IL17A (200ng|24hr);IL1a (10ng|24hr);IL1a (25ng|48hr)
	23	LCE1C; LCE2B;LCE2C; LCE2D;C1orf68	EDC; keratinization	HR		IL36g (5Kng|24hr);IL36RA (5Kng|24hr);IL36b (5Kng|24hr)
	26	CNFN; TREX2;SPRR2B;SPRR2D;SULT2B1	keratinization; EDC;cornified envelope	GLIS1	IL17A+TNF (10ng|24hr);IL17A+TNF (1ng|24hr);TNF (1ng|24hr)	IL13 (20ng|24hr)

Only the top 3 most significant results are shown.

Hub genes were defined as the genes having the highest similarity in the module, as assessed by the strength of correlation and the degree of shared neighbors based on the gene expression data, with other genes in the same module.

IL37 (Nold ) was the most down-regulated gene (FC = 0.10) in psoriatic module P23, which is strongly enriched for genes expressed in the stratum granulosum and regulated by TINCR and STAU1 (Figure 3). We performed QRT-PCR to measure IL37 mRNA and immunostaining to localize IL-37 protein with respect to the known keratinocyte differentiation marker, loricrin. As shown in Figure 4, IL37 mRNA was expressed at much lower levels in lesional psoriatic skin, compared to non-lesional psoriatic skin and normal skin. IL-37 protein was readily detectable and confined to the stratum granulosum of normal and non-lesional psoriatic skin, where it co-localizes with loricrin, whereas it was undetectable in lesional psoriatic skin.

Figure 4

IL-37 is co-expressed with loricrin in the granular layer of non-lesional skin but significantly reduced in lesional psoriasis skin

(a) Immunohistochemistry revealed that while IL-37 protein (DAB, brown) is expressed in the granular layer of the epidermis of normal appearing skin it is not detected in lesional psoriasis skin. Scale bar indicates 25 µm. (b) QRT-PCR results demonstrating down-regulation of IL-37 and loricrin mRNAs in lesional psoriatic skin compared with non-lesional psoriatic skin. (c) Fluorescence immunohistochemistry showing that IL-37 (green) co-localized with loricrin (red) in the granular layer of non-lesional skin and this was undetectable in lesional skin. DAPI counterstaining of nuclei is shown in blue. Scale bar indicates 200 µm. Statistical significance as assessed by 2-tailed t-test: ** p<0.01, *** p<0.001,.

Discussion

While not the first study to apply RNA-seq to the study of psoriasis (Jabbari ; Joyce ), the availability of a much larger sample provides a more comprehensive view of the processes influencing the psoriatic transcriptome. Many of the differentially-expressed genes that we have detected by RNA-seq were not represented or detected in previous microarray studies on psoriasis (Gudjonsson ; Reischl ; Suarez-Farinas ; Suarez-Farinas ; Yao ; Zaba ; Zhou ), including a recent meta-analysis (Tian ). Comparison of similarly powered experiments using the same statistical criteria revealed a substantial increase in the number of DEGs detected by RNA-seq, with little effect on the percentage of overlapping DEGs across studies (Table S2).

Inspired by earlier observations of coordinate gene expression in psoriasis (Elder ; Elder and Zhao, 2002; Hardas ; Kojima ), a major feature of this study is the application of gene co-expression network analysis to the RNA-seq data. While similar analyses have been fruitfully applied to the study of mouse skin carcinogenesis (Quigley ; Quigley ), to our knowledge this is its first application to human skin.

Links to immunological networks

WGCNA identified three co-expression modules in psoriatic skin that were significantly annotated for immune function and enriched for signature genes of myeloid cells (P21), or both myeloid and T-cells (P13 and P19) (Figure 2). Consistent with the known inflammatory / immune cell infiltration characteristic of psoriasis lesions, each of these modules manifested increased proportions of up-regulated genes in psoriatic skin (Figure 1, Table S6).

Among modules significantly enriched for epidermis or keratinocyte signature genes, modules N15, P20, and P26 were enriched for transcripts induced by IL-17 treatment of keratinocytes in vitro (Table 1). Module P20 was also enriched for genes regulated by activation of STAT1, STAT3, and IRF3 (Table 1). Together, these findings suggest that psoriatic keratinocytes are responding to a proinflammatory cytokine milieu in which IL-17 is upregulated (Table S1) and biologically available (Di Meglio ; Shi ; Zhang ). This is of clinical relevance because biologicals targeting IL17, its receptor, and both subunits of IL-23 are highly effective treatments for psoriasis (Lowes ), and IL-23 is a major signal supporting the expansion and survival of T-cells expressing IL-17 (Kryczek ). Because module P26 is highly upregulated in psoriasis lesions (Figure 1), keratinocyte differentiation-related (Figure 3) and IL-17 responsive (Table 1), it seems likely that IL-17 may promote innate host defense in the skin by signaling to the outer layers comprising the epidermal barrier.

Links to epidermal lineage commitment

During development, lineage commitment occurs when multipotent stem cells are specified into 3 lineages―the epidermis, the hair follicle, and the sebaceous gland (Hsu ). In this study, we have identified gene co-expression modules that can be attributed to each lineage. Modules N10, N15, P23 and P26 were highly expressed in epidermis (Figure 2). There was significant (p<1×10−4) overlap between genes mapping to these modules and genes expressed in the stratum granulosum (Mattiuzzo ; Toulza ) (Figure 3). HOPX, which encodes a homeodomain protein involved in regulating the balance between epidermal proliferation and differentiation (Obarzanek-Fojt ), maps to modules N15 and P23. HOPX is a target gene of the master epidermal development transcription factor p63 (TP63), as are several other downstream transcriptional regulators of epidermal differentiation including GRHL3, KLF4, and ZNF750 (Sen ; Zarnegar ), all of which mapped to module N10. Notably, each of these genes maps to reported psoriasis susceptibility loci (Birnbaum ; Tsoi ; Yang ). Modules N15 and P26 are significantly enriched (p < 1×10−4) with genes whose expression was altered by the knockdown of KLF4 and/or ZNF50 (Zarnegar ) (Figure 3). We also found that modules N15 and P23 were significantly (p< 1×10−4) enriched with genes whose expression was altered by deficiency of the long non-coding RNA TINCR or the STAU1 gene that is regulated by TINCR. Both of these genes are specifically involved in regulating the terminal phases of epidermal differentiation (Kretz ). The mapping of these transcription factors and their target genes to epidermal co-expression modules demonstrates that WGCNA is capable of detecting, in vivo, the activity of a p63-initiated and progressively focused transcription factor network that is engaged in the control of epidermal differentiation.

Regarding the hair follicle, the two modules most significantly enriched for function ("keratin filament", Table S6), N14 and P18, were strongly overlapping and contained many hair keratins (Schweizer ) and keratin-associated proteins that are specifically expressed in the hair follicle (Rogers ). Both clusters also contain LHX2, which encodes a hair follicle stem cell master transcription factor (Mardaryev ), and module N14 contains MSX2, another such factor. Both of these transcription factors, and most of the genes encoding hair keratins and keratin-associated proteins have also been identified as part of a hair follicle gene co-expression network in mouse skin (Quigley and Balmain, 2009; Quigley ). This conclusion is further reinforced by the pronounced down-regulation of many of the genes mapping to these modules in bald scalp (Garza ).

Regarding the sebaceous gland, module N17 was significantly enriched for signature genes expressed by adipose tissue (Plaisier ) (Figure 2). We suspect that module N17 is reflective of sebaceous glands because all visible subcutaneous fat was trimmed from our biopsies prior to analysis, and because this module was significantly enriched for PPAR-regulated genes (Table S6). PPAR-γ stimulates sebocyte lipid production (Trivedi ) and protects sebocytes from apoptosis (Schuster ). Intimate anatomic and molecular relationships exist between the hair follicle, the sebaceous gland, and the arrector pili muscle, and (Fujiwara ; Lever and Schaumburg-Lever, 1990), and a clear connection has been observed between muscle and hair follicle-related gene co-expression networks in mouse skin (Quigley ). Consistent with these observations, we found that module N16 was significantly enriched for genes involved in muscle contraction (Table S6).

Tissue architecture as a determinant of gene co-expression

Unlike cell cultures, tissue biopsies consist of multiple substructures comprised of various cell types, which may differ considerably in normal vs. pathological states (Swindell ; Swindell ). Upon annotation of the modules identified by WGCNA, we noticed that many of the down-regulated modules were enriched for processes involving the extracellular matrix, which in skin is mainly found in the dermis. Noting that a much greater fraction of the cellular mass of a psoriatic biopsy is epidermal in origin than in a normal skin biopsy (Figure S7e), we reasoned that a larger fraction of the RNA in a psoriatic biopsy would be epidermal, and hence the dermal fraction would be lower. To test this, we identified two sets of genes from a dataset obtained through LCM of epidermal vs. appendage-free dermal tissue from normal vs. uninvolved psoriatic skin (Mitsui ). These gene sets shared two properties, including (a) similar expression in LCM-dissected lesions and uninvolved skin, and (b) a trend towards dermis- or epidermis-specific expression. Querying our RNA-seq dataset with the resulting two gene lists, we found an approximately 50% reduction in dermis-specific transcripts, (Figure S7). These results suggest that fold-change estimates for dermis-derived transcripts are driven downward by epidermal expansion in psoriasis lesions. This process is independent of intracellular shifts in gene transcription, but may nonetheless by a key factor influencing gene expression in full-thickness skin biopsies. Such changes would be predicted to involve not only dermis-derived mRNAs but also those from other cell types, including hair follicles, sweat glands, and immune/inflammatory cells, and indeed this is what we observe for hair follicle-derived clusters N14 and P18 (Table S7). Our findings thus have general implications for transcriptome studies of inflammatory skin disease, or any context involving tissue samples that feature a complex mixture of minority and majority cell types.

In summary, these results bring unprecedented depth and resolution to the psoriatic transcriptome. They also illuminate an intricate gene co-expression network driven by not only by inflammatory cells and transcription factors, but also by disease-related differences in tissue architecture. Applied in tandem with histologically-targeted approaches such as microdissection-coupled transcript profiling, more detailed examination of coordinate gene expression offers great promise for an integrated understanding of the genetics, immunology, and cell biology of psoriasis in the coming years.

Methods

All subjects involved in this study provided written informed consent under a protocol adherent to the Helsinki Guidelines and approved by the Institutional Review Board of the University of Michigan Medical School. Libraries for RNA-seq were generated from polyadenylated skin punch biopsy RNA and sequenced at one library per lane on the Illumina Genome Analyzer IIx. Reads were aligned to the reference genome NCBI build 37 using TopHat (Garber ). and expression was normalized to the number of reads per kilobase per million mapped reads (RPKM). We used the Wilcoxon rank-sum test to identify differentially expressed genes. Significant differentially expressed genes were detected based p < 10−6 (corresponding to Family-wise error rate (FWER) < 0.025) with a fold change greater than 2. Functional annotation was performed using Gene Ontology (Ashburner ), KEGG (Kanehisa ), and Reactome (Matthews ). Transcription factor analyses were performed using Ingenuity Pathway Analysis software (www.ingenuity.com). Significance thresholds were determined by Bonferroni correction. Coordinate gene expression analysis was performed using the weighted gene co-expression network (WGCNA) package (Horvath and Dong, 2008), with normal and lesional psoriatic skin samples being analyzed separately. Gene expression signatures were identified as described (Swindell ; Swindell ). Additional details are provided in the Supplementary Information.