Identification of Long Noncoding RNA Associated ceRNA Networks in Rosacea

Rosacea is a chronic and relapsing inflammatory cutaneous disorder with highly variable prevalence worldwide that adversely affects the health of patients and their quality of life. However, the molecular characterization of each rosacea subtype is still unclear. Furthermore, little is known about the role of long noncoding RNAs (lncRNAs) in the pathogenesis or regulatory processes of this disorder. In the current study, we established lncRNA-mRNA coexpression networks for three rosacea subtypes (erythematotelangiectatic, papulopustular, and phymatous) and performed their functional enrichment analyses using Gene Onotology, KEGG, GSEA, and WGCNA. Compared to the control group, 13 differentially expressed lncRNAs and 525 differentially expressed mRNAs were identified in the three rosacea subtypes. The differentially expressed genes identified were enriched in four signaling pathways and the GO terms found were associated with leukocyte migration. In addition, we found nine differentially expressed lncRNAs in all three rosacea subtype-related networks, including NEAT1 and HOTAIR, which may play important roles in the pathology of rosacea. Our study provided novel insights into lncRNA-mRNA coexpression networks to discover the molecular mechanisms involved in rosacea development that can be used as future targets of rosacea diagnosis, prevention, and treatment.


Introduction
Rosacea is a complex chronic and recurrent dermatological condition characterized by flushing, transient/persistent erythema, telangiectasia, inflammatory papules and pustules on the central face, phymatous changes, and ocular manifestations [1,2]. Rosacea is generally classified into four major subtypes: erythematotelangiectatic (ETR), papulopustular (PPR), phymatous (PhR), and ocular (OR) [2]. Among these subtypes, ETR is the most common one, followed by PPR, PhR, and OR (not discussed in the current work).
Although clinical subtypes are clearly described, the various clinical presentations make the pathophysiology of rosacea elusive; therefore, its precise pathogenesis remains unclear. Several studies suggest that dysregulation of the immune pathways and neurovascular changes are found to varying degrees in different rosacea subtypes [3]. e most studied pathway in rosacea is the cathelicidin activation pathway, followed by inflammasome-associated pathways [3,4]. T cells response involved in rosacea pathogenesis is dominated by 1/ 17-polarized immune cells [5].
Long noncoding RNAs (lncRNAs), a class of over 200 nucleotides (nt), are a significant category of ncRNAs involved in a series of biological functions [6]. Abnormal expression of ncRNAs is often related to various human diseases, including cancer, inflammation, and autoimmune diseases. For example, the roles of lncRNAs in cutaneous squamous cell carcinoma [7], atopic dermatitis [8], psoriasis [9,10], and chronic actinic dermatitis [11] have recently been demonstrated, including immune response, epidermis development, and regulation of leukocyte-mediated cytotoxicity [10,12]. MicroRNAs (miRNAs), the small ncRNAs with 20-22 nt length, are key controllers of gene expression by targeting messenger RNAs (mRNAs) [13]. Advanced evidence illustrates that lncRNAs and miRNAs not only play important roles in the progression of cutaneous tumor or inflammatory dermatosis separately but also work together to modulate pathways [14][15][16]. Competing endogenous RNAs (ceRNAs) are noncoding/coding RNAs that compete for the same miRNA to reduce the amount of miRNAs available for interaction via the binding of the miRNA recognition/response elements (MREs); LncRNAs can act as ceRNAs to target and recruit miRNAs [17,18]. For example, in ovarian cancer, lncRNA (NEAT1) functions as a ceRNA for miRNA506 (miR-506) to promote cell proliferation and migration [19].
To our knowledge, the role of lncRNAs in rosacea, specifically in immune response, is poorly understood. We evaluated the changes in mRNA and lncRNA expression in each rosacea subtype after reannotating microarray data. Our results revealed that both mRNAs and lncRNAs expression profiles varied in each rosacea subtype. Furthermore, we established lncRNA-related ceRNA networks in rosacea to explore the possible correlations among differentially expressed (DE) lncRNAs, miRNAs, and mRNAs.
e present study provides a novel and comprehensive overview of the functions and mechanisms of differentially expressed lncRNAs in rosacea pathophysiology.
ere was a total of 58 samples, including 38 rosacea samples (14 ETR, 12 PPR, 12 PhR) and 20 normal samples and each sample consisted of the corresponding RNA-sequence data and mRNA-sequence data. ese data were analyzed to identify all gene features in the human genome using a combination of lncRNA reannotation and bioinformatics analysis. Since there were around 1,000 lncRNAs annotated in the library files of Affymetrix HG-U133 plus 2.0 platform [20], we performed the lncRNA reannotation of probes based on the protocols reported previously [21][22][23]. e retained probes should match the chromosomal position of corresponding lncRNAs retrieved from the GENCODE project (human GENCODE Release 29, GRCh38.p12). e annotated transcripts longer than 200 nt and categorized as "lincRNA," "antisense," "processed_transcript," and "non_coding" were identified as lncRNAs. After removing the multiple annotated probes of both lncRNAs and coding genes, and the lncRNA transcripts with less than three independent probes, we obtained 4,508 lncRNAs for Affymetrix HG-U133 plus 2.0 microarray. Before processing the differential expression analysis by limma package in R, all the expression values were normalized and log 2 transformed [24].

