The raw dataset of the gene expression profile GSE90594 generated by Michel et al. (5) was downloaded from the GEO database (6). The data was obtained from scalp biopsies taken from the vertex region of 14 young males with premature alopecia (age 29.4 ± 3.4 years, stage V–VII as per Hamilton-Norwood classification) and 14 healthy volunteers with less than 2% white hairs (age 26.1 ± 3.6 years, Stage I or II according to Hamilton- Norwood classification). Both the alopecia and healthy group did not have any other skin involvement, autoimmune disorders, and systemic diseases (5).

limma v3.50.3 (Linear Models for Microarray Data) package, a R/Bioconductor software package, which provides an integrated solution for analyzing gene expression data from microarray technologies was utilized for data analysis (7). The Data pre-processing included background correction using normexp method and quantile normalization. Boxplot and cluster analyses were performed to identify and remove outliers in the samples. Then the control probes and the unexpressed probes are filtered out while the probes that are expressed above background are retained for further analysis. In addition, for multiple probes corresponding to the same genes in the arrays their average expression value was computed by avereps function in limma. Then the DEGs for the alopecia samples compared to the healthy samples were mined using the single-channel design matrix provided in the limma package. Benjamini and Hochberg’s method was utilized to compute the adjusted p-values (False Discovery Rate, FDR) (8). The probes with adjusted p-value (FDR) < 0.05 were selected as differentially expressed.

ToppGene Suite¹ (updated: Mar 2021) was employed to perform gene ontology (GO) functional and pathway analysis to identify items in gene lists that may have relevance to the biological question being investigated (9, 10). The ToppFun function in the ToppGene Suite was utilized to carry out GO (biological process and molecular function), gene family (source: genenames.org), and pathway enrichment (source: Biosystems-Reactome) analyses for the DEGs. All genes that are detected in the microarray analysis were used as background gene set in the ToppGene Suite for these analyses. The probability density function which is the default method for p-value and FDR calculation was selected. Gene count >2 and FDR B&H q-value < 0.05 were chosen as the cut-off criteria for the analyses.

Windows of 50 kb, 100 kb and 500 kb flanking the 107 lead SNPs associated from 8 genome wide association studies [Study Accession IDs: GCST000250 (11), GCST000251 (12), GCST001548 (13), GCST001297 (14), GCST005116 (15), GCST90043616 (16), GCST003983 (17), and GCST90043619 (16)] for the trait androgenetic alopecia indexed in the NHGRI-EBI GWAS catalog were prepared by calculating coordinates of 50 kb, 100 kb and 500 kb distance on either sides from the SNP position. Gene coordinates of DEGs transcript(s) were annotated using RefSeq Identifiers (Hg38). The flanking coordinates of SNPs were overlapped with the coordinates of DEGs utilizing the intersect function in Bedtools v2.30.0 (18). An overlap is only considered when there is a minimum of 1 bp overlap between the coordinates of DEG transcripts and the flanking coordinates of the lead SNPs (19).

The promoter regions of up and down-regulated DEGs were separately subjected to motif analysis utilizing the gene-based analysis method in HOMER v4.11 software² (20). 2,000 bp upstream and 200 bp downstream relative to the transcriptional start site of the genes were considered as promoter regions (19) and the promoter sets for the DEGs were constructed based on RefSeq genes (Hg38). Motifs of length up to 12 bases were probed with Benjamini-Hochberg-corrected p-value ≤ 0.05 as cut-off value.

The protein-protein interaction (PPI) interaction network for the DEGs were computed through the STRING database. The online web resource STRING v11.5³ is a biological database that includes direct (physical) and indirect (functional) protein-protein association data which are both specific and biologically meaningful (21). The PPI interaction network for the DEGs were computed through the stringApp plugin v1.7.1 in Cytoscape v3.9.1 (22). An interaction score of 0.900 (highest confidence) was used as the cut off criterion for constructing the PPI network.

Reactome functional interaction (FI) network was constructed for the DEGs utilizing the Cytoscape application ReactomeFIViz v8.0.4 which probe for disease-related pathways and network patterns using the Reactome functional interaction (FI) network (23, 24) created based on the well-known biological pathway database Reactome⁴ (25, 26). Reactome FI network 2021 version was used to construct the FI network for the DEGs. Gene ontology biological process and pathway enrichment analysis for the nodes (genes) mapped in the network was carried out through the inbuilt Reactome FI network analysis tool.

