Identification of Critical Biomarkers and Mechanisms of Fructus Ligustri Lucidi on Vitiligo Using Integrated Bioinformatics Analysis

Introduction

Vitiligo is a chronic autoimmune disease of the skin characterized by the appearance of patches of depigmentation that are visible as white spots, accounting for approximately 0.5–2% of the world’s population.¹ Recent studies have yielded a robust mechanistic understanding of vitiligo, including environmental stimuli, oxidative stress, and interactions between innate and adaptive immunity.² Although many mechanisms have been proposed, the exact pathogenesis of vitiligo remains unclear, posing a difficult challenge in its treatment. A significant advance in understanding the pathological basis has led to the discovery of potential biomarkers for vitiligo,³ which are used in clinical trials and daily practice.⁴ However, critical biomarkers of vitiligo have not been fully clarified.

The treatment methods used for vitiligo mainly include glucocorticosteroids, immunosuppressive agents, and calcineurin inhibitors.⁵ However, such long-term treatment methods often lead to serious side effects such as local skin irritation, atrophy, and inflammation. Since this complicated mechanism induces the onset of vitiligo, definitive safe and efficacious treatments for vitiligo are still lacking. Therefore, the development of a clinically effective and safe drug with fewer side effects for treating vitiligo is urgently needed. It is a feasible method to develop drugs based on critical biomarkers of vitiligo.

As a Chinese medicine, Fructus Ligustri Lucidi (FLL), the fruit of Ligustrum lucidum Ait. (Oleaceae), has been used for 2000 years in China.⁶ It maintains energy levels, nourishes the liver and kidneys, and strengthens bones.⁷ The antioxidant and anti-inflammatory activities of FLL are associated with the amelioration of vitiligo.^6,8 FLL promotes the adhesion and migration of human epidermal melanocytes, and decreases skin depigmentation, suggesting that FLL is a potential therapy for vitiligo.^9,10 Various extracts and compounds derived from FLL were detected using liquid chromatography mass spectroscopy.¹¹ However, the possible active compounds and mechanisms of FLL-induced melanogenesis remain unclear.

Microarrays based on high-throughput platforms have been widely used to explore and identify promising biomarkers for the diagnosis and prognosis of diseases at the genomic level because of their general and comprehensive analysis of genes.^12,13 Numerous studies have demonstrated that the pathophysiological process underlying the development of vitiligo is closely associated with mutations and abnormal gene expression.^14,15 Although integrating bioinformatics analyses have been attempted before, critical biomarkers and mechanisms of FLL on vitiligo remain unclear. In this study, two microarray datasets were used to screen the differentially expressed genes (DEGs) between the vitiligo and healthy samples. The mechanisms of vitiligo were explored using an enrichment analysis of functions and pathways. Based on the analysis above, network pharmacology techniques and molecular docking analysis were used to explore the underlying mechanism of FLL in the treatment of vitiligo.

Materials and Methods

Access to Gene Expression Omnibus (GEO) Datasets

