Pathogenesis of Alopecia Areata Based on Bioinformatics Analysis.

BACKGROUND: Alopecia areata (AA) describes a sudden localized patchy alopecia. The cause of AA is not completely clear and its incidence may be related to genetic, autoimmune, and environmental factors. AIM: To explore the possible mechanisms of AA and to provide a basis for the early diagnosis and treatment of AA.
METHODS: Gene microarray data from 122 scalp skin biopsy tissue samples from patients with AA or healthy controls from the Gene-Cloud of Biotechnology Information database were analyzed using bioinformatics analysis methods. Molecular network analysis of the differentially expressed genes (DEGs) was conducted by Cytocluster using the Molecular Complex Detection (MCODE) algorithm.
RESULTS: The gene expression profile of skin lesions from patients with AA was significantly altered, with 111 DEGs found in the skin lesions of AA, compared with that of the healthy skin. The DEGs were mainly related to biological processes such as the development of the epidermis and inflammatory reaction. The protein-protein interaction network analysis of DEGs revealed bone morphogenetic protein 2 (BMP2) as a core protein interaction network. BMP2 acted not only via the inflammatory response but also via the signaling pathways in epithelial cell development and epidermal cell differentiation to affect the epidermal development. MCODE analysis further showed that keratins (KRTs) and keratin-associated proteins (KRTAPs) can affect the epidermal development via the epidermal development pathway.
CONCLUSIONS: The abnormal development of the epidermis and inflammatory reactions in skin tissue play important roles in the pathogenesis of AA and are closely related to BMP2, KRTs, and KRTAPs genes. LIMITATIONS: Our study was limited by experimental verification.

Introduction

Alopecia areata (AA) is a sudden localized patchy alopecia, which can occur on any part of the body and at any age. It is common in young and middle-aged adults, with a prevalence of 0.9%–6.9%.[12] Although AA does not endanger life, it can bring about significant psychological burden to patients and can seriously affect their quality of life.[34] At present, the cause of AA is not completely clear and its incidence may be related to genetic, autoimmune, and environmental factors.[5678] Twin studies show that genetic predisposition is an important factor in the development of AA and that there is a significant association with family history.[910] In the present study, we analyzed gene microarray data from the Gene-Cloud of Biotechnology Information (GCBI) database from patients with AA and from healthy controls to explore the molecular level changes in the skin of patients with AA and the possible underlying molecular causes of AA in order to identify the putative candidates for the early diagnosis and intervention of AA.

Methods

Data source

Data from gene chip GSE68801 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68801) were obtained from the GCBI high-throughput data analysis platform using the search term “alopecia areata”, and the data which were submitted by Jabbari et al., were analyzed on May 12, 2015 using Affymetrix Human Genome U133 Plus 2.0 Array Gene Expression Chip Platform GPL570, and updated on September 15, 2017. A total of 122 samples were collected for this analysis, including 60 scalp skin lesional samples from patients with AA, 26 scalp skin nonlesional samples from patients with AA, and 36 healthy control samples. The GSE68801 dataset contained a total of 54,675 gene expression data points.

Screening of differentially expressed genes

The GSE68801 dataset was imported into the GCBI analysis laboratory and normalized according to Log2. In this study, we identified the differentially expressed genes (DEGs) according to the standard: P< 0.05, false discovery rate < 0.05, and fold change >2, and independent-samples t-test was performed. The list of DEGs between scalp skin lesional samples and nonlesional samples from patients with AA were marked as A, and the list of DEGs between scalp skin nonlesional samples from patients with AA and healthy control samples were marked as B. We removed the coincidence gene in A and B from A, and then the DEGs were marked as C. The DEGs in group C could be considered as genes having specific etiology to AA.

Identification pathway analysis

All the genes in Group C were uploaded to the GCBI analysis laboratory for gene ontology analysis (GO analysis) and Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway enrichment analysis.

Protein–protein interaction networks

