Canonical and Interior Circular RNAs Function as Competing Endogenous RNAs in Psoriatic Skin.

(1) Background: Understanding the function of circular RNAs (circRNAs), a class of noncoding RNA, in psoriatic skin can provide important insights into the complex regulation of genes contributing to the pathogenesis of psoriasis. (2)
Methods: A novel method was applied to RNA-seq datasets from 93 skin biopsy samples to comprehensively identify circRNAs of all types, i.e., canonical circRNAs from the intron-exon junctions of mRNAs and interior circRNAs (i-circRNAs) from the interior regions of exons, introns, and intergenic regions. Selected circRNAs were experimentally validated by qRT-PCR and Sanger sequencing. CircRNAs with abundant and differential expression were identified and their putative function as competing endogenous RNAs (ceRNAs) was analyzed by an integrated analysis of circRNAs, microRNAs, and mRNAs. (3)
Results: With a comprehensive search using no information of splicing signals, we systematically identified 179 highly abundant circRNAs in psoriatic skin. Many of these were reported for the first time and many were differentially expressed in involved versus normal or uninvolved skin. Validation based on three additional RNA-seq datasets confirmed most of the identified circRNAs in psoriatic skin. Experimental analyses confirmed the expression of the well-known circRNA CDR1as, a canonical circRNA, and a novel i-circRNA in psoriasis. We also identified many circRNAs that may act as ceRNAs to regulate the expression of mRNA genes in psoriasis-related signaling pathways in psoriasis. (4) Conclusions: The result of the study suggested that circRNAs are abundant in psoriatic skin, have distinct characteristics, and contribute to psoriatic pathogenesis.

1. Introduction

Psoriasis (PS) is a chronic, inflammatory, and immune-mediated skin disease, characterized by raised, red scaly plaques [1]. Its prevalence ranges from between 0.09% and 11.4%, affecting at least 100 million individuals worldwide [2]. Besides the long-lasting and high recurrence rate, PS may increase the risk of stroke [3], myocardial infarction [4], type 2 diabetes [5], metabolic syndrome [6,7], and cancer [8,9]. The etiology and pathogenesis of PS remain poorly understood. Studies in the past decades have identified many genetic risk factors [10,11] and aberrant expression of many transcripts including non-coding RNAs such as microRNAs (miRNAs) [12,13,14,15] and long non-coding RNAs (lncRNAs) [16]. The expression of many PS related mRNAs may be regulated by miRNAs [14,15]. For example, miR-21 contributes to T-cell derived psoriatic skin inflammation [17], and miR-31 overexpression enhances the production of inflammatory cytokines and chemokines [18]. Many lncRNAs also function in PS, e.g., PRINS [19,20], lnc-IL7R [21], and LincR-Ccr2-5′AS [22].

Circular RNAs (circRNAs), a new class of non-coding RNAs, have been identified in nearly all eukaryotic species [23,24,25,26]. CircRNAs have many molecular functions [27]. For example, they act as sponges of miRNAs [28], promote transcription [29], alter the expression of mRNAs [30,31] and proteins [32], and can even turn into functional polypeptides [33,34,35]. Many circRNAs exhibit aberrant expression in complex diseases such as Alzheimer’s disease (AD) [36,37], osteoarthritis [38], and cardiovascular disease [39,40] as well as cancer [41,42].

However, the expression and potential function of circRNAs in PS are understudied and poorly understood. Three previous studies have been reported: A microarray-based analysis of psoriatic lesions identified 4,956 differentially expressed (DE) circRNAs [43]; a study of mesenchymal stem cells (MSC) of psoriatic skin described 129 DE circRNAs [44]; and a recent study of six paired lesional and non-lesional skin samples reported 148 DE circRNAs [45]. The microarray study missed novel circRNAs, the MSC study was unlikely to have captured many circRNAs in psoriatic skin, and the skin study compared circRNAs between lesional and non-lesional skin but not those of healthy controls.

