Advertisement

Landscape of Long Noncoding RNAs in Psoriatic and Healthy Skin

Open ArchivePublished:December 17, 2015DOI:https://doi.org/10.1016/j.jid.2015.12.009
      We used RNA sequencing to study and characterize the long noncoding RNA (lncRNA) transcriptome in lesional skin from psoriasis patients before (PP) and after treatment (PT) with adalimumab and in normal skin from healthy individuals (NN). To this end, we sequenced total RNA from 18 psoriasis patients and 16 healthy controls. We merged three lncRNA reference datasets to create a single combined reference of 67,157 lncRNA transcripts with no overlaps. We identified differential expression of 971 lncRNAs between PP and NN, 157 between PP and PT, and 377 between PT and NN. Using differentially expressed lncRNAs between PP and NN, we identified a molecular lncRNA signature that distinguishes psoriatic skin from healthy skin. Furthermore, we performed an unsupervised hierarchical analysis that revealed distinct clustering of PP samples from NN. A coding noncoding network analysis revealed a large network of highly correlated lncRNA and protein coding transcripts that provided insight into the potential functions of unannotated lncRNAs. To the best of our knowledge, this description of both polyadenylated as well as nonpolyadenylated lncRNA transcripts in psoriasis has not been previously reported. Our findings highlight the potential importance of lncRNAs in the biology of psoriasis and response to treatment.

      Abbreviations:

      CNC (coding noncoding), DE (differentially expressed), FPKM (fragments per kilobase per million reads mapped), lncRNA (long noncoding RNA), NN (normal skin from healthy controls), PP (psoriatic plaque skin before treatment), PT (psoriatic plaque skin post-treatment), qPCR (quantitative PCR), TNF-α (tumor necrosis factor-α)

      Introduction

      Psoriasis is a common, immune-mediated inflammatory disease of the skin that affects 1–3% of the world’s population. It is a complex, multifactorial disorder that is influenced by both genetic and environmental factors. Transcriptome analysis using RNA seq, an approach that uses deep sequencing technologies to sequence RNAs, has revealed approximately 3,500 protein coding genes that are differentially expressed (DE) in lesional skin (
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      ).
      The sequencing of the human genome showed that there are only 20,000–25,000 protein coding genes, which represent <2% of the total genomic sequence (
      International Human Genome Sequencing Consortium
      Finishing the euchromatic sequence of the human genome.
      ). Furthermore, it was shown that transcription is not limited to protein coding regions but instead is pervasive throughout the mammalian genome (
      Encode Project Consortium
      Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project.
      ), and many novel non-protein coding transcripts were identified. Recently, it was also shown that approximately 80% of the genome is transcribed from one or both strands (
      Encode Project Consortium
      An integrated encyclopedia of DNA elements in the human genome.
      ), which could result in a large number of overlapping transcripts that include thousands of long RNAs with little or no protein coding potential.
      Long noncoding RNAs (lncRNAs) are defined as sequences longer than 200 nucleotides and are among the most highly transcribed of the noncoding RNAs (
      • Mattick J.S.
      • Makunin I.V.
      Non-coding RNA.
      ). Although the function of many of these lncRNAs remains to be identified, an increasing number of targeted functional studies suggest that they have diverse and important biological roles, ranging from regulation of development and differentiation to regulation of epigenetic processes by guiding chromatin-modifying enzymes to their sites of action and/or by acting as scaffolds for chromosomal organization, and playing a role in RNA modification, evolution, and inheritance.
      However, little is known about lncRNAs in psoriasis.
      • 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 noncoding RNA gene, PRINS.
      identified and characterized a complementary DNA, PRINS (Psoriasis susceptibility related RNA gene Induced by Stress), with unknown function. This complementary DNA was predicted to function as a non-coding RNA based on an in silico approach. In a recently published study, 1,214 lncRNAs (709 annotated and 505 novel) were found to be DE in psoriatic plaque skin before treatment (PP) versus normal skin from healthy controls (NN) (

      Tsoi LC, Iyer MK, Stuart PE, Swindell WR, Gudjonsson JE, Tejasvi T, et al. Analysis of long non-coding RNAs highlights tissue-specific expression patterns and epigenetic profiles in normal and psoriatic skin. Genome Biol 2015;16:1–15.

      ). However, this study focused on polyadenylated lncRNAs and did not examine nonpolyadenylated lncRNAs.
      To decipher the role of lncRNAs in psoriasis, we first used a candidate approach to screen 90 common lncRNAs in PP (n = 9) and NN (n = 11) using real-time quantitative PCR (qPCR). After quality control and filtering, relative expression of 33 lncRNAs was analyzed using glyceraldehyde-3-phosphate dehydrogenase as a reference gene (see Supplementary Table S1 online). We identified ST7OT, LOC285194, and Car Intergenic 10 as showing significantly decreased expression in psoriasis patients (P < 0.05). Meg9 and NEAT1 showed a trend toward decreased expression, but the results were not statistically significant (P > 0.05, see Supplementary Figure S1 online).
      Encouraged by these initial findings, we decided to perform an in-depth analysis of expression of lncRNAs in psoriasis skin. For this, we used the sensitive and robust RNA-seq approach to comprehensively examine the expression of lncRNAs in psoriatic plaque skin before treatment (PP) and after treatment (PT) with a systemic tumor necrosis factor (TNF)-α inhibitor, as well as in the skin of healthy controls (NN). For this analysis, we used total RNA to detect both polyadenylated as well as nonpolyadenylated lncRNA. Here we report the identification of a large number of lncRNAs that are DE in psoriatic skin compared to healthy skin. We also observed differential expression of lncRNAs in PP versus PT psoriatic skin. We were able to validate our findings through qPCR and by reanalyzing recently published RNA-seq data on psoriasis skin. We further performed coexpression analysis of our lncRNA data with mRNA data to help infer the functions of unannotated lncRNAs.

      Results

      Differential expression of lncRNAs in pretreatment lesional skin in comparison to healthy controls

      To profile DE lncRNAs in PP and NN, we performed RNA seq on 18 psoriasis lesional samples and 16 healthy control samples using total RNA. We used an Illumina HiSeq 2500 (Illumina, San Diego, CA) for sequencing and obtained an average of 52 million 101 bp paired-end reads per sample. Using the RefSeq database for noncoding transcripts, we identified 263 DE lncRNAs with a q value or false discovery rate < 0.05. Of these, 109 were overexpressed and 154 were underexpressed in PP relative to NN. However, because the RefSeq database does not comprehensively annotate lncRNAs, we used two additional reference datasets, GENCODE v19 and
      • Hangauer M.J.
      • Vaughn I.W.
      • McManus M.T.
      Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs.
      , to annotate lncRNAs from our data. Using these reference databases, we created a combined reference set of 67,157 lncRNAs using the three individual reference sets (see Materials and Methods). Using the combined reference, we observed 62,826 expressed lncRNAs (fragments per kilobase per million reads mapped [FPKM] > 0.2) in our samples. In the PP versus NN comparison, 971 lncRNAs were DE (q value < 0.05). Of these 971 lncRNAs, 572 were overexpressed whereas 399 were underexpressed in PP. Figure 1 shows a heatmap of 971 DE lncRNA transcripts between PP and NN skin. This revealed distinct separation of PP versus NN samples based on lncRNA expression only, with only one outlier. The top 30 up-regulated and top 30 down-regulated lncRNAs observed using the combined reference are given in Supplementary Tables S2a (online) and S2b (online). A detailed list of all DE transcripts is provided in Supplementary Table S3 (online). Our results demonstrate that a large number of lncRNAs are DE in psoriatic skin compared to healthy skin and that these DE lncRNAs form a genomic signature that can distinguish psoriatic from healthy skin.
      Figure 1
      Figure 1Heatmap of differentially expressed transcripts between PP and NN skin. Differentially regulated transcripts between PP and NN skin (n = 971) are shown in a heatmap image. The green color indicates high expression levels, and the red color indicates low expression levels. The black bar above the heatmap indicates NN skin, and the purple bar above the heatmap indicates PP skin. Note that the number of up-regulated transcripts is more than the number of down-regulated transcripts in PP. lncRNA, long noncoding RNA; NN, normal skin from healthy controls; PP, psoriatic plaque skin before treatment.

      Validation of DE lncRNAs

      To validate our RNA-seq findings, we performed reverse transcriptase qPCR on RNA derived from 17 PP and 14 NN skin biopsies. We chose four DE RefSeq lncRNAs (TRHDE-AS1, CYP4Z2P, HINT1, and RPSAP58) in PP versus NN skin for validation. Our qPCR results showed that three of the four tested lncRNAs (TRHDE-AS1, CYP4Z2P, and HINT1) were DE in psoriasis lesional skin in the same direction as observed in our RNA-seq data (see Supplementary Figure S2 online).
      To further validate our findings, we performed reanalysis of a previously published psoriasis RNA-seq dataset (
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      ). Although Li et al. mostly reported mRNA expression of protein coding genes in their study, we reanalyzed their raw data, aligning reads to the human genome (hg19) and examining lncRNA differential expression of psoriasis and healthy skin using our combined lncRNA reference. We observed that of the top 30 overexpressed lncRNAs in our dataset (see Supplementary Table S4a online), 28 were present in the data of Li et al. Of these 28, 26 showed overexpression and 20 were statistically significant (P < 0.05). Moreover, 5 were overexpressed but no P value was calculated by the CuffDiff program because <10 reads were present. Of the top 30 underexpressed lncRNAs from our data (see Supplementary Table S4b online), 27 were present in the data from Li et al. All 27 of them showed underexpression, with 12 having P < 0.05 and 15 having no P value calculated because <10 reads were present.

      Unsupervised hierarchical clustering of lncRNAs distinguishes psoriatic skin from healthy skin

      Although heatmap analysis (Figure 1) showed that lncRNA expression is different in PP and NN skin and that lncRNAs show a distinct molecular signature in psoriasis, we further asked if we could observe similar separation of samples using the whole lncRNA dataset (DE or not). For this, we performed unsupervised hierarchical clustering of all the lncRNAs expressed (n = 61991) in PP and NN and observed a near-complete separation of PP samples from NN with some minor overlaps (Figure 2). Two controls (NN14, NN15) clustered separately from the other controls, and one of the PP samples (PP8) clustered with the NN samples. All other samples showed complete separation based on disease status. These results indicate that, as a whole, lncRNA expression in PP skin differs markedly from that of NN skin.
      Figure 2
      Figure 2Unsupervised hierarchical clustering of transcripts expressed in PP and NN skin. Unsupervised hierarchical clustering was performed on all the lncRNAs expressed (n = 61991) in PP and NN. NN, normal skin from healthy controls; PP, psoriatic plaque skin before treatment.

      Differential expression of lncRNAs in psoriatic skin before and after anti-TNF treatment

      We also compared differential expression of lncRNAs in 18 pairs of psoriatic skin before (PP) and after 1 month of treatment (PT) with adalimumab, a humanized monoclonal antibody against TNF-α. The number of DE lncRNAs between PP and PT samples was 157 for the combined lncRNA. We also assessed whether PT samples showed differential lncRNA expression compared to NN samples. A large number of lncRNAs were found to be DE, and the number of DE lncRNAs between PT and NN (377 transcripts) was higher than that between PP and PT (157 transcripts) (Figure 3a). A paired analysis of lncRNAs between PP and PT samples using EdgeR, a software package used to analyze differential expression of gene expression data, is shown in Supplementary Table S5 (online). Furthermore, we calculated an lncRNA molecular distance to health metric for PP and PT patients using gene expression data from NN patients as baseline expression. This approach involves computation of a score that represents the “molecular distance” of a given sample relative to a baseline (e.g., healthy controls; see Materials and Methods for details). We observed PT patients to be closer to baseline than PP patients (Figure 3b).
      Figure 3
      Figure 3Differential expression of lncRNAs in psoriatic skin before and after treatment. (a) Number of differentially expressed lncRNA transcripts in three comparisons: PP versus NN, PP versus PT, and PT versus NN. (b) lncRNA molecular distance to health analysis of PP and PT using NN as baseline. Molecular distance to health for PP and PT is shown. Expression levels of NN were used as baseline. Note that PP skin has higher molecular distance to health than PT. lncRNA, long noncoding RNA; NN, normal skin from healthy controls; PP, psoriatic plaque skin before treatment; PT, psoriatic plaque skin post-treatment.

      Correlation of lncRNA expression with mRNA expression

      To assess the biological functions of expressed lncRNAs, we first performed gene ontology enrichment analysis of the 971 DE lncRNAs. Of these 971 lncRNAs, only 188 were annotated in the gene ontology database. Analysis for molecular functions of these 188 lncRNAs revealed the following three gene ontology terms to be significantly enriched: “molecular functions,” “protein binding,” and “binding.” Because the majority of lncRNAs were not annotated, we used a different method to assess potential functions of these lncRNAs. Because previous studies have indicated that genes with similar coexpression patterns show similar functions, we used coding noncoding (CNC) network analysis to correlate lncRNA expression with that of protein coding genes in our dataset. For this analysis, we used GENCODE lncRNAs and compared their expression with RefSeq mRNA transcripts. After filtering (see Materials and Methods), we observed 933 lncRNAs to show high correlation (correlation coefficient ≥ 0.95 and P ≤ 0.05) with 1,143 protein coding genes. Of these 933 lncRNAs, 11 were DE in our dataset. A network of the candidate genes associated with those 11 lncRNAs is shown in Supplementary Figure S3a (online). FPKM, fold change, correlations, and q values of the lncRNA-mRNA pair shown in Supplementary Figure S3a are listed in Supplementary Table S6a (online). Gene ontology term enrichment analysis of the protein coding genes associated with those 11 lncRNAs revealed sensory perception of smell (P = 1.77E-06), sensory perception of chemical stimulus (P = 1.04E-05), and G-protein–coupled receptor signaling pathway (P = 5.58E-03) as some of the significant biological processes. Furthermore, we asked whether any of the top 100 DE protein coding genes were also present in our lncRNA-mRNA association matrix. We identified seven protein coding genes in our lncRNA-mRNA association matrix. These seven protein coding genes with their corresponding associated lncRNAs are shown in Supplementary Figure S3b (online). FPKM, fold change, correlations, and q values of the lncRNA-mRNA pair shown in Supplementary Figure S3b are listed in Supplementary Table S6b (online).

      Proximity of lncRNAs to genome-wide association study signals in psoriasis

      Next, we examined whether DE lncRNAs were located near known psoriasis genome-wide association study loci. For this we used 36 known and newly identified psoriasis loci as reported by
      • Tsoi L.C.
      • Spain S.L.
      • Knight J.
      • Ellinghaus E.
      • Stuart P.E.
      • Capon F.
      • et al.
      Identification of 15 new psoriasis susceptibility loci highlights the role of innate immunity.
      . We observed CARD14 lncRNA (location: chr17:78143790-78183130) to overlap with CARD14 (chr17:78,161,599-78,172,723), a gene involved in the activation of NF-κB and known to be genetically associated with both common and autosomal dominant forms of psoriasis (
      • Jordan C.T.
      • Cao L.
      • Roberson E.D.
      • Pierson K.C.
      • Yang C.F.
      • Joyce C.E.
      • et al.
      PSORS2 is due to mutations in CARD14.
      ). The CARD14 lncRNA can be distinguished from the CARD14 mRNA by the presence of an additional 5′ transcribed exon. It was observed that CARD14 lncRNA was overexpressed in psoriatic lesional skin compared to healthy control (false discovery rate = 0.001, 2.9-fold increase). Although we were not able to confirm this using qPCR, most likely because of very low expression, we were able to see the alignment of reads specific for CARD14 lncRNA compared to CARD14 mRNA in Integrated Genomics Viewer (IGV), a tool to visualize large genomic data sets, (
      • Robinson J.T.
      • Thorvaldsdóttir H.
      • Winckler W.
      • Guttman M.
      • Lander E.S.
      • Getz G.
      • et al.
      Integrative Genomics Viewer.
      ) for PP and NN skin (see Supplementary Figure S4 online). Additionally, we found that the CARD14 lncRNA was present in 84% of the psoriasis and healthy controls samples in the dataset of
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      , and that CARD14 lncRNA was significantly overexpressed in psoriatic lesional skin compared to healthy skin (false discovery rate = 0.003, 2.3-fold increase), thus replicating our results. Previous studies have shown that mutations in CARD14 genes are associated with susceptibility to psoriasis. Thus, we used the web-based program RNAsnp to predict the effect of all the reported CARD14 SNPs associated with psoriasis on CARD14 lncRNA secondary structure. The program predicted that SNP rs9902358 can cause significant changes in secondary structure of the CARD14 lncRNA (see Supplementary Figure S5 online).
      The other two DE lncRNAs that were found to be located closely to psoriasis genome-wide association study signal were LINC00302 and IL12RB2, located <50 kb to single nucleotide polymorphisms rs6677595 (LCE3B, LCE3C) and rs9988642 (IL23R), respectively.

      Discussion

      The pathogenesis of psoriasis has been linked to pathways affecting innate and adaptive immunity, as well as epidermal differentiation. Microarray studies have reported that a large number of genes are DE in psoriatic versus healthy control skin (
      • Bowcock A.M.
      • 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.
      ,
      • Haider A.S.
      • Peters S.B.
      • Kaporis H.
      • Cardinale I.
      • Fei J.
      • Ott J.
      • et al.
      Genomic analysis defines a cancer-specific gene expression signature for human squamous cell carcinoma and distinguishes malignant hyperproliferation from benign hyperplasia.
      ,
      • Oestreicher J.L.
      • Walters I.B.
      • Kikuchi T.
      • Gilleaudeau P.
      • Surette J.
      • Schwertschlag U.
      • et al.
      Molecular classification of psoriasis disease-associated genes through pharmacogenomic expression profiling.
      ,
      • Suarez-Farinas M.
      • Lowes M.A.
      • Zaba L.C.
      • Krueger J.G.
      Evaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (GSEA).
      ,
      • Swindell W.R.
      • Johnston A.
      • Voorhees J.J.
      • Elder J.T.
      • Gudjonsson J.E.
      Dissecting the psoriasis transcriptome: inflammatory- and cytokine-driven gene expression in lesions from 163 patients.
      ,
      • Zaba L.C.
      • Suarez-Farinas M.
      • Fuentes-Duculan J.
      • Nograles K.E.
      • 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.
      ,
      • Zhou X.
      • Krueger J.G.
      • Kao M.C.
      • 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.
      ). More recently RNA seq, a technique that involves direct sequencing RNA/complementary DNA using high-throughput next-generation sequencing technologies, has been used to analyze psoriasis transcriptome (
      • Jabbari A.
      • Suárez-Fariñas M.
      • Dewell S.
      • Krueger J.G.
      Transcriptional profiling of psoriasis using RNA-seq reveals previously unidentified differentially expressed genes.
      ,
      • Joyce C.E.
      • 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.
      ,
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      ). In one such study (
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      ), the authors performed RNA-seq analysis of skin biopsies from 174 individuals and observed that >80% of genes identified by microarray were also identified by RNA seq but only 22% of genes identified by RNA seq were identified by microarray. These studies show that RNA seq is capable of capturing the full dynamic range of the psoriasis transcriptome and can better assist in identifying key genes involved in pathogenesis of psoriasis. Furthermore, RNA seq can also help in finding novel transcripts, which cannot be identified using microarray.
      Previous studies on lncRNAs have indicated that lncRNAs can be polyadenylated as well as nonpolyadenylated. In view of this, we undertook this study to examine both polyadenylated and nonpolyadenylated lncRNA profiles in psoriatic skin using ribosomally depleted total RNA. Furthermore, there currently exist different lncRNA reference datasets with different annotations. Previous studies on lncRNA profiling in various diseases have mostly used one or two datasets. This practice limits the identification of all the lncRNAs; as a result, many lncRNAs remain unidentified. Here, we assembled and combined three different lncRNA reference sets: RefSeq, GENCODE, and
      • Hendriks A.G.
      • van der Velden H.M.
      • Wolberink E.A.
      • Seyger M.M.
      • Schalkwijk J.
      • Zeeuwen P.L.
      • et al.
      The effect of adalimumab on key-drivers in the pathogenesis of psoriasis.
      (see Materials and Methods), to create a more comprehensive reference dataset for annotation. Using this combined reference, we observed 971 lncRNAs to be DE in PP versus NN skin. This number is likely to be a conservative estimate of the true number of DE lncRNAs in psoriasis as we applied a heavy false discovery rate correction to our results based on >60,000 lncRNA transcripts. When we created a heatmap of those 971 transcripts, we observed a distinct molecular signature of psoriasis based on lncRNAs alone. Furthermore, hierarchical clustering analysis using the entire lncRNA data also revealed that the lncRNA expression profile of PP skin was distinct from that of NN skin.
      We validated our findings using qPCR and through reproduction of our findings using an independent RNA-seq dataset (
      • Li B.
      • Tsoi L.C.
      • Swindell W.R.
      • Gudjonsson J.E.
      • Tejasvi T.
      • Johnston A.
      • et al.
      Transcriptome analysis of psoriasis in a large case-control sample: RNA-seq provides insights into disease mechanisms.
      ). We observed significant congruence between our top DE lncRNAs and the same transcripts in their dataset, with the slight differences primarily due to lower transcript abundance in several of the lncRNAs examined in their data resulting in the absence of a P-value calculation.
      In an effort to understand if any of the lncRNAs observed in our study are related in expression with the protein coding genes, we generated a CNC network. For this we generated a correlation matrix of lncRNAs with that of protein coding transcripts. Of the 933 lncRNAs present in the matrix, we found 11 to be DE in PP versus NN skin. We observed these 11 lncRNAs to be coexpressed along with many protein coding genes. Functional analysis of the functions of those protein coding genes showed perception to smell and G-protein–coupled receptor signaling pathway as some of the significant processes.
      We also analyzed expression of lncRNAs in skin derived from patients treated with the TNF-α inhibitor adalimumab. A comparison of DE lncRNAs in PP versus PT and in PT versus NN revealed that the number of DE lncRNAs in PP versus NN was higher than that observed in PP versus PT. Furthermore, a molecular distance to health analysis revealed that molecular distance of lncRNAs expressed in PT was closer to healthy controls than PP. Because TNF-α plays an important role in the inflammatory process and has been found in elevated levels in psoriasis patients, the targeting and inhibition of TNF-α has become a key approach to treating psoriasis. It has been shown that treatment with adalimumab causes a decrease in expression of various cytokines and transcription factors, such as IL-17 and beta-defensins, in psoriatic skin (
      • Hendriks A.G.
      • van der Velden H.M.
      • Wolberink E.A.
      • Seyger M.M.
      • Schalkwijk J.
      • Zeeuwen P.L.
      • et al.
      The effect of adalimumab on key-drivers in the pathogenesis of psoriasis.
      ). Our results indicated that the lncRNA profile of posttreated skin was significantly different from that of pretreated skin and closer to that of normal skin.
      In this study, we identified three DE lncRNAs in proximity to known psoriasis susceptibility loci at CARD14, LCE3B/LCE3C, and IL23R. Whether these lncRNAs potentially regulate the function of these psoriasis genes is an intriguing question that will require additional investigation.
      Overall, this study identified a large number of lncRNAs present in human skin (>60,000), of which approximately 1,000 are DE in psoriasis. Many of these DE lncRNAs are novel, and future studies will be aimed at deciphering the function of these important biological mediators.

      Materials and Methods

      Patient enrollment and library preparation

      We recruited 18 adult subjects who had chronic plaque psoriasis with affected body surface area >10% and who were not taking systemic medications at the time of enrollment at the UCSF Psoriasis and Skin Treatment Center. The mean baseline psoriasis area and severity index of the patients was 14.6 (standard deviation 3.9). Punch biopsies (5 mm) were taken from the edge of a psoriatic plaque before treatment and after 1 month of treatment with adalimumab. The mean improvement in the psoriasis area and severity index over this time period was 53.1%. We obtained 16 normal skin samples from healthy control surgical discard specimens. We obtained written, informed consent for study participation from all the subjects under the approval of the local Institutional Review Boards. We extracted total RNA from biopsy samples using the Qiagen RNeasy mini kit (Qiagen, Valencia, CA). After ribosomal RNA depletion, we prepared RNA-seq libraries using ScriptSeq complete kits from Epicenter (Madison, WI).

      RNASeq workflow

      We obtained an average of 52.3 million 101 bp paired-end reads per sample (minimum 37.1 million reads, maximum 89.4 million reads). Quality was validated using FastQC (version 0.11.1). Quality and adapter trimming was performed with Trimmomatic (version 0.3). Reads were aligned to the human genome (hg19) using Tophat2 (version 2.0.9). An average of 34.8 million paired-end reads were aligned for each library (minimum 25.3 million, maximum 61.4 million). For differential gene expression, the RefSeq gene set (downloaded from the UCSC genome browser, January 2013), GENCODE lncRNA set (GENCODE version 19), and the
      • Hangauer M.J.
      • Vaughn I.W.
      • McManus M.T.
      Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs.
      lncRNA set (Dataset_S3) were used for annotation. To create a combined lncRNA reference, a script was used to extract RefSeq entries that were noncoding and >200 bp (parameters for lncRNAs) and then merge GENCODE transcripts, RefSeq lncRNAs transcripts, and
      • Hangauer M.J.
      • Vaughn I.W.
      • McManus M.T.
      Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs.
      transcripts that were nonoverlapping. After annotation, expression levels were first normalized by estimating the FPKM values for each gene as follows:
      FPKMi=Xil˜iN×109,


      where Xi is the count of gene i, li˜ is the effective transcript length of gene i, and N is the number of reads sequenced (
      • Trapnell C.
      • Hendrickson D.G.
      • Sauvageau M.
      • Goff L.
      • Rinn J.L.
      • Pachter L.
      Differential analysis of gene regulation at transcript resolution with RNA-seq.
      ). Cuffdiff (version 2.1.1), a program from the Cufflinks package, was used to calculate differential gene expression between each of the two conditions and was run with default parameters. Cuffdiff uses a simple T test as follows:
      T=E[log(y)]Var[log(y)],


      where y is the ratio of the normalized counts between conditions. The RNA-seq data have been deposited in NCBI’s Gene Expression Omnibus, and the GEO accession number is GSE74697 (available at http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74697).

      Validation of findings using Li et al. data (GSE54456)

      We obtained the .sra files from the study by Li et al. and converted each .sra file into a .fastq file using sratoolkit (version 2.4.2). Reads were aligned to the human genome (hg19) using Tophat2 (version 2.0.9). Samples were split into psoriasis and controls based on the provided sraruntable from the study. Cuffdiff (version 2.1.1), a program from the Cufflinks package, was used to calculate the differential gene expression between the two conditions and was run with default parameters. RefSeq (from January 2013) and the combined lncRNA dataset were used as references.

      CNC network analysis

      Cufflinks (version 2.0.2) was used to calculate the RefSeq gene expression. The FPKMs of protein coding genes were extracted using a Perl script. CNC analysis was done using R software (version 2.10.0). The results were filtered using the following criteria: correlation coefficient ≥ |0.95|and P ≤ 0.005. We performed CNC analysis on the PP and NN groups separately, and then filtered the candidate genes in the control group from the PP group. KEGG pathway information of the candidate protein coding genes was then obtained using Perl scripts. Finally, Cytoscape (version 3.1.1) was used to visualize the network of lncRNAs and genes.

      Molecular distance to health

      We performed molecular distance to health analysis (
      • Pankla R.
      • Buddhisa S.
      • Berry M.
      • Blankenship D.M.
      • Bancroft G.J.
      • Banchereau J.
      • et al.
      Genomic transcriptional profiling identifies a candidate blood biomarker signature for the diagnosis of septicemic melioidosis.
      ) on all the DE lncRNAs on PP versus NN comparison. Briefly, we first established the baseline by calculating average expression level and standard deviation of each gene for the healthy skin group. Next, we calculated the “distance” of an individual gene from the baseline, which is the difference in raw expression level from the baseline average of a gene for a given sample, and then we calculated the number of standard deviations from baseline levels that the difference in expression represents. In step 3 we applied filters to establish qualify genes. The qualifying genes must differ from the average baseline expression by at least 200 and 2 standard deviations. In step 4 we calculated a global distance from baseline, which is the total of the number of standard deviations for all qualifying genes.

      ORCID

      Conflict of Interest

      The authors state no conflict of interest.

      Acknowledgments

      This work was supported in part by National Institutes of Health Grants R01AR065174 and K08AR057763 to WL.

      Supplementary Material

      References

        • Bowcock A.M.
        • 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-1805
        • Encode Project Consortium
        Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project.
        Nature. 2007; 447: 799-816
        • Encode Project Consortium
        An integrated encyclopedia of DNA elements in the human genome.
        Nature. 2012; 489: 57-74
        • Haider A.S.
        • Peters S.B.
        • Kaporis H.
        • Cardinale I.
        • Fei J.
        • Ott J.
        • et al.
        Genomic analysis defines a cancer-specific gene expression signature for human squamous cell carcinoma and distinguishes malignant hyperproliferation from benign hyperplasia.
        J Invest Dermatol. 2006; 126: 869-881
        • Hangauer M.J.
        • Vaughn I.W.
        • McManus M.T.
        Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs.
        PLoS Genet. 2013; 9: e1003569
        • Hendriks A.G.
        • van der Velden H.M.
        • Wolberink E.A.
        • Seyger M.M.
        • Schalkwijk J.
        • Zeeuwen P.L.
        • et al.
        The effect of adalimumab on key-drivers in the pathogenesis of psoriasis.
        Br J Dermatol. 2013; 170: 571-580
        • International Human Genome Sequencing Consortium
        Finishing the euchromatic sequence of the human genome.
        Nature. 2004; 431: 931-945
        • Jabbari A.
        • Suárez-Fariñas M.
        • Dewell S.
        • Krueger J.G.
        Transcriptional profiling of psoriasis using RNA-seq reveals previously unidentified differentially expressed genes.
        J Invest Dermatol. 2012; 132: 246-249
        • Jordan C.T.
        • Cao L.
        • Roberson E.D.
        • Pierson K.C.
        • Yang C.F.
        • Joyce C.E.
        • et al.
        PSORS2 is due to mutations in CARD14.
        Am J Hum Genet. 2012; 90: 784-795
        • Joyce C.E.
        • 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: 4025-4040
        • Li B.
        • Tsoi L.C.
        • Swindell W.R.
        • Gudjonsson J.E.
        • 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: 1828-1838
        • Mattick J.S.
        • Makunin I.V.
        Non-coding RNA.
        Hum Mol Genet. 2006; 15: R17-29
        • Oestreicher J.L.
        • Walters I.B.
        • Kikuchi T.
        • Gilleaudeau P.
        • Surette J.
        • Schwertschlag U.
        • et al.
        Molecular classification of psoriasis disease-associated genes through pharmacogenomic expression profiling.
        Pharmacogenomics J. 2001; 1: 272-287
        • Pankla R.
        • Buddhisa S.
        • Berry M.
        • Blankenship D.M.
        • Bancroft G.J.
        • Banchereau J.
        • et al.
        Genomic transcriptional profiling identifies a candidate blood biomarker signature for the diagnosis of septicemic melioidosis.
        Genome Biol. 2009; 10: R127
        • Robinson J.T.
        • Thorvaldsdóttir H.
        • Winckler W.
        • Guttman M.
        • Lander E.S.
        • Getz G.
        • et al.
        Integrative Genomics Viewer.
        Nat Biotechnol. 2011; 29: 24-26
        • 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 noncoding RNA gene, PRINS.
        J Biol Chem. 2005; 280: 24159-24167
        • Suarez-Farinas M.
        • Lowes M.A.
        • Zaba L.C.
        • Krueger J.G.
        Evaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (GSEA).
        PLoS One. 2010; 5: e10247
        • Swindell W.R.
        • Johnston A.
        • Voorhees J.J.
        • Elder J.T.
        • Gudjonsson J.E.
        Dissecting the psoriasis transcriptome: inflammatory- and cytokine-driven gene expression in lesions from 163 patients.
        BMC Genomics. 2013; 14: 527
        • Trapnell C.
        • Hendrickson D.G.
        • Sauvageau M.
        • Goff L.
        • Rinn J.L.
        • Pachter L.
        Differential analysis of gene regulation at transcript resolution with RNA-seq.
        Nat Biotechnol. 2013; 31: 46-53
      1. Tsoi LC, Iyer MK, Stuart PE, Swindell WR, Gudjonsson JE, Tejasvi T, et al. Analysis of long non-coding RNAs highlights tissue-specific expression patterns and epigenetic profiles in normal and psoriatic skin. Genome Biol 2015;16:1–15.

        • Tsoi L.C.
        • Spain S.L.
        • Knight J.
        • Ellinghaus E.
        • Stuart P.E.
        • Capon F.
        • et al.
        Identification of 15 new psoriasis susceptibility loci highlights the role of innate immunity.
        Nat Genet. 2012; 44: 1341-1348
        • Zaba L.C.
        • Suarez-Farinas M.
        • Fuentes-Duculan J.
        • Nograles K.E.
        • 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-1395
        • Zhou X.
        • Krueger J.G.
        • Kao M.C.
        • 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