Identification of a Potentially Functional circRNA-miRNA-mRNA Regulatory Network in Melanocytes for Investigating Pathogenesis of Vitiligo

CircRNAs have been reported to play essential roles in regulating immunity and inflammation, which may be an important regulatory factor in the development of vitiligo. However, the expression profile of circRNAs and their potential biological functions in vitiligo have not been reported so far. In our study we found there are 64 dysregulated circRNAs and 14 dysregulated miRNAs in the patients with vitiligo. Through the correlation analysis, we obtained 12 dysregulated circRNAs and 5 dysregulated miRNAs, forming 48 relationships in the circRNA-miRNA-mRNA regulatory network. Gene Ontology analysis indicated dysregulated circRNAs in vitiligo is closely related to the disorder of the metabolic pathway. The KEGG pathway of dysregulation of circRNAs mainly enriched in the biological processes such as ubiquitin mediated proteolysis, endocytosis and RNA degradation, and in Jak-STAT signaling pathway. Therefore, we found the circRNA-miRNA-mRNA regulatory network are involved in the regulation of numerous melanocyte functions, and these dysregulated circRNAs may closely related to the melanocyte metabolism. Our study provides a theoretical basis for studying the vitiligo pathogenesis from the perspective of circRNA-miRNA-mRNA network.


INTRODUCTION
Vitiligo is characterized by the progressive disappearance of melanocytes, resulting in depigmentation of the skin and/or hair (Malhotra and Dytoc, 2013). The exact pathogenesis remains unknown. Several hypotheses about the pathogenesis mechanisms have been implicated, including genetic, autoimmune, biochemical factors, oxidative stress, and neural or viral causes (Iannella et al., 2016).
Circular RNAs (circRNAs) are a type of covalently closed, single-stranded circular non-coding RNA (Chen and Yang, 2015). Accumulating evidence has shown that circRNAs are widespread and diverse throughout eukaryotic cells and that they have multiple biological functions and play potentially important roles in various diseases (Lei et al., 2018). CircRNAs play an important role in regulating cell proliferation, inflammation and immunity (Chen et al., 2019b;Liu et al., 2019). Therefore, circRNAs are considered to be an important regulator of the onset and development of vitiligo and may be a new therapeutic target for vitiligo. Circular RNAs are rich in miRNA binding sites (miRNA response elements, MREs), serving as miRNA sponges or competitive endogenous RNAs (ceRNAs). CircRNA regulates mRNA at the post-transcriptional level by completely binding to shared miRNAs (Kristensen et al., 2019). However, circRNAs as ceRNAs mediating pathological processes has not been reported in vitiligo.
Studies on circRNA regulation in vitiligo are limited, the expression profile of circRNAs and their potential biological functions in vitiligo have not been reported so far. In our study, we used RNA-seq to investigated comprehensive circRNA expression profile between vitiligo lesions and normal healthy skin tissue, further to explore the possible regulated mechanism of circRNAs in the pathogenesis of vitiligo.

Patients and Tissue Specimens
Fresh tissues were obtained from patients with vitiligo and normal donors, from April 2019 to April 2020. The study was approved by the committee of People's Hospital of Guangxi Zhuang Autonomous Region and all subjects signed an informed consent form under the code KY-ZGR-2019-013. This study included a total of 6 adult patients and 6 control subjects, male and female adults between the ages of 30-60. The patients with vitiligo were diagnosed with vitiligo vulgaris, defined as a condition with depigmentation between 10% and 80% of the body surface. The detail information was shown in Supplementary  Table 1.
Skin tissue obtained by skin punch biopsy (4 mm diameter) of the normal areas from normal donors and the lesions from vitiligo patients. Two skin biopsies were obtained from each patient with vitiligo and control subjects. The biopsy was obtained from the central part of the affected skin from patients with vitiligo. Control skin biopsy was taken from healthy control subjects at the same part vs. patient with vitiligo.

RNA Extraction and Qualification?
Total RNA was extracted from each skin biopsies by RNAprep Pure Tissue Kit (TIANGEN BIOTECH, Beijing, China) in accordance with the manufacturer's instructions.
We monitored RNA degradation and contamination used 1% agarose gels. And using the NanoPhotometer R spectrophotometer (IMPLEN, CA, United States) to check RNA purity. Then we measured RNA concentration by Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer (Life Technologies, CA, United States) and assessed RNA integrity by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States).

RNA-Seq
Details of the circRNA-seq and miRNA-seq methods are described in Supplementary Materials. Briefly describe as that:

CircRNA-Seq
A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext R UltraTM RNA Library Prep Kit for Illumina R (NEB, United States) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. In order to select cDNA fragments of preferentially 150∼200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, United States). Then 3 µl USER Enzyme (NEB, United States) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated. Index of the reference genome was built using STAR and paired-end clean reads were aligned to the reference genome using STAR (v2.5.1b). The circRNA were detected and identified using find_circ and CIRI2 (Memczak et al., 2013;Gao et al., 2018). HTSeq v0.6.0 was used to count the reads numbers mapped to each gene. Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.10.1). Genes with an adjusted P-value < 0.05 found by DESeq2 were assigned as differentially expressed.

MiRNA-Seq
A total amount of 3 µg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB, United States) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. The clustering of the indexcoded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50 bp single-end reads were generated. The small RNA tags were mapped to reference sequence by Bowtie without mismatch to analyze their expression and distribution on the reference (Langmead et al., 2009). miRBase20.0 was used as reference, modified software mirdeep2 and srna-tools-cli were used to obtain the potential miRNA and draw the secondary structures. To remove tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker, Rfam database or those types of data from the specified species itself. The characteristics of hairpin structure of miRNA precursor can be used to predict novel miRNA. The available software miREvo and mirdeep2 were integrated to predict novel miRNA through exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps. miRNA expression levels were estimated by TPM (transcript per million) through the following criteria: Normalization formula: Normalized expression = mapped readcount/Total reads * 1000000. Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.8.3). The P-values was adjusted using the Benjamini   and Hochberg method. Corrected P-value of 0.05 was set as the threshold for significantly differential expression by default.

Real-Time qPCR
To validate the RNA-Seq data, we randomly selected 2 of circRNA and 2 of miRNA for qRT-PCR analysis. Total RNA was extracted from each skin tissue sample, and then reversetranscribed into cDNA using PrimeScript TM RT reagent Kit with gDNA Eraser (Takara, Dalian, China) according to the manufacturer's instruction. Quantitative PCR (qPCR) was performed using the SYBR R Premix Ex TaqTM II (Tli RNase H Plus) Kit with a Bio-Rad CFX Manager 3.1 real-time PCR system (CFX96 TM Real-Time PCR, Bio-Rad, United States). The circRNA expression were normalized to GAPDH as an endogenous reference transcript, the miRNA expression were normalized to that of U6 (Zhong et al., 2019). The relative expression levels of circRNA and miRNA were calculated using the 2 − Ct method. The specific primers for each gene are listed in Supplementary Table 2. Data shown represent the means of 3 experiments.

GO Annotations and KEGG Pathway Analyses
Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the clusterProfiler R package, in which gene length bias was corrected. GO terms with corrected P-value less than 0.05 were considered significantly enriched by differential expressed genes.
KEGG is a database resource for understanding highlevel functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecularlevel information, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies 1 . We used clusterProfiler R package to test the statistical enrichment of differential expression genes in KEGG pathways. 1 http://www.genome.jp/kegg/   Annotation for circRNA-miRNA-mRNA Interaction Predicting circRNA-miRNA interactions and the target gene of miRNA were performed by TargetScan 2 and miRanda 3 . Upregulation and downregulation of circRNA and miRNA were identified. Then, the circRNA-ceRNA regulatory networks were constructed through the Cytoscape software V2.7.0 (San Diego, CA, United States).

Differential Expression of circRNAs and miRNAs in Vitiligo Patients and Control Subjects
We identified 64 significantly dysregulated circRNA transcripts between the two groups in skin tissue of subjects: 34 circRNAs were upregulated and 30 were downregulated in vitiligo patients relative to the control subjects, respectively ( Figure 1A and Table 1). We conducted the cluster analysis of the circRNAs expression, and the heatmap results were shown in Figure 1B.
We obtained 14 significantly dysregulated miRNAs: inducing 6 miRNAs were upregulated and 8 were downregulated in skin tissue of vitiligo patients relative to the control subjects ( Figure 1C and Table 2). We conducted the cluster analysis of the miRNAs expression, and the heatmap results were shown in Figure 1D.
The data showed that these differentially expressed circRNAs were scattered different chromosomes: 34 upregulated circRNAs

Confirmation the Differential Expression of circRNAs and miRNAs
We sought to confirm the differential expression identified in our RNA-seq experiments using qPCR. The differentially expressed transcripts of circRNAs and miRNAs were randomly selected. As shown in Figure 2, all selected transcripts were detected in skin tissue of vitiligo patients and the control subjects and exhibited significant differential expressions.