We set forth to identify novel circRNAs and study their potential functions in PS. We studied a cohort of 93 skin tissue samples from which a large collection of circRNA-enriched RNA-seq libraries from 28 psoriatic-involved (PP) lesions, 38 normal (NN) controls, and 27 psoriatic-uninvolved (PN) skin samples were derived. Besides studying circRNAs arising from intron-exon boundaries, we also studied interior circRNAs (i-circRNAs) arising from the interior regions of introns, exons, and intergenic regions [46]. We confirmed our discoveries with three validation sets of RNA-seq data from psoriatic skin samples and cell lines. We also experimentally validated and studied three circRNAs in psoriatic skin with qRT-PCR and Sanger sequencing. Furthermore, we studied the potential function of circRNAs as competing endogenous RNAs (ceRNAs) in PS.

2. Results

2.1. Detection and Profiling of circRNAs in PS

A new method for finding circRNAs of all types, termed CAT, was developed to comprehensively identify circRNAs and i-circRNAs from RNA-seq data (Figure 1a, see Methods and an early version of CAT in [46]). CAT uses no information of splicing signals or genome annotation for circRNA detection but uses genome annotation to classify circRNAs into the categories of canonical circRNAs with two back-fusion (BF) points at intron-exon boundaries, complete i-circRNAs with both BF points not residing at intron-exon boundaries, and partial i-circRNAs with one BF point from an intron-exon boundary [46]. When applied to paired-end RNA-seq data, it detects a circRNA if one end of a paired-end read is mapped to the genome in an orientation-reversed splitting form (see Methods); the splitting points define the back-splicing or BF point of the circRNA (Figure 1a). Although a circRNA may span more than one exon or intron, it is unreliable to infer its actual structure based on RNA-seq data alone, particularly if RNA-seq data were not enriched for circRNAs. Therefore, we focused on candidate circRNAs from single exons, single introns, pairs of adjacent exons and introns, and intergenic transcripts. We compared CAT with the four most popular existing methods for circRNA detection, i.e., CIRI2 [47], CIRCexplorer2 [48], DCC [49], and find_circ [23], on a single-end stranded HeLa dataset [50] following the scheme proposed in [51]. CAT had a comparable performance with the other four methods, i.e., the estimated error rate of CAT was ~3% lower than that of DCC and ~2% higher than that of find_circ (Figure 1b). In summary, CAT is an exclusively RNA-seq data-driven approach for identifying circRNAs from any genomic locus without using any information of splicing signals or genome annotation.

Figure 1

Workflow of CAT and analysis of circRNAs in PS. (a) Illustration of split mapping of paired-end RNA-seq reads to candidate circRNAs. If Read1 is completely mapped and Read2 is split mapped to the reference genome, a circular structured RNA, circRNA, could be formed from this genomic locus with the back-fusion point identified from the split mapping. (b) Stacked bar plots of the circRNAs in HeLa cells (RNase R vs. control) that were predicted by five circular RNA identification methods, with results separated into RNase R resistant (≥ fivefold enrichment, green), unaffected (one to fivefold enrichment, gray) and RNase R sensitive (depleted in RNase R treated samples, red). Percentages refer to the fraction of RNase R sensitive circRNAs defined as false positives. The minimal number of supporting reads was set to 3 for all five methods compared. (c) Workflow for identification and analysis of circRNAs in the ribo-zero, stranded paired-end RNA-seq data from 93 PS samples. In total, 8856 unique circRNAs were identified and 2863 circRNAs were annotated and classified by single exons, single introns, pairs of adjacent exons and introns, and intergenic non-coding transcripts. A total of 179 highly expressed circRNAs were selected by three criteria, i.e., expressed in more than 10 samples, supported by more than five reads in at least one sample and filtered based on a mapping quality threshold of 40 for both anchors. This resulted in 51 differentially expressed circRNAs when PP vs. NN skin was compared by Wilcoxon rank sum test. (d) The number of BF reads supporting 8856 unique circRNAs, and (e) the number of samples expressing 8856 unique circRNAs, where the y-axes are logarithmic (log10).

