Transcriptome wide analysis of long non‐coding RNA‐associated ceRNA regulatory circuits in psoriasis

Abstract Long non‐coding RNAs (lncRNAs) play critical roles in regulating immune‐associated diseases and chronic inflammatory disorders. Here, we found that lncRNAs involve in the pathogenesis of psoriasis through integrative analysis of RNA‐seq data sets from a psoriasis cohort. Then, lncRNA‐protein‐coding genes (PCGs) co‐expression network analysis demonstrated that lncRNAs extensively interact with IFN‐γ signalling pathway‐associated genes. Further, we validated 3 lncRNAs associate with IFN‐γ signalling pathway activation upon IFN‐γ stimulated in HaCaT cells, and loss of function experiments indicate their functional roles in the activation of inflammatory cytokine genes. Additionally, microRNA target screening analysis showed that lncRNAs may regulate JAK/STAT pathway activity through complete endogenous RNA (ceRNA) mechanism. Further experimental validation of PRKCQ‐AS1/STAT1/miR‐545‐5p regulatory circuitry showed that lncRNAs regulate the expression of JAK/STAT signalling pathway genes through competing for miR‐545‐5p. In summary, our results demonstrated that dysregulation of lncRNA‐JAK/STAT pathway axis promotes the inflammation level in psoriasis and thus provide potential therapeutic targets for psoriasis treatments.

Therefore, the functional roles of lncRNAs in psoriasis require further investigations.
Interferons (IFNs) play an important role in innate immune responses through IFNγ signalling and activation of intracellular JAK/STAT signalling pathway. 10 Disruption of the JAK/STAT signalling pathway might lead to various autoimmune diseases and inflammatory diseases. 11,12 Recent studies showed that the JAK/STAT pathway might involve inflammatory and neoplastic skin diseases, like psoriasis, atopic dermatitis, vitiligo and melanoma. 11 Besides, the JAK/STAT signalling pathway is essential for the activation of a wide range of cytokines and growth factors, leading to critical cellular events, such as cell proliferation, differentiation and apoptosis. 13 Psoriasis is a prototypic Th17 cell-mediated autoimmune disease with high expression of IL-17, TNF, IFNγ, IL-22 and IL-23. 14 Numerous studies have shown that critical pathogenic mediators of psoriasis in the JAK/STAT signalling pathway. 8,15 For instance, Th1 cytokine IFNγ induces phosphorylation of STAT1 by Janus kinases (JAK). 13 IFNγ induced the expression of numerous genes in the skin contributing to chronic inflammation and implicated in the pathogenesis of psoriasis. 16 Additionally, STAT1 expression is elevated in lesion psoriatic skin, suggesting that the JAK/STAT signalling pathway is associated with psoriasis. 17 Furthermore, IL-23 produced by dendritic cells and macrophages promotes Th17 cell expansion and survival within psoriatic skin, engagement with its cognate receptors, resulting in activation of STAT3 and STAT4. 18 However, whether lncRNAs participate in the activation of the JAK/STAT pathway in the onset of psoriasis remains mostly unknown.
In this study, we found that lncRNAs involve in the pathogenesis of psoriasis through integrative analysis of RNA-seq data from a psoriasis cohort. Then, lncRNA-protein-coding genes (PCGs) co-expression network analysis demonstrated that lncRNAs extensively interact with IFNγ signalling pathway-associated genes. Further expression dynamic analysis and loss of function analysis demonstrated that 3 ln-cRNAs associate with the activation of JAK/STAT signalling pathway upon IFNγ stimulated in HaCaT cells, indicating their functional roles in promoting inflammatory cytokines production in psoriatic keratinocytes. Moreover, microRNA target screening analysis demonstrated that lncRNAs and JAK/STAT pathway-associated genes co-regulated by the same set of miRNAs, suggesting lncRNAs may regulate the JAK/STAT pathway through complete endogenous RNA (ceRNA) mechanism. Further experimental validation of PRKCQ-AS1/STAT1/ miR-545-5p regulatory circuitry showed that lncRNA might regulate the activity of the JAK/STAT signal pathway through the ceRNA mechanism. In conclusion, our results demonstrated the functional mechanism of lncRNAs in psoriasis, which will benefit our further understanding of the pathogenesis of various skin diseases.

| Data collection
RNA sequencing (RNA-seq) and small RNA sequencing (small RNAseq) data sets derived from psoriasis and healthy cohort were downloaded from the NCBI Gene Expression Omnibus database (GEO) 19 with accession number GSE54456, GSE63979 (RNA-seq) 4,20 and GSE31037 (small RNA-seq). 21 Briefly, we obtained 189 RNA-seq data sets from 99 psoriatic and 90 healthy punch biopsies, while 44 small RNA-seq data sets from 24 psoriasis biopsies and 20 healthy controls. Raw sequencing reads were extracted from SRA files by SRAtoolkits. 22 All data sets used in this study were listed in Table S1.

