Analysis of long non-coding RNAs highlights tissue-specific expression patterns and epigenetic profiles in normal and psoriatic skin
- Lam C Tsoi1,
- Matthew K Iyer2,
- Philip E Stuart3,
- William R Swindell3,
- Johann E Gudjonsson3,
- Trilokraj Tejasvi3, 4,
- Mrinal K Sarkar3,
- Bingshan Li1, 5,
- Jun Ding1, 6,
- John J Voorhees3,
- Hyun M Kang1,
- Rajan P Nair3,
- Arul M Chinnaiyan2, 7, 8,
- Goncalo R Abecasis1Email author and
- James T Elder3, 4, 9Email author
© Tsoi et al.; licensee BioMed Central. 2015
Received: 9 September 2014
Accepted: 11 December 2014
Published: 30 January 2015
Although analysis pipelines have been developed to use RNA-seq to identify long non-coding RNAs (lncRNAs), inference of their biological and pathological relevance remains a challenge. As a result, most transcriptome studies of autoimmune disease have only assessed protein-coding transcripts.
We used RNA-seq data from 99 lesional psoriatic, 27 uninvolved psoriatic, and 90 normal skin biopsies, and applied computational approaches to identify and characterize expressed lncRNAs. We detect 2,942 previously annotated and 1,080 novel lncRNAs which are expected to be skin specific. Notably, over 40% of the novel lncRNAs are differentially expressed and the proportions of differentially expressed transcripts among protein-coding mRNAs and previously-annotated lncRNAs are lower in psoriasis lesions versus uninvolved or normal skin. We find that many lncRNAs, in particular those that are differentially expressed, are co-expressed with genes involved in immune related functions, and that novel lncRNAs are enriched for localization in the epidermal differentiation complex. We also identify distinct tissue-specific expression patterns and epigenetic profiles for novel lncRNAs, some of which are shown to be regulated by cytokine treatment in cultured human keratinocytes.
Together, our results implicate many lncRNAs in the immunopathogenesis of psoriasis, and our results provide a resource for lncRNA studies in other autoimmune diseases.
Long non-coding RNAs (lncRNAs) have received much attention in the past several years. Coincident with improved annotation of functional elements [1,2], it has been appreciated that a large portion of the genome is transcribed during the course of development, much of which represents lncRNA . LncRNAs resemble mRNAs because they are typically transcribed from active chromatin , polyadenylated, and capped; however, they do not direct protein synthesis . LncRNAs play important functional roles in epigenetic regulation, by forming networks of ribonucleoprotein complexes with chromatin regulators , and targeting their action to appropriate genomic regions both in cis and in trans [5,7,8]. Initially recognized for their role in X-chromosome inactivation , lncRNAs are increasingly being implicated in a variety of disease states, including susceptibility to infection , neurodegenerative diseases , and cancer [12-15].
Recently, discovery and analysis of non-coding RNAs have been enhanced by RNA-seq technology, and different pipelines have been developed to identify novel lncRNAs using RNA-seq data [1,13,16,17]. However, despite successful studies highlighting important roles of lncRNAs in different tissues and diseases [12,13,15,18], little is known about the roles of lncRNAs in human autoimmune diseases . Furthermore, biological inference of lncRNA function remains a challenging task, given their currently-limited annotation status and low expression levels .
Psoriasis is a chronic immune-mediated inflammatory and hyperproliferative disease of skin and joints, affecting around 2% of the population [20,21]. Recently, we  and others  have applied RNA-seq technology to the analysis of protein-coding genes in psoriatic lesions, compared to uninvolved skin from the same individual , or to normal skin . Aided by the availability of large numbers of samples, we identified multiple networks of coordinately expressed genes in normal versus psoriatic skin utilizing weighted gene co-expression network analysis . Moreover, we found strong enrichment for genes that are known targets of lncRNA-mediated terminal differentiation in both normal and psoriatic skin . In the skin, several lncRNAs have been shown to play key roles in both epidermal specification during development  and epidermal terminal differentiation . A lncRNA named PRINS (Psoriasis susceptibility-related RNA Gene Induced by Stress) has also found to be essential in the survival of keratinocytes under stress condition and may contribute to psoriasis susceptibility . These observations motivated us to consider psoriasis as an exemplar for the roles of lncRNAs in autoimmune disease.
To accomplish this, we first applied a computational approach  and stepwise filtering procedures to identify high-confidence lncRNAs expressed in the RNA-seq cohort. We then developed an analytical pipeline to globally characterize skin-expressed lncRNAs, with a particular focus on novel lncRNAs that were not previously annotated. We tailored our analysis to ask: (i) whether newly identified lncRNAs would have different expression behaviors compared to mRNAs and previously-annotated lncRNAs; and (ii) whether we could use existing biological information and data to infer the functional roles of the identified lncRNAs. Applying the aforementioned tools, we identified tissue-specific expression patterns and epigenetic profiles for novel lncRNAs, which are more pronounced than for previously-annotated (known) lncRNAs, with a significantly higher proportion of novel lncRNAs differentially expressed in lesional psoriatic skin. By examining patterns of co-expression with mRNAs of known function, we found strong enrichment for immune-related functions among all identified lncRNAs, particularly those that are differentially expressed in psoriatic skin. Overall, our study highlights the importance of lncRNAs in the pathogenesis of psoriasis, and provides a valuable resource for lncRNA studies in other autoimmune diseases.
Identification of lncRNAs in normal, uninvolved, and lesional psoriatic skin
As the identification of unannotated transcripts could be due to artifacts involving mappability or artifactual nomination of immature mRNA fragments, we applied several filtering steps to remove potential artifacts (Materials and methods; Figure 1). First, we identified transcripts with coverage in at least 5% (that is, ≥11) of the samples in our full dataset. Next, we removed false lncRNA predictions due to premature mRNA fragments by calculating, for each annotated transcript, the genomic distance from its closest annotated gene exon (Additional file 3: Figure S1) and removing unannotated transcripts which map within 2 kb of the nearby annotated exons. Next, we imposed mappability and length filters to remove unannotated transcripts in genomic regions of low complexity and with short lengths (<200 bp) . These filters are described in more detail in the Materials and methods and in Additional file 3: Figures S2-5. In subsequent analyses, we only considered transcripts with ≥1 read per sample on average. An expressed transcript identified in this study is thus defined as a transcript that passes the above quality control filters.
Transcripts remaining after application of various filtering steps
≥5% (11) samples
Distance (≥2 kb)
Length (≥200 bp)
≥216 mapped reads
To evaluate the efficiency of the lncRNA discovery pipeline we used a subset of the data (that is, the 174 samples described in ) to identify novel lncRNAs using the procedures described above, and asked whether the identified transcripts are expressed in the remaining 42 independent samples. We found that over 95% of the novel lncRNAs identified in these 174 samples were also expressed in the 42 independent samples (Additional file 3: Figure S9). The result illustrates the robustness of the pipeline and the filtering procedures in the identification of previously unannotated lncRNAs.
Differences in gene expression patterns between novel lncRNAs and annotated transcripts
Numbers and proportions (in percentage) of differentially expressed genes for different gene categories under three different comparisons: normal vs. lesional psoriatic skin (NN vs. PP), uninvolved vs. lesional psoriatic skin (PN vs. PP), and normal vs. uninvolved skin (NN vs. PN)
No. DEGs (%)
Protein coding (n = 14,011)
Antisense (n = 2,294)
Pseudogene (n = 1,476)
Annotated lncRNA (n = 2,942)
Total (n = 1,080)
Intronic (n = 196)
Intergenic (n = 840)
Interleaving (n = 15)
Encompassing (n = 29)
NN vs. PP
PN vs. PP
NN vs. PN
Previous studies have noted that lncRNAs can be classified as enhancer-associated (elncRNAs) or promoter-associated (plncRNAs) [36-38]. In an effort to understand the potential impact or effect of these elncRNAs/plncRNAs on neighboring genes, we first determined the candidate elncRNAs or plncRNAs using predicted enhancers/promoters in NHEK as a reference. The elncRNAs were identified as having close proximity (<5 kb) with enhancers but not with promoters. Conversely, the plncRNAs were identified as being in close proximity (<5 kb) to promoters but not to enhancers. This analysis identified 764 elncRNAs (483 annotated and 281 novel) and 369 plncRNAs (342 annotated and 27 novel), resulting in a substantially higher proportion of novel lncRNAs among the elncRNAs than among the plncRNAs (full list in Additional file 6). Additionally, the skin specificity index/FC value (in the PP vs. NN comparison) obtained from the elncRNAs is significantly correlated (by Spearman’s correlation test) with that obtained from the protein-coding genes neighboring the elncRNAs (GelncRNAs); similar results were obtained for the plncRNAs and the protein-coding genes neighboring the plncRNAs (GplncRNAs). Specifically, significance of the correlation in skin specificity index between elncRNAs and GelncRNAs was P = 1.5 × 10−3 (ρ = 0.4), and between plncRNAs and GplncRNAs was P = 5.3 × 10−5; ρ = 0.3); and the significance of correlation in FC in the PP/NN comparisons was P = 1.5 × 10−6 (ρ = 0.6) for elncRNAs vs. GelncRNAs and P = 3 × 10−8 (ρ = 0.4) for plncRNAs vs. GplncRNAs. In contrast, we did not identify any difference when comparing the GelncRNAs to GplncRNAs (P >0.05). Furthermore, we observed a significant (P = 2.9 × 10−4) negative correlation between the skin specificity index (T s ) of the annotated gene with its distance to the closest novel lncRNA (distance restricted to ≤1 Mb). Taken together, our results show that while both elncRNA and plncRNAs exhibit tissue specificity and differential expression profiles that are similar to their corresponding neighboring protein-coding genes, novel lncRNAs are more likely to be enhancer-associated than promoter-associated. These results further highlight the potential biological and functional roles of the novel lncRNAs, which constitute more than 37% of the elncRNAs we identified.
Functional characterization of the identified lncRNAs
To better understand the functional ramifications of the expressed lncRNAs in our dataset, we deployed an analytical pipeline assessing the relationship of the identified lncRNAs to: (i) known psoriasis susceptibility loci; (ii) co-expressed annotated mRNAs; and (iii) lncRNAs expressed by cytokine-stimulated keratinocytes. We first asked whether the expressed lncRNAs are enriched in regions of biological interest in the cutaneous context. Notably, two of the regions of highest lncRNA densities (Figure 2) were the Epidermal Differentiation Complex (EDC, chromosome 1q21.3, 150–155 Mb) [39,40] and the Major Histocompatibility Complex (MHC, chromosome 6p21.3, 26–34 Mb) , both of which contain psoriasis-associated genes . In the EDC, we identified 16 annotated and 12 novel lncRNAs, yielding a significant enrichment for novel lncRNAs mapping to the EDC when using annotated lncRNAs as background (P = 3 × 10−2, hypergeometric test). In contrast, for the MHC region, the enrichment of novel lncRNAs is not significant (20 annotated vs. 6 novel lncRNAs, P = 0.92) (see Discussion).
We next estimated the density of lncRNAs within known psoriasis susceptibility regions , and provided a comprehensive catalog of expressed lncRNAs in the psoriasis loci in Additional file 7 (Table S6). We identified 103 expressed lncRNAs (31 of them novel), and 26 of them were differentially expressed in PP vs. NN skin. Moreover, we found that the psoriasis locus 9q31.2, which has no expressed protein-coding gene, contains two expressed lncRNAs (see Discussion).
Enriched (FDR ≤ 0.1) inferred functions among all the differentially expressed lncRNAs (DE lncRNAs) and differentially expressed novel lncRNAs (DE novel lncRNAs) in psoriatic skin
No. of genes with the function
No. of DEGs with the function
DE lncRNAs (total)
Extracellular region part
Cellular lipid metabolic process
Fatty acid metabolic process
Lipid metabolic process
Monocarboxylic acid metabolic process
Lipid biosynthetic process
Regulation of inflammatory response
Cytokine receptor interaction
DE novel lncRNAs
Extracellular region part
Cellular lipid metabolic process
Lipid metabolic process
Biosynthesis of unsaturated fatty acid
Interleukin-17 (IL-17) and tumor necrosis factor (TNF) are key proinflammatory cytokines in psoriasis , and we have recently shown that the RNA-seq based psoriatic skin transcriptome manifests a significant IL-17 stimulation signature . To explore the effects of these key cytokines on lncRNA expression profiles, we utilized an RNA-seq dataset generated to characterize the transcriptome of keratinocytes stimulated for 24 h with IL-17 and TNF. We evaluated whether the cytokine-enhanced or cytokine-repressed lncRNAs would be enriched among lncRNAs that are up- or downregulated in the psoriatic skin (cytokine-enhanced and cytokine-repressed lncRNAs are shown in Additional file 10: Table S9). Consistent with the above findings showing enrichment of immune-related functions for the identified lncRNAs, the upregulated lncRNAs in psoriatic skin were significantly enriched for cytokine-stimulated lncRNAs induced by IL17 + TNF treatment of keratinocytes (P <3.3 × 10−2 after multiple testing correction, and observed-to-expected ratio >2; Additional file 11: Table S10). We did not identify any enrichment for downregulated lncRNAs among cytokine-repressed transcripts, whether annotated or novel.
Previous transcriptomic studies have emphasized the analysis of protein-coding transcripts, using mRNA expression levels to characterize the patterns and potential functional roles of their translated proteins . The development of next generation sequencing technology has greatly accelerated the discovery and characterization of a new class of biologically-significant RNA transcripts - the lncRNAs [46,47]. As a class, lncRNAs tend to be under less stringent evolutionary constraint, to be expressed at lower levels, and are preferentially enriched in the nucleus . Moreover, the structural rules governing lncRNA function are just now beginning to come to light .
Using a computational approach , we first enumerated known and novel lncRNAs in biopsies of normal (NN) skin from healthy controls, from lesional psoriatic (PP) skin, and in a subset of 27 affected individuals, from paired biopsies of uninvolved psoriatic (PN) skin. The ENCODE project identified 9,000 lncRNA gene in the human genome [1,3], and we could detect expression of approximately 40% of these in our skin samples. We further identified over 1,000 novel lncRNAs in our data, many of which were not well-expressed in other tissue types (Figure 4). Based on this finding, we would speculate that many more tissue-specific lncRNAs remain to be identified in other tissue/cell types.
A notable feature of the novel lncRNAs we identified is their tissue specificity. We assessed this parameter by comparison to lncRNAs we identified in data generated by the BodyMap 2.0 project  as well as in an independent skin RNA-seq sample  (Figure 4). In addition to protein-coding transcripts, antisense and pseudogene transcripts manifested substantial overlap between expression in skin and other tissues. While this overlap was somewhat less for known lncRNAs, these four transcript classes were clearly distinguished from the novel lncRNAs. Taken together with the observations that lncRNAs play crucial roles in both the development (ANCR)  and the terminal differentiation of skin (TINCR) , these results further suggest important, tissue-specific roles for the identified lncRNAs in skin development and differentiation. Indeed, seven lncRNAs (three of them novel) were shown to be highly correlated in their expression (that is, ρ ≥0.7) with genes in the differentiation-associated clusters identified in our previous study of protein-coding genes (annotated as N15 and P23 in ). To our knowledge, this is the first study to show that novel lncRNAs identified in a differentiated tissue behave substantially differently than do annotated lncRNAs in terms of tissue specificity in both the transcriptomic and epigenetic scales. More complete analysis of other tissues will be needed to determine whether this conclusion can be extended to tissue-specific versus more widely expressed lncRNAs.
In an effort to better understand the biological significance of the novel lncRNAs, we took a clue from a recent study comparing intergenic lncRNAs arising from enhancer-associated elements (elncRNAs) to those arising from promoter-associated elements (plncRNAs) in mouse erythroblasts . Utilizing a panel of nine cell types in which promoters and enhancers have been well-mapped in ENCODE [34,35], we found that novel lncRNAs mapped significantly closer to enhancer sites in the two ectodermally-derived cell lines (namely human mammary epithelial cells (HMEC) and normal human epidermal keratinocytes (NHEK), relative to the other seven non-ectodermally-derived lines. In contrast, there was no significant difference in the relative distance to promoter elements in these comparisons. Both enhancer-associated and promoter-associated lncRNAs were highly correlated with the expression of nearby protein-coding genes. Based on these findings, we would speculate that the novel lncRNAs identified in this study may participate in the control of tissue-specific gene expression.
We explored the potential functions of the identified lncRNAs by correlating their co-expression patterns across combined NN and PP samples with those of known, biologically-annotated mRNA transcripts. This analysis identified enriched functions consistent with those identified in previous transcriptome analyses of psoriasis that were focused on protein-coding genes [22,49-55]. Furthermore, the novel lncRNAs that are differentially expressed yielded higher observed to expected ratios for the enrichment of inferred immunological functions, compared to all of the differentially expressed lncRNAs as a whole (Table 3). When taken together with the higher percentage of differentially-expressed novel lncRNAs compared to differentially expressed known lncRNAs (Table 2), these observations suggest that the novel lncRNAs we have identified may exert important biological functions in psoriatic skin.
Compared to known lncRNAs, significantly more of the novel lncRNAs were differentially expressed in psoriasis (Table 2). While intriguing, this result needs to be interpreted with caution for several reasons, including potential differences in RNA preparation and purity across studies, tissue-dependent expression, as well as differential overall expression of known vs. novel lncRNAs. However, we evaluated the latter possibility and found no significant difference in the distribution of absolute expression levels compared to extent of differential expression in PP versus NN skin for either annotated or novel lncRNAs. When we imposed a minimum median RPKM (among cases or controls) of 0.1 (a value comparable with those used by previous studies [56,57]), we also observed similar proportions of differential expression as we observed in this study. Moreover, genes identified as differentially expressed in our initial set of 174 samples also tended to be differentially-expressed in an independent set of 42 samples. These results suggest the identified novel lncRNAs in this study are true positives with high confidence. Furthermore, using genes with high skin specificity (T s >0.4), we observed higher proportions of differential expression in psoriatic skin for protein-coding genes (46%) and annotated lncRNAs (38%), when comparing with results from Table 2: 17% for protein-coding genes and 24% for annotated lncRNAs. In concordance with the observation that novel lncRNA has high skin-specificity and proportion of differential expression, these findings suggest that skin specificity is an important factor contributing to the dysregulation of transcription in psoriatic skin.
In our previous study , we showed the decrease in expression of dermal-specific genes could be due to the relative decrease in the dermal compartment in psoriatic skin when compared to normal skin . We acknowledge the fact that difference in cellular compositions (for example, decrease in relative abundance of dermal cells or infiltration of immune cells) could play a role in yielding some differentially expressed genes when comparing the psoriatic and normal skin. The laser-capture microdissection experiments we used in our previous study to identify the epidermal- and dermal-specific genes are microarray-based. This limited our ability to evaluate whether the differentially expressed novel lncRNAs are cell intrinsic or due to the change in cellular proportions within skin tissue. However, our work has generated important hypotheses and questions regarding the study of transcriptomic architecture in complex tissues. Future studies measuring gene expression for specific cell types captured by fluorescence-activated cell sorting and/or laser capture microdissection will be able to facilitate the assessment and evaluate whether the differentially expressed lncRNAs in skin tissues are cell intrinsic or due to the change in cellular proportions.
One genomic region known to contain many genes that are selectively expressed in skin is the EDC [31,39,40,59]. Indeed, we found that the EDC was among one of the genomic regions of highest lncRNA densities (Figure 2), with significant enrichment for novel lncRNAs (12 novel out of 28 total, P = 3 × 10−2). In contrast, the other high lncRNA density region, the MHC, showed no significant enrichment (6 novel out of 26 total, P = 0.92). While the MHC does contain the corneodesmosin (CDSN) gene, which is relatively specifically expressed in skin and hair, it is not generally enriched in genes specifically expressed in skin. We acknowledge the fact that the MHC contains several hundred genes [1,41] and it is more challenging to identify novel lncRNAs in this area after imposing the ‘distance to known annotated genes’ filter (Figure 1). However, abandoning this filter would increase the identification of artefactual novel lncRNAs arising to immature transcripts from previously annotated genes.
As observed for many other complex diseases [60,61], the majority of genetic susceptibility loci identified for psoriasis fall into non-coding regions [42,62,63]. These findings challenge us to identify the non-coding elements in these regions which might be responsible for conferring susceptibility. One possible scenario is that genetic variants (expression quantitative trait loci, or eQTLs) mapping to the enhancer/promoter regions could play important roles in regulating gene expression levels in different tissues . Indeed, we have shown enrichment for psoriasis susceptibility signals in psoriasis eQTLs, relative to non-eQTLs . However, variations in lncRNAs should also be taken into account. The fact that we observed similar ratios of expressed lncRNAs (103 out of 4,022: 2.6%) and protein-coding mRNAs (450 out of 14,011: 3.2%) within the susceptibility loci suggest the identification and interpretation for causal elements for disease association should not be restricted to protein-coding transcripts. Indeed, single-nucleotide variations in lncRNAs (‘riboSNitches’) have been shown to map to eQTLs in disease susceptibility regions, suggesting that they may directly confer risk for complex traits , and are strong candidates for further investigation.
In conclusion, we have identified a substantial number of skin-specific lncRNAs in this study, and we have imputed potential immunological functions for them in the pathogenesis of psoriasis. Moreover, our results provide interesting potential clues into the mechanisms of tissue-specific gene regulation. As the roles of lncRNAs in other human autoimmune diseases have not yet been fully identified and understood, this analysis should provide valuable resource and information for the future studies.
Materials and methods
The preparation methods and quality control procedures used for the initial set of 174 RNA-seq samples (92 PP, 82 NN) have been described , and the same protocol was used for the additional 42 RNA-seq samples (7 PP, 8 NN, and 27 PN). In order to limit the variability of expression caused by treatment, we required a washout period prior to biopsy: at least 1 week for topical medications and 2 weeks for phototherapy/systematic medications. Informed consent was obtained from all subjects under protocols approved by the University of Michigan Institutional Review Board (HUM00037994) and adheres to the Declaration of Helsinki Principles.
Identification of unannotated transcripts
We used Tophat  (version 1.3.3) to perform alignment and Cufflink  version 2.1.1 to perform ab initio transcript assembly. We then estimated read counts for each identified genes using the ReadCount software . We then employed an enhanced version of a computational approach  which combined the information across the samples in the dataset to detect any unannotated transcripts using gene models from Ensembl version 74 as reference . The identified transcripts were classified in eight different categories (Additional file 2: Table S2). We used the presence of AT/GU splice site sequences to predict the strand orientation for transcripts with at least 2 exons.
We applied different filtering procedures to remove potential artifacts for novel transcripts identified. We first required that transcripts have coverage in at least 5% of the samples in our dataset (that is, ≥11 samples). To remove the unannotated transcripts that may be fragments of premature mRNA from annotated genes, we estimated for each annotated gene i the distance (d i ) to the exon of another annotated gene. We then used the median distance of annotated lncRNA to determine the distance threshold (<2 kb) for removing potential premature mRNA fragments. We acknowledge that this criterion will remove true positive novel lncRNAs (assuming the same distribution of d i for annotated and novel lncRNAs), but in this study we imposed stringent criteria to reduce false positive results and as shown in Table 1 this distance filter removes more than half of the unannotated transcripts in some categories of unannotated transcripts. Sequence reads mapped to regions of low mappability could generate potential false positive transcripts; we obtained the uniqueness (35-mers) and alignability (75-mers) tracks from the UCSC Genome Browser  and computed the gene-wise mappability measure (that is, uniqueness and alignability scores) using the bigWigSummary tool from the UCSC Genome Browser. We retained novel transcripts with score of at least 0.9 in both mappability measures; for comparison, 80% and 88% of annotated ncRNAs transcripts have scores of ≥0.9 for the uniqueness and alignability measures, respectively. Next, we applied a minimum transcript length criterion (≥200 bp) to remove short transcripts.
Evaluation of novel lncRNAs
We evaluated the coding potential of the identified novel lncRNAs using TransDecoder  and txCdsPredict from UCSC Genome Browser. For the prediction using txCdsPredict, a score greater than 800 was used as a criterion (90% predictive of protein coding genes ). Additional file 3: Figure S7 shows the percentage of genes predicted to be candidates of coding transcripts by different approaches. Only two of the identified novel lncRNAs were predicted to have coding potential by both approaches.
We used DESeq, which implements a read count model based on negative binomial distribution, to perform the expression normalization and differential expression analysis  for three different comparisons: (i) lesional psoriatic (PP) versus normal skin (NN); (ii) PP versus paired uninvolved (PN) skin from 27 psoriatic patients; and (iii) uninvolved (PN) versus NN skin. Significantly differentially expressed genes were declared to have FDR ≤0.1 and |log2 fold change| ≥1.
qRT-PCR analysis of selected genes was performed using six lesional skin samples from psoriasis patients (PP), six uninvolved skin from psoriasis patients (PN), and six normal skin from control subjects (NN). These samples were independent of the samples we used in the RNA-seq experiments. Skin biopsies were flash-frozen in liquid nitrogen and stored at 80°C. RNA extractions were performed using RNeasy columns (Qiagen, Cat # 74136). RNA quantity and quality were measured on an Agilent 2100 Bioanalyzer (Agilent Technologies), and only samples yielding intact 18S and 28S ribosomal RNA profiles were used. Reversed transcription was performed using High Capacity cDNA Transcription kit (Applied Biosystems, Cat # 4368813). Transcripts were quantified by SYBR green fluorescence (Applied Biosystems, Cat # 4367659) using 7300 Real-Time PCR system (Applied Biosystems). Relative expression was quantified using large ribosomal protein P0 (RPLP0) as an internal reference. In the qPCR process, primers for G36220, 5′-AGG ATG TTC CCC TGC TTT TT-3′ and 5′-CAC TCT TGC GAT GAA GTG ATG-3′; for G25746, 5′-CCC CTG AGA CAT TTC TTC CA-3′ and 5′-AGC CTT GGA GGG TTT CAA AT-3′; for G2608, 5′-GGC CTT ATC TTT TGC ACC TG-3′ and 5′-CAA CCA GCC AAA TTC CTG TT-3′; and for RPLP0, 5′- GCT GAT CCA TCT GCC TTT GT-3′ and 5′- AAG TTG GTT GCT TTT TGG TGA-3′. All custom primers were purchased from Sigma-Aldrich.
Evaluation of tissue specificity
We downloaded RNA-seq sequence data from two independent RNA-seq cohorts [23,33] from the NCBI Gene Expression Omnibus , and performed alignment by Tophat using the same parameters and arguments as we described above. The Body Map 2.0 project  consists of 16 different tissues, and the other study  consists of three pairs of PN/PP skin samples. We then measured the expression level in these RNA-seq datasets for each gene we identified in this study. The tissue specificity (T s ) of a gene in tissue s was calculated as the fraction of expression (in RPKM) relative to the sum of its expression in all 17 tissues. We averaged the RPKM values between the three PN samples in the Jabbari et al. study to estimate the skin’s RPKM values.
We downloaded the chromatin state segmentation [34,35] files for nine different cell lines (GM12878, H1-hESC, K562, HepG2, HUVEC, HMEC, HSMM, NHEK, and NHLF) from the UCSC Genome Browser, and retrieved the predicted strong enhancer and active promoter elements for each cell line. For each gene, we computed the distance to the closest enhancer/promoter from the starting position of the gene. To evaluate whether the enhancer/promoter regions are closer to novel lncRNAs in ectodermally-related cell types (NHEK and HMEC) than in other cell types, we computed the relative distance (D ecto /D average ), where D ecto is the closest distance to the enhancer (or promoter) in ectodermally-related cell type (NHEK or HMEC), and D average is the average closest distance to the enhancer (or promoter) in the other cell types. The means and standard error bars from Figure 5 were computed after removing the outlier distances (that is, 1.5 times the inter-quartile range).
Functional characterization of the identified lncRNAs
We first examined if any of the expressed transcripts are within the previously identified regions of psoriasis susceptibility loci . The associated regions were defined by ±500 kb intervals (±3 Mb for MHC) with respect to the best genome-wide significant signal.
Preparation and analysis of keratinocyte RNA-seq libraries
Normal human epidermal keratinocytes prepared from adult skin as described  were grown to post-confluence as described , prior to addition of recombinant human IL-17 and/or TNF (each at 10 ng/mL). After 24 h of treatment, total RNA was isolated using Qiagen RNeasy Minikits (Valencia, CA, USA) and RNA quality and quantity were assessed using an Agilent Bioanalyzer. Libraries for high throughput sequencing were prepared using the Illumina mRNA-Seq kit according to the manufacturer’s description (Illumina, San Diego, CA, USA) and libraries were sequenced on the Illumina Genome Analyzer IIx. The alignment and expression quantification procedures were identical to those described above. We performed differential expression analysis using DESeq. Among genes with a differential expression FDR ≤0.1, those with a FC >2 were declared as significantly enhanced and those with a FC <0.5 as significantly repressed.
The RNA-seq data used for this analysis are accessible through GSE63980 (superseries of GSE54456 and GSE63979).
Differentially expressed genes
Enhancer-associated long non-coding RNA
Long non-coding RNA
Promoter-associated long non-coding RNA
High-throughput cDNA sequencing
We thank the many volunteers who provided skin biopsies for this study. This research was supported by NIH R01 grants AR042742, AR050511, AR054966, AR062382, and AR065183 to JTE, as well as K08 AR060802 to JEG. GRA is supported by research grants R01HG007022 from the National Human Genome Research Institute and R01EY022005 from the National Eye Institute. JTE and TT are supported by the Ann Arbor Veterans Affairs Hospital, and JEG is supported by the Doris Duke Charitable Foundation Grant No. 2013106 and Taubman Medical Institute (as Kenneth and Frances Eisenberg Emerging Scholar). AMC is supported by the Howard Hughes Medical Institute, and is also supported as a Taubman Scholar. We also acknowledge generous support from the Dawn and Dudley Holmes Memorial Fund and the Babcock Endowment Fund to the Department of Dermatology at the University of Michigan.
- Bernstein BE, Birney E, Dunham I, Green ED, Gunter C, Snyder M. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489:57–74.View ArticleGoogle Scholar
- Esteller M. Non-coding RNAs in human disease. Nat Rev Genet. 2011;12:861–74.PubMedView ArticleGoogle Scholar
- Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, et al. The GENCODE v7 catalog of human long non-coding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 2012;22:1775–89.PubMed CentralPubMedView ArticleGoogle Scholar
- Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458:223–7.PubMed CentralPubMedView ArticleGoogle Scholar
- Mercer TR, Mattick JS. Structure and function of long non-coding RNAs in epigenetic regulation. Nat Struct Mol Biol. 2013;20:300–7.PubMedView ArticleGoogle Scholar
- Kretz M, Webster DE, Flockhart RJ, Lee CS, Zehnder A, Lopez-Pajares V, et al. Suppression of progenitor differentiation requires the long non-coding RNA ANCR. Genes Dev. 2012;26:338–43.PubMed CentralPubMedView ArticleGoogle Scholar
- Rinn JL, Chang HY. Genome regulation by long non-coding RNAs. Annu Rev Biochem. 2012;81:145–66.PubMedView ArticleGoogle Scholar
- Nakagawa S, Kageyama Y. Nuclear lncRNAs as epigenetic regulators–Beyond skepticism. Biochim Biophys Acta. 1839;2013:215–22.Google Scholar
- Lee JT, Bartolomei MS. X-inactivation, imprinting, and long non-coding RNAs in health and disease. Cell. 2013;152:1308–23.PubMedView ArticleGoogle Scholar
- Gomez JA, Wapinski OL, Yang YW, Bureau JF, Gopinath S, Monack DM, et al. The NeST long ncRNA controls microbial susceptibility and epigenetic activation of the interferon-gamma locus. Cell. 2013;152:743–54.PubMed CentralPubMedView ArticleGoogle Scholar
- Wapinski O, Chang HY. Long non-coding RNAs and human disease. Trends Cell Biol. 2011;21:354–61.PubMedView ArticleGoogle Scholar
- Maruyama R, Suzuki H. Long non-coding RNA involvement in cancer. BMB Rep. 2012;45:604–11.PubMed CentralPubMedView ArticleGoogle Scholar
- Prensner JR, Iyer MK, Balbin OA, Dhanasekaran SM, Cao Q, Brenner JC, et al. Transcriptome sequencing across a prostate cancer cohort identifies PCAT-1, an unannotated lincRNA implicated in disease progression. Nat Biotechnol. 2011;29:742–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Huarte M, Guttman M, Feldser D, Garber M, Koziol MJ, Kenzelmann-Broz D, et al. A large intergenic non-coding RNA induced by p53 mediates global gene repression in the p53 response. Cell. 2010;142:409–19.PubMed CentralPubMedView ArticleGoogle Scholar
- Prensner JR, Iyer MK, Sahu A, Asangani IA, Cao Q, Patel L, et al. The long non-coding RNA SChLAP1 promotes aggressive prostate cancer and antagonizes the SWI/SNF complex. Nat Genet. 2013;45:1392–8.PubMedView ArticleGoogle Scholar
- Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Integrative annotation of human large intergenic non-coding RNAs reveals global properties and specific subclasses. Genes Dev. 2011;25:1915–27.PubMed CentralPubMedView ArticleGoogle Scholar
- Weikard R, Hadlich F, Kuehn C. Identification of novel transcripts and non-coding RNAs in bovine skin by deep next generation sequencing. BMC Genomics. 2013;14:789.PubMed CentralPubMedView ArticleGoogle Scholar
- Kretz M, Siprashvili Z, Chu C, Webster DE, Zehnder A, Qu K, et al. Control of somatic tissue differentiation by the long non-coding RNA TINCR. Nature. 2013;493:231–5.PubMed CentralPubMedView ArticleGoogle Scholar
- van Roosbroeck K, Pollet J, Calin GA. miRNAs and long non-coding RNAs as biomarkers in human diseases. Expert Rev Mol Diagn. 2013;13:183–204.PubMedView ArticleGoogle Scholar
- Gudjonsson JE, Elder JT. Psoriasis. In: Goldsmith L, Katz SI, Gilchrest BA, Paller AS, Woolf K, Leffell DJ, editors. Dermatology in General Medicine, 1. 8th ed. New York: McGraw-Hill; 2012. p. 197–231.Google Scholar
- Perera GK, Di Meglio P, Nestle FO. Psoriasis. Annu Rev Pathol. 2012;7:385–422.PubMedView ArticleGoogle Scholar
- 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 Investig Dermatol. 2014;134:1828–38.PubMed CentralPubMedView ArticleGoogle Scholar
- Jabbari A, Suarez-Farinas M, Dewell S, Krueger JG. Transcriptional profiling of psoriasis using RNA-seq reveals previously unidentified differentially expressed genes. J Investig Dermatol. 2012;132:246–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Sonkoly E, Bata-Csorgo Z, Pivarcsi A, Polyanka H, Kenderessy-Szabo A, Molnar G, et al. Identification and characterization of a novel, psoriasis susceptibility-related non-coding RNA gene. PRINS J Biol Chem. 2005;280:24159–67.View ArticleGoogle Scholar
- Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11.PubMed CentralPubMedView ArticleGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.PubMed CentralPubMedView ArticleGoogle Scholar
- Kersey PJ, Allen JE, Christensen M, Davis P, Falin LJ, Grabmueller C, et al. Ensembl Genomes 2013: scaling up access to genome-wide data. Nucleic Acids Res. 2013;42:D546–52.PubMed CentralPubMedView ArticleGoogle Scholar
- Xie C, Yuan J, Li H, Li M, Zhao G, Bu D, et al. NONCODEv4: exploring the world of long non-coding RNA genes. Nucleic Acids Res. 2014;42:D98–103.PubMed CentralPubMedView ArticleGoogle Scholar
- Gonzalez-Porta M, Calvo M, Sammeth M, Guigo R. Estimation of alternative splicing variability in human populations. Genome Res. 2011;22:528–38.PubMedView ArticleGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.PubMed CentralPubMedView ArticleGoogle Scholar
- Zhao XP, Elder JT. Positional cloning of skin-specific genes from the human epidermal differentiation complex. Genomics. 1997;45:250–8.PubMedView ArticleGoogle Scholar
- 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 Investig Dermatol. 2009;129:2795–804.PubMed CentralPubMedView ArticleGoogle Scholar
- Farrell CM, O’Leary NA, Harte RA, Loveland JE, Wilming LG, Wallin C, et al. Current status and new features of the consensus coding sequence database. Nucleic Acids Res. 2014;42:D865–72.PubMed CentralPubMedView ArticleGoogle Scholar
- Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, et al. Mapping and analysis of chromatin state dynamics in nine human cell types. Nature. 2011;473:43–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Ernst J, Kellis M. Discovery and characterization of chromatin states for systematic annotation of the human genome. Nat Biotechnol. 2010;28:817–25.PubMed CentralPubMedView ArticleGoogle Scholar
- Wang D, Garcia-Bassets I, Benner C, Li W, Su X, Zhou Y, et al. Reprogramming transcription by distinct classes of enhancers functionally defined by eRNA. Nature. 2011;474:390–4.PubMed CentralPubMedView ArticleGoogle Scholar
- Ren B. Transcription: Enhancers make non-coding RNA. Nature. 2010;465:173–4.PubMedView ArticleGoogle Scholar
- Marques AC, Hughes J, Graham B, Kowalczyk MS, Higgs DR, Ponting CP. Chromatin signatures at transcriptional start sites separate two equally populated yet distinct classes of intergenic long non-coding RNAs. Genome Biol. 2013;14:R131.PubMed CentralPubMedView ArticleGoogle Scholar
- Kypriotou M, Huber M, Hohl D. The human epidermal differentiation complex: cornified envelope precursors, S100 proteins and the ‘fused genes’ family. Exp Dermatol. 2012;21:643–9.PubMedView ArticleGoogle Scholar
- Elder JT, Zhao X. Evidence for local control of gene expression in the epidermal differentiation complex. Exp Dermatol. 2002;11:406–12.PubMedView ArticleGoogle Scholar
- Trowsdale J, Knight JC. Major histocompatibility complex genomics and human disease. Annu Rev Genomics Hum Genet. 2013;14:301–23.PubMedView ArticleGoogle Scholar
- Tsoi LC, Spain SL, Knight J, Ellinghaus E, Stuart PE, Capon F, et al. Identification of 15 new psoriasis susceptibility loci highlights the role of innate immunity. Nat Genet. 2012;44:1341–8.PubMed CentralPubMedView ArticleGoogle Scholar
- Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40:D109–14.PubMed CentralPubMedView ArticleGoogle Scholar
- Lee HK, Hsu AK, Sajdak J, Qin J, Pavlidis P. Coexpression analysis of human genes across many microarray data sets. Genome Res. 2004;14:1085–94.PubMed CentralPubMedView ArticleGoogle Scholar
- Mattick JS. The genetic signatures of non-coding RNAs. PLoS Genet. 2009;5:e1000459.PubMed CentralPubMedView ArticleGoogle Scholar
- Ponting CP, Oliver PL, Reik W. Evolution and functions of long non-coding RNAs. Cell. 2009;136:629–41.PubMedView ArticleGoogle Scholar
- Wan Y, Qu K, Ouyang Z, Chang HY. Genome-wide mapping of RNA structure using nuclease digestion and high-throughput sequencing. Nat Protoc. 2013;8:849–69.PubMedView ArticleGoogle Scholar
- 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:1793–805.PubMedView ArticleGoogle Scholar
- 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 Investig Dermatol. 2010;130:1829–40.PubMed CentralPubMedView ArticleGoogle Scholar
- Reischl J, Schwenke S, Beekman JM, Mrowietz U, Sturzebecher S, Heubach JF. Increased expression of Wnt5a in psoriatic plaques. J Investig Dermatol. 2007;127:163–9.PubMedView ArticleGoogle Scholar
- 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:e10247.PubMed CentralPubMedView ArticleGoogle Scholar
- 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:e2737.PubMed CentralPubMedView ArticleGoogle Scholar
- 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:1022–10.e1-1395.PubMed CentralPubMedView ArticleGoogle Scholar
- 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:69–78.PubMedGoogle Scholar
- Toung JM, Morley M, Li M, Cheung VG. RNA-sequence analysis of human B-cells. Genome Res. 2011;21:991–8.PubMed CentralPubMedView ArticleGoogle Scholar
- Hackett NR, Butler MW, Shaykhiev R, Salit J, Omberg L, Rodriguez-Flores JL, et al. RNA-Seq quantification of the human small airway epithelium transcriptome. BMC Genomics. 2012;13:82.PubMed CentralPubMedView ArticleGoogle Scholar
- Quigley D. RNA-seq permits a closer look at normal skin and psoriasis gene networks. J Investig Dermatol. 2014;134:1789–91.PubMed CentralPubMedView ArticleGoogle Scholar
- Jackson B, Tilli CM, Hardman MJ, Avilion AA, MacLeod MC, Ashcroft GS, et al. Late cornified envelope family in differentiating epithelia–response to calcium and ultraviolet irradiation. J Investig Dermatol. 2005;124:1062–70.PubMedView ArticleGoogle Scholar
- Welter D, MacArthur J, Morales J, Burdett T, Hall P, Junkins H, et al. The NHGRI GWAS Catalog, a curated resource of SNP-trait associations. Nucleic Acids Res. 2013;42:D1001–6.PubMed CentralPubMedView ArticleGoogle Scholar
- Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature. 2014;506:376–81.PubMed CentralPubMedView ArticleGoogle Scholar
- Nair RP, Duffin KC, Helms C, Ding J, Stuart PE, Goldgar D, et al. Genome-wide scan reveals association of psoriasis with IL-23 and NF-kappaB pathways. Nat Genet. 2009;41:199–204.PubMed CentralPubMedView ArticleGoogle Scholar
- Stuart PE, Nair RP, Ellinghaus E, Ding J, Tejasvi T, Gudjonsson JE, et al. Genome-wide association analysis identifies three psoriasis susceptibility loci. Nat Genet. 2010;42:1000–4.PubMed CentralPubMedView ArticleGoogle Scholar
- Lonsdale J, Thomas J, Salvatore M, Phillips R, Lo E, Shad S, et al. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45:580–5.View ArticleGoogle Scholar
- Ding J, Gudjonsson JE, Liang L, Stuart PE, Li Y, Chen W, et al. Gene expression in skin and lymphoblastoid cells: Refined statistical method reveals extensive overlap in cis-eQTL signals. Am J Hum Genet. 2010;87:779–89.PubMed CentralPubMedView ArticleGoogle Scholar
- Center for Statistical Genetics. http://genome.sph.umich.edu/wiki/Bam_read_count. Accessed 22 Dec 2014.
- USCS Genome Bioinformatics. http://genome.ucsc.edu/index.html. Accessed 22 Dec 2014
- TranscriptDecoder. http://transdecoder.sourceforge.net/. Accessed 22 Dec 2014
- Gene Expression Omnibus. http://www.ncbi.nlm.nih.gov/geo. Accessed 22 Dec 2014
- Croft D, Mundo AF, Haw R, Milacic M, Weiser J, Wu G, et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 2013;42:D472–7.PubMed CentralPubMedView ArticleGoogle Scholar
- Stoll SW, Johnson JL, Bhasin A, Johnston A, Gudjonsson JE, Rittie L, et al. Metalloproteinase-mediated, context-dependent function of amphiregulin and HB-EGF in human keratinocytes and skin. J Investig Dermatol. 2010;130:295–304.PubMed CentralPubMedView ArticleGoogle Scholar
- Johnston A, Gudjonsson JE, Aphale A, Guzman AM, Stoll SW, Elder JT. EGFR and IL-1 signaling synergistically promote keratinocyte antimicrobial defenses in a differentiation-dependent manner. J Investig Dermatol. 2011;131:329–37.PubMed CentralPubMedView ArticleGoogle Scholar
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.