Go and KEGG Analyses for the Differential Expressed circRNAs
We performed GO analyses on these dysregulated circRNAs in the networks and the results of the top highly enriched GO terms of cellular component (CC), biological process (BP), and molecular function (MF) are shown in Figure 3.
The top terms were molecular function (GO:0003674), cellular process (GO:0009987), and cell (GO:0005623). Several metabolism-associated terms were also observed, such as metabolic process (GO:0008152), organic substance metabolic process (GO:0071704), cellular metabolic process (GO:0044237), regulation of biological process (GO:0050789), cellular macromolecule metabolic process (GO:0044260), Frontiers in Genetics | www.frontiersin.org and regulation of metabolic process (GO:0019222). Thus, we hypothesis that the pathogenesis of vitiligo may be closely related to many metabolic-related pathways regulated by circRNAs. We further conducted KEGG pathway analysis to determine the signaling cascades related to dysregulated circRNAs by using the Q-value scale from 0 to 1, and the top 5 significantly enriched pathways were identified: synthesis and degradation of ketone bodies, terpenoid backbone biosynthesis, butanoate metabolism, propanoate metabolism, and citrate cycle (TCA cycle) (Figure 4).

Construction of a circRNA-miRNA-mRNA Regulatory Network
According to ceRNA hypothesis, RNA transcripts effectively communicate with one another. The members of ceRNA compete for the same miRNA response elements (MREs) to regulate one another. We have conducted the regulatory network analysis of circRNA, miRNA, and their target genes. These data show that there are significant differential expression between vitiligo patients and the control subjects (corrected p < 0.05). Analysis based on the RNA-seq data, we selected 12 circRNAs and 5 miRNAs that were differentially expressed, containing 48 relationships in the circRNA-miRNA-mRNA regulatory network ( Figure 5A and Table 3).
The ceRNA network covered two cases: one was circRNA (upregulated in vitiligo)-miRNA (downregulated in vitiligo)-mRNA (upregulated in vitiligo), and the other one was circRNA (downregulated)-miRNA (upregulated)-mRNA (upregulated). Among them, there are 46 relationships in the manner of up-down-up, occupying the leading position, and only 2 relationships in the manner of down-up-down.
These RNA interactions may supply a novel perspective for the pathogenesis of vitiligo. We observed that one miRNA could be regulated by multiple circRNAs, further regulated a variety of downstream genes that are closely related to melanocyte proliferation or melanin synthesis, such as SEMA4D (ENSG00000187764), RPL22 (ENSG00000116251), ATG13 (ENSG00000175224), and ANKRD6 (ENSG00000135299). For example, hsa-miR-149-5p was co-related with hsa_circ_0007716, novel_circ_0001407, hsa_circ_0030716 and hsa_circ_0001377, further to regulate SEMA4D, ATG13, and ANKRD6.
Furthermore, we detected the expression of target genes such ANKRD6, ATG13, SEMA4D, SMAD3, and PAXILLIN to verify that these circRNA-ceRNA networks are involved in the pathological process of vitiligo. As the results shown that the expression of ANKRD6, ATG13, SEMA4D, SMAD3, and PAXILLIN were all significantly downregulated in vitiligo ( Figure 5B). And, only the expression change of PAXILLIN conforms to the profile of circRNA-miRNA-mRNA. Therefore, we hypothesize that circ_0087961-miR-27a-3p-PAXILLIN might play an important regulatory role in the pathology of vitiligo.