| RNA-seq data analysis
We adopt a transcriptome analysis pipeline established in our previous studies. 23,24 Briefly, low-quality sequencing reads were filtered using adapter Removal 25 if read length <50 bps and quality value <3. 26 The remaining reads were then aligned to the human reference genome (hg19) using STAR aligner. 27 We then calculated gene expression levels of all of the transcripts as fragments per kilobase of exon per million fragments mapped (FPKM) by cufflinks. 28 Differentially expressed genes were detected by comparing psoriasis and healthy group using DEseq2. 29 Finally, genes were identified as differentially expressed genes using Benjamini-Hochberg adjusted P-value < .05 and | Log 2 (Foldchange) | ≥ 2 as the cut-off.

| Small RNA-seq data analysis
Small RNA-seq from psoriasis and healthy groups was analysed using miRDeep2. 30 Briefly, low-quality reads were eliminated using Cutadapt (version 2.10) 31 ; the filtered reads were then aligned to miRbase and the human reference genome (hg19). The expression level of each microRNA is quantified as Tag Per Millions reads (TPM).
The differentially expressed small RNA and microRNAs were identified if | Log 2 (Foldchange) | of TPM is higher than 2.0 and the adjusted P-values are less than P < .05. 32

| Construction of the lncRNA-PCG coexpression network
A gene expression matrix was constructed by integrating gene expression profiles from 189 RNA-seq data sets of 99 psoriatic and 90 healthy punch biopsies. We then applied WGCNA to construct gene co-expression network using the established gene expression matrix as input. 33 The weak co-expression gene pairs were eliminated with adjacency threshold higher than 0.02 and interacting with at least one lncRNA molecule, resulting in a reliable lncRNA-PCG coexpression network.

| Identification of ceRNA pairs
The differentially expressed lncRNAs, protein-coding genes (PCGs) and miRNAs were used to establish ceRNA interacting pairs. Briefly, we first screen potential miRNA target binding sites at 3' UTR of

| Functional enrichment analysis
The differentially expressed genes were subjected to DAVID functional annotation pipeline for gene function enrichment analysis. 39 KEGG pathways 40 were identified if the adjusted P-value was <.05.
The pre-ranked gene list with Log2 (Foldchange) was analysed by GSEA package against MSigDB. 41

| Statistical analysis
Data were analysed using GraphPad Prism (version 8; GraphPad Software). Data were represented as mean ± SD. All tests were two sided, and P < .05 was considered statistically significant.   Table S2.

