Abstract Objectives This study aimed to identify significant mechanisms and potential treatments for temporomandibular joint internal derangement (TMJD) through integrated bioinformatics analysis. Materials and Methods Gene expression data sets ([24]GSE66864) from the Gene Expression Omnibus (GEO) database were downloaded. Differentially expressed genes (DEGs) were identified both in the treatment groups and in controls by R packages. Network analysis of protein–protein interaction (PPI) and Human Protein Atlas was used to explore DEGs' potential function. DGIdb database was utilized to gain potential drug targets. Results In conclusion, 126 DEGs were selected for TMJD through bioinformatics analysis. Both GO and Kyoto Encyclopedia of Genes and Genomes analyses combined showed the pathways involved in TMJD. A PPI network was constructed to select the top 10 hub genes, of which five hub genes were chosen for further investigation. Moreover, the microenvironment of immune cells related to hub genes was evaluated by R packages. Macrophages play an important role in inflammation and oral‐related tumors. The Human Protein Atlas analysis indicated that the five hub genes are highly related to head and neck cancer. Finally, eight potential drugs were selected for two genes using the DGIdb database. Conclusion Our integrated bioinformatics analysis identified DEGs in TMJD and provided potential ideas for further research and treatment approaches. However, experimental validation of the hub genes and potential drug targets is still needed. The key mechanisms of the identified genes and their potential roles as biomarkers or drug targets in TMJD require further investigation. Keywords: hub genes, immune cells, integrated bioinformatics analysis, temporomandibular joint internal derangement 1. INTRODUCTION Temporomandibular joint Internal derangement (TMJD) is one kind of the most complex human diseases (Granados, [25]1979); it is a common disease affecting more than 15% of adults, peaking in incidence among those aged 20–40 years (Gauer & Semidey, [26]2015). The main etiology of TMJD is complex, including biological, social, and emotional factors. TMJD is associated with primary headaches, neurological disorders, dental occlusion confusion, and a series of oral–maxillofacial diseases. Previous studies focus on the treatment of TMJD, such as Arthrocentesis versus glucocorticosteroid injection (AbdulRazzak et al., [27]2021). Few studies discuss the biological pathway of TMJD. As a result, currently, the molecular process associated with pathways in the TMJ is still unclear. Therefore, we aim to elucidate the biological progression of TMJD by bioinformatic analysis. As a complex multifactorial disease, multiple pathways of hub genes and biological procession are involved in TMJD. High‐throughput sequencing technologies are increasingly used in tumor genomics and can also provide insights into nontumor diseases. Bioinformatics analysis has been shown to discover potential early diagnostic biomarkers for osteoarthritis (Cao et al., [28]2021). However, there is little research about TMJD on the gene expression profile and the biological pathways to explore critical signaling pathways and some essential hub genes for the process, complication, and further development of TMJD. Therefore, this study aims to figure out the hub genes and key biological pathways of TMJD, thus opening a novel perspective on this problem. Elucidating the molecular basis of TMJD can guide the development of targeted therapies and diagnostic tests for this complex disease. A better understanding of the biological mechanisms underlying TMJD may help clinicians manage the condition more effectively. In this study, we identified differentially expressed genes (DEGs) from the Gene Expression Omnibus (GEO) database in published transcriptomic data sets of TMJD in the condylar cartilage of rabbits from the GEO database ([29]GSE66864). DEGs that were highly concerning TMJD underwent Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and the weighted gene coexpression network analysis (WGCNA). Protein–protein interaction (PPI) network analysis was utilized to figure out the hub genes, potential key biological procession, and signaling pathways of TMJD. Rabbit genes share similarities with human genes, including the class I major histocompatibility complex gene (Marche et al., [30]1985). Hence, this study utilized the Human Protein Atlas to analyze the key protein function in the human tissue, hoping to find a similar biological progression in humans. Altogether, an integrated analysis was performed in this study to identify the hub genes of TMJD and illustrate the potential biological procession and signaling pathways for the progression of TMJD. Some information can be gathered for the diagnosis, monitoring, and treatment of TMJD. Some insights may aid the diagnosis, monitoring, and treatment of TMJD, providing a novel perspective for future TMJD research. 2. MATERIALS AND METHODS 2.1. Data acquisition The [31]GSE66864 gene expression series matrix was downloaded from the GEO database ([32]https://www.ncbi.nlm.nih.gov/geo/). Raw database [33]GSE66864 was presented by the [34]GPL13288 platform, including one control group and four treatment groups. Four rabbits were selected in one group as control (without any treatment); 16 rabbits were divided randomly into four treatment groups, which were used for histological analysis, RNA microarray analysis, and real‐time PCR at 1 week, 2 weeks, 4 weeks, and 8 weeks, respectively (Wang et al., [35]2015). 2.2. DEG identification Background calibration, normalization, and data filter were presented to DEGs through the R software (version 4.1.3) and GEO2R ([36]https://www.ncbi.nlm.nih.gov/geo/geo2r). DEGs with a statistically significant difference were identified through standard adjustment p < .05 and log fold‐change (log FC)| ≥ 2. R package “ggplot2” and “ggrepel,” were used to plot the volcano map. A heatmap was generated using TBtools (Chen et al., [37]2020). WGCNA was utilized to confirm the reliability of DEGs in the different samples. 2.3. Analysis of biological functions and signaling pathways The GO terms and KEGG pathways of DEGs concerning the analysis of biological functions were enriched via the R package “clusterProfiler” (Yu et al., [38]2012). Three types of biological analyses were performed in the GO terms: biological process (BP), cellular component (CC), and molecular function (MF) in GO terms. The p [adj] < .05 and q< .05 were set for the enriched DEGs and their further data visualization as a threshold. Gene enrichment in KEGG pathways shows similar results. 2.4. Construction of PPI network In our study, we established the PPI network that contained the hub genes via the STRING database ([39]https://cn.string-db.org/) (version 11.5) (Szklarczyk et al., [40]2019), and the minimum required interaction score was set at 0.4. The differentially expressed proteins were all utilized to construct PPI visualization networks and submit the DEGs by Cytoscape (Shannon et al., [41]2003). The top 10 DEGs were selected as nodes based on the Degree algorithm (Chin et al., [42]2014). 2.5. Analysis of immune cell infiltration The immune cell infiltration was analyzed using cell type identification by estimating related subsets of RNA transcripts (CIBERSORT) (Newman et al., [43]2015). The relationship between the hub gene and the 547 genes, which distinguish 22 human immune cell phenotypes, was visualized and performed using a heatmap. The hub genes were also explored between the immune cell infiltration and human tissue. The hub genes were submitted to the Human Protein Atlas ([44]https://www.proteinatlas.org/) to identify similarities in human tissue. 2.6. Drug interaction of hub genes The DGIdb database ([45]http://dgidb.genome.wustl.edu/) was utilized to figure out potential drug targets and potential therapeutic candidates for the upregulated hub genes. The DGldb database was established with various sources of drugs, including immunotherapy drugs and antineoplastic drugs (Cotto et al., [46]2018). 3. RESULTS 3.1. Identification of DEGs One control group and four TMJD groups were contained in the [47]GSE66864 database. We plotted a heatmap based on the expression of genes in each sample (Figure [48]1a). A total of 104 upregulated and 22 downregulated genes were shown in Figure [49]1b. Five hub genes were selected including ITGA6, PTHLH, FPR1, ANXA8, and ATP6V0A2, which are highly upregulated in the DEGs. The WGCNA of sample clustering showed no significant difference between each sample (Figure [50]1c). The soft threshold was set to 7 (R ^2 = 0.8) to make the construction of the scale‐free network (Figure [51]1d). Several modules were identified based on average hierarchical clustering and dynamic tree clipping (Figure [52]1e,f). The brown‐red module had a high correlation with pathological grades. Figure 1. (a) Identification of differentially expressed genes (DEGs) in temporomandibular joint internal derangement (TMJD). Heatmap of top 500 DEGS from [53]GSM1633547 to [54]GSM1633566 in [55]GSE66864. Red past represented the genes that are upregulated (up), while green parts represented the genes that are downregulated (down). (b) Volcano plot of DEGs. Red represents upregulated DEGs; blue represents downregulated DEGs; gray represents genes with no significant difference. (c) Clustering dendrogram of 20 samples. (d) Analysis of the scale‐free index with various soft‐threshold power. R ^2 is equal to 0.8, which is mapped out by the red line. (e) Clustering dendrogram of all samples. (f) Network heatmap plot of the selected key genes modules and correlation between each module. NA, not available; not sig., not significant. graphic file with name CRE2-9-641-g003.jpg graphic file with name CRE2-9-641-g005.jpg [56]Open in a new tab 3.2. Analysis of biological functions and signaling pathways The GO enrichment analysis via dot plot showed that muscle system processes were highly associated with the BP group genes, regulation of ion transmembrane, and procession of muscle contraction (Figure [57]2a). In addition, the CC group genes were especially associated with the cation channel complex, ion channel complex, sarcomere, myofibril, contractile fiber, transmembrane transporter complex, and transporter complex (Figure [58]2b). MF groups were highly associated with voltage‐gated ion channel activity, voltage‐gated channel activity, voltage‐gated cation channel activity, metal ion transmembrane transport, gated channel activity, and passive transmembrane transport (Figure [59]2c). The KEGG pathway enrichment analysis shows a similar result of hub genes between human and rabbits (Figure [60]2d), which inspire us to associate rabbit genes obtained in this study with human genes. Figure 2. Figure 2 [61]Open in a new tab (a) Cellular component (CC), (b) biological process (BP), (c) molecular function (MF), and (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GO, Gene Ontology. 3.3. Construction of PPI network The 126 key genes were submitted to the String database, constructing PPI networks under the threshold of 0.4 interaction score (Figure [62]3a). Cytoscape software was used to illustrate the interaction result. The top 10 genes were selected by Degree, which is one of the algorithms: ATP6V0A2, PTHLH, FPR1, ANXA8, ITGA6, PLG, SERPINF2, ATP6V0A2, ITGB4, and PTHLH1R (Figure [63]3b). Figure 3. Figure 3 [64]Open in a new tab (a) Protein–protein interaction (PPI) network was constructed based on all the 126 differentially expressed genes (DEGs) by using the STRING database under the threshold of 0.4 interaction score. (b) The PPI network of DEGs in temporomandibular joint internal derangement (TMJD) (STRING). The top 10 genes included ATP6V1A, ATP6V0D1, ATP6V1B2, ATP6V1D, ATPV1E1, ATP6V0A2, ATP6V1B1, ATP6V0D2, ATP6V0C, and ITGAV, with a high degree of connections. 3.4. Immune cell infiltration in TMJD The bars depict the abundance of 22 immune cell infiltrates in the temporomandibular joint tissue samples (Figure [65]4a). Obviously, tumor‐associated macrophages and inflammation were detected at high levels, including T cells CD4 memory activated, B cell memory, and dendritic cells activated. The Human Protein Atlas was used to find the information about the hub genes, including ITGA6 (Figure [66]4b), FPR1 (Figure [67]4c), PTHLH (Figure [68]4d), ATP6V0A2 (Figure [69]4e), and ANXA8 (Figure [70]4f). The database shows that the ITGA6 and PTHLH are highly relevant to head and neck cancer, which are specifically expressed in the tissue. Meanwhile, the FPR1 is enriched in Immune cells (neutrophils). Figure 4. (a) Immune cell infiltration and abundance percentage of immune infiltrating in temporomandibular joint internal derangement (TMJD) tissues. (b) Percentage of ITGA6 expression of cancer tissues from high to low. Head and neck cancer displayed the highest expression of ITGA6. (c) FPR1 expression of immune cell types from high to low. Neutrophils displayed the highest expression of FPR1. (d) PTHLH expression of cancer tissues. Head and neck cancer displayed the highest expression of PTHLH. (e) ATP6V0A2 expression of immune cell types from high to low. Eosinophil displayed the highest expression of ATP6V0A2. (f) Survival probability of ANXA8 in patients between high‐expression group and low‐expression group. DC, dendritic cell; NK, natural killer; PBMC, peripheral blood mononuclear cell; TCGA, The Cancer Genome Atlas. graphic file with name CRE2-9-641-g004.jpg graphic file with name CRE2-9-641-g001.jpg [71]Open in a new tab 3.5. Drug interaction of hub genes In our study, highly expressed hub genes were submitted to the DGIdb database ([72]https://dgidb.org/) and eight potential drug interactions were obtained (Table [73]1). Two drugs interact with PTHLH, while five other drugs interact with FPR1. No drug interactions were found with ITGA6, ANXA8, and ATP6V0A2. The drug CAL functions as an inhibitor for PTHLH, and the drug CHEMBL1290251 functions as an agonist for FPR1. Table 1. Eight potential drug interactions between DEGs. Gene Drug Interaction types Sources PTHLH Vinblastine NCI PTHLH CAL Inhibitor ChemblInteractions FPR1 Penicillin G potassium DTC FPR1 CHEMBL1290251 Agonist DTC FPR1 Phenylbutazone DTC FPR1 Sulfinpyrazone DTC FPR1 CHEMBL1290365 DTC FPR1 Rebamipide TTD [74]Open in a new tab Abbreviation: DEG, differentially expressed gene. 4. DISCUSSION TMJD is a common joint disease often associated with TMJ synovitis (Ibi, [75]2019). Besides, TMJD also has a noninflammatory origin. Open surgical intervention is a traditional treatment for TMJD (Buckley et al., [76]1993). However, the open surgery did not offer a pleasing performance in patients' experience of postoperative recovery. Given that TMJD has complex and sometimes controversial etiologies, elucidating the potential biological mechanisms of TMJD is essential to develop novel diagnostic strategies for TMJD patients. In this study, biological analysis of a gene expression series matrix was conducted and identified 104 genes that are highly expressed and 22 genes that are downregulated between TMJDs and normal tissues. DEGs were submitted to the GO and KEGG analysis to explore the potential pathogenesis and developmental mechanisms of TMJDs. Five hub genes were selected including ITGA6, PTHLH, FPR1, ANXA8, and ATP6V0A2, which are upregulated in the DEGs. The upregulated genes are highly related to a muscle system process, regulation of ion transmembrane, voltage‐gated channel activity including PI3K‐Akt signaling pathway, immune responses, and inflammatory responses. Many researchers have confirmed that the PI3K‐Akt signaling pathway may be involved with TMJD caused by pressure and inflammation by increasing the cadherin‐11 in synovial fibroblasts(Wu et al., [77]2013). Meanwhile, some research also shows that the activation of neural pathways, for instance, the agents such as serotonin and nitric oxide mediated the pain of TMJD (Patil & Kirkwood, [78]2007). A PPI network was constructed by submitting the hub genes to the STRING database and Cytoscape to identify the main proteins related to TMJD. To better illustrate the mechanism underlying TMJD pathogenesis, light was shed on the immune cell infiltration of the TMJD tissue via the R package and discovered the biological functions of these hub genes in the inflammation procession. The immune cell infiltration is intricately linked to the tumor microenvironment (Gajewski et al., [79]2013). Our results demonstrate that TMJD tissue had higher immune cell infiltration compared to normal tissue. The level of T‐cell CD4 memory activated was greatly enhanced. Some researchers have found that the infiltration of memory CD4+ T cells is directly implicated in rheumatoid arthritis, which is a typical arthritis with inflammation and a high amount of CD4+ T cells always releases a high level of interleukin‐17 (IL‐17) and tumor necrosis factor‐β (TNF‐β), which are the most important inflammatory factors (Chemin et al., [80]2019). The high frequencies of the expression of CD4+ T cells always relate to chronic inflammatory syndromes. Our study has discovered the relationship between TMJD and the infiltration of CD4+ T cells. We hypothesize that the CD4+ T cells may have a similar function in the TMJD tissues. Some research studies indicate that the TNF‐α has an impact on inhibiting the CD28 gene expression (Bryl et al., [81]2001). The fusion protein CTLA4 can prevent the T cells from activating by binding to CD80 and CD86 on antigen‐presenting cells and blocking cohesive the process in the engagement of CD28 (Kremer et al., [82]2003). According to research conducted by Joel M. Kremer et al. ([83]2003), the upregulated hub genes may provide a new vision of biological treatment. So, the hub genes were conducted to the DGIdb database, and two and six drugs interacting with PTHLH and FPR1, respectively, were uncovered through the DGIdb database, which gives the possibility of inventing potential therapeutic candidates for TMJD. ITGA6 encodes a member of the integrin ⍺‐chain family of proteins, which are involved in cell surface adhesion and signaling. The α6β1 integrin may promote tumorigenesis. Some research studies show that the high level of transcriptions of the ITGA6 can be specifically found in the tumor tissues, displaying an oral‐cancer‐related biomarker (Lo et al., [84]2012). On the contrary, the low expression of ITGA6 can impair the activity of the PI3K/AKT pathway (Chen & Zhang, [85]2022) to exert a tumor‐suppressive function. ITGA6 can also promote cell invasion by changing the microenvironment of the extracellular matrix (ECM) (You et al., [86]2021), suggesting a new way for oral cancer treatment by silencing the molecule. Meanwhile, ITGA6 is involved in macrophage function, particularly migration and M2 activation (Sima et al., [87]2016), and the mechanisms of proinflammatory cytokines, which are derived from M1 and M2 macrophages. Macrophages' activation can be inhibited by IL‐37 through the NLRP3 pathway and by suppressing the lipopolysaccharide and interferon‐γ‐induced extracellular signal‐regulated kinase (ERK) and nuclear factor‐κB (NF‐κB) activation in human M1 macrophages (Luo et al., [88]2020). The PTHLH gene is related to the development of endochondral bone. Studies showed that PTHLH is highly expressed in oral cancer (Lv et al., [89]2014) and tongue cancer (Suwa et al., [90]2014). PTHrP, which is the protein encoded by PTHLH, is mostly produced and expressed by many tumor tissues and just a few normal cells. The normal tissues have a low concentration of PTHrP (Park & McCauley, [91]2012). Cell proliferation was reported to be highly associated with the upregulation of the PTHLH gene (Hameetman et al., [92]2005). Estrogen can also affect PTHLH expression in the tissue which may enhance the tumor growth in the cartilaginous tissue. Some report also shows that Emodin can suppress the PTHLH expression in related cancer by suppressing the TCF4/TWIST1‐induced complex (Fang et al., [93]2022). In this study, the PTHLH is upregulated in the DEGs; we can infer that the high expression of this gene is linked to oral‐related cancer. The FPR1 gene found in mammalian phagocytic cells encodes a G‐protein‐coupled receptor. The protein plays an important role in host defense and inflammation. FPRs are expressed in nearly all kinds of immune cells, including neutrophils, macrophages, and natural killers. Among those, the neutrophils were mostly influenced by the FPR1 gene, attracting more leukocytes to migrate to the sites of inflammation (Leslie et al., [94]2020). Annexin A1 (ANXA1) is a ligand of FPR1 relevant for cancer immunosurveillance (Vacchelli et al., [95]2020). In this study, FPR1 is highly expressed in the temporal–mandibular joint tissue, which may cause inflammation by leukocyte migration. The high expression of FPR1 has also been proven to activate biological procession that contributed to chronic inflammation such as ERK, mitogen‐activated protein kinase, NF‐κB, and signal transducer and activator of transcription 3 (Cao & Zhang, [96]2018). The FPR1 expression can predict early cancer and function as a pharmacologic target for innate immune responses (Liu et al., [97]2012); it plays a significant role in phagocyte accumulation and promotes innate immunity. High ANXA8 expression suggests that it may act as a receptor for FPR1. Further research is needed to test this hypothesis. ANXA8, the gene encodes a member of the annexin protein family that binds to cell membranes. Rosenbaum et al. ([98]2011) discover that the ANXA8 lacks the ability to bind to phosphatidylserine (PS) and apoptotic cells. If more apoptotic cells are cleared and identified, it may suppress the inflammation reaction by preventing the release of detrimental contents in cells. The loss of PS in the cells is the biomarker that the early apoptotic cells contained. The receptors in the phagocyte cell surface consist of the family members of scavenger receptors; moreover, TIM1 and TIM4 can also bind to PS straight on dying cells. ANXA5 especially interacted with PS, which is a specific partner (van Genderen et al., [99]2006), while ANXA8 lacks the ability, which may explain why the high expression of ANAX8 is the factor of inflammation and cell apoptosis. The high expression of ANXA8 is also associated with oral squamous cell carcinoma (OSCC) by metastasizing the cervical lymph node, suggesting potential treatment avenues. Some reports detected the expression of ANXA8 expression in the different groups of cervical lymph nodes (Oka et al., [100]2016) and found that the ANAX8 messenger RNA expression was rapidly detected in the metastasis‐positive group and rarely detected in the normal group. In some reports, ANXA8 was frequently overexpressed in pancreatic cancer (Karanjawala et al., [101]2008), but its role in other cancers, including OSCC, is unclear, including OSCC and other oral‐related cancer. In our study, we also found the survival probability between high‐ and low‐expression patients who have oral cancers. The high‐expression group has a lower survival probability. All of these pieces of evidences may indicate that the ANXA8 may play a crucial role in inflammation, cell apoptosis, and oral cancer. Further research is needed to validate these mechanisms. ATP6V0A2, the gene encodes an enzyme that contributed to the formation of acid microenvironment in cells. Studies suggest ATP6V0A2 upregulation triggers macrophage transport, initiating an inflammatory state in cells (Jaiswal et al., [102]2012). During the window of implantation, the expression of V‐ATPase has an important function in the communication between cells and trophoblast invasion (Satokata et al., [103]1995). Moreover, ATP6V0A2 may also relate to cancer cell ECM and control tumor metastasis (Katara et al., [104]2018) by damaging the ECM proteins by way of defective glycosylation and producing the compromised ECM. PHY34, a synthetic small molecule can target and inhibit the ATP6V0A2 subunit of V‐ATPase (Salvi et al., [105]2022) and significantly reduce tumor burden. ATP6V0A2 also expresses in numerous immune cells, especially highly expresses in eosinophils according to our research through the Human Protein Atlas. Several limitations still exist in this study. First, we obtained the data from the GEO database in this study, while there is a lack of clinical experimental research certifications for hub genes. A small number of samples with a total of one control group and four treatment groups are presented in our study, for the relative biological process of the disease still remained undiscovered. Second, our integrative bioinformatics approach requires further experimental validation to fully explore TMJD pathogenesis. However, our study provides insights into potential molecular mechanisms, biomarkers, and drugs for TMJD. We plan experimental confirmation of gene expression changes and further clinical investigation of DEGs to develop new treatments and strategies for TMJD. Despite these caveats, our study gives out further exploration into the significant molecular mechanisms of TMJD, as well as searches for potential biomarkers and therapeutic candidate drugs, providing novel treatment and biological targets for TMJD patients. We planned to perform more experiments to confirm the process of gene expression. Further clinical investigations on the DEGs may help develop the treatment and discover new diagnostic and therapeutic strategies for TMJD. 5. CONCLUSION This study provides new insights into the biological factors underlying TMJD. Using integrative bioinformatics analysis, we identified five hub genes (ITGA6, PTHLH, FPR1, ANXA8, and ATP6V0A2) that likely contribute to TMJD pathogenesis. These genes were found to be involved in processes including inflammation, ECM regulation, immune responses, and tumorigenesis, shedding light on potential mechanisms of TMJD. We also identified potential drugs that target two of the hub genes, representing candidate therapeutics for TMJD. However, further experimental validation is needed to confirm the roles of these genes in TMJD and evaluate the efficacy of potential drug treatments. In summary, this study represents an important first step in understanding the biological processes and molecular mechanisms involved in TMJD. Future research should experimentally test the findings and advance the development of effective TMJD treatments and diagnostic strategies. With additional investigation and confirmation, our work may help to improve TMJD management and patient outcomes. AUTHOR CONTRIBUTIONS Junda Yang: Conceptualization; methodology; software; data curation; writing—original draft preparation. CONFLICT OF INTEREST STATEMENT The author declares no conflict of interest. Yang, J. (2023). Integrated bioinformatics analysis of differentially expressed genes in the temporomandibular joint internal derangement. Clinical and Experimental Dental Research, 9, 641–652. 10.1002/cre2.768 DATA AVAILABILITY STATEMENT Data are openly available in a public repository that issues data sets with DOIs. REFERENCES