DISCUSSION
Vitiligo is a common chronic acquired disease of pigmentation. Recent clinical and experimental studies suggest that the main mechanism leading to vitiligo involves in autoimmune dysfunction, neurotoxic factor destruction, oxidative stress and genetic factors (Iannella et al., 2016). These factors further lead to loss of melanocyte adhesion and destruction of melanocytes, thereby causing pigmentation loss (Speeckaert and van Geel, 2017;Delmas and Larue, 2019). Vitiligo is a complex genetic disease. Fifty genes at least have already been evaluated in order to identify a link with vitiligo (Sharma et al., 2017;Yuan et al., 2019). Recently studies have suggested that the occurrence and progression of vitiligo are due to multiple factors and gene interactions in which non-coding RNAs contribute to an individual's susceptibility to vitiligo (Mansuri et al., 2016;Vaish et al., 2019).
CircRNAs are non-coding RNAs characterized by a unique covalent closed-loop structure. It is widely known that circRNAs are rich in miRNA binding sites (miRNA response elements, MREs), serving as miRNA sponges or competitive endogenous RNAs (ceRNAs), which leading the inhibitory effect of miRNAs on their target genes and further increased expression of the target genes (Salzman, 2016;Kristensen et al., 2019). In this mechanism the circRNAs acting as miRNA sponges, play an important role in the regulation of diseases via their interaction with disease-associated miRNAs.
In our study we found there are 64 dysregulated circRNAs in the patients with vitiligo, including 34 circRNAs were upregulated and 30 were downregulated. And also there are 14 dysregulated miRNAs identified in vitiligo, 6 miRNAs were upregulated and 8 were downregulated. Through the correlation analysis, we obtained 12 dysregulated circRNAs and 5 dysregulated miRNAs, forming 48 relationships in the circRNA-miRNA-mRNA regulatory network. Our results showed that 12 dysregulated circRNAs had miRNAs binding sites, and were thus predicted to play a regulatory role via the ceRNA mechanism. And these circRNA-miRNA-mRNA regulatory networks may play an important role in the pathological mechanism of vitiligo. For instance, hsa-miR-149-5p targets genes SEMA4D (ENSG00000187764) and ANKRD6 (ENSG00000135299) are involved in regulating melanocyte survival and growth (Soong et al., 2012;Prabhakar et al., 2019). Sema4D is a paracrine factor that protects normal human melanocytes to survive under ultraviolet (UV) radiation. Sema4D could stimulate cell proliferation and regulates the activity of c-Met receptors, which was determined promotes melanocyte migration (Gotte et al., 2007;Soong et al., 2012). Soong et al. (2012) reported that Sema4D as the ligand for Plexin B1 could suppresses c-Met activation and migration and promotes melanocyte survival and growth. Hsa-miR-320a-3p regulated the target gene SMAD3 (ENSG00000166949). Smad3, a signal mediator of the activin/TGF-beta pathway, which associated with microphthalmia-associated transcription factor (MITF) (Funaba et al., 2003;Wehbe et al., 2012). Previous study indicated MITF is a crucial transcription factor expressed in melanocytes (Mizutani et al., 2010;Vachtenheim and Borovansky, 2010;Hu et al., 2020). Hsa-miR-27a-3p targets gene Paxillin (ENSG00000089159). Paxillin, a 69-kDa vinculin binding protein, is significant higher in neonatal melanocytes than in fetal melanocytes (Scott et al., 1995). It indicates that paxillin is closely related to the adhesion ability of melanocytes (Skoczynska et al., 2016).
GO analysis was performed to further annotate the biological functions of the host genes of the differentially expressed circRNAs. Interestingly, the top 20 GO terms of dysregulated circRNAs in vitiligo compared with normal are most enrich in the regulation of primary metabolic process, regulation of cellular metabolic process, negative regulation of cellular metabolic process, primary metabolic process, organic substance metabolic process, cellular metabolic process, organic cyclic compound metabolic process and cellular macromolecule metabolic process (Pietrzak et al., 2012;Atas and Gonul, 2017;Sahoo et al., 2017). This indicates that the dysregulation of circRNAs is closely related to the disorder of the metabolic pathway. In addition, the GO analysis of miRNA also showed that dysregulated miRNAs are closely related to biological or cellular metabolism (Supplementary Figure 1A).
The KEGG pathway of dysregulation of circRNAs mainly enriched in the biological processes such as ubiquitin mediated proteolysis, endocytosis and RNA degradation, and in Jak-STAT signaling pathway. Some studies have indicated that the Jak-STAT signaling pathway is involved in the immune-mediated inflammatory skin diseases such as vitiligo (Chen et al., 2019a;Gomez-Garcia et al., 2019;Montilla et al., 2019;Samaka et al., 2019). Furthermore, through KEGG analysis of miRNAs, some miRNA dysregulated lead to the disorder of certain signaling pathways (Supplementary Figure 1B), such as p53 signaling pathway, mTOR signaling pathway, Insulin signaling pathway and HIF-1 signaling pathway. Studies have shown that p53 signaling pathway and mTOR signaling pathway all played an important role in the synthesis or metabolism of melanin (Kadekaro et al., 2012;Chang et al., 2017;Yang et al., 2018;.
We elucidated the circRNA-associated-ceRNA profiles of vitiligo using deep RNA-seq analysis. From the analysis results, we found the circRNA-miRNA-mRNA regulatory network are involved in the regulation of numerous melanocyte functions, further GO and KEGG analysis show that the regulatory functions of these dysregulated circRNAs may closely related to the melanocyte metabolism. Through the analysis of the circRNA-miRNA-mRNA network, we found circ_0087961-miR-27a-3p-PAXILLIN might play an important regulatory role in the pathology of vitiligo. Our study provides a theoretical basis for studying the vitiligo pathogenesis from the perspective of circRNA-miRNA-mRNA network.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the NCBI repository (https://www.ncbi.nlm.nih.gov/sra/PRJNA712982), accession number PRJNA712982.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by People's Hospital of Guangxi Zhuang Autonomous Region. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
LL conceived and wrote the article. LL and ZX designed and performed the research. XQ, TW, MJ, and JQ analyzed the data. CW, RW, and CS contributed essential reagents or tools. All authors read and approved the manuscript.