Integrated bioinformatics analysis of differentially expressed genes in the temporomandibular joint internal derangement

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 (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.


| INTRODUCTION
Temporomandibular joint Internal derangement (TMJD) is one kind of the most complex human diseases (Granados, 1979); it is a common disease affecting more than 15% of adults, peaking in incidence among those aged 20-40 years (Gauer & Semidey, 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., 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. Highthroughput 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., 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 (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., 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.

| Data acquisition
The GSE66864 gene expression series matrix was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Raw database GSE66864 was presented by the 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., 2015).

| DEG identification
Background calibration, normalization, and data filter were presented to DEGs through the R software (version 4.1.3) and GEO2R (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., 2020). WGCNA was utilized to confirm the reliability of DEGs in the different samples.

| 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., 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.

| Construction of PPI network
In our study, we established the PPI network that contained the hub genes via the STRING database (https://cn.string-db.org/) (version 11.5) (Szklarczyk et al., 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., 2003). The top 10 DEGs were selected as nodes based on the Degree algorithm (Chin et al., 2014).

| 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., 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 (https:// www.proteinatlas.org/) to identify similarities in human tissue.

| Drug interaction of hub genes
The DGIdb database (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., 2018).

| Identification of DEGs
One control group and four TMJD groups were contained in the GSE66864 database. We plotted a heatmap based on the expression of

| 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 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 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 2c). The KEGG pathway enrichment analysis shows a similar result of hub genes between human and rabbits (Figure 2d), which inspire us to associate rabbit genes obtained in this study with human genes.

| Drug interaction of hub genes
In our study, highly expressed hub genes were submitted to the DGIdb database (https://dgidb.org/) and eight potential drug interactions were obtained (Table 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.

| DISCUSSION
TMJD is a common joint disease often associated with TMJ synovitis (Ibi, 2019). Besides, TMJD also has a noninflammatory origin. Open surgical intervention is a traditional treatment for TMJD (Buckley  (Wu et al., 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, 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., 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., 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., 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., 2003). According to research conducted by Joel M. Kremer et al. (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., 2012). On the contrary, the low expression of ITGA6 can impair the activity of the PI3K/AKT pathway (Chen & Zhang, 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., 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., 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., 2020).
The PTHLH gene is related to the development of endochondral bone. Studies showed that PTHLH is highly expressed in oral cancer F I G U R E 4 (Continued). (Lv et al., 2014) and tongue cancer (Suwa et al., 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, 2012). Cell proliferation was reported to be highly associated with the upregulation of the PTHLH gene (Hameetman et al., 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., 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., 2020). Annexin A1 (ANXA1) is a ligand of FPR1 relevant for cancer immunosurveillance (Vacchelli et al., 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, 2018). The FPR1 expression can predict early cancer and function as a pharmacologic target for innate immune responses (Liu et al., 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. (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., 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., 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., 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., 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., 1995). Moreover, ATP6V0A2 may also relate to cancer cell ECM and control tumor metastasis (Katara et al., 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., 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.
T A B L E 1 Eight potential drug interactions between DEGs.

| 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.

CONFLICT OF INTEREST STATEMENT
The author declares no conflict of interest.

DATA AVAILABILITY STATEMENT
Data are openly available in a public repository that issues data sets with DOIs.