The topological properties of the PPI and FI network were analyzed through the Cytoscape pre-installed network analyzer v4.4.8 tool (27). Cytohubba v0.1 plugin was used to identify hub proteins in the PPI and FI network and rank them based on topological algorithms and centralities such as Maximal Clique Centrality (MCC), Maximum Neighborhood component (MNC), Density of Maximum Neighborhood Component (DMNC), Degree, Closeness, and betweenness (28). The clusters in the networks were determined using the MCODE plugin with specific parameters including a degree cut-off of 2, fluff node density cut-off of 0.1, node score cut-off of 0.2, K-core of 2, and max depth of 100 to determine the highly interconnected nodes (29).

The GSE90594 dataset contained 28 samples of which 14 samples are from men with premature AGA and 14 samples from healthy men without hair loss. Cluster analysis of the samples after background correction and normalization of the arrays revealed 9 samples (5 alopecia and 4 healthy samples) as outliers (Supplementary I-1, 2). The outlier samples were removed, and differential gene expression analysis was carried out between 9 alopecia and 10 healthy samples. The probes were annotated with Entrez Gene ID, Gene Symbol, and Gene names using the clusterProfiler 4.0 v4.4.3 package in R by querying the Reference Seq ID (30). From this list, the probes that have a valid Entrez Gene ID are selected for further analysis. Subsequently probes with similar Probe IDs (Probe Names) are averaged using avereps function in limma and the probes with different probe ID for same genes are kept as such.

The analysis returned a total of 289 DEGs (33 up-regulated and 256 down-regulated DEGs) for a threshold cut off of value log₂FC > | 1| and q-value < 0.05 (Supplementary II-6). Further to construct a big and detailed STRING PPI and reactome FI networks a total of 2,439 unique DEGs (1,261 up-regulated, 1,171 down-regulated, and 7 genes with probes expressed in both directions) that falls within a cut off value of log₂FC > | 0.3| and q-value < 0.05 were mined (Figure 1). The gene family enrichment analysis of these 2,439 DEGs are given in Figure 2 and Supplementary II-7. GO functional analyses revealed that the up-regulated genes enriched for immune system mediated GO terms implying a heightened immune response in hairless scalp, while the down-regulated genes enriched for hair growth related GO terms as expected (Table 1 and Supplementary II-8). The Reactome pathway enrichment analysis also enriched pathways such as keratinization, developmental biology G2/M DNA replication checkpoint for down-regulated genes, wherein for up-regulated genes innate immune system, Cytokine signaling, interferon signaling, adaptive immune system, and antigen processing cross presentation pathways were enriched (Table 2 and Supplementary II-9). The DEG list was inspected for genes known to be involved in various signaling pathways such as Wnt, NF-κB, TGF-β, BMP, and Vitamin D metabolism and the mapped DEGs for these signaling pathways were provided in the Supplementary I-3.

Mapping of SNPs identified through GWAS to DEGs annotates genetic variants located in or near the gene regions that are differentially expressed, helps to understand the functional role of DEGs and their association with disease. To identify genes that potentially contribute to AGA pathology, we annotated the coordinates of 107 genomic loci associated with AGA risk in men identified through GWAS with our DEGs. The analysis identified 51 DEGs within the window of 500 kb of 73 AGA risk SNPs and 14 DEGs within the window of 50 kb of 16 lead SNPs (Table 3). Some of the DEGs were mapped to reported AGA risk SNPs including MEMO1 at loci 2p22.3, SRD5A2 at 2p23.1, FOXL2NB at 3q23, FGF5 at 4q21.21, DKK2 at 4q25, EBF1 at 5q33.3, IRF4 at 6p25.3, CENPW at 6q22.32, PAGE2 at Xp11.21, highlighting their association with AGA pathology. Our analysis also identified other DEGs, such as HOXD9 at 2q31.1, LHPP at 10q26.13, CRHR1 at 17q21.31, STH at 17q21.31, and PAGE2B at Xp11.21, as more likely to be candidate gene risks for AGA than the mapped genes in GWAS studies. MEMO1 was down-regulated, and it plays a crucial role in regulating cell proliferation, survival, and differentiation in the hair follicle (30). SRD5A2, whose product inhibits hair growth, was up-regulated (31). HOXD9, a member of the HOX family of genes that plays a crucial role in the development and patterning of various tissues and organs in the body, was up-regulated in our analysis, although its role in hair growth is unknown (32). FGF5, which inhibits hair growth and is involved in the transition of hair follicles from anagen to catagen phase, was down-regulated (33). DKK2, a Wnt inhibitor that leads to hair growth inhibition was up-regulated (1). FOXL2NB, IRF4, CENPW, EBF1, LHPP, CRHR1, and STH located within the AGA risk loci warrant further investigation. The mapped DEGs within 500 kb of lead SNPs may also be considered for future investigation of their association with hair growth and AGA (Table 3).