| Identification of differentially expressed lncRNAs in psoriasis
Numerous studies have been shown that thousands of lncRNAs are differentially expressed in various disease conditions, 42 but the biological function of most lncRNAs in skin diseases remains unknown. To gain more insights into the functional role of lncRNA in psoriasis, we analysed a set of RNA-seq data derived from a psoriasis cohort including 90 normal skin and 99 psoriatic skin. 4 Then, we applied a customized bioinformatics pipeline to identify differentially expressed genes (DEGs) and retrieve differentially expressed lncRNAs (DE-lncRNAs) ( Figure 1A). As a result, 1319 up-regulated genes and 2173 down-regulated genes were observed in psoriatic skin compared with normal skin (Table S3). Then, we carried out principal component analyses (PCA) among normal and psoriatic skin using all expressed lncRNAs in healthy and psoriatic skin.
Expectedly, PCA analysis showed that the psoriasis group is clearly separated from the normal group, suggesting that the expression of lncRNA is significantly different between the healthy and psoriatic groups ( Figure 1B

| Long non-coding RNAs may involve in IFNγ signaling pathway in psoriasis
Although extensive long non-coding RNAs (lncRNAs) are elevated in psoriatic skin compared with normal skin, how those lncRNAs act their functional role still warrants further investigation. 42 Towards this end, we sought to establish a lncRNA-PCG interaction network including 503 lncRNAs and 12,263 PCGs in psoriasis group using WGCNA. Then, a subnet network with 50 DE-lncRNAs and 480 PCGs was extracted (Table S4)  . All experiments were performed in triplicate and the average expression level was reported. Significant differences were determined using the Mann-Whitney unpaired two tailed t test. (*indicates significance, P < .05; **indicates significance, P < .01; ***indicates significance, P < .001) stimulation have been used to mimic the innate immune response in the progression of numerous skin diseases. 44 We sought to validate the expression and functional role of lncRNA in RNA-seq data set and external clinical samples. As a result, we found PRKCQ-AS1, SH3PXD2A-AS1 and CERNA2 are significantly up-regulated in psoriasis skin compared to healthy skin ( Figure 3A-C). Expectedly, the expression of PRKCQ-AS1 and SH3PXD2A-AS1 was markedly elevated in psoriatic skin tissues (n = 10) compared with non-psoriatic skin tissues (n = 10), except CERNA2 ( Figure 3D-F). We observed the up-regulation of CERNA2 in psoriatic skin compared with normal skin but not significantly due to the heterogenous of clinical samples ( Figure 3F

| LncRNAs may regulate JAK/STAT pathway through ceRNA mechanism
Previous studies have been demonstrated that lncRNAs act as ceR-NAs to fine-tune the expression of target genes through completing for shared miRNAs in diverse biological processes. 45,46 Thus, we hypothesized that lncRNAs might participate in the regulation of inflammation through the ceRNA mechanism. To this end, we analysed small RNA-seq data sets from 24 psoriatic and 20 healthy skin. Interestingly, we found that 362 miRNAs were differentially expressed in psoriatic skin compared with healthy skin, including 205 up-regulated and 157 down-regulated miRNAs (Table S6).
Then the differentially expressed microRNAs were identified with the P-values less than 0.05 ( Figure 4A). We then established putative miRNA-lncRNAs, miRNA-PCGs interactions using miranda, and further KEGG enrichment analysis showed that the predicted target genes of miRNAs were enriched in the JAK/STAT signalling pathway ( Figure 4B; Table S7 (Table S8). To gain more insights into the regulatory mechanism of lncRNAs, we retrieved a ceRNA sub-network mediated by STAT1, which composited by 5 lncRNAs, 29 miRNAs and 1 PCG ( Figure 4C; Table S9). Particularly, we observed potential hsa-miR-545-5p target sites in PRKCQ-AS1 and STAT1, while potential hsa-miR-665 target sites in SH3PXD2A-AS1, STAT1 and LIFR-AS1, suggesting that lncRNA may regulate the activity of JAK/STAT signalling pathway through lncRNA-miRNA-PCG circuitries. In summary, our analyses demonstrated that lncRNA may regulate JAK/STAT pathway through the ceRNA mechanism.

| Experimental validation of PRKCQ-AS1/miR-545-5p/STAT1 regulatory circuit
Richard et al have demonstrated that thousands of lncRNAs are upregulated in psoriatic skin compared with healthy skin. 47 However, only a small proportion of them have been experimentally tested, lots of them warranting further investigations. Therefore, we sought to validate the lncRNA-miRNA-JAK/STAT pathway regulatory circuit comprised of PRKCQ-AS1, has-miR-545-5p and STAT1 using HaCaT cell in vitro ( Figure 5A). Miranda target site analyses showed that STAT1 and PRKCQ-AS1 contain the binding sites of has-miR-545-5p, indicating that PRKCQ-AS1 may regulate STAT1 expression through competing for has-miR-545-5p ( Figure 5B-C). Furthermore, we performed luciferases assay to validate the binding events between miR-545-5p and PRKCQ and STAT1 ( Figure 5D). The relative luciferase activity of STAT1 mRNA and PRKCQ-AS1 was obviously reduced after co-transfection of miR-545-5p in HEK293T cells, indicating that PRKCQ-AS1 functions as a ceRNA sponge for miR-545-5p to fine-tune the gene expression of targeted genes. Further expression analysis again confirmed a strong correlation between STAT1 and PRKCQ-AS1 ( Figure 5E). To test this, we checked the expression dynamic of PRKCQ-AS1, STAT1 and has-miR-545-5p during the differentiation process of HaCaT upon IFNγ stimulation. As a result, we observed that the expression of STAT1 and PRKCQ-AS1 was activated in response to IFNγ treatment but has-miR-545-5p remained constant ( Figure 5F-H), suggesting that PRKCQ-AS1 and STAT1 are strongly associated with the activity of JAK/STAT signalling pathway and has-miR-545-5p act as a brake to modulate this process. Interestingly, the expression of PRKCQ-AS1 dramatically increased in 12 hour after IFNγ treatment, while STAT1 was activated as early as 6h ( Figure 5F-H). It means that has-miR-545-5p act as a feedback regulator to finetune the expression of STAT1, while PRKCQ-AS1 act as a sponge to F I G U R E 5 Experimental validation of PRKCQ-AS1/miR-545-5p/STAT1 regulatory circuit. A, A 3-node network motif composed of STAT1, PRKCQ-AS1 and has-miR-545-5p. B-C, Predicted has-miR-545-5p binding sites within 3'UTR of STAT1 (B) and the whole gene body of PRKCQ-AS1 (C). D, Luciferase reporter assay for the interaction between miR-545-5p and PRKCQ and STAT1. The pmirGLO luciferase reporter plasmids with STAT1 mRNA coding regions and PRKCQ-AS1 were transiently transfected into 293T cells with the miR-545-5p precursor and the negative control. And luciferase activities were measured after 72 hours. E, Scatter plot showed the expression of STAT1, and PRKCQ-AS1 is positively correlated. F-H, lncRNAs regulate the activity of JAK/STAT signalling pathway through miRNA-mediated ceRNA mechanism. Quantitative RT-PCR analysis showed the similar expressed pattern of STAT1 (F) and PRKCQ-AS1 (G) upon IFNγ stimulation in HaCaT cells, and the expression of has-miR-545-5p (H) is up-regulated. I, A schematic mechanistic model of lncRNA functions in psoriasis. In healthy skin, the expression of lncRNAs is low and the activity of JAK/STAT signalling pathway is maintained in a low level in keratinocytes; when the skin suffers chronic inflammatory stimulation, the expression of lncRNAs is elevated and lncRNAs fine-tuning the expression level of the key factors in JAK/STAT signalling pathway by competing for miRNAs, resulting the decreased suppression of JAK/ STAT signalling pathway-associated genes. Then, it then results the activation of JAK/STAT signalling pathway and lead to the secretion of cytokines in keratinocytes and exacerbate the inflammation level in psoriatic skin control the abundance of miRNAs, suggesting PRKCQ-AS1 acts as competing endogenous RNA for regulating the activity of JAK/STAT signalling pathway. Taken together, our experimental analysis showed that dysregulation of the ceRNA circuit comprised of PRKCQ-AS1, has-miR-545-5p and STAT1 act as a vital regulatory role in JAK/STAT signalling pathway activation, which involve in the hyper-inflammation of keratinocytes in psoriasis ( Figure 5I). LncRNA serve important roles in regulating various functions in the immune system. 45 Our result confirmed that lncRNA dysregulation is associated with a variety of cell signalling pathways and plays a vital role in the pathogenesis of psoriasis. Previous studies have showed that signalling pathways regulated by lncRNA may modulate multiple biological processes. 28 JAK/STAT signalling was a critical pathway for progression of inflammatory and autoimmune diseases, such as rheumatoid arthritis, psoriasis and inflammatory bowel disease. 11,12 In the present study, the co-expression results support a network model that lncRNA and protein together to trigger dysregulation of JAK/STAT signalling pathway in psoriasis. Those findings suggested that lncRNA potential acts its function through interacting with RNA binding proteins. We then investigated the regulatory associations among lncRNA and the JAK/STAT pathway by integrative analysis of the lncRNA-miRNA-PCG network. Interestingly, we found that the differentially expressed lncRNA, miRNAs and PCGs form 3-node network motifs to modulate the JAK/STAT pathway activity through ceRNA mechanism. Notably, most of the differentially expressed lncRNA and PCG pairs are regulated by the same set of miRNAs in psoriasis, indicating that lncRNAs regulate the targeted PCGs abundance through completing interaction with miRNAs. The established ceRNA network and network motif analysis methods will benefit further mechanistic investigations on psoriasis and other inflammatory disorders.

| D ISCUSS I ON
Current study determined the function of lncRNAs in the regulation of JAK/STAT pathway in psoriasis. To validate the network analysis results and investigate the potential role of lncRNA in psoriasis, 3 dysregulated lncRNA associated with JAK/STAT pathway were randomly selected for RT-qPCR analysis. Furthermore, we found that the expression of PRKCQ-AS1 and SH3PXD2A-AS1 is elevated in HaCaT cells stimulated by interferon gamma (IFNγ). Those results conferred that lncRNAs are activated in inflammatory skin, while the increased expression of lncRNA modulates the activities of JAK/ STAT pathways. The up-regulation of lncRNAs amplify the expression of key regulators in JAK/STAT signalling pathway and activate autoimmune responses, eventually lead to hyper-inflammation in psoriatic skin ( Figure 5I). In summary, the present study enhances our understanding of lncRNA function in psoriasis. Moreover, we validated that dysregulation of lncRNA-JAK/STAT pathway axis plays an important role in the pathogenesis of psoriasis and thus provide potential therapeutic targets for psoriasis treatments.

ACK N OWLED G M ENTS
We thank the high computational resource provide by biomedical computing centre in Dermatology Hospital of Southern Medical University.

CO N FLI C T O F I NTE R E S T
The authors confirm that there are no conflicts of interest.