The protein names corresponding to 111 DEGs between scalp skin lesional samples and nonlesional samples were uploaded to the protein–protein interaction on-line analysis software STRING 10.0 (https://string-db.org/). The reliability and additional node parameters were properly adjusted according to the specific analysis results, and the DEGs’ protein–protein interaction map was drawn.

Kyoto Encyclopedia of Genes and Genomes signaling pathway network

To further analyze the potential interaction between the key gene and the KEGG signaling pathway in the process of AA, Cytoscape's CLUEGO platform was used to establish a gene–pathway interaction network.

The Molecular Complex Detection analysis

The specific biological network of AA was analyzed based on the potential physiological effects of the 111 DEGs on the background of GeneMania human genetic biological network; we identified the important gene cluster in this network using the Molecular Complex Detection (MCODE) method and predicted the potential genes involved in the development of AA.

Results

Gene chip data preprocessing

After we normalized the gene chip data from GSE68801, it was evident that the entire dataset was centered around zero with small fluctuations [Figure 1], which indicated that the data quality was eligible and highly stable and that the data could be used for the follow-up analysis.

Figure 1

Data distribution of scalp skin samples gene transcriptional chip of alopecia areata and healthy controls. The gene chip data from GSE68801 was centered around zero with small fluctuations, which indicated that the data quality was eligible and highly stable

Analysis of differentially expressed genes among the three groups of scalp skin samples

There were no significant DEGs between the nonlesional samples from patients with AA and the healthy control samples, so the change of transcriptional profile of the two groups was not obvious. According to the same screening criteria, there were 111 DEGs between scalp skin lesional samples and nonlesional samples from patients with AA. As shown in Figure 2, 107 DEGs expressed downward accounted for 96.40%, while 4 DEGs expressed upward accounted for 3.60%; this overall change in the DEG expression profile indicated that the gene expression profile of scalp skin lesional samples from patients with AA was significantly changed when compared to nonlesional samples.

Figure 2

Gene expression profiles of scalp skin lesional samples and nonlesional samples from patients with alopecia areata. The red dots represent differentially expressed genes, and the blue dots indicated no statistically significant genes. And, there were 111 differentially expressed genes between scalp skin lesional samples and nonlesional samples from patients with alopecia areata

Differentially expressed genes, signaling pathway, and gene ontology analysis

The 111 DEGs were involved in six signaling pathways, which GO analysis aggregated to three signaling pathways, among which the epidermal development pathway was the most influential [Figure 3 and Supplementary Table 1]. Among the DEGs, there were seven genes (KRTAP5-9, KRT34, KRT85, KRT32, KRT16, KRT83, and KRT31) related to epidermal development [Table 1], seven genes (FGF18, BMP2, FGF5, CCL18, GJA3, CXCL10, and CXCL9) involved in cell–cell signaling [Table 1], and six genes (CHST2, CXCL10, CXCL9, BMP2, SERPINA3, and CCL18) related to inflammatory response. These indicated that epidermal development, inflammatory response, and other biological processes played important role in the pathogenesis of AA.

Figure 3

Differentially expressed genes analyzed by Kyoto Encyclopedia of Genes and Genomes pathway and gene ontology analysis. The circle meant the counts of differentially expressed genes in pathway. The color from green to red indicated that the difference was from small to large, while the deeper the color, the greater was the difference

Table 1

sDifferentially expressed genes related to Gene ontology

Differentially expressed genes related to KEGG Pathways

Click here for additional data file.

Protein–protein interaction network of the differentially expressed genes

The corresponding proteins of DEGs were imported into the NetworkAnalyst (http://www.networkanalyst.ca/) to establish the protein–protein interaction network. We found that bone morphogenetic protein 2 (BMP2) interacts with other proteins in the network [Figure 4], and it was the core of the entire network. Therefore, it was likely that BMP2 played a key role in AA.

Figure 4

Protein–protein interaction network of differentially expressed genes. Bone morphogenetic protein 2 interacts with other proteins in the network; it is the core of the entire network

Differentially expressed genes and signaling pathway interaction network

The coexpression of interaction of DEGs and signaling pathways revealed that BMP2 directly interacts with signaling pathways in epithelium development and epidermal cell differentiation [Figure 5] and then affected the signaling pathway in epidermal development, leading to epidermal dysplasia that caused AA.

Figure 5

Coexpression network of differentially expressed genes and signaling pathways. BMP2 directly interacts with the signaling pathways in epithelium development and epidermal cell differentiation and then affects the signaling pathways in epidermal development

Molecular Complex Detection analysis and candidate gene prediction

The physiological function-specific network of AA was established for the 111 DEGs using MCODE analysis, and the top three clusters are shown in Figure 6. There were a total of 38 genes in cluster 1 including 12 (31.58%) predicted genes [Figure 6a], 21 genes were in cluster 2 including 2 (9.52%) predicted genes [Figure 6b], 4 genes were in cluster 3 with no predicted genes found [Figure 6c], and there were 14 potential genes associated with AA [Table 2]. Interestingly, keratins (KRTs) and keratin-associated proteins (KRTAPs) might play important roles in the development of AA.

Figure 6

Molecular Complex Detection analysis of alopecia areata-related genes. The genes marked red were differentially expressed genes listed in C while the green were predicted genes based on its biological function. a: The cluster 1 (cluster score = 38) showed role of keratin associated protein; b: The cluster 2 (cluster score = 21) indicated role of keratins in the development of alopecia areata; c: The cluster 3 (cluster score = 4) with no predicted genes

Table 2

Candidate genes which had potential biological associations with alopecia areata

Discussion

AA typically manifests as the sudden appearance of round or oval patches of hair loss, with a diameter of 1–10 cm and varying in number. The lesional area is smooth, with no outward signs of inflammation, scaling or scarring. Patients may continue to develop more severe hair loss, including complete loss of eyebrows, eyelashes, armpit hair, pubic hair, and body hair. Depending on the degree of hair loss, AA is grouped into patchy AA, alopecia totalis (AT), and alopecia universalis (AU); the course of AT and AU marked with bad prognosis. Younger the age at onset, worse is the prognosis.

In this study, we carefully explored and analyzed the gene chip data of scalp skin samples from patients with AA using genomic analysis. Gene expression profiling showed that there were 111 significant DEGs between scalp skin lesional samples and nonlesional samples from patients with AA, while there were no significant DEGs between scalp skin nonlesional samples from patients with AA and healthy control samples, illustrating that gene expression profile of scalp skin lesional samples from patients with AA was significantly changed.

Gene ontology (GO) analysis showed that biological processes, such as epidermal development and inflammatory response, played important roles in the morbidity of AA. Among the DEGs, KRTAP5-9, KRT34, KRT85, KRT32, KRT16, KRT83, and KRT31 were related to epidermal development, while CHST2, CXCL10, CXCL9, BMP2, SERPINA3, and CCL18 were related to inflammatory response.

Research had shown that lesional skin had a high level of immune activity, and hair loss in the lesional area was related to infiltration of CD4+ T cells around the hair follicle and CD8+ T cell infiltration in the hair follicle.[11] Thein et al.[12] demonstrated that T cells in the lesional area could release interferon γ and tumor necrosis factor α (TNF-α) to inhibit the growth of keratinocyte and the proliferation of hair follicle epithelial cells which then interfered with hair growth using T cell cloning technology. The hair follicles in anagen were attacked by T cell-mediated immune responses in patients with genetic susceptibility to AA who were missing immune privilege in hair follicles.[1314] Another study confirmed[15] that Th1 cytokine was downregulated in AA while interleukin-10, which was secreted by Th2 cells, was upregulated after treatment with diphenylcyclopropenone.

The protein–protein interaction networks of all the DEGs show an obvious protein–protein interaction network with BMP2 at the core. BMP2, from the BMP family, regulates the cycle of hair follicles, inhibits the regeneration of hair follicles, and also plays an important role in the transformation of hair follicles from the growth phase to the degenerative phase and maintenance of hair follicles.[1617] Our study showed that the level of BMP2 in the lesional skin of patients with AA was decreased, indicating that as the inhibition of hair follicle regeneration in the lesional skin was weak, hair follicles could enter the growth phase from the resting phase and promote the regeneration of hair follicles, in accordance with the Rui et al.

's report,[18] which showed that BMP2 was highly expressed in hair follicle in resting phase and lowly expressed in cashmere goats. BMPs can inhibit the proliferation, recombination, and migration of horn cells in epithelial cells of regenerated skin and play certain effects on hair follicle formation at the embryonic stage, the cycle of hair follicles after birth, the differentiation, development and growth of hair follicles.[19] TNF-α induced epithelial–mesenchymal transition during wound healing process is achieved through BMP2.[20] GO analysis showed that BMP2 was related to the inflammatory response in the lesional area of the patients with AA. The coexpression of DEGs and signaling pathways revealed that BMP2 could affect the epidermal development by interfering with epidermal cell differentiation. The decrease of BMP2 expression in lesional skin of patients with AA interfered with the normal development and differentiation of epithelial cells and led to epidermal dysplasia that caused AA.

KRTs and KRTAPs play important roles in the development of AA, as determined from the physiological function-specific network generated using MCODE analysis of AA based on the 111 DEGs. There were 19 KRTs and 26 KRTAPs in the DEGs, and another 14 potential genes related to AA were found (two KRTs and 12 KRTAPs). KRTs are distributed in human and animal skin tissues and are important components of hoof, hair, horns, etc., and are expressed in hair follicles to maintain the structure of hair follicles, form hair follicle cells, and partake in cell signal transduction and apoptosis at the same time. At present, more than 80 KRTAPs genes and 25 gene families have been identified in humans, and all of the found KRTAPs genes are expressed in human hair follicle with the exclusion of KRTAP16, KRTAP22, KRTAP25, and KRTAP27.[212223] The KRT and KRTAP families are a class of genes that are specifically expressed in hair follicle tissue. The present study showed that KRTAP5-9, KRT34, KRT85, KRT32, KRT16, KRT83, and KRT31 were related to epidermal development, while BMP2 could regulate the hair follicle cycle and also affected epidermal development by interfering with epithelial cell development and epithelial differentiation pathway. This indicated that BMP2 might interact with KRTAPs and KRTs and affected the development of AA via the epidermal development pathway.

Conclusion

In conclusion, the gene expression profile of lesional area of scalp skin of patients with AA changes significantly during the development of AA. Based on the biological function of DEGs (scalp skin lesional samples from patients with AA vs. scalp skin nonlesional samples from patients with AA), we found that the epidermal developmental injury and inflammatory response in the skin tissue played a key role in the development of AA; KRTs and KRTAPs might play important role in the epidermal developmental injury and might have direct or indirect interaction with BMP2. The above biomarkers might be of use in the diagnosis, treatment, and prognosis of AA, but the mechanism and clinical value still needed future studies. Therefore, we believe that an increase in the expressions of the genes, such as BMP2, KRT16, and KRTAP5, in the lesional skin tissue may suggest a better prognosis; while a low value may indicate a relapse or worsening of the prognosis. In our study, we found that there were no significant changes in gene expression between the nonlesional samples from AA patients and the healthy control samples, while the gene expression profiles of lesional samples were significantly altered in comparison with nonlesional samples from AA patients. The abnormal development of the epidermis and inflammatory reactions in skin tissue played important roles in the pathogenesis of AA and are closely related to BMP2, KRTs, and KRTAPs genes.

Financial support and sponsorship

This study was funded by the Scientific and Technological Achievements Project of Xianyang City (2017KT-16).

Conflicts of interest

There are no conflicts of interest.