Identification of Long Noncoding RNA Associated ceRNA Networks in Rosacea

Wang, Lian; Lu, Ruifeng; Wang, Yujia; Wang, Xiaoyun; Hao, Dan; Wen, Xiang; Li, Yanmei; Zeng, Minghui; Jiang, Xian

doi:https://doi.org/10.1155/2020/9705950

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Disclosure Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2020 | Article ID 9705950 | https://doi.org/10.1155/2020/9705950

Identification of Long Noncoding RNA Associated ceRNA Networks in Rosacea

Lian Wang,¹Ruifeng Lu,²Yujia Wang,¹Xiaoyun Wang,¹Dan Hao,¹Xiang Wen,¹Yanmei Li,¹Minghui Zeng,³and Xian Jiang¹

Academic Editor: Nobuo Kanazawa

Received12 Nov 2019

Revised24 Jan 2020

Accepted27 Jan 2020

Published24 Feb 2020

Abstract

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.

1. 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]. The 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 Th1/Th17-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–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. The present study provides a novel and comprehensive overview of the functions and mechanisms of differentially expressed lncRNAs in rosacea pathophysiology.

2. Materials and Methods

2.1. Microarray Data Process and LncRNA Reannotation

Rosacea gene expression data were obtained from Gene Expression Omnibus (GEO, GSE65914) [5].

There 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. These 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–23]. The retained probes should match the chromosomal position of corresponding lncRNAs retrieved from the GENCODE project (human GENCODE Release 29, GRCh38.p12). The 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].

2.2. 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 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].

2.3. 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. The “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.

2.4. Construction of lncRNA-miRNA-mRNA Networks (ceRNA Networks)

The miRcode (http://www.mircode.org/) database was used to predict lncRNA-miRNA interactions, and the miRNAs sequences were modified by StarBase v2.0 database (http://starbase.sysu.edu.cn/) and MicroRNA ENrichment TURned NETwork (MIENTURNET, http://userver.bio.uniroma1.it/apps/mienturnet/) [28]. miRNA targeted mRNAs were obtained from miRDB (http://mirdb.org), miRTarBase (http://miRTarBase.mbc.nctu.edu.tw/), and TargetScan databases (http://www.targetscan.org/vert_71/) [29, 30]. Only mRNAs recognized by all three databases were considered as candidate mRNAs and intersected with the DE mRNAs to screen out the DE mRNAs targeted by the DE miRNAs. Significantly correlated pairs of DE lncRNA-miRNA and miRNA-mRNA interactions were included to create the coexpression network by Cytoscape 3.6.1.

2.5. Patients and Specimens

This 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.

2.6. Reverse Transcription-Quantitative PCR (RT-qPCR) Assay

Total RNA was extracted from rosacea patients using TRIzol reagent (Invitrogen, USA) and synthesized into cDNA using Reverse Transcription Kits (Applied Biosystems, USA). Quantitative PCR was performed by SYBR PCR kits (Bio-Rad) on Real-Time PCR Detection System (Bio-Rad). RNA concentration was measured by Nanodrop (Invitrogen, USA). Relative lncRNA (HOTAIR, ZNF667-AS1, and NEAT1) levels were normalized against Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and calculated using the 2^–ΔΔCt methods. Each sample was analyzed in triplicate. The primer sequences were as follows: GAPDH-F: CGCTGAGTACGTCGTGGAGTC GAPDH-R: GCTGATGATCTTGAGGCTGTTGTC NEAT1-F: CACAGGCAGGGGAAATGTCT NEAT1-R: TGCTGCGTATGCAAGTCTGA HOTAIR-F: AGAGAAAAGGCTGAAATGGAGGAC HOTAIR-R: AGCTTGGGTGTAATTGCTGGTTT ZNF667-AS1-F: TCTCACCTTGTTTTGTCCTGT ZNF667-AS1-R: TCGTTTAATGGGTAGAGTTTC

2.7. Fluorescence In Situ Hybridization (FISH)

HOTAIR, ZNF667-AS1 and NEAT1 probes were designed and synthesized by RiboBio (Guangzhou, China), and the probe sequences are available upon request

The probe signals were detected with a FISH Kit (RiboBio, Guangzhou, China) according to the manufacturer’s instructions. Briefly, rearrangements were analyzed on 5 μm formalin fixed and paraffin embedded (FFPE) sections. Images were captured using a fluorescence microscope (Leica, Germany). At least 50 nuclei were scored for each probe and case.

2.8. Statistical Analysis

All statistical analyses were conducted on Graph Prism 7. The 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. was considered as significant.

3. Results

3.1. 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1

LncRNA and mRNA expression profile of normal group, ETR, PhR, and PPR groups. (a) Heat map of lncRNAs expression in normal group, ETR, PhR, and PPR groups. Red represents upregulated lncRNAs and green represents downregulated lncRNAs. (b) Heat map of mRNAs expression in normal group, ETR, PhR, and PPR groups. Red represents upregulated mRNAs and green represents downregulated mRNAs. (c) Volcano plots of lncRNAs for normal group, ETR, PhR, and PPR groups. The horizontal axis represents fold change (log 2) and the vertical axis is value (–log 10). Red points (fold change > 1) indicate upregulated lncRNAs and blue points (fold change < –1) indicate downregulated lncRNAs. (d) Hierarchical clustered heat maps showing the log 10 transformed expression values for differently expressed lncRNAs among normal group, ETR, PhR, and PPR groups. The horizontal axis represents samples, while the vertical axis represents the biological elements studied. Red represents higher expression and blue represents lower expression. (e) Volcano plots of mRNAs for normal group, ETR, PhR, and PPR groups. The horizontal axis is fold change (log 2) and the vertical axis is value (–log 10). Red points (fold change > 1) indicate upregulated mRNAs and blue points (fold change < –1) indicate downregulated mRNAs. (f) Hierarchical clustered heat maps showing the log10 transformed expression values for differently expressed mRNAs among normal group, ETR, PhR, and PPR groups. Red represents higher expression and blue represents lower expression.

The 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)). The 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. The 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.

