Gene Dysregulation in Peripheral Blood of Young Females With Temporomandibular Joint Osteoarthritis

Early onset of the disease and female preponderance are the unique features of the temporomandibular joint (TMJ) osteoarthritis (OA). The immune modulation mechanisms related to etiology of OA from other joints such as knee or hip have been suggested, but the immune-associated pathophysiology of TMJ OA, especially in young females, has not been elucidated. The present study aimed to investigate the immune-related pathophysiology of TMJ OA by analyzing transcriptional proles of peripheral blood mononuclear cells which identify the differentially expressed genes (DEGs) in young females with TMJ OA. Methods


Introduction
Osteoarthritis (OA) is characterized by the degradation of components of the extracellular matrix within the articular cartilage and the simultaneous remodeling of the underlying subchondral bone, in association with low in ammatory changes [1]. The epidemiology of OA of the TMJ is different from OA of the hand, knee, or hip joint, the incidence of which is associated with aging. Previous reports have already described female preponderance and early onset of the disease, especially from the pubertal phases to the early 20 s [2,3]. The pathophysiology of the TMJ OA is multifactorial and complex, which includes diverse etiological factors such as prolonged parafunctional habits, abnormal occlusal relationship, sustained masticatory muscle tension, and hormonal imbalance [4,5]. Given the early onset of the condition, TMJ OA cannot be considered to be just a simple degenerative disease such as the OA of other joints associated with the aging process, so other etiological factors could be assumed.
Nevertheless, the clear pathophysiology of the TMJ OA, particularly in young patients has not yet been elucidated so far.
Many emerging evidences have demonstrated the role of immune modulation mechanisms in the development and progression of OA [6][7][8][9][10][11][12][13]. These processes involved immune-modulating agents, in both innate and adaptive compartment such as cytokines, chemokines, T cells, and B cells [6][7][8][9][10][11][12][13]. TMJ OA has been considered as a low in ammatory arthritic condition and mainly depends on in ammation and elevated levels of in ammatory mediators and cytokines in the TMJ synovial uid [14][15][16][17]. Even though one report suggested the possibilities of involvement of systemic immune dysfunction in occurrence of TMJ OA [18], the in uences of systemic immune function and composition of immune cells in peripheral blood on incidence and progression of TMJ OA have not been fully clari ed.
RNA sequencing (RNA-seq) technology has been utilized as a powerful tool to discover potential molecular mechanisms or therapeutic targets in a variety of diseases [19]. However, to the best of the knowledge, there are yet no reported studies which adopted RNA-seq technology to the patients with TMJ OA. Therefore, the molecular and genetic background of TMJ OA remains obscure and clinical treatment effects for TMJ OA are limited. Understanding comprehensive molecular pro ling of TMJ OA is an essential step in discovering new candidate target molecules that are potentially involved in the pathogenesis of TMJ OA. Hence, the aim of the present study was to investigate the role of the immunerelated pathophysiology of TMJ OA by analyzing the transcriptional pro les of RNA from peripheral blood mononuclear cells (PBMC) which identify the differentially expressed genes (DEGs) of young females with TMJ OA.

Participants
In the present study, a total of 35 young females (mean age 19.7 ± 3.1 years; age range 15-25 years) were enrolled. Twenty-four female patients (19.3 ± 3.1 years; age range 15-25 years) with TMJ OA on at least one side of the condyles were consequently recruited from those who attended the TMD·Orofacial Pain Clinic at a university hospital from January, 2019 to November, 2019 (RNAOA). Eleven young females (mean age 20.5 ± 3.7 years; age range 15-25 years) without any signs of the TMD and/or TMJ OA who voluntarily participated in the present study served as control (CON). Patients with the following conditions were excluded from the study: history of head and neck trauma prior to at least 6 months prior to study entry; autoimmune diseases; craniofacial anomalies; and neurodegenerative disorder. All participants in the RNAOA did not show any sign of capsule or myofascial pain in the orofacial area for at least 3 months before enrollment; hence the effect of pain condition on RNA transcription pro les could be excluded. Clinical parameters such as the degree of maximum unassisted opening as well as the duration of TMD symptoms including TMJ noise and di culties in opening and/or closing the mouth were evaluated. Participants were diagnosed following the Diagnostic Criteria for Temporomandibular Disorder Axis I [20,21]. All female participants in the RNAOA showed erosive osseous bony changes on computed tomography images indicating the area of discontinuation of the cortical lining and adjacent bone. The research protocol was approved by the Institutional Review Board of the University Hospital (AJIRB-MED-GEN-18-449). Informed consents were obtained from all participants.