Using CAT, we detected and characterized circRNAs in 93 biopsy skin samples for circRNA discovery (Figure 1c and Supplementary Figure S1a, Table S1). We identified 8856 circRNAs supported by at least two sequencing reads (Figure 1c–e). To increase the confidence of reporting genuine circRNAs, we focused on 2863 circRNAs originating from single exons, single introns, pair of adjacent exons and introns, and intergenic transcripts. As expected, splicing signals were observed on most canonical circRNAs arising from intron-exon boundaries but not on complete i-circRNAs (Figure S1b). In total, 91.69% (739/806) and 72.65% (263/362) of the canonical and partial boundary circRNAs were accompanied by a splicing signal, whereas this ratio decreased to 11.98% (203/1,695) for complete i-circRNAs. Combined, the results showed that i-circRNAs existed broadly in skin tissues and that most i-circRNAs were not accompanied by the splicing signal.

2.2. circRNAs in PS and their Characteristics

We identified 179 circRNAs expressed in at least 10 samples and supported by at least five reads in one of these samples (Figure 1c, Table S2). Among the 179 circRNAs, 64 (35.8%) were identified for the first time and 65 (36.3%) were i-circRNAs, 49 of which were novel. As a validation of the result, we examined the expression of the 179 circRNAs in three additional RNA-seq datasets (Figure 2a). In total, 170 circRNAs were expressed in at least one of the three validation datasets. In particular, 166 (92.7%) of the 179 circRNAs were expressed in another large RNA-seq dataset of 34 psoriatic and normal skin biopsy samples in dataset hsa_skin2, and 37 circRNAs were highly expressed in all three validation datasets. Among these, three were exon or intron i-circRNAs with no associated splicing signal (Figure 2a). A substantial portion of the boundary and partial i-circRNAs, but very few of the complete i-circRNAs, of the 179 circRNAs carried the splicing signal (Figure 2b).

Figure 2

Distribution and characteristics of abundant circRNAs. (a) Names and number of samples in the discovery and validation datasets. In total, 179 abundant circRNAs were identified in the discovery dataset, among which 166, 83, and 51 were independently validated by the three validation datasets. (b) Distribution of the 179 circRNAs in the discovery datasets. Boundary circRNAs are canonical circRNAs, and i-circRNAs are further classified into partial and complete i-circRNAs. All circRNAs are further annotated as to whether they are adjacent to splicing signals. (c) An example of RNA folding structure and complementary sequences flanking an exon i-circRNA, hsa_skin_038904. (d) Distribution of the lengths of short homologous sequences (SHS). (e) An example of two circRNAs, hsa_skin_175896 and hsa_skin_226345, arising from FLG (filaggrin). (f) Average expression (spliced reads per billion mappings [SRPBM]; see Methods) of the 179 highly expressed circRNAs in PP and NN skin. (g) Volcano plot visualizing differential expression of circRNAs when PP and NN skin samples were compared. The red and black dots in the plot represent significantly differentially expressed (p-value < 0.01, see Methods) and not significantly expressed circRNAs, respectively. Circle and cross marks represent interior and boundary circRNAs, respectively.

Many circRNAs have two features: (1) they have complementary sequences flanking their BF points. Twenty-four circRNAs had paired flanking sequences, among which 17 were i-circRNAs including 15 with no splicing signal (see Figure 2c for an example). This suggested that the production of these i-circRNAs may be assisted by fold-back structures of complementary sequences; and (2) the BF points of many circRNAs, particularly i-circRNAs, resided within short homologous sequences (SHS), with lengths of 2- to 56-nt and peaked at 2-nt (Figure 2d). Among the 179 abundant circRNAs, 101 (56.4%) were associated with SHS that were longer than 2-nt, 29 of which had no splicing signal (Table S2), indicating a potential role of SHS in circRNA biogenesis [46].