The transcription factor motif enrichment analysis on the promoter regions of the differentially expressed genes was carried to identify potential transcription factors involved in the AGA pathology. The top transcription factor motifs enriched for the down-regulated genes are LEF1 (Lymphoid Enhancer binding Factor 1), HOXB13 (Homeobox B13), NEUROD1 (Neuronal Differentiation 1), ZNF189 (Zinc Finger protein 189), and MEF2C (MADS Box Transcription factor 2, Polypeptide C) (Figure 3 and Supplementary II-10). The transcription factor LEF1 actively participates in the Wnt signaling pathway by activating the transcription of target genes in the presence of β-catenin. Wnt/β-catenin Signaling plays a crucial role in hair follicle differentiation and morphogenesis (31). The transcription factor HOXB13 belongs to HOX gene family which plays a crucial role in regulating embryonic development including hair formation. HOXB13 is implicated in skin development and low level of its expression is associated with telogen hair follicle (32, 34). The transcription factor NEUROD1 is primarily involved in the development and differentiation of the nervous system. NEUROD1 acts by controlling the expression of genes involved in neuronal development and in the formation of axons and dendrites (35). ZNF189 belongs to the zinc finger protein family which play important roles in various biological processes including transcriptional regulation, DNA repair, and cellular signaling. MEF2C belongs to the MADS box transcription factor 2 (MEF2) family of transcription factors and is involved in myogenesis (32). Many transcription factor motifs belonging to the SMAD, HOX, STAT, ZNF, NEURO, FOX, and FOS gene families (Figure 3) are enriched for the down-regulated genes indicating their role in hair growth which has to be studied further.

The motifs for IRF3 (Interferon Regulatory Factor 3), PRDM1 (PR/SET Domain 1), IRF8 (Interferon Regulatory Factor 8), SPI1 (Spi-1 Proto-Oncogene), SPI1:IRF8, ISRE (Interferon-sensitive response element), IRF2 (Interferon Regulatory Factor 2), IRF1 (Interferon Regulatory Factor 1) and SF1 transcription factors were enriched as the top motifs for the up-regulated genes (Figure 3 and Supplementary I-11). The Interferon regulatory factors (IRFs) are a family of transcription factors that regulate various aspects of the immune system from promoting immune cell development to immune cell differentiation. They play a central role in controlling the innate and adaptive immune responses to pathogens (33). IRF1 and IRF2 are important in regulating dendritic cells which participates in antigen presentation and bridge the innate and adaptive immune system. IRF3 involves in type I interferon production and IRF8 regulate myeloid cell development (33). PRDM1 coordinates several important functions in the adaptive immune system that support the key effector functions of B and T lymphocytes (36). SPI-1 encodes an ETS-domain transcription factor that control gene expression involving in the development of myeloid and B-lymphoid immune cells (37). The enrichment of these transcription factor motifs as the top motifs in the up-regulated genes of bald scalp implies a state of heightened immune response in AGA.