GSE65127 and GSE75819 were selected using the GEO (http://www.ncbi.nlm.nih.gov/geo) database,^14,15 in which all identify genes and pathways involved in the formation of vitiligo compare with normal individuals. As the datasets used are from publicly available data of NCBI, the data accessed complied with relevant data protection and privacy regulations. The GSE65127 series on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), and the GSE75819 series on the GPL6244 platform (Illumina HumanWG-6 v3.0 expression beadchip). The basic information of GSE65127 and GSE75819 is shown in Table 1. Both the series matrix and platform files were obtained from the GEO. The density plots of the applied dataset were performed by GEO2R (Supplementary Figure 1). Log2 transformation was performed on each expression matrix using R language commands. After that, the gene probe IDs were transformed into gene symbols by R language commands.

Table 1 Data Set from GEO Database

DEGs Identified by R

Experimental groups for the two GEO series were established. Data processing was performed using R software (https://www.r-project.org/) and Bioconductor packages (http://www.bioconductor.org/). The differentially expressed genes were screened, and the differentially expressed genes were screened on the condition of p-value ≤ 0.05 and |Fold change|≥ 1. Genes without a corresponding gene symbol and genes with more than one probe set were separately removed. In order to identify common significant targets from DEGs, the Venn online tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html) was used to draw a Venn map, and overlapping DEGs were retained for further analysis. The volcano maps, heatmap, and pie charts were drawn using R.

Functional Annotation and Pathway Enrichment Analysis

To functionally annotate DEGs identified by the comparison mentioned above groups, annotation and visualization of gene ontology (GO) enrichment were mainly analyzed the biological process (BP), cellular composition (CC), and molecular function (MF) by R. BP, CC, and MF were performed by R. The DEGs were then introduced into the R for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Identified the DEGs for vitiligo were provided as input in the STRING database (http://string.embl.de/), was used to determine the protein-protein interactions (PPI) and included data on interacting proteins or genes within the human species.¹⁶ Interactions among the list of putative targets were retrieved using a threshold score of 0.7.

Validation of the Key Biomarkers in Vitiligo

The expression levels of intersection targets were extracted and divided into the normal group and vitiligo group. The statistical analysis was performed by GraphPad Prism 8.0.2. The correlation between key biomarkers for the applied datasets were performed by Prism 8.0 (Supplementary Figure 2).

Target Prediction of FLL Components

The BATMAN-TCM database (http://bionet.ncpsb.org/batman-tcm), an online bioinformatics analysis tool, was used to identify the main active compounds of FLL.¹⁷ The putative targets of FLL were uniformly transferred into a list of ingredients. In addition, the potential targets of these ingredients were predicted by the similarity-based target prediction method. The candidate targets were predicted (score cutoff >20), and used for further bioinformatics analyses to obtain an accurate prediction performance. To mine the essential compounds from FLL, a compounds-targets network of FLL was constructed. The compounds-targets relation was calculated via a network topological analysis using Cytoscape 3.7.1 using the NetworkAnalyzer plugin.¹⁸

Molecular Docking

Molecular docking was performed by using the Auto Dock Vina software. The structure of disease-targets used for docking was obtained from Protein Data Bank (http://www.rcsb.org). The compound’s structure was downloaded from the PubChem database (https://www.ncbi.nlm.nih.gov/pccompound/?term=). PyMol software was used to delete the water molecules and small molecule ligands from the protein structure. Then the compounds and proteins were imported into the Auto Dock Tools to dock separately by format files. Finally, the lowest binding energy data was obtained from the molecular docking. The binding force between the active ingredient and the target protein was calculated. A drug molecule with binding energy ≤ −5.0 kJ/mol was considered to have a good binding activity to the targets.

Cell Viability Assay

B16F10 cells were cultured and preserved in DMEM/F12 medium, and supplemented with 10% FBS and 1% penicillin-streptomycin in an incubator (37°C, 5% CO₂). MTT (BioFroxx, Guangzhou, China) assay was used to assess cell viability after treat with FLL.¹⁹ The B16F10 cells were seeded at a density of 5×10³ cells per well in 96-well plates, treated with FLL (0.315ng/mL~315000 ng/mL). After 24h, 10μL MTT solution was added, incubated at 37°C for 4h, and absorbance was determined at 490nm.

H₂O₂ Injury Model and Treatment

The B16F10 cells were seeded at a density of 5×10³ cells per well in 96-well plates and 1×10⁴ cells per well in 6 -well plates, respectively. Then the cells were treated with the H₂O₂ (800 μM) for 24 h as model group and H₂O₂-free as blank control group. The cells were treated with various concentrations of FLL as experimental group, and treated with the edaravone (100 μM) as a positive group. The survival rate of the cells was measured by using a MTT assay as described in 2.7.

Western Blot

Western blot was performed according to the methods of previous published paper.²⁰ The following primary antibodies were used: TYRP1 (ab235447, 1:1000, Abcam) and GAPDH (60004-1-Ig, 1:8000, PROTEINTECH GROUP). After four washes with Tris buffered saline Tween, the membranes were incubated with Horseradish Peroxidase (HRP)-conjugated goat anti-rabbit IgG secondary antibodies (SA00001-2, 1:5000, PROTEINTECH GROUP) or goat anti-mouse IgG secondary antibodies (SA00001-1, 1:5000, PROTEINTECH GROUP) and were visualized by MiniChem610 (SAGECREATION, Beijing, China). The intensity of protein bands was quantified with Image J 1.53t software.

Results

Identification and Analysis of DEGs in Datasets

Based on high-throughput analysis, DEGs in the two microarray datasets (GSE65127 and GSE75819) were screened after the chip results were normalized (Table 1). All microarray datasets were analyzed using R language, with P<0.05, and |logFC|≥1 as the screening conditions. GSE65127 was obtained with 45 DEGs, including 2 up-regulated genes and 43 down-regulated genes; GSE75819 was obtained with 911 DEGs, including 857 upregulated and 54 downregulated genes (Figure 1).

Figure 1 Identification and analysis of DEGdds in GSE65127 and GSE75819 datasets. (A) The volcano plot illustrates DEGs in GSE65127; (B) the heatmap illustrates DEGs in GSE65127; (C) the pie chart illustrates DEGs in GSE65127; (D) the volcano plot illustrates DEGs in GSE75819; (E) the heatmap illustrates DEGs in GSE75819; (F) the pie chart illustrates DEGs in GSE75819. In the volcano plot and heatmap, red represents high expression, green represents low expression, and black represents no difference. The heat map is divided into the normal group and the vitiligo group. The red indicates that the gene is highly expressed in this group, and the green indicates the low expression. The pie chart makes a statistic for the differential genes.

Screening of Common Targets, Functional Enrichment, and Target Interaction Analysis

As shown in the Venn map, 13 common targets overlapped among the differential targets of GSE65127 and GSE75819 (Figure 2A). Thirteen common targets include Transient receptor potential cation channel subfamily M member 1 (TRPM1), Melanoma antigen recognized by T-cells 1 (MLANA), Membrane-associated transporter protein (SLC45A2), Tyrosinase (TYR), Neuronal membrane glycoprotein M6-b (GPM6B), Dipeptidyl aminopeptidase-like protein 6 (DPP6), L-dopachrome tautomerase (DCT), Leucine-rich repeat LGI family member 3 (LGI3), Transcription factor SOX-10 (SOX10), Carbonic anhydrase 14 (CA14), Cellular retinoic acid-binding protein 1 (CRABP1), Tubulin beta-2B chain (TUBB2B), and 5.6-dihydroxyindole-2-carboxylic acid oxidase (TYRP1). GO enrichment analysis was used to evaluate the potential mechanism of common targets from the biological process and cellular component categories based on these common targets. BP is involved in development of pigmentation, melanin biosynthetic processes, and melanin metabolic processes; CC is involved in melanosome membrane, pigment granule membrane, and melanosome; and MF is involved in copper ion binding, and monooxygenase activity (Figure 2B and C). In the KEGG enrichment analysis, the pathway is involved in tyrosine metabolism, melanogenesis, and nitrogen metabolism (Figure 2D). The interaction relationships between the targets were obtained using the STRING database (Figure 2E).

Figure 2 Screening of common targets, functional enrichment, and target interaction analysis. (A) The common targets in the Venn map; (B) BP part for GO enrichment analysis; (C) CC part for GO enrichment analysis; (D) KEGG enrichment analysis; (E) targets interaction analysis.

Validation of the Key Biomarkers in Vitiligo

To better understand the common targets as critical biomarkers in vitiligo, differences in the protein expression of 13 targets were analyzed (Figures 3 and 4). Through GSE65127, the results showed that the expression of 13 proteins such as MLANA, DPP6, DCT, SOX10, and TYRP1, was significantly different between the normal and vitiligo groups (Figure 3). Through GSE75819, the results showed that the expression of eight proteins including MLANA, TYR, DCT, LGI3, SOX10, CRABP1, TUBB2B, and TYRP1, was significantly different (Figure 4). Eight intersecting proteins were obtained, which are the key targets of MLANA, TYR, LGI3, CRABP1, TUBB2B, DCT, SOX10, and TYRP1.

Figure 3 The differences in protein expression of 13 targets on GSE65127. “N” represent non-lesional skin, and “L” represent lesional skin. Compared with the non-lesional skin group at ***P < 0.001, and ****P < 0.0001.

Figure 4 The differences in protein expression of 13 targets on GSE75819. “N” represent non-lesional skin, and “L” represent lesional skin. Compared with the non-lesional skin group at **P < 0.01, ***P < 0.001, and ****P < 0.0001.

Target Prediction of FLL Components

To determine the mechanisms of the effects of FLL in vitiligo, for the query of composite compounds, only the predicted candidate target proteins with scores ≥ 20 were presented by BATMAN-TCM. Eight compounds were identified in FLL. Cytoscape was used to analyze the components of the protein interaction network (Figure 5). Furthermore, to determine the critical compounds in FLL, Cytoscape was used to analyze the network, and it was found that the top vital compounds included ursolic acid and oleanolic acid (Figure 5).

Figure 5 Target prediction of FLL using bioinformatic tools. (A) Construction of compounds and targets of vitiligo; (B) the degree of compounds of FLL; (C) the betweenness centrality of compounds of FLL; (D) the closeness centrality of compounds of FLL.

Verification by Molecular Docking

To further explore the regulatory effect of FLL on vitiligo, molecular docking verification was performed on the above-mentioned main chemical components and predicted key targets through AutoDock Vina molecular docking. The virtual map of energy and molecular docking are shown in Table 2 and Figure 6, respectively. These results show that the selected active ingredients can bind to critical targets, further proving the reliability of network pharmacology for predicting active ingredients and their targets.

Table 2 Binding Energy of Ursolic Acid and Oleanolic Acid

Figure 6 Ursolic acid and oleanolic acid are bound to residues outside the active pockets of DCT, MLANA, and TYPR1. (A) Ursolic Acid is bound to residues outside the active pockets of DCT, MLANA, and TYPR1; (B) oleanolic acid is bound to residues outside the active pockets of DCT, MLANA, and TYPR1.

Verification by MTT and Western Blot

The viability of B16F10 cells was few changes after treated by FLL range from 0.315 to 31,500  ng/mL (Figure 7A). The result showed that FLL significantly increased the cell viability in the H₂O₂-induced B16F10 cell death model in the concentration range of 0.315–315 ng/mL (Figure 7B). The verification data found that, the expression of TYRP1 was significantly increased in the H₂O₂-treated B16F10 cells with a FLL treatment at 315 ng/mL (Figure 7C).

Figure 7 FLL increased the cell viability, and expression of TYRP1 in H2O2-induced B16F10 cell death. (A) B16F10 cells were treated with different concentrations of FLL; (B) the cell viability of B16F10 cells induced by H2O2 was determined by MTT assay; (C) the protein expression of TYRP1 was measured by Western blot. Data are presented as mean ± SEM. Compared with the control group at **P < 0.01, and compared with the H2O2-induced model group at #P < 0.05, and ##P < 0.01.

Discussion

Vitiligo is an autoimmune skin disease that affects people worldwide.²¹ Modern clinical treatment often combines cell transplantation,²² laser therapy,²³ and other methods; however, the mechanism of vitiligo remains unclear. In this study, two microarray datasets from GEO to screen for DEGs between the vitiligo and healthy samples were analyzed. Based on bioinformatics and network pharmacology, we found that the critical targets of MLANA, DCT, and TYRP1 may be related to the mechanism of action of FLL in the treatment of vitiligo.

Compare with traditional methods such as HitPredict, Traditional Chinese Medicine Integrated Database, and Traditional Chinese Medicine Systems Pharmacology, BATMAN-TCM can provide target prediction online for each ingredient. In addition, the advantage of this method is that it integrates the data of targets from the samples of human species, compounds, and molecular docking.¹⁷ First, we identified and analyzed vitiligo-related DEGs using the GSE65127 and GSE75819 datasets, which contain information of samples of human species. As a result, 13 common targets overlapped among the differential targets of GSE65127 and GSE75819. Moreover, GO analysis revealed that the vitiligo-related DEGs were mainly enriched in development of pigmentation, melanin biosynthetic processes, and melanin metabolic processes. Previous studies have shown that vitiligo may be closely related to melanin metabolism,²⁴ supporting our findings. KEGG enrichment analysis revealed that the pathways associated with DEGs were significantly involved in tyrosine metabolism, melanogenesis, and nitrogen metabolism. According to relevant studies, tyrosine metabolism is enriched in patients with vitiligo.²⁴ Therefore, these DEGs may play essential roles as potential drug targets for vitiligo via the pathways described above.

Our data indicated that MLANA, DCT, and TYRP1 might be related to the mechanism of action of FLL in the treatment of vitiligo. MLANA is involved in melanosome biogenesis and the expression, stability, trafficking, and processing of the melanocyte protein PMEL, which is critical for the formation of stage II melanosomes.^25,26 DCT can catalyze the conversion of L-dopachrome into 5.6-dihydroxyindole-2-carboxylic acid, which is related to melanin and pigment biosynthesis.²⁷ TYRP1 plays a vital role in melanin biosynthesis.²⁸ Some studies have promoted the potential mechanism of melanin production by activating the CREB-MITF signaling pathway by increasing the expression levels of genes such as TYRP1 and DCT.^29,30 In addition, memory T cells induced after tumor resection recognize the melanocyte differentiation antigens TRP-2/DCT and gp100 to promote the formation of autoimmune vitiligo, suggesting that DCT is a critical target of vitiligo.^31,32 Studies have identified 12 genes (including MLANA and TYRP1) that affect melanocyte melanogenesis and cellular oxidative stress by integrating transcriptome and differential methylation analyses.³³ The significance of vitiligo treatment will help to better understand the molecular mechanisms of vitiligo, and explore new treatment strategies.³³ Thus, DCT, TYRP1, and MLANA may be key potential targets in the treatment of vitiligo. According to the results of molecular docking, the expression of TYRP1 was verified in vitro experiments by using Western blot method. The data found that FLL significantly increased the expression of TYRP1. However, further in vitro or in vivo studies are required to confirm our findings.

Another vital finding of our study is that ursolic acid and oleanolic acid from FLL may be active compounds corresponding to MLANA, DCT, and TYRP1 in vitiligo. Ursolic acid and oleanolic acid are pentacyclic triterpenoids that have long been used in cosmetics and healthcare products because of their anti-inflammatory and antioxidation.^10,34,35 Both the FLL extract and ursolic acid can significantly improve the adhesion of human epidermal melanocytes. In addition, the FLL extract also increased the activation rate of tyrosinase, while ursolic acid and oleanolic acid are the active components of FLL that promote the proliferation of melanocytes and the synthesis of melanin.^10,36 Moreover, its therapeutic potential in vitiligo warrants further investigation. We could use the experiments in vivo in the follow-up study to explore the role of the active components of FLL.

In conclusion, TYRP1, as a melanocyte molecular biomarker, may be closely related to the underlying mechanism of FLL in the treatment of vitiligo via the inhibition of melanocyte death.