RNA extraction from PBMC and RNA seq
To minimize the effect of circadian variation and menstruation cycle, peripheral blood was collected from all participants between 9:00 a.m. and 11:00 a.m. during their mid-luteal phases. For each individual, 3 ml peripheral blood was put into a PAXgene Blood RNA tube vacutainer tube (Qiagen, New York, USA). PBMC from the venous blood from all participants were isolated using Ficolle-Paque PLUS (Sigma-Aldrich). Trizol reagent (Thermo Fisher Scienti c Invitrogen Inc., MA, USA) was utilized to isolate the total RNA from PBMC. Genomic DNA contamination was removed using RNase-free DNase I. The amount of total RNA was measured using NanoDrop 2000 (Thermo Fisher Scienti c Invitrogen Inc., MA, USA) and the integrity and quality of RNA samples were analyzed using an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies, Inc., CA, USA). All samples passed quality control RNA integrity analysis (RIN ≥ 7). cDNA libraries were constructed using the TruSeq Stranded Total RNA Sample Preparation Kit (Illumina Inc, USA) according to the manufacturer's protocol. To produce 100 bp paired-end reads, the total RNA was sequenced using the Illumina NovaSeq 6000 system (Macrogen, Seoul, Korea). After sequencing, the indexed samples were demultiplexed before the generation of FASTQ les for analysis and assessed by FastQC version 0.11.7.
Using HISAT2 version 2.1.0 with the best score matches reported for each read, all libraries were aligned to hg19 assembly of human genome. The mapped reads were assembled using STRING Tie version 1.3.4d. Inter gene expression comparisons were based on calculated fragments per kilobase of transcript per million (FPKM) mapped reads.

Differential expression analysis of RNA-seq data
The expression level was normalized by calculating FPKM mapped reads. For DEG analysis, the value of log 2 (fold change) were calculated. Fold change values were calculated by dividing FPKMs from RNAOA by those from CON. The DEGs with an adjusted P ≤ 0.05 and log 2 (fold changes) ≤ -1 or ≥ 1 were determined. Using heatmap function, the hierarchical clustering of the expression pro les of detected DEGs was conducted.