More than one circRNA may arise from a genomic locus. Among the 161 host genes of the 179 circRNAs, 13 (including FLG, FLG2, KRT10, KRT2, and LPP; Table S2) produced multiple circRNAs. For example, hsa_skin_175896 contained hsa_skin_226345 (Figure 2e).

2.3. Aberrantly Expressed circRNAs in PS

The digital expression level of a circRNA was quantified by the number of reads splitting mapped to its BF point. The digital expression was normalized to the reads per billion mappings (SRPBM; see Methods). Using SRPBM, 47 (26.3%) of the 179 circRNAs were expressed at higher levels between PP vs. NN (Figure 2f) and 51 (28.5%) circRNAs were differentially expressed (DE) (see Methods) in PP vs. NN skin (Figure 2g).

A putative function of circRNAs is acting as competing endogenous RNAs (ceRNAs) [52] for regulating mRNA genes that share common regulating miRNAs with circRNAs. We named these circRNAs ce-circRNAs and the associated mRNAs circRNA-associated genes. There were 8096 mRNA genes significantly DE between PP vs. NN skin (Table S3), among which 4034 were associated with 51 DE circRNAs.

From these circRNAs, we selected two canonical circRNAs, hsa_skin_194345 (CDR1as) and hsa_skin_088763 (hsa_circ_0109327), and three i-circRNAs, hsa_skin_100269, hsa_skin_143837, and hsa_skin_052271 for experimental validation using psoriatic skin biopsy samples (PP and PN) and HaCaT keratinocyte cells (see Methods). Using sequencing reads across BF points as surrogates to circRNAs, the five circRNAs had 1,440, 405, 41, 24, and 22 reads in the discovery dataset, respectively. Note that these five circRNAs were all detected in the validation datasets, providing the first validation of these circRNAs in psoriasis. For experimental validation, divergent and convergent PCR primers (Table S4) were applied separately to the RNA and DNA of these circRNAs. Three circRNAs, CDR1as, hsa_skin_052271, and hsa_skin_088763, were experimentally validated (Figure 3a) and their BF points were confirmed with Sanger sequencing (Figure 3b). Sanger sequencing successfully recovered the full-length of hsa_skin_052271. Furthermore, the three validated circRNAs were subjected to RNase R treatment (Figure S2a) and the changes in their expression in PP vs. PN were confirmed. The fold changes from the RNase R experiments were consistent with and even more pronounced than those of the untreated ones (Figure 3c and Supplementary Figure S2b). qRT-PCR was also applied to validate their RNA-seq based DE in PP vs. PN skin and consistent fold change directions were confirmed (Figure 3d and Supplementary Figure S2c). We also determined that these circRNAs were DE in PP vs. NN skin by an analysis of sequencing data from the hsa_skin2 validation dataset (Figure 3e).

Figure 3

Experimental validation of three identified circRNAs in psoriatic skin and HaCaT keratinocyte cells. (a) Validation of a well-known intergenic i-circRNA (CDR1as), a novel exon i-circRNA (hsa_skin_052271), and a canonical circRNA (hsa_skin_088763) by PCR. The divergent and convergent arrows above the gel image represent, respectively, the divergent and convergent PCR primers used in RNA (cDNA) and DNA (gDNA), and the white arrows in the gel image point to circRNAs. Numbers on the left ranging from 100 to 1kb+ represent sizes in DNA ladder. GAPDH was used as an internal control. (b) Validation of CDR1as, hsa_skin_052271, and hsa_skin_088763 by Sanger sequencing. (c) Box and whisker plot for RNase R validation of differential expression of CDR1as, hsa_skin_052271, and hsa_skin_088763 when PP vs. PN skin was compared. (d) Quantitative real-time PCR (qRT-PCR) validation of differentially expressed circRNAs hsa_skin_088763, CDR1as, and hsa_skin_052271 in PP vs. PN skin. (e) RNA-seq (hsa_skin2 of validation set) validation of differentially expressed circRNAs CDR1as, hsa_skin_052271, and hsa_skin_088763 in PP vs. NN skin.

