Bioinformatics analysis of genes associated with the patchy-type alopecia areata: CD2 may be a new therapeutic target

Background. Alopecia areata (AA) is mainly a T cell-medicated autoimmune disease with non-scarring hair loss and limited treatment options. Of these, the patchy-type alopecia areata (AAP) is the most common and relatively easy to treat due to smaller areas of the scalp affected. To understand the pathogenesis of AAP and explore the therapeutic target, we focus on the molecular signatures by comparing AAP and normal subjects. Methods. The gene expression profile (GSE68801) was obtained from Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) were identified between AAP patients and normal controls using the GEO2R. Then the Gene ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and Protein-Protein interaction (PPI) network analysis were performed for DEGs. Results. A total of 185 DEGs were identified, including 45 up-regulated genes and 140 down-regulated genes. The upregulated DEGs were related to the immune response and chemokine signaling pathway. Meanwhile, down-regulated DEGs were enriched in keratin filament and intermediate filament. Subsequently, the top 10 hub genes were picked out in the PPI network, among them, CD2 showed the highest connectivity degree and central roles. Conclusion. Our data suggest that the CD2 may be a new therapeutic target for AAP. Further study is needed to explore the value of CD2 in the treatment of alopecia areata.


INTRODUCTION
Alopecia areata (AA) is a multifactorial disease with nonscarring hair loss resulting from disrupted immune privilege, autoimmune-mediated destruction, and up-regulated cytokine pathways of hair follicle, which was a great mental burden on patients 1 . The prevalence of alopecia varies by race and region, ranging from 0.1% to 0.2% of the general population 2 . Due to the difference of clinical manifestations, AA was classified into three major phenotypic variants patchy-type AA (AAP), alopecia totalis (AT) and alopecia universalis (AU) (ref. 3 ).
Histological examination shows that hair follicles of AA are infiltrated by CD4 + and CD8 + T cells, natural killer (NK) cells and macrophages, etc 4 , but no essential differences have been established between AAP and AT/ AU samples. Besides this, the genome wide association study (GWAS) and subsequent Genome-wide meta-analysis found that ULBP are risk genes for activating ligands of the NK cell receptor NKG2D, which interacts with mainly CD8 + NKG2D + T cells, accelerating the pathogenesis of AA (ref. [5][6][7]. Other recent studies also showed that cytotoxic CD8 + NKG2D + T cells were critical for AA in mouse models of disease and inhibition of Janus kinases (JAKs) has potential clinical utility, based on the AA signature connected with interferon-γ, IL-2, TH1, TH2, IL-23, and IL-9/TH9 cytokine activation [8][9][10] .
Although much research on AA has been done, the precise pathogenesis and etiology of AA remains unknown, especially the specific subtype. Of three major phenotypic variants, AAP is limited to well-demarcated patches and relatively easy to treat early before progressing to the entire scalp or body hairs. Thus, it is necessary to gain insight into the pathogenesis of AAP, using gene expression analysis to explore therapeutic targets. In the present study, we analyzed the gene expression profile of AAP patients versus normal controls to screen for differentially expressed genes (DEGs). Next, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for DEGs. In the end, we created a PPI network to explore the hub genes related to AAP.

Data resource
The gene expression of GSE68801 was obtained from Gene Expression Omnibus database (GEO, http://www. ncbi.nlm.nih.gov/geo/) of NCBI, which was performed on Affymetrix Human Genome U133 Plus 2.0 Array platform. Our dataset included 22 AAP patients samples and 36 normal controls.

Screening of DEGs
The GEO2R tool was performed to identify differently expressed genes (DEGs) between AAP patients and normal controls. GEO2R employed GEO query of BioConductor and package limma of R to perform the statistical analysis 11 . Adjusted P-value (adj. P Val)<0.05 and a fold-change (FC)<1.0 were set as the thresholds for DEGs.

GO term and KEGG pathways analysis
GO annotation and KEGG pathway enrichment analysis of DEGs were performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) 12 , which is a bioinformatics library that integrates biological data and analysis tools to provide systematic functional annotations for large-scale gene or protein lists. P<0.05, FDR (false discovery rate) <0.05 and the count ≥ 5 were set as thresholds for the pathway enrichment analysis. In the volcano plot, X-axis is fold change (log2) and Y-axis is adjust P value (-log10). Red points (fold change>1) indicate upregulated genes, whereas blue points (fold change<-1) indicate downregulated genes. Moreover, the darker red indicates a stronger upregulation in expression and the darker green indicates a stronger downregulation.

Construction of the protein-protein interaction (PPI) network
PPI network of DEGs was constructed using STRING database (http://string-db.org/), which maps a set of genes to be studied into network of relevant protein interac-tions to analyze hub proteins and screen out key genes associated with the specific disease. Interactions with Combined Score>0.4 were retained and then visualized using Cytoscape software (www.cytoscape.org/).

DEGs Identification
According to the predetermined threshold (adj. P Val<0.05 and |log 2 FC|>1.0), a total of 185 DEGs, including 45 up-regulated genes and 140 down-regulated genes, were screened out. The top 10 up-and down-regulated DEGs are listed in Table 1. Subsequently, a heatmap of the top 35 up-regulated and down-regulated DEGs and volcano plot of all DEGs was created to visualize expression variation between AAP patients and normal controls (Fig. 1a, b).

Function enrichment analysis of DEGs related to AAP
GO terms and KEGG pathways enrichment analysis was performed using DAVID to demonstrate the function of screened DEGs (Table 2). GO enrichment analysis indicated that the up-regulated genes were mainly enriched in biological processes associated with the immune response, signal transduction and inflammatory response ( Table 2 and Fig. 2a). Down-regulated genes were mainly enriched in the hair cycle and aging (Table 2 and Fig. 2b). In terms of the cellular component, up-regulated genes were mainly enriched in the plasma membrane, including the integral component of plasma membrane and external side of plasma membrane, while down-regulated DEGs were primarily enriched in keratin filament and intermediate filament. In terms of molecular function, up-regulated genes were significantly enriched in chemokine activity, whereas down-regulated genes were mainly enriched in structural molecule activity. Subsequently, the results of KEGG pathway analysis indicated that the up-regulated DEGs were primarily enriched in Chemokine signaling pathway. However, there is no pathway for down-regulated genes to be enriched in.

PPI network construction and hub genes identification
By mapping separately the up-and down-regulated DEGs into STRING database, the PPI networks of DEGs were visualized using Cytoscape (Combined Score ≥ 0.4). A total of 93 nodes and 184 edges were contained in the PPI network (Fig. 3). Degrees ≥10 is used as a cutoff criterion for judging of hub protein. The top 10 hub genes were screened out, including CD2, CCL5, CXCL9, CD3D,   (Table  3). Among them, CD2 showed the highest connectivity degree, and all the genes except KRTAP17-1 and BMP2 were up-regulated.

DISCUSSION
Alopecia areata (AA) is a multifatorial, non-scarring hair loss disease in which common clinical manifestations is a clear round or oval patchy hair loss on the scalp. Among three major phenotypic variants, AAP is the most common and relatively easy to treat due to smaller areas of the scalp area or beard area. Although the key role of autoimmunity in the pathogenesis of AA has been proved, the core genes and therapeutic targets of AAP still need further identification.
Hence, in our study, GSE68801 was utilized to screen the DEGs between 22 AAP patients and 36 normal controls from the GEO database. We identified a total of 185 DEGs, in which 45 genes were up-regulated and 140 genes were down-regulated. Based on the results of GO and KEGG pathway enrichment analysis, the up-regulated DEGs were mainly enriched in immune response, which are consistent with previous reports 9 . The down-regulated DEGs were mainly involved in cell component associated with keratin filament and intermediate filament. Further, it was reported that only these four functional keratins genes (KRT81, KRT83, KRT85 and KRT86) are related to human genetic disorders of the hair and/or nails [13][14][15] . Besides this, KRT85 plays a more critically important role in the structure of hair and nail than other keratins, therefore, mutations in KRT85 may cause alopecia universalis and exceptionally abnormal nails, due to recessive ectodermal dysplasia of hair and nails 14 . According to our research, although the keratins DEGs of KRT81, KRT83, and KRT85 were all significantly down-regulated, these patients still suffer from AAP, implying the close connection between the genes encoding four keratins and pathogenesis of AA. Further, the PPI analysis showed Fig. 3. Protein-protein interaction network constituted by up-and down-regulated DEGs and visualized using Cytoscape (Combined Score>0.4). Red nodes stand for up-regulated genes, while green nodes stand for down-regulated genes. The lines represent interaction relationship between nodes. Besides, the size of nodes represents the connectivity degree and the larger the size, the larger the degree. that desmoglein 4 (Dsg4) and plakophilins 1 (PKP1) were connected with the KRT83, which was in accordance with relevant researches on the importance of PKP1 and Dsg4 for the hair shaft to maintain integrity [15][16][17][18] . Therefore, down-regulated Dsg4 and up-regulated PKP1 may play critical roles in the pathogenesis of AAP.
CD2, also known as LFA-2, acts as an adhesion molecule expressed on T cells and natural killer cells that binds CD58 on APCs, facilitating TCR binding and signal transduction. Meanwhile, various studies have demonstrated the important role of CD2 in T-cell activation, the production of cytokines by T-cells and T-or NK-mediated effector responses [19][20][21] . Besides this, the protein encoded by CD3D is part of the T-cell receptor/CD3 complex (TCR/ CD3 complex) and is involved in T-cell development and signal transduction. More importantly, an increasing number of studies show that anti-CD2 treatment that block the CD2-CD58 co-stimulatory signals could be the therapeutic strategies for dampening autoimmune diseases, since CD2-CD58 interactions promote the synthesis of IL-2 and interferon-γ (IFN-γ), as well as T and NK cells-mediated cytotoxicity 20,22,23 . Hence, although there is no direct research on anti-CD2 treatment to reverse the alopecia areata, the prospect is still cheerful.
In addition to CD2 and CD3D, another T cells surface antigen, CD69, an early T cells activation marker, has been reported to distinguish tissue-resident memory T cells from circulating subsets, and thus, the up-regulation of CD69 may trigger retention of memory CD4 + and CD8 + T cells in relevant tissues, resulting in autoimmune diseases such as rheumatoid arthritis, psoriasis as well as alopecia areata [24][25][26][27] . Therefore, our findings are consistent with the report by Suárez-Fariñas et al. (ref. 8 ). and Krueger et al.(ref. 28 ).
Also, it is interesting that CXCL9/CXCL10, known as the monokine and interferon γ-induced protein 10 (IP-10), respectively, are strongly induced by IFN-γ, leading to migration of immune cells to focal sites, such as scalp and hair follicles, through the CXCL9, CXCL10, CXCL11/ CXCR3 axis 29,30 . Besides this, Xing et al. (ref. 7 ) pointed out the potential clinical application of JAK inhibition in human AA by downregulating effectors of a series of cells surface cytokine receptors, including mainly IFN-γ and IL-2, which are the two most abundant cytokines produced by human CD4 + and CD8 + T cells. CCL5, one of the subfamilies of chemokines (C), plays a similar function to CXCL9 and CXCL10 in the recruitment and retention of T cells or other lymphocytes, contributing to this mainly T cell-mediated autoimmune hair loss disease. Collectively, our results support the view that CXCL9, CXCL10 and CCL5 are associated with up-regulation of inflammation of alopecia areata, rheumatoid arthritis and psoriasis, which was also involved with the JAK/STAT signaling pathway [31][32][33][34] .
Combined with the roles that the above 6 hub genes play, we verified that the CD2/CD3D/CD69 promote T cell activation, and subsequently CD4 + and CD8 + T cells produce cytokines IFN-γ and IL-2 to induce up-regulation of CXCL9, CXCL10 and CCL5. Our research confirms further the effect of T cells-mediated autoimmune and the upregulated cytokine pathways on hair follicle damage in alopecia areata.
In addition, the Keratin-associated protein 17-1 (KRTAP17-1) is a member of the ultrahigh sulfur keratinassociated protein (KAP) family 35 . The KAP proteins are involved in forming a matrix of keratin intermediate filaments that support the structure of hair fibers, in turn, the down-regulation of KRTAP17-1 may accelerate the pathogenesis of AAP. Granzyme A (GZMA) is an abundant serine protease in the cytolytic granules of cytolytic T cells and NK cells, inducing caspase-independent cell death when delivered into target cells by perforin 36 . Bone morphogenetic protein 2 (BMP2), an extra-follicular macro-environment modulator, inhibited the transition of hair cycle from telogen to anagen and kept the hair stem cells relatively dormant 37 . It means that the downregulated BMP2 could promote the hair regeneration. CD163, one member of the scavenger receptors, is exclusively expressed in monocyte/macrophage lineage, and increased soluble CD163 was significantly involved with the pathogenesis of diverse autoimmune diseases, including rheumatoid arthritis, multiple sclerosis and systemic lupus erythematosus 38 .

CONCLUSION
In summary, the present study identified DEG differences between AAP and normal healthy control with bioinformatics data mining and found the potential key genes and pathway of AAP. Of ten hub genes, CD2 may be selected as a new therapeutic target for AAP. In addition, the role of keratins in the pathogenesis of AAP should not be ignored. However, further validation of these related genes and pathways is needed to investigate the pathogenic mechanism of AAP.