(a)

(b)

(c)

(d)

(e)

Figure 2

(a) The overlapping Venn diagram of lncRNAs among ETR, PhR, and PPR groups. There are 6 upregulated lncRNAs and 3 kinds of downregulated lncRNAs. (b) The overlapping Venn diagram of mRNAs among ETR, PhR, and PPR groups. (c) Gene ontology (GO) analysis of differentially expressed lncRNAs in rosacea patients and normal group. The 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 value. (d) The annotation terms are displayed as an interaction network using the BinGO plug-in for Cytoscape. Bubble scale represents number of genes; depth of bubble color represents value. (e) The enriched GO biological process terms of differentially expressed mRNAs involved in the lncRNAs network.

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.

3.2. 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)). The 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.

(a)

(b)

(c)

(d)

(e)

Figure 3

(a) Clustering dendrogram of genes. The dissimilarity of genes is based on the topological overlap. The genes are assigned to different modules and are identified using different colors. (b) Interactions between genes in the coexpression modules. The different colors on the horizontal and vertical axis represent different groups and modules. The colors in the middle represent the relativity among each module. (c) Gene set enrichment analysis (GSEA) in differentially expressed mRNAs linked to lncRNAs in the ETR group. Ferroptosis regulated by GPX4, MAP1LC3B, VDAC3 and PPAR signaling pathway regulated by PPARG, CD36, MMP1 were involved. The horizontal axis represents the rank in all the ordered dataset, and the vertical axes represent enrichment score and ranked list metric. (d) Ferroptosis regulated by GPX4, MAP1LC3B, VDAC3, NOD-like receptor signaling pathway regulated by STAT1, IL1B, NLRX1, IL6, VDAC3, and PPAR signaling pathway regulated by PPARG, CD36, MMP1 were involved in PhR group. (e) NF-kappa B signaling pathway regulated by IL1B, TNFSF13B, VCAM1, NOD-like receptor signaling pathway regulated by STAT1, IL1B, IL6, NLRX1, and PPAR signaling pathway regulated by MMP1, PPARG, CD36 were involved in PPR group.

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. The core regulators on ferroptosis included GPX4, MAP1LC3B, and VDAC3 genes. On the other hand, the core regulators on PPAR 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. The 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, NOD-like, and PPAR signaling pathways. IL1B, TNFSF13B, and VCAM1 genes were regulators of NF-κB signaling pathway.

3.3. Construction of the ceRNA Network and GO Analyses and KEGG Pathway