The stringApp generated 1967 PPI pairs for the submitted DEGs. The main PPI network, which consisted of 749 nodes (447 up-regulated genes, 273 down-regulated genes, and 29 linker genes) and 1,856 edges, was selected for further analysis while disconnected nodes and small isolated PPI pairs were discarded (Figure 4 and Supplementary II-12). The PPI network had a clustering coefficient of 0.334, a characteristic path length of 5.683, a network diameter of 19, a network density of 0.007, and an average of 4.956 neighbors. The functional enrichment analysis performed using the inbuilt STRING tool on the Reactome and Wikipathway databases revealed that the PPI network was enriched for several immune response-related pathways (Supplementary II-13). The pathway terms related to Cytokine Signaling in the Immune System, Interferon Signaling, T cell receptor signaling, Signaling by Interleukins, Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell, Adaptive and Innate immune systems, and TCF-dependent signaling in response to WNT were enriched from the Reactome database. Additionally, the Wikipathways database identified significant pathway terms related to the Inflammatory response pathway, Development of pulmonary dendritic cells and macrophage subsets, B cell receptor signaling pathway, and the Vitamin D receptor pathway (Supplementary II-13). These results confirm the credibility of the PPI network and reinforce the observation of immune response-related and hair follicle-related pathways.

The top 20 ranking hub nodes (genes) in the PPI network were identified using the Cytoscape plugin Cytohubba based on four topological analysis methods and two centralities (MCC, DMNC, MNC, Degree, Closeness, and Betweenness) and are listed in Table 4. Out of these, a total of 15 hub genes that appeared in at least three of these categories were considered as significant hub genes, and the frequently appeared genes are highlighted in the Table 4. The MCODE cluster analysis performed on the String PPI network revealed 8 clusters when using the 15 hub genes as roots for clustering. The 8 clusters were comprised of 53, 63, 101, 59, 37, 41, 53, and 9 nodes, respectively (Supplementary II-14). The top 3 highly interconnected clusters were selected for further analysis (Supplementary 1–4). Cluster 1 had 10 hub genes, cluster 2 contained 5 hub genes, and cluster 3 had 6 hub genes. Our analysis revealed that two hub genes LCK and STAT5A appeared in all 3 clusters strongly suggesting their putative role in AGA.

The ReactomeFIViz tool was utilized to construct the FI network for the DEGs resulting in an initial network of 1,092 connected nodes, 1,340 unconnected nodes, and 4,047 edges (Supplementary II-15). The unconnected nodes were discarded from the analysis and the final FI network consisted of 1,014 nodes (581 up-regulated and 433 down-regulated genes) with 3,980 edges as shown in Figure 5. The FI network had a clustering coefficient of 0.280, a network diameter of 11, a network density of 0.008, and an average number of neighbors of 7.850. The pathway enrichment and GO Biological process analyses were conducted using the inbuilt ReactomeFIViz – analysis network function tool (Supplementary II-16). Reactome pathway terms such as Extracellular matrix organization, Keratinization, Interferon gamma signaling, Regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins, Response to elevated platelet cytosolic Ca2 +, Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell, and WNT ligand biogenesis and trafficking were enriched for the genes present in the FI network. GO Biological process terms such as immune response, transmembrane receptor protein tyrosine kinase signaling pathway, canonical Wnt signaling pathway, and inflammatory response were also enriched.

The top 20 ranked hub genes in the FI network identified based on the six algorithms including MCC, DMNC, MNC, Degree, Closeness, and Betweenness are presented in the Table 5. Out of these, a total of 19 hub genes that appeared in at least three of the categories were considered significant hub genes, and the frequently appeared genes are highlighted in the Table 5. The MCODE cluster analysis of the FI network revealed 11 clusters when using the 19 hub genes as roots for clustering. The 11 clusters had node numbers of 189, 55, 216, 47, 34, 59, 153, 180, 29, 69, and 6, respectively (Supplementary II-17). The top 3 clusters ranked based on their cluster score were selected for further analysis (Supplementary I-5). Cluster 1 consisted of 14 hub genes, cluster 2 contained 1 hub genes, and cluster 3 had 9 hub genes. Our results showed that seven hub genes (HCK, GNAI3, RAC2, PDGFRB, EGF, NRAS, and STAT5A) were present in two of the selected clusters indicating their potential role in AGA.