2.4. Putative circRNA Functions in PS

The three experimentally validated circRNAs may potentially function as ce-circRNAs in psoriatic skin. CDR1as may regulate 149 out of 655 targets of miR-7-5p via 67 binding sites and 242 out of 962 targets of miR-135b-5p via four binding sites (Figure 4a, Table S5). CDR1as was sharply down-regulated 3.7-fold between PP vs. NN (Table S6), whereas both miR-7-5p and miR-135b-5p are up-regulated in psoriasis [12,15,53]. The target genes that were down-regulated in PP skin and positively correlated with the expression of CDR1as (Spearman’s rank-order correlation analysis) were called circRNA-associated genes of CDR1as (see Methods, Figure 4a). Among the associated genes were EGR3, GATA6, GATA3, and FOXN3, which play important roles in psoriasis. EGR3 can regulate late epidermal differentiation and contribute to the keratinocyte differentiation-related module in a skin-specific network [54]. GATA6 has significantly lower expression in psoriatic dermal MSCs [55]. FOXN3, regulating cell differentiation and cell cycle, is down-regulated along with GATA3 in psoriasis [56].

Figure 4

Putative functions of three circRNAs as ceRNAs. The genomic origin of (a) intergenic i-circRNA CDR1as, (b) exon boundary circRNA hsa_skin_088763, and (c) exon i-circRNA hsa_skin_052271, their psoriasis-related binding miRNAs predicted by miRanda, and their circRNA-associated genes. Volcano plots visualize the differential expression of target genes between PP and NN skin. The red and black dots in the plot represent significantly down-regulated (p-value < 0.01, see Methods) and other target genes, respectively. Among the significantly down-regulated target genes, those that are positively correlated with the expression of circRNAs are defined as circRNA-associated genes (see Methods). Scatter plots show Spearman’s rank-order correlation between three circRNAs and one of their associated genes, where r is the Spearman correlation coefficient.

CircRNA hsa_skin_088763 from pseudogene RP11-255H23.2 may regulate 109 out of 560 targets of miR-338-3p via four binding sites and 340 out of 1,508 targets of miR-23a/b-3p via 10/11 binding sites (Figure 4b, Table S5). hsa_skin_088763 was down-regulated 2.6-fold when PP vs. NN was compared (Table S6), and miR-338-3p and miR-23a/b-3p were up-regulated in psoriasis, consistent with this trend in their putative regulatory circRNA [57]. Several psoriasis-related genes were among the associated genes of hsa_skin_088763, including GATA6, SIK2 (Figure 4b), IL17RD [58], EGR3, FAS, LRIG1, and PPARGC1A [59]. LRIG1 negatively regulates growth factor signaling to regulate epidermal stem cell quiescence [60]; SIK2 modulates cytokine responses during innate immune activation [61]; FAS signaling is essential for inducing key inflammatory cytokines in psoriasis [62]. All these results indicated that hsa_skin_088763 functions as a ceRNA in psoriasis.

hsa_skin_052271, an i-circRNA from the third exon of gene FLG2, may mediate three target genes of miR-135b-5p, two target genes of miR-205-5p, and nine target genes of miR-27a-3p (Figure 4c, Table S5). Hsa_skin_052271 was down-regulated 2.8-fold between PP and NN (Table S6), and miR-135b-5p [12,15,53], miR-205-5p [12] and miR-27a-3p [57] were up-regulated in psoriasis. Psoriasis-related gene GATA6 was an associated gene of hsa_skin_052271 (Figure 4c).

We performed GO and KEGG pathway analyses separately on the 1,558 up-regulated and 2,476 down-regulated circRNA-associated transcripts (see Methods). The results revealed that most of the top 20 KEGG pathways and many of the top 20 GO terms were psoriasis-related (Figures S3–S6). The above pathway analysis added additional pieces of evidence that CDR1as, hsa_skin_088763, and hsa_skin_052271 may play essential roles in psoriatic pathogenesis by acting as ceRNAs.