LncRNA-miRNA interactions and miRNA-mRNA interactions were combined to establish a complete lncRNA-miRNA-mRNA network, which consisted of 5 lncRNAs, 5 miRNAs, and 56 mRNAs, totaling to 241 interactions (Figure 5(a)). LncRNAs (NEAT1, HOTAIR, and ZNF667-AS1) acted as ceRNAs to target miRNAs (hsa-miR-148b-3p, hsa-miR-148a-3p, hsa-miR-296-3p, hsa-miR-378g, and hsa-miR-152-3p), which play important roles in this network. Upregulation of mRNAs (e.g., CCL19, IL21R, WNT2B, STAT3, ICAM1, GLUL) and downregulation of mRNAs (e.g., IL20RB, KLF4, SH3PXD2A) had strong connections with these miRNAs in the network. The GO analyses revealed T cell activation, cell response to an external stimulus, and response to alcohol (Figure 5(b)). The KEGG pathway revealed JAK-STAT signaling pathway (Figure 5(c)). The interaction network on GO terms illustrated major GO terms were interacted with various mRNAs (Figure 5(d)).

(a)

(b)

(c)

(d)

3.4. RT-qPCR and FISH

The 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)).

(a)

(b)

4. Discussion

The 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]. The 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–39]. Recently, lncRNAs broad expression in cutaneous disorders including systemic lupus erythematosus (SLE), sensitive skin, psoriasis, and aging process has been documented [40–43]. These 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.

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. This indicates an increase in innate immunity and inflammatory reactions in rosacea pathophysiology. The 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). These 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]. These 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]. This 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 and is associated with cell proliferation, migration, and invasion [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.

5. 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. The 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

The 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

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Lian Wang and Ruifeng Lu contributed equally to this work.

Acknowledgments

The authors would like to sincerely thank their colleagues for their comments on this article. This research was funded by the National Natural Science Foundation of China, grants nos. 81872535 and 21772131.

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. The 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 value. (b) The 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. The 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) The 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. The 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) The three kinds of enriched GO biological process terms of differentially expressed mRNAs involved in the lncRNAs network in the PPR group. (Supplementary Materials)

References