The 25 hub genes identified from the analyses of the PPI and FI networks constructed based on the DEGs were considered key genes in the pathology of AGA (Table 6). Out of these 25 hub genes, 21 genes (BTK, ESR1, HCK, ITGB7, LCK, LCP2, LYN, PDGFRB, PIK3CD, PTPN6, RAC2, SPI1, STAT3, STAT5A, VAV1, PSMB8, HLA-A, HLA-F, HLA-E, IRF4, and ITGAM) were found to be up-regulated, while 4 genes (CTNNB1, EGF, GNAI3, and NRAS) were downregulated. The results of the GO biological process and pathway enrichment analysis, conducted using the Toppgene suite, revealed that the hub genes were associated with immune and inflammatory processes (Table 7). The significant biological terms enriched for the hub genes included regulation of immune system process, T cell activation, immune response-regulating cell surface receptor signaling pathways. Furthermore, the significant pathway terms enriched for the hub genes included cytokine signaling in the immune system, signaling by interleukins, signaling by the B cell receptor (BCR), and signaling by SCF-KIT, interleukin-3, 5, and GM-CSF signaling, and the innate immune system.

CTNNB1 (Catenin Beta 1) is a crucial downstream component of the Canonical Wnt Signaling Pathway. In the presence of Wnt ligand, β-catenin accumulates in the nucleus and functions as a coactivator for the transcription factors TCF/LEF, leading to the activation of Wnt responsive genes (35). The Wnt/β-catenin signaling pathway is essential for hair growth and its inhibition, driven by 5α-dihydrotestosterone through the androgen receptor, can result in hair loss in AGA (1). GNAI3 (G Protein Subunit Alpha I3) functions as a downstream transducer of G protein-coupled receptors (GPCRs) in various signaling pathways (38). GPCRs play a role in regulating skin homeostasis and maintaining hair growth (39–41). NRAS (NRAS Proto-Oncogene, GTPase) is a membrane protein that travels between the plasma membrane and Golgi apparatus (42). EGF (Epidermal Growth Factor) acts as a switch in the hair growth cycle (43). It regulates the expression of hair follicle regulatory genes through Wnt//β-catenin signaling (44). Thus, these 4 downregulated hub genes which are involved in hair growth mechanisms are crucial and their downregulation in AGA is expected.