3. Discussion

CircRNA has not been well-investigated in psoriasis. The only previous study on this species of RNA in psoriatic lesional skin used a small cohort of six patients but no healthy controls [45]. That study revealed that circRNAs in psoriatic skin exhibited a lower abundance in lesional skin than in non-lesional skin, which was also observed in our current study (Figure 2f,g). That study also reported that miRNAs were in large part not responsive to the differential expression of circRNAs because the miRNAs targeted by circRNAs were profiled as a group so that the regulatory relationships of pairs of circRNAs and miRNAs were lost.

In the study described here, we leveraged a large collection of RNA-seq data from psoriatic and normal skin biopsies to detect and profile the expression of circRNAs and identify their potential targets contributing to PS. Among the 179 highly abundant circRNAs detected, 64 (35.8%) were identified the first time, and 51 (28.5%) exhibited aberrant expression in PS. Interior circRNAs (i-circRNAs), a newly characterized type of circRNA [46], were also detected in PS. Specifically, 65 (36.3%) of the 179 circRNAs were i-circRNAs and 49 of the 65 i-circRNAs were novel.

We validated most of the identified circRNAs and i-circRNAs with three new validation RNA-seq datasets and experimentally confirmed and studied the expression of three circRNAs (CDR1as, hsa_skin_088763, and hsa_skin_052271) in PS by PCR and Sanger sequencing (Figure 3, Figure S2).

The significance of this experimental confirmation is threefold. First, it demonstrated for the first time that the most well-studied circRNA CDR1as can potentially function in psoriasis. CDR1as is known to play a role in cancer [63,64], neurological diseases such as Alzheimer’s disease [36,37], and many other diseases, and during development [65]. CDR1as was the fourth most abundant and the topmost DE circRNA in psoriatic skin (Tables S2 and S6). It was not only aberrantly expressed but also functioned as a ce-RNA [31,52] through a circRNA-miRNA-mRNA regulatory cascade by regulating miR-7-5p and miR-135b-5p in psoriasis. Psoriasis is an inflammatory skin disorder where hyperproliferation of keratinocytes is accompanied by abnormal differentiation of immune cells [66]. The theme of abnormal cell proliferation that is common in psoriasis and cancer explains, in part, the aberrant expression of CDR1as in these different diseases. Indeed, our results suggest that many aberrantly expressed circRNAs regulated many genes function in cancer-related pathways (Figures S3–S5). miR-135b-5p is also reported to target KLF4 [67], an important transcription factor in skin barrier formation [68]; a pathway that is disrupted in psoriasis.

The second significance of this experimental confirmation is that it showed, for the first time, that interior circRNA (i-circRNA) exist and potentially function in human cells [46]. The novel i-circRNA hsa_skin_052271 emerged from an exon and its BF points were not adjacent to the splicing signal. Importantly, hsa_skin_052271 may function as a ce-RNA, mechanistically acting as a sponge of miR-27a-3p and thus positively regulating the expression of GATA6, which is a target of miR-27a-3p and down-regulated significantly in psoriatic dermal MSCs [55].

The third is that the result revealed a possible connection between a circRNA and a pseudogene. CircRNA hsa_skin_088763 arose from pseudogene RP11-255H23.2. Many pseudogenes are known to function as ceRNAs as they are prone to become the targets of miRNAs [52]. Therefore, circRNAs and pseudogenes may crosstalk by competing for miRNA binding. This is likely the case in our study since both hsa_skin_088763 and pseudogene RP11-255H23.2 are down-regulated in psoriatic lesions (Tables S2 and S3).

The i-circRNAs reported in the current study were genuine circRNAs, which was supported by the validation datasets. i-circRNAs were previously identified from intergenic transcripts in rice [69] and oncogenic chimeric transcripts in cancer [70]. We previously reported a large number of i-circRNAs in humans, mice, and rice [46]. The results from the current study expanded our understanding of i-circRNAs in complex disease.