E. J. van Zuuren, “Rosacea,” New England Journal of Medicine, vol. 377, no. 18, pp. 1754–1764, 2017.
View at: Publisher Site | Google Scholar
R. L. Gallo, R. D. Granstein, S. Kang et al., “Standard classification and pathophysiology of rosacea: the 2017 update by the national rosacea society expert committee,” Journal of the American Academy of Dermatology, vol. 78, no. 1, pp. 148–155, 2018.
View at: Publisher Site | Google Scholar
A. D. Holmes and M. Steinhoff, “Integrative concepts of rosacea pathophysiology, clinical presentation and new therapeutics,” Experimental Dermatology, vol. 26, no. 8, pp. 659–667, 2017.
View at: Publisher Site | Google Scholar
M. Steinhoff, J. Buddenkotte, J. Aubert et al., “Clinical, cellular, and molecular aspects in the pathophysiology of rosacea,” Journal of Investigative Dermatology Symposium Proceedings, vol. 15, no. 1, pp. 2–11, 2011.
View at: Publisher Site | Google Scholar
T. Buhl, M. Sulk, P. Nowak et al., “Molecular and morphological characterization of inflammatory infiltrate in rosacea reveals activation of Th1/Th17 pathways,” Journal of Investigative Dermatology, vol. 135, no. 9, pp. 2198–2208, 2015.
View at: Publisher Site | Google Scholar
Y. G. Chen, A. T. Satpathy, and H. Y. Chang, “Gene regulation in the immune system by long noncoding RNAs,” Nature Immunology, vol. 18, no. 9, pp. 962–972, 2017.
View at: Publisher Site | Google Scholar
M. Piipponen, J. Heino, V. M. Kahari, and L. Nissinen, “Long non-coding RNA PICSAR decreases adhesion and promotes migration of squamous carcinoma cells by downregulating alpha2beta1 and alpha5beta1 integrin expression,” Biology Open, vol. 7, no. 11, 2018.
View at: Publisher Site | Google Scholar
Y. Ding, X. Shao, X. Li et al., “Identification of candidate genes in atopic dermatitis based on bioinformatic methods,” International Journal of Dermatology, vol. 55, no. 7, pp. 791–800, 2016.
View at: Publisher Site | Google Scholar
R. Gupta, R. Ahn, K. Lai et al., “Landscape of long noncoding RNAs in psoriatic and healthy skin,” Journal of Investigative Dermatology, vol. 136, no. 3, pp. 603–609, 2016.
View at: Publisher Site | Google Scholar
R. Ahn, R. Gupta, K. Lai et al., “Network analysis of psoriasis reveals biological pathways and roles for coding and long non-coding RNAs,” BMC Genomics, vol. 17, no. 1, p. 841, 2016.
View at: Publisher Site | Google Scholar
D. Lei, L. Lv, L. Yang et al., “Genome-wide analysis of mRNA and long noncoding RNA profiles in chronic actinic dermatitis,” BioMed Research International, vol. 2017, Article ID 7479523, 15 pages, 2017.
View at: Publisher Site | Google Scholar
M. Huarte, “The emerging role of lncRNAs in cancer,” Nature Medicine, vol. 21, no. 11, pp. 1253–1261, 2015.
View at: Publisher Site | Google Scholar
F. Russo, G. Fiscon, F. Conte, M. Rizzo, P. Paci, and M. Pellegrini, “Interplay between long noncoding RNAs and MicroRNAs in cancer,” Methods in Molecular Biology, vol. 1819, pp. 75–92, 2018.
View at: Publisher Site | Google Scholar
G. J. Yu, Y. Sun, D. W. Zhang, and P. Zhang, “Long non-coding RNA HOTAIR functions as a competitive endogenous RNA to regulate PRAF2 expression by sponging miR-326 in cutaneous squamous cell carcinoma,” Cancer Cell International, vol. 19, p. 270, 2019.
View at: Publisher Site | Google Scholar
Q. Liu, D.-H. Wu, L. Han et al., “Roles of microRNAs in psoriasis: immunological functions and potential biomarkers,” Experimental Dermatology, vol. 26, no. 4, pp. 359–367, 2017.
View at: Publisher Site | Google Scholar
Q. Zhou, Q. Yu, Y. Gong et al., “Construction of a lncRNA-miRNA-mRNA network to determine the regulatory roles of lncRNAs in psoriasis,” Experimental and Therapeutic Medicine, vol. 18, no. 5, pp. 4011–4021, 2019.
View at: Google Scholar
L. Poliseno, L. Salmena, J. Zhang, B. carver, W. J. Haveman, and P. P. Pandolfi, “A coding-independent function of gene and pseudogene mRNAs regulates tumour biology,” Nature, vol. 465, no. 7301, pp. 1033–1038, 2010.
View at: Publisher Site | Google Scholar
L. Salmena, L. Poliseno, Y. Tay, L. Kats, and P. P. Pandolfi, “A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?” Cell, vol. 146, no. 3, pp. 353–358, 2011.
View at: Publisher Site | Google Scholar
W. Yong, D. Yu, Z. Jun et al., “Long noncoding RNA NEAT1, regulated by LIN28B, promotes cell proliferation and migration through sponging miR-506 in high-grade serous ovarian cancer,” Cell Death & Disease, vol. 9, no. 9, p. 861, 2018.
View at: Publisher Site | Google Scholar
Q. Liao, C. Liu, X. Yuan et al., “Large-scale prediction of long non-coding RNA functions in a coding-non-coding gene co-expression network,” Nucleic Acids Research, vol. 39, no. 9, pp. 3864–3878, 2011.
View at: Publisher Site | Google Scholar
Z. Du, T. Fei, R. G. W. Verhaak et al., “Integrative genomic analyses reveal clinically relevant long noncoding RNAs in human cancer,” Nature Structural & Molecular Biology, vol. 20, no. 7, pp. 908–913, 2013.
View at: Publisher Site | Google Scholar
C. Qian, H. Li, D. Chang, B. Wei, and Y. Wang, “Identification of functional lncRNAs in atrial fibrillation by integrative analysis of the lncRNA‐mRNA network based on competing endogenous RNAs hypothesis,” Journal of Cellular Physiology, vol. 234, no. 7, pp. 11620–11630, 2019.
View at: Publisher Site | Google Scholar
Y. Zhang, L. Jiao, L. Sun et al., “LncRNA ZFAS1 as a SERCA2a inhibitor to cause intracellular Ca²⁺ overload and contractile dysfunction in a mouse model of myocardial infarction,” Circulation Research, vol. 122, no. 10, pp. 1354–1368, 2018.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, p. e47, 2015.
View at: Publisher Site | Google Scholar
B. Zhu, M. Xu, H. Shi, X. Gao, and P. Liang, “Genome-wide identification of lncRNAs associated with chlorantraniliprole resistance in diamondback moth Plutella xylostella (L.),” BMC Genomics, vol. 18, no. 1, p. 380, 2017.
View at: Publisher Site | Google Scholar
G. Yu, L.-G. Wang, Y. Han, and Q.-Y. He, “Clusterprofiler: an R package for comparing biological themes among gene clusters,” OMICS: A Journal of Integrative Biology, vol. 16, no. 5, pp. 284–287, 2012.
View at: Publisher Site | Google Scholar
M. Kanehisa, M. Araki, S. Goto et al., “KEGG for linking genomes to life and the environment,” Nucleic Acids Research, vol. 36, no. 36, pp. D480–D484, 2008.
View at: Publisher Site | Google Scholar
V. Licursi, F. Conte, G. Fiscon, and P. Paci, “MIENTURNET: an interactive web tool for microRNA-target enrichment and network-based analysis,” BMC Bioinformatics, vol. 20, no. 1, p. 545, 2019.
View at: Publisher Site | Google Scholar
Y. Chen and X. Wang, “miRDB: an online database for prediction of functional microRNA targets,” Nucleic Acids Research, vol. 48, no. D1, 2020.
View at: Publisher Site | Google Scholar
C.-H. Chou, S. Shrestha, C.-D. Yang et al., “miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions,” Nucleic Acids Research, vol. 46, no. D1, pp. D296–d302, 2018.
View at: Publisher Site | Google Scholar
B. M. Rainer, S. Kang, and A. L. Chien, “Rosacea: epidemiology, pathogenesis, and treatment,” Dermato-endocrinology, vol. 9, no. 1, Article ID e1361574, 2017.
View at: Publisher Site | Google Scholar
K. Yamasaki, K. Kanada, D. T. Macleod et al., “TLR2 expression is increased in rosacea and stimulates enhanced serine protease production by keratinocytes,” Journal of Investigative Dermatology, vol. 131, no. 3, pp. 688–697, 2011.
View at: Publisher Site | Google Scholar
K. Yamasaki, A. DI Nardo, A. Bardan et al., “Increased serine protease activity and cathelicidin promotes skin inflammation in rosacea,” Nature Medicine, vol. 13, no. 8, pp. 975–980, 2007.
View at: Publisher Site | Google Scholar
J. Tschopp and K. Schroder, “NLRP3 inflammasome activation: the convergence of multiple signalling pathways on ROS production?” Nature Reviews Immunology, vol. 10, no. 3, pp. 210–215, 2010.
View at: Publisher Site | Google Scholar
A. Weber, P. Wasiliew, and M. Kracht, “Interleukin-1 (IL-1 ) processing pathway,” Science Signaling, vol. 3, no. 105, p. cm2, 2010.
View at: Publisher Site | Google Scholar
L. Tang, W. Zhang, B. Su, and B. Yu, “Long noncoding RNA HOTAIR is associated with motility, invasion, and metastatic potential of metastatic melanoma,” BioMed Research International, vol. 2013, Article ID 251098, 7 pages, 2013.
View at: Publisher Site | Google Scholar
M. Li, J. Wang, D. W. Liu, and H. Huang, “High-throughput sequencing reveals differentially expressed lncRNAs and circRNAs, and their associated functional network, in human hypertrophic scars,” Molecular Medicine Reports, vol. 18, no. 6, pp. 5669–5682, 2018.
View at: Publisher Site | Google Scholar
Y. Xia, Z. Q. Yan, Y. Wan et al., “Knockdown of long noncoding RNA GHET1 inhibits cell-cycle progression and invasion of gastric cancer cells,” Molecular Medicine Reports, vol. 18, no. 3, pp. 3375–3381, 2018.
View at: Publisher Site | Google Scholar
Y. H. Zhang, H. M. Li, W. Y. Zhang et al., “LASSO-based Cox-PH model identifies an 11-lncRNA signature for prognosis prediction in gastric cancer,” Molecular Medicine Reports, vol. 18, pp. 5579–5593, 2018.
View at: Google Scholar
F. Zhang, L. Wu, J. Qian et al., “Identification of the long noncoding RNA NEAT1 as a novel inflammatory regulator acting through MAPK pathway in human lupus,” Journal of Autoimmunity, vol. 75, pp. 96–104, 2016.
View at: Publisher Site | Google Scholar
L. Yang, L. Lyu, W. Wu et al., “Genome-wide identification of long non-coding RNA and mRNA profiling using RNA sequencing in subjects with sensitive skin,” Oncotarget, vol. 8, no. 70, pp. 114894–114910, 2017.
View at: Publisher Site | Google Scholar
K. Yo and T. M. Rünger, “The long non-coding RNA FLJ46906 binds to the transcription factors NF-κB and AP-1 and regulates expression of aging-associated genes,” Aging, vol. 10, no. 8, pp. 2037–2050, 2018.
View at: Publisher Site | Google Scholar
B. Ke, M. Tian, J. Li, B. Liu, and G. He, “Targeting programmed cell death using small-molecule compounds to improve potential cancer therapy,” Medicinal Research Reviews, vol. 36, no. 6, pp. 983–1035, 2016.
View at: Publisher Site | Google Scholar
S. J. Dixon, K. M. Lemberg, M. R. Lamprecht et al., “Ferroptosis: an iron-dependent form of nonapoptotic cell death,” Cell, vol. 149, no. 5, pp. 1060–1072, 2012.
View at: Publisher Site | Google Scholar
B. Grygiel-Gorniak, “Peroxisome proliferator-activated receptors and their ligands: nutritional and clinical implications--a review,” Nutrition Journal, vol. 13, p. 17, 2014.
View at: Publisher Site | Google Scholar
Y. K. Kim, J.-S. Shin, and M. H. Nahm, “NOD-like receptors in infection, immunity, and diseases,” Yonsei Medical Journal, vol. 57, no. 1, pp. 5–14, 2016.
View at: Publisher Site | Google Scholar
S. Banerjee, A. Biehl, M. Gadina, S. Hasni, and D. M. Schwartz, “JAK-STAT signaling as a target for inflammatory and autoimmune diseases: current and future prospects,” Drugs, vol. 77, no. 5, pp. 521–546, 2017.
View at: Publisher Site | Google Scholar
T. Li, Q. Zeng, X. Chen et al., “The therapeutic effect of artesunate on rosacea through the inhibition of the JAK/STAT signaling pathway,” Molecular Medicine Reports, vol. 17, pp. 8385–8390, 2018.
View at: Publisher Site | Google Scholar
R. Liu, A. Tang, X. Wang et al., “Inhibition of lncRNA NEAT1 suppresses the inflammatory response in IBD by modulating the intestinal epithelial barrier and by exosome-mediated polarization of macrophages,” International Journal of Molecular Medicine, vol. 42, no. 5, pp. 2903–2913, 2018.
View at: Google Scholar
K. Wu, F. Liu, W. Wu, Y. Chen, H. Wu, and W. Zhang, “Long non-coding RNA HOX transcript antisense RNA (HOTAIR) suppresses the angiogenesis of human placentation by inhibiting vascular endothelial growth factor A expression,” Reproduction, Fertility, and Development, vol. 31, no. 2, p. 377, 2019.
View at: Publisher Site | Google Scholar
M. Obaid, S. M. N. Udden, P. Deb et al., “LncRNA HOTAIR regulates lipopolysaccharide-induced cytokine expression and inflammatory response in macrophages,” Scientific Reports, vol. 8, no. 1, p. 15670, 2018.
View at: Publisher Site | Google Scholar
J. L. Pang, J. W. Wang, P. Y. Hu, J. S. Jiang, and C. Yu, “HOTAIR alleviates ox-LDL-induced inflammatory response in Raw264.7 cells via inhibiting NF-kappaB pathway,” European Review for Medical and Pharmacological Sciences, vol. 22, no. 20, pp. 6991–6998, 2018.
View at: Google Scholar
W. Meng, W. Cui, L. Zhao et al., “Aberrant methylation and downregulation of ZNF667-AS1 and ZNF667 promote the malignant progression of laryngeal squamous cell carcinoma,” Journal of Biomedical Science, vol. 26, no. 1, p. 13, 2019.
View at: Publisher Site | Google Scholar
L. Vrba, J. C. Garbe, M. R. Stampfer, and B. W. Futscher, “A lincRNA connected to cell mortality and epigenetically-silenced in most common human cancers,” Epigenetics, vol. 10, no. 11, pp. 1074–1083, 2015.
View at: Publisher Site | Google Scholar
J. W. Li, Y. Kuang, L. Chen, and J. F. Wang, “LncRNA ZNF667-AS1 inhibits inflammatory response and promotes recovery of spinal cord injury via suppressing JAK-STAT pathway,” European Review for Medical and Pharmacological Sciences, vol. 22, no. 22, pp. 7614–7620, 2018.
View at: Google Scholar

Copyright

Copyright © 2020 Lian Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1186

Downloads

1373

Citations