Putative Target Gene Prediction and Enrichment
Analysis. To predict the putative target genes of significantly DE lncRNAs, we applied coexpression and colocation methods. Coexpression was predicted on the basis of the expression correlation coefficient (Pearson's correlation >0.95 or <-0.95). In colocation, all the protein-coding genes 10 kb upstream or downstream of the differently expressed lncRNAs were searched [25]. Gene ontology (GO) enrichment analysis was carried out using the clusterProfiler package [26], and GO terms with corrected p values <0.05 were regarded as significantly enriched terms. Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology based annotation system (KOBAS) software was applied to test the significance of DE lncRNAs putative target genes in KEGG pathways (http://www.genome.jp/kegg) [27].

Module
Analysis and GSEA. Different modules were identified by a weighted gene coexpression network analysis module (WGCNA) algorithm. A heat map was used to depict the strength of interactions between different modules embedded in R language. All DE genes (DEGs) were included to establish coexpression modules. e "flashclusttool" package was applied to execute the cluster and trait analysis. Clustered Spearman correlation matrix for different RNA-seq replicates for all the groups. Gene set enrichment analysis (GSEA) was used to function annotation of all gene expression data.

Patients and Specimens.
is study was conducted in accordance with the ethical guidelines of the Declaration of Helsinki and approved by the Medical Ethics and Human Research Committee of West China Hospital (No. 2017.163). A total of 7 samples and matched normal adjacent specimens (NAS) used to analyze lncRNAs expression were obtained from patients who were diagnosed with rosacea at West China Hospital, Sichuan University, China. All patients were informed with consent.

Statistical
Analysis. All statistical analyses were conducted on Graph Prism 7.
e data measurements were presented as mean ± standard deviation (SD). Statistical analyses were carried out by Student's t-test and one-way ANOVA as appropriate. p < 0.05 was considered as significant.

Reannotation of lncRNAs and Bioinformatics Analysis of DE RNAs.
Expression levels of lncRNAs and mRNAs were used to produce heat maps Both lncRNAs and mRNAs were self-segregated into three subtypes of rosacea and normal groups, as shown in Figures 1(a) and 1(b). Volcano plots estimated the statistical significance of differentially expressed (DE) lncRNAs and mRNAs (Figures 1(c) and 1(e)). Subsequently, DE lncRNAs and mRNAs were screened between different subtypes of rosacea and normal groups by unsupervised hierarchical clustering analysis (Figures 1(d) and 1(f )). We identified 13 DE lncRNAs between the normal and ETR groups that were statistically significant. Five out of 13 were upregulated and eight downregulated. Moreover, there were 25 statistically significant DE lncRNAs between the normal and PhR groups, of which 15 were upregulated and 10 downregulated. In addition, the number of significant DE lncRNAs between the normal and PPR groups was 25, where 18 lncRNAs were upregulated and 7 downregulated. Regarding mRNAs, we identified 777 DE mRNAs with statistical significance between the normal and ETR groups, of which 505 genes were upregulated and 272 downregulated. Furthermore, a total of 1051 DE mRNAs showed statistical significance between the normal and PhR groups, with 532 upregulated genes and 494 downregulated genes. Finally, there were 1069 DE mRNAs with statistical significance between the normal and PPR groups, among which 643 genes were upregulated and 426 were downregulated. e overlapping Venn diagram showed that six DE lncRNAs identified were significantly upregulated, while three were significantly downregulated in all rosacea subtypes (Figure 2(a)). Besides, 525 differentially annotated mRNAs were obtained in all subtypes of rosacea ( Figure 2(b)). e functions of differentially expressed lncRNAs were predicted through the GO analysis. Trans targets were enriched in ten biological processes (BP), in eight cellular components (CC), and in nine molecular functions (MF) (Figure 2(c)). GO analysis in different groups are shown in Figures S1a, S2a, and S3a. Most significant enriched BP identified in ETR group were leukocyte migration, epidermis development, and skin development. e top three enriched BP identified in PhR group were T cell activation, leukocyte migration, and neutrophil activation. Finally, leukocyte migration, cell chemotaxis, and leukocyte chemotaxis were the most significant enriched BP identified in PPR. All enriched GO terms were displayed as an interaction network (Figure 2(d)), which were divided into three categories: leukocyte migration, neutrophil activation, and T cell activation.
Enriched GO biological processes of DE mRNAs involved in the lncRNAs network are shown in Figure 2(e). Forty and ten genes were up-and downregulated, respectively. Most of these genes are associated with chemotaxis and migration. Different candidate genes identified in the three most significant BP coexpression network in the three rosacea subtypes were represented in Figures S1b, S2b, and S3b.

Construction of Cluster Dendrogram and WGCNA.
We identified thirteen distinct coexpressing gene modules in rosacea groups. Next, we constructed modules with different colors and clustered the DEGs (Figure 3(a)). Interactions between genes by weighted gene coexpression network analysis (WGCNA) module illustrated that there were notable differences in the correlation between normal and the three rosacea subtypes, especially in the PPR group ( Figure 3(b)). e association between clinical phenotypes and genome-wide expression profiles was studied using "Eigengene Networks" tool ( Figure 4). Here, the "green" and "red" modules were highly associated with ETR. On the other hand, the "blue," "magenta" and "grey" modules were highly associated with PPR and the "pink" module was most highly associated with PhR.
Gene set enrichment analysis (GSEA) was conducted on the genes in the thirteen constructed modules (Figures 3(c)-3(e)). Genes in ETR group were mainly enriched in ferroptosis and PPAR signaling pathway. e core regulators on ferroptosis included GPX4, MAP1LC3B, and VDAC3 genes. On the other hand, the core regulators on PPAR Gene expression signaling pathway were PPARG, CD36, and MMP1 genes. Similarly, genes in PhR group were mainly enriched in ferroptosis, NOD-like receptor, and PPAR signaling pathways. e core regulators genes on NOD-like receptor signaling pathway were STAT1, IL1B, NLRX1, IL6, and VDAC. Genes in PPR group were enriched in NF-κB, NODlike, and PPAR signaling pathways. IL1B, TNFSF13B, and VCAM1 genes were regulators of NF-κB signaling pathway.

RT-qPCR and FISH.
e results by RT-qPCR revealed the increased expression of NEAT1 and decreased expression of HOTAIR, ZNF667-AS1 in rosacea tissues, which was consistent with microarray ( Figure 6(a)). FISH staining (Merge) confirmed the increased expression of NEAT1 and decreased expression of HOTAIR and ZNF667-AS1 in rosacea compared to normal adjacent specimens (Figure 6(b)).

Discussion
e exact molecular mechanisms of rosacea are unclear; however, this disorder is expected to be the result of a combination of several causes with genetic predisposition [3]. Accumulating evidence points to activated cellular pattern recognition receptors, such as toll-like receptor 2 (TLR2), transient receptor potential ion channels, and inflammatory mediators as actors of critical steps in the clinical diagnosis of rosacea [31]. e cathelicidin activation pathway, initiated by activated TLR2, leads to upregulation and release of matrix metalloproteases (MMPs) and kallikrein 5 (KLK5) [32]. MMPs and KLK5 participate in a proteolytic cascade and activate the cathelicidin peptide LL-37 from its inactive precursor, inducing a sequence of proinflammatory events, such as leucocyte chemotaxis, activation of NF-κB, and angiogenesis [33]. Besides, the inflammasome components NOD-like receptor 3 and caspase-1 can activate interleukin-1β (IL-1β). Both NF-κB responsive and inducible genes enhance the expression of tumor necrosis factor (TNF), IL-8, and cyclooxygenase-2, which mediate neutrophil chemotaxis, inflammation amplification, prostaglandin E2 synthesis, and downstream vascular effects [34,35].
Dysregulation of lncRNAs has been reported to be associated with several human diseases, including autoimmune diseases and different types of tumors [7,8,21,[36][37][38][39]. Recently, lncRNAs broad expression in cutaneous disorders including systemic lupus erythematosus (SLE), sensitive skin, psoriasis, and aging process has been documented [40][41][42][43]. ese results raise the possibility that lncRNAs play central functional roles in the pathology of inflammatory and immune disorders. Moreover, these effects would be mainly mediated by the NF-κB, AP-1, focal adhesion, and PI3K-Akt signaling pathways [41,42]. However, the relationships between the lncRNA-mRNA network and rosacea have not been elucidated so far. We applied an integrated bioinformatics approach to establish the ceRNA network, RT-qPCR, and FISH to identify dysregulated lncRNAs that may play important roles in rosacea from the database. logFC S100A9 S100A12 CXC L9 MM P1 CX CL 13 C XC L1 0 C C L 18 C X C L 1 1 JC H A IN S 1 0 0 A 7 C C L 1 9 C D 2 M M P 9 I L 3 6 G S E L L C C L 8 C C L 2 0 C C L 4 C C R 5 S1 00 A 8 miRNAs function in RNA silencing and posttranscriptional regulation of gene expression. Recent advances reveal miRNAs are associated with inflammatory skin disease (e.g., psoriasis) [15]. In the current study, a total of 5 DE-miRNAs were screened out, including hsa-miR-148b-3p, hsa-miR-148a-3p, and hsa-miR-296-3p, which targeted the most  genes. Furthermore, the three miRNAs were significantly enriched in T cell activation and JAK-STAT signaling pathway. Besides, in the present study, we observed an enrichment of lncRNAs in leukocyte migration and cell chemotaxis in the three rosacea subtypes. is indicates an increase in innate immunity and inflammatory reactions in rosacea pathophysiology. e interaction network obtained suggests that T cell activation plays an important role in rosacea development. Similarly, enrichment analyses showed that genes mainly involved in chemotaxis and migration of cells were upregulated (including S100A9, S100A12, CXCL9, MMP1, and CXCL13) and downregulated (including IL37, HOXA7, CCL21, LEP, and AGTR1). ese data indicate the importance of these processes in the pathogenesis of rosacea. In addition, KEGG pathway enrichment analyses reveal that ferroptosis, PPAR, NOD-like, and NF-κB and JAK-STAT signaling pathways are involved in rosacea. Ferroptosis is a distinct form of programmed cell death that is associated with the generation of iron-dependent reactive oxygen species (ROS) [44]. Peroxisome proliferator-activated receptors (PPARs) are ligand-activated transcription factors, expressed primarily in adipose tissue and the immune system that regulate genes in cell differentiation and in multiple metabolic processes [45]. Once activated the PPAR signaling pathway participates in the amplification of the inflammatory process. NOD-like receptors are cytoplasmic receptors that play a vital role in the innate immune system, recognizing molecular patterns associated with pathogens and damage. Activation of the NOD-like signaling pathway induces the innate immune response [46]. NF-κB is a protein complex that performs transcription factor functions, among them acting as a mediator of the inflammatory response regulation. Activation of NF-κB signaling pathway results in proinflammatory cytokine secretion [3]. Finally, the JAK-STAT signaling pathway that transduces intracellular signals is associated with inflammatory and autoimmune diseases [47]. In vitro, Li et al. showed artesunate could treat rosacea via the inhibition of the JAK/STAT signaling pathway [48]. ese evidences strongly indicate lncRNAs play a role in the pathology of rosacea through the five pathways (ferroptosis, PPAR, NOD-like, NF-Κb, and JAK-STAT). In addition, lncRNAs and miRNAs may work together via JAK-STAT pathway.
Next, we showed that DE lncRNA NEAT1 was upregulated, while HOTAIR and ZNF667-AS1 were downregulated. NEAT1 is located on chromosome 11q13.1 and is expressed in various tissues and cell types. NEAT1 can dysregulate proinflammatory chemokines and cytokines in monocytes by mediating the TLR4-mediated inflammatory pathway in patients with SLE [40]. Further, NEAT1 is overexpressed in inflammatory bowel disease patients and participates in the inflammatory response by regulating the intestinal epithelial barrier and exosome-mediated polarization of macrophages [49]. Here, for the first time, it was identified a relationship between the NEAT1 and rosacea. Our data suggest that NEAT1 may cause damage in rosacea by regulating genes associated with immune dysregulation and inflammation, as indicated by NEAT1-regulated genes that are enriched in chemokine activity, which play crucial roles in inflammatory infiltration. HOTAIR is located in the Homeobox C gene cluster on chromosome 12q13.13 and is important for suppressing the proliferation, invasion, migration, and tube formation of vein endothelial cells by downregulating the expression of vascular endothelial growth factor A [50]. Furthermore, HOTAIR is a critical participant in lipopolysaccharide-induced cytokine expression, immune and inflammatory response in macrophages [51].
is lncRNA also relieves inflammation by regulating the expression levels of inflammatory factors and inhibiting the activation of NF-κB pathway [52]. In our study, the downregulated expression of HOTAIR in all subtypes of rosacea indicates that HOTAIR may improve rosacea by alleviating inflammation and inhibiting the NF-kappa B pathway.
ZNF667-AS1, namely Zinc finger protein 667-antisense RNA 1, is located on human chromosome 19q13. 43   10 BioMed Research International [53]. In human normal cells, ZNF667-AS1 is highly expressed while reduced in immortalized cells and multiple tumor tissues [54]. Li et al. illustrated that ZNF667-AS1 could inhibit the inflammatory response and promote spinal cord injury recovery via suppressing the JAK-STAT pathway [55]. Although there is no research about the relationship between ZNF667-AS1 and rosacea, our data indicate that the downregulated ZNF667-AS1 may play an important role in rosacea. Whether ZNF667-AS1 prevents the development of rosacea via JAK-STAT pathway needs deep research.

Conclusions
In conclusion, this is the first study to identify the ceRNA network in rosacea. Pathways including ferroptosis, PPAR signaling pathway, NOD-like signaling pathway, NF-kappa B signaling pathway, and JAK-STAT signaling pathway can be considered important therapeutic targets for rosacea. e lncRNAs and mRNAs identified in our study may play key roles in rosacea development. Upregulated NEAT1 can be a possible biomarker for rosacea diagnosis as well as HOTAIR and ZNF667-AS1 can be a target for rosacea prevention and treatment. More in-depth researches are needed to better characterize the lncRNAs and mRNAs in rosacea.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
Ruifeng Lu is a co-first author.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Supplementary Materials
See Figures S1-S3 in the Supplementary Material for comprehensive image analysis. Figure S1: (a) the enriched bubble diagrams of the GO analysis with differentially expressed genes of the lncRNAs in ETR group. e horizontal axis represents the proportion of those genes accounted for in all the GO annotated genes, the left side of the vertical axis represents the annotation terms, and the right side of the vertical axis represents biological process (BP) terms, cellular component (CC) terms, and molecular function (MF) terms. Bubble scale represents the number of genes in each GO term; depth of bubble color represents p value. (b) e three kinds of enriched GO biological process terms of differentially expressed mRNAs involved in the lncRNAs network in ETR group. Figure S2: (a) the enriched bubble diagrams of the GO analysis with differentially expressed genes of the lncRNAs in PhR group. e horizontal axis represents the proportion of those genes accounted for in all the GO annotated genes, the left side of the vertical axis represents the annotation terms, and the right side of the vertical axis represents biological process (BP) terms, cellular component (CC) terms, and molecular function (MF) terms. Bubble scale represents the number of genes in each GO term; depth of bubble color represents p value. (b) e three kinds of enriched GO biological process terms of differentially expressed mRNAs involved in the lncRNAs network in the PhR group. Figure S3: (a) the enriched bubble diagrams of the GO analysis with differentially expressed genes of the lncRNAs in the PPR group. e horizontal axis represents the proportion of those genes accounted for in all the GO annotated genes, the left side of the vertical axis represents the annotation terms, and the right side of the vertical axis represents biological process (BP) terms, cellular component (CC) terms, and molecular function (MF) terms. Bubble scale represents the number of genes in each GO term; depth of bubble color represents p value. (b) e three kinds of enriched GO biological process terms of differentially expressed mRNAs involved in the lncRNAs network in the PPR group. (Supplementary Materials)