The LCK (LCK Proto-Oncogene, Src Family Tyrosine Kinase) gene, which encodes a non-receptor protein-tyrosine kinase, is a crucial signaling molecule in the selection and maturation of developing T cells and plays a key role in T cell receptor signal transduction pathways (25, 26). The up-regulation of the LCK gene is also associated with alopecia areata (27). The LYN (LYN Proto-Oncogene, Src Family Tyrosine Kinase) gene encodes a non-receptor tyrosine-protein kinase and is crucial for regulating innate and adaptive immune responses, integrin signaling, growth factor and cytokine responses, and hematopoiesis (24). BTK (Bruton Tyrosine Kinase) and plays a key role in B lymphocyte development and is a target for inflammatory diseases (45). Inhibition of BTK by inhibitors leads to changes in hair and nails texture (38). PIK3CD (Phosphatidylinositol-4,5-Bisphosphate 3-Kinase Catalytic Subunit Delta) is involved in immune system response (46). PTPN6 (Protein Tyrosine Phosphatase Non-Receptor Type 6) is critical for the function of lymphoid and myeloid cells (47). SPI1 (Spi-1 Proto-Oncogene) encodes a transcriptional activator specifically involved in the development of macrophages and B cells. This protein also regulates pre-mRNA splicing (23). STAT3 (Signal Transducer and Activator of Transcription 3) is activated by cytokines and growth factors. This gene plays an important role in maintaining the homeosis of skin (48). STAT5A (Signal Transducer and Activator of Transcription 5A) protein serves a dual function of signal transduction and activation of transcription in cells exposed to cytokine and other growth factors. This protein also mediates cellular responses to activated FGFR1, FGFR2, FGFR3 and FGFR4 (31). Also, STAT5 activation is important for hair growth phase induction in hair dermal papilla cells (DPCs) (34). VAV1 (Vav Guanine Nucleotide Exchange Factor 1) encoded protein is important in hematopoiesis and plays a role in the development and activation of T-cell and B-cell (32). PSMB8 (Proteasome 20S Subunit Beta 8) plays an important role in cellular homeostasis through selective destruction of ubiquitinated proteins. Mutations in this gene are associated with autoinflammatory responses (49). ESR1 (Estrogen Receptor 1) is a nuclear sex steroid hormone receptor which regulates many genes responsible for growth, metabolism and reproductive functions. This gene is known to express in hair follicle cells (50). HCK (HCK Proto-Oncogene, Src Family Tyrosine Kinase) participates in the regulation of innate immune responses by inducing monocyte, neutrophil, macrophage and mast cell functions. This gene is recently reported to play a role in hair regenerative potential of stem cells (51). ITGB7 (Integrin Subunit Beta 7) is an adhesion receptor which mediates signaling from the extra cellular matrix to the cell. They also function as a homing receptor for lymphocytes migration (46). LCP2 (Lymphocyte Cytosolic Protein 2) acts as a substrate for the T cell antigen receptor mediated intracellular tyrosine kinase pathway (46). PDGFRB (Platelet Derived Growth Factor Receptor Beta) gene encodes a cell surface tyrosine-protein kinase receptor for the members of the platelet-derived growth factor family. It plays an essential role in cell proliferation, differentiation, survival, chemotaxis, and migration (52). RAC2 (Rac Family Small GTPase 2) involve in phagocytosis of apoptotic cells and epithelial cell polarization (46). IRF4 (Interferon Regulatory Factor 4) regulates interferon signaling and negatively regulates Toll like receptor in the induction of innate and adaptive immune systems (42). ITGAM (Integrin Subunit Alpha M) functions as macrophage receptor and plays a key role in the adherence of monocytes and neutrophils (42). HLA-A (Major Histocompatibility Complex, Class I, A), HLA-F (Major Histocompatibility Complex, Class I, F) and HLA-E (Major Histocompatibility Complex, Class I, E) plays a central role in immune system by participating in cell presentation for recognition by T cell receptor (42). A majority of the hub genes namely PTPN6, LCK, LCP2, LYN, HCK, VAV1, STAT3, STAT5A, and BTK belongs to the Src homology 2 (SH2) domain containing tyrosine kinases and participate in the immune system process.

We conducted ClueGO reactome pathway enrichment analysis for the genes that were identified by at least two algorithms of the Cytohubba analysis of the biological networks (PPI and FI) as well as 289 DEGs that met the cut-off value of log₂FC > | 1| using the ClueGo plugin v2.5.9 in Cytoscape (53). The results were presented as a network of pathways with genes participating in the pathways, which are illustrated in Figure 6. The analysis revealed pathways such as keratinization, formation of the cornified envelope, developmental biology, interferon alpha/beta signaling, cytokine signaling in the immune system, receptor tyrosine kinase signaling, PI5P, PP2A, and IER3 regulation of PI3K/AKT signaling, immunoregulatory interactions between lymphoid and non-lymphoid cells, costimulation by the CD28 family, and the GPVI-mediated activation cascade are the predominant pathways for our input genes. The enrichment of pathways involved in immune system function are consistent with our findings suggesting that immune system dysregulation plays a role in AGA pathology.

To validate the results of our analysis of the GEO dataset GSE90594, we compared the DEGs obtained with other datasets available in the GEO database. As of November 1, 2022, we found that no profile in the database contained samples from men with AGA and from normal haired men, except for the profile we analyzed in this study. The few available datasets related to AGA lacked control samples from normal men and the quality of the microarray and RNA-Seq data was questionable. Despite these limitations, we selected two datasets, GSE66663 (which includes hTERT-immortalized DPCs derived from balding frontal and non-balding occipital scalp samples from men with AGA) and GSE212301 (which contains RNA-Seq data from balding vertex and non-balding occipital scalp samples of 10 men with AGA), and performed differential gene expression analyses. The common DEGs between the datasets are presented in the Supplementary II-18. We discovered 490 DEGs were common between GSE90594 and GSE66663 dataset in which 190 genes were differentially regulated in same directions. Whereas 180 DEGs were common between GSE90594 and GSE212301 dataset in which 44 genes were differentially regulated in same directions.