Bioinformatics analysis
Using online tool STRING analysis (http://string-db.org, version 11.0), gene ontology (GO) pathway enrichment analysis was conducted to compare gene transcription patterns between RNAOA and CON and assess the functional association between encoded proteins. False discovery rate (FDR) adjusted P values were calculated for each enriched biological pathway and the threshold was set to an adjusted P < 0.01. Gene set enrichment analysis (GSEA, version 4.0.3) was used to compare gene dysregulation patterns in RNAOA with those in CON.
The threshold was set to an unadjusted P < 0.05 to fully explore the results by ingenuity pathway analysis (IPA; Ingenuity System Inc., Redwood City, CA, USA). Using IPA, biological processes, canonical pathways, and networks of analysis were analyzed. An enrichment score measures the overlap of observed and predicted gene sets. A z-score assesses the match of observed and predicted regulation patterns, serving as a predictor for the activation state of the identi ed regulator molecule.

Quantitative real time polymerase chain reaction validation
To validate the differential expression, gene expression levels were examined using quantitative real time polymerase chain reaction (qRT-PCR). qRT-PCR reactions were conducted via the Step ONE Plus (ABI, Life Technologies, CA, USA) using SYBR premix EX Taq II (Applied Biosystems, Foster City, CA, USA), and cDNA was synthesized using 1 µg of mRNA. The gene expression results were obtained using the formula 2^-(ΔCt), and the fold change was calculated by the formula 2^-(ΔΔCt). The qRT-PCR values were normalized using the average of the expression of the reference gene, GAPDH. Finally, the averaged fold ratios from the reference housekeeping gene were used as the relative mRNA level. Each experiment was conducted in triplicate. Two reactions, one without template and one without reverse transcriptase were also performed.

Statistical analysis
The differences in demographic features and TMD characteristics between CON and RNAOA were evaluated using Mann-Whitney U tests. A t-test identi ed signi cantly changed gene transcripts between RNAOA and CON. The gene expressions were de ned as fold changes. Statistical analysis was achieved by setting the change in gene expression threshold to an unadjusted P < 0.01 for encoded protein functional network analysis (FNA). The threshold was set to an unadjusted P < 0.05 to fully explore the results by IPA. To compare gene expression levels obtained by qRT-PCR between RNA OA and CON, Mann-Whitney U tests were utilized. The level of signi cance was set as P < 0.05.

Demographic features and TMD symptoms in participants
The differences in age (P = 0.268) and body mass index (P = 0.740) did not show statistical signi cance between CON and RNAOA. Participants in RNAOA showed 16.4 ± 28.0 months of TMD symptom durations and the amount of maximum unassisted opening showed a signi cant difference between CON and RNAOA (P = 0.036) ( Table 1). Gene dysregulation in peripheral blood of patients with TMJ OA There were 27685 transcripts in total, and transcripts with an FPKM value less than 5 were excluded, leaving 12788 transcripts to be analyzed. After excluding the following FPKM baseline gene, 440 genes were differentially expressed in the RNAOA compared with those in the CON (|fold change| > 1, P < 0.05). A total of 41 genes were expressed signi cantly differentially between RNAOA and CON with 17 genes upregulated and 24 down regulated in the RNAOA (log 2 (fold change) ≥ 1 or ≤ -1, P < 0.05) ( Table 2). A heatmap using hierarchical clustering analysis showing the expression levels of each of these genes per individual are provided (Fig. 1). The genes with the highest fold change in RNAOA include FGFR2, EREG, WTH3DI, CXCL8, and LINC02458 and those with the lowest fold change were SYNGR1, CYP27A1, IL11RA, MAVS, and MIR941-2 ( Table 2).
Pathway analysis using IPA IPA was carried out on the transcriptome dataset with signi cance set at P < 0.05 and |log 2 (fold change)| > 1. IPA identi ed 20 annotated categories of diseases and functions that were signi cantly upregulated (z-score > 2.0) and included a total of 13 of the DEGs from the study. The 5 top canonical pathways were associated with role of natural killer cells, communication between innate and adaptive immune cells, and T cell signaling (Table 3). Nine genes were involved in ve top canonical pathways namely, FCGR3B, HLA-C, HLA-F, HSPA-6, CXCL8, ITGB8, FCER1A, IL11RA, and IL13RA1. Hence, these pathways reinforce signi cant immune and in ammatory dysregulation occurring in the pathogenesis of TMJ OA. qRT-PCR The gene expression patterns of 6 hub genes which showed signi cant results in both GO enrichment analysis and IPA including HLA-C, HLA-F, CXCL8, IL11RA, IL13RA1, and FCGR3B transcripts were tested via qRT-PCR. There were signi cant differences in HLA-C (P = 0.030), CXCL8 (P = 0.022), and IL11RA transcripts (P = 0.018) (Fig. 4).

Discussion
Early onset of diseases, especially during pubertal phases and early 20 s is the unique feature of TMJ OA [2]. The progression of TMJ OA may accompany compromised masticatory function and altered craniofacial morphology, which could affect an individual's quality of life [22][23][24]. However, to the date, exact pathophysiology of TMJ OA particularly in young patients and clear molecular mechanisms of the development of TMJ OA have not yet been elucidated, so far. Immune dysfunction has been considered as one of the main etiological factors which may result in bony destruction in patients with OA from other joints [6-9, 12, 13], but studies which clearly elucidate the immunological background of TMJ OA especially in young patients are sparse. RNA-seq is an effective tool for clarifying molecular mechanisms in genetic levels and could provide the novel therapeutic targets involved in a certain condition [19]. To the best of the knowledge, no study ever attempted to reveal the immune related etiology of TMJ OA using RNA-seq technology. Consequently, the aim of the present study was to investigate the role of immune dysfunction by analyzing transcriptional pro les of PBMC and search for new therapeutic targets for the management of TMJ OA in young patients.
Declined acquired immune responses accompanied with increased auto-reactivity have been detected in the elderly [25] and relationships in this altered innate immune function, T cell and B cell responses, and cartilage and bone degradation have been reported [26,27]. Aforementioned results from IPA demonstrated the associations among various autoimmune disorders, including autoimmune thyroid disease and systemic lupus erythematosus (SLE), and occurrence of TMJ OA. Previous studies already mentioned the high prevalence of temporomandibular disorders in patients with autoimmune thyroid disorders or SLE [28][29][30]. Even though severe condylar resorption in patients with rheumatoid arthritis or juvenile idiopathic arthritis have been detected [31], few reports ever mentioned the bony changes of condyles in other autoimmune diseases such as autoimmune thyroid disorders or SLE. Because patients with autoimmune disorders were excluded from the present study, clear associations between autoimmune thyroid disorders or SLE and TMJ OA could not be clari ed. Even though, previous studies which dealt with the associations between innate immunity and OA progression focused on increased auto-reactivity related with aging process, the role of altered innate immune function in bony destruction of TMJ condyles even in young females could be assumed through this study.
Accumulating evidence suggested the involvement of in ammatory and immune responses in the pathogenesis of OA [32]. The results from the present study demonstrated the increased expression levels of chemokine(C-X-C motif) ligand 8 (CXCL8), the involvement of signal transducers and activators of transcription 3 (STAT3) pathway, and cytokine mediated immune response in the development of TMJ OA in young females. Interleukin-8/CXCL8 (IL-8/CXCL8) has been found to be an attractant for neutrophils and a population of lymphocytes [33]. The increased level of IL-8 in synovial uid from the knee OA and TMD patients with disc displacement have been reported [34][35][36][37], but the mechanisms of IL-8 in subchondral bone changes in TMJ have not been clearly revealed. One study has reported the elevated serum levels of CXCL8 and expression levels of STAT3 in knee OA patients [38]. This study suggested that CXCL8 may inactivate the mitosis of chondrocytes and further aggravate OA progression and indirectly induce the Janus kinase (JAK)/STAT3 signaling in chondrocytes. Although the speci c molecular mechanisms of CXCL8 and STAT3 could not be revealed in the present results, the associations between CXCL8, STAT3 signaling, and subchondral bony destruction of TMJ condyles in young females could be suspected.
The peripheral blood of OA patients have been analyzed and revealed that patients with OA have shown altered levels of CD8 + T cells and more cytotoxic pro les in comparison with healthy controls [11,[39][40][41]. One animal study suggested the molecular mechanisms of the role of CD8 + T cells on the OA process that CD8 + T cells were activated once OA had been initiated and cartilage degeneration occurred more slowly in CD8 + T cell knockout mice than in wild-type [42]. The results from this study also showed the relationships between the occurrence of TMJ OA and T cell activity through increased expression levels of genes related to Tec kinase signaling, nuclear factor of activated T cells (NFAT) regulation, and OX40 signaling. Furthermore, cytotoxic T lymphocyte mediated apoptosis and CTLA4 signaling pathway which have associations with cytotoxic T cell function, seem to have roles in the incidence of TMJ OA. A previous study has focused on the role of T helper cell in synovial uid on the subchondral bony changes in TMJ OA patients [16], but no study ever attempted to clarify altered cytotoxic T cell activities in blood in patients with TMJ OA.
To the best of our knowledge, the present study is the rst study which attempted to reveal the genetic and molecular background of TMJ OA in young patients using RNA-seq technology. Several previous reports which dealt with the systemic immune dysfunction and the development of OA have focused on the role of the aging process and the senescence of immune cells. In addition, most of the studies which analyzed the relationships between levels of the immune-modulating agents in blood and OA dealt with weight bearing axial joints such as the knee and hip not with peripheral joints. Aforementioned results showed that alteration systemic immune function could have a role in the development of arthritic conditions even in young patients with OA from small peripheral joints such as TMJ. However, the present study still has several limitations, rst of all, RNA-seq was conducted using only PBMC, not using the synovial tissue. Secondly, only female participants were included. Even though previous studies con rmed the female preponderance of TMJ OA, the absence of male participants would limit the understanding of the genetic etiology of the TMJ OA. Thirdly, owing to cross-sectional characteristics of the study, the information regarding the changes in the transcriptome pro les accordance with the progression or recovery of TMJ OA could not be provided. Future prospective RNA-seq studies with large samples including both male and female participants, using both PBMC and synovial tissue would be required.

Conclusion
In conclusion, in the present study, young female TMJ OA patients showed alterations in the blood immune cell expression and some of the changes may re ect in ammation, auto-reactivity, and altered T cell functions. For proper management and successful treatment of TMJ OA in young females, future research regarding immune-modulation based therapy would be warranted.

Consent for publication
Not applicable Availability of data and materials The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The author certify that no con icts of interests were involved in this paper.

Funding
This study was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (2018R1C1B6007671).