The discovery of i-circRNAs was thought-provoking and shed new light on the biogenesis of circRNAs. While the splicesome-based biogenesis model is well-supported, it cannot explain all canonical circRNAs. In particular, 19,444 (21.05%) of the 92,369 circRNAs in circBase have no splicing signal adjacent to their back splicing points. Interestingly, the back fusion points of most i-circRNAs and some canonical circRNAs were embedded in short homolog sequences (SHS, Figure 2d) [46]. We hypothesized earlier that i-circRNAs and canonical circRNAs with SHS are products of RNA polymerase template switching during RNA processing [46]. This needs to be tested in future studies.

4. Materials and Methods

4.1. RNA-seq Data

Four paired-end RNA-seq datasets from human skin were used to investigate circRNAs in PS (Figure S1a, Table S1). One was used for discovery, the other three were used for validation. The discovery dataset, hsa_skin1 (GSE121212), was derived from human skin biopsies where the TruSeq Stranded Total RNA Protocol was used in combination with the RiboZero rRNA removal Kit. It consisted of 28 PP, 38 NN, and 27 PN samples. One of the validation datasets, hsa_skin2 (GSE74697), was generated from human skin samples with the ScriptSeq complete kit from Epicenter and contained 18 PP and 16 NN samples. The other two datasets, hsa_MSC1 (GSE81106) and hsa_MSC2 (GSE89725), were from skin-derived Mesenchymal stem cells (MSC) and derived with the TruSeq Stranded Total RNA Library Prep kit using the CircRNA Enrichment kit and ribo-zero rRNA Removal kit, respectively. Each of the MSC datasets was from three PP and three NN samples.

4.2. Psoriatic Skin Biopsy Samples and Cells for Validation

Psoriatic samples (PP and PN) from skin biopsies of 30 PS patients and HaCaT keratinocyte cells were used for validation. These psoriatic samples were from patients with untreated psoriasis of moderate-to-severe intensity (PASI > 12, BSA at least 10%) and they were all collected with ethical approval obtained from The Rockefeller University Institutional Review Board (RU IRB, New York city, NY, USA). The RNA was extracted with the miRNeasy kit for purification of total RNA (Qiagen, Germantown, MD, USA).

4.3. RNase R Digestion

RNA was digested with RNase R (BioVision, Milpitas, CA, USA) (a 5′ to 3′ exonuclease) to degrade linear RNA and enrich for circular RNA. RNase R digested RNA was prepared with a reaction containing 2 μg of prepared RNA, 1 μL (40U) RiboLock (ThermoFisher-Scientific, Waltham, MA, USA), 2 μL 10× RNase R reaction buffer, and 1 μL (10U) RNase R; adjusting the volume to 20 μL with nuclease-free water [71]. RNA was then converted to cDNA with a Quantitect Reverse Transcription kit (Qiagen, Germantown, MD, USA).

4.4. PCR and Sanger Sequencing

Divergent primers (FisherScientific, Waltham, MA, USA) were designed with the CircBase IDs of circRNAs in the CircInteractome webtool, if available; otherwise, primers were manually designed across the back-fusion (BF) points. The primers used for validation are listed in Table S4. Divergent primers are outward-facing primers that should only amplify the circular RNA. Convergent primers were also designed as a control. Convergent primers can generate products from genomic DNA, whereas divergent primers cannot. The correctness of the BF sites of these amplified products was verified with Sanger sequencing (Psomagen, Brooklyn, NY, USA).

4.5. Quantitative Real-Time PCR

Divergent primers were used to amplify cDNA in a qRT-PCR reaction with SYBR Green (Kapa Biosystems, Wilmington, MA, USA). GAPDH was used as a reference. Both target and reference were amplified in triplicate and the relative level of each circRNA was calculated with the 2−∆∆Ct method.

4.6. Identification of circRNAs of All Types—The CAT Method

A new method was developed to comprehensively search for circular RNAs of all types (i.e., the CAT method) using paired-end RNA-seq data without information of splicing signals. CAT first maps sequencing reads to the reference genome with Bowtie2 [72]. The mappable reads are discarded since they are from linear RNAs. CAT takes a left and a right x-mer on the 5′-end and 3′-end of an unmapped read, called the left and right anchors, respectively (Figure 1a). It attempts to split-map the read to the genome—the left anchor maps to a locus that is downstream from the locus to which the right anchor maps (Figure 1a). If successful, it then searches for the split point within the sequence between the two anchors such that when the read is split at the split point, the two segments of the read can be aligned to the genomic loci that the anchors determine (Figure 1a). In the current CAT implementation for sequencing read length of 100 bp, x was initially set to 20, and no more than two mismatches were allowed in the mapping. Once a candidate BF point was found from one read of a pair, an additional criterion was adopted to ensure the other read was completely mapped to the region within the two BF points.

CAT retains those BF points that are supported by at least k (e.g., k > 2 in most of our analyses) split-mapped reads in an RNA-seq library. Furthermore, only unambiguous candidate BF points are considered. Specifically, while there certainly exist circRNAs spanning more than one exon, it is difficult to infer such circRNAs with certainty based on RNA-seq data alone, particularly on RNA libraries not enriched for circRNAs. Therefore, circRNAs from single introns, single exons, pairs of adjacent exons and introns, and intergenic non-coding transcripts are considered if RNA-seq libraries are not RNase R treated (Table S2). The CAT software in python is available in the github repository at https://github.com/xiaoxin8712/find_circ_strand (accessed on 27 July 2019).

4.7. Identification of Differentially Expressed circRNAs

Two criteria were introduced to identify abundantly expressed circRNAs. First, the minimal number of reads mapped to a BF point in at least one sample is no less than k = 5. Second, the total number of samples in which the circRNA appears is no less than m, and a parameter is adjusted for different sample sizes; in our experiment, m = 10. The circRNAs that did not meet these two criteria were considered as having low expression. The expression levels of highly expressed circRNAs were normalized by SRPBM (spliced reads per billion mappings), i.e., the number of circular reads/number of mapped reads (units in billion), to ameliorate the effects of sequencing depth and batch effects. Differentially expressed (DE) circRNAs were identified by the RankSum method [73] with a threshold of p-value ≤ 0.01 based on SPRBM.

4.8. Identification of Differentially Expressed Genes

DE genes of PP versus NN and PN were obtained with the Tophat and Cufflinks pipelines [74], with the threshold being adjusted to q < 0.05 and fold change >2.

4.9. GO and KEGG Pathway Analyses

Analyses of GO and KEGG pathway enrichment were carried out with Metascape [75] with the default parameters: minimum overlap of 3, p-value cutoff of 0.01, and minimum enrichment of 1.5, with GO biological processes and KEGG pathway for enrichment generated separately.

4.10. CircRNA-Associated Genes

CircRNA-associated mRNAs were defined as DE mRNAs that were positively regulated by circRNAs via circRNA-miRNA-mRNA regulatory cascades. That is, a DE mRNA was associated with a circRNA if both were bound and regulated by a common miRNA. mRNA and circRNA targets of mature miRNAs [76] were predicted with miRanda [77] using the default parameters of a score of 140.0 and energy of −1.0. CircRNA-associated genes were selected from the DE genes that were positively correlated with their corresponding circRNAs, i.e., Spearman’s rank correlation coefficient > 0.0 and p-value < 0.01.

5. Conclusions

Overall, our results reveal and demonstrate the diversity of circRNAs and the functional complexity of circRNAs operating in psoriasis. There exist not only circRNAs originating from exon-intron boundaries, but also i-circRNAs from interior regions of exons, introns, and intergenic regions. Importantly, circRNAs may function as ce-RNAs and contribute to psoriasis pathogenesis. Our results suggest that CDR1as, hsa_skin_088763, and hsa_skin_052271 should be analyzed in future studies as novel biomarkers for PS diagnosis.