Transcriptional Expression in Human Periodontal Ligament Cells Subjected to Orthodontic Force: An RNA-Sequencing Study

This study was performed to investigate the changes in gene expression in periodontal ligament (PDL) cells following mechanical stimulus through RNA sequencing. In this study, premolars extracted for orthodontic treatment were used. To stimulate the PDL cells, an orthodontic force of 100× g was applied to the premolar (experimental group; n = 11), whereas the tooth on the other side was left untreated (control group; n = 11). After the PDL cells were isolated from the extracted teeth, gene set enrichment analysis (GSEA), differentially expressed gene (DEG) analysis, and real-time PCR were performed to compare the two groups. GSEA demonstrated that gene sets related to the cell cycle pathway were upregulated in PDL. Thirteen upregulated and twenty downregulated genes were found through DEG analysis. Real-time PCR results confirmed that five upregulated genes (CC2D1B, CPNE3, OPHN1, TANGO2, and UAP-1) and six downregulated genes (MYOM2, PPM1F, PCDP1, ATP2A1, GPR171, and RP1-34H18.1-1) were consistent with RNA sequencing results. We suggest that, from among these eleven genes, two upregulated genes, CPNE3 and OPHN1, and one downregulated gene, PPM1F, play an important role in PDL regeneration in humans when orthodontic force is applied.


Introduction
Periodontal ligament (PDL), which is a group of connective tissue fibers, connects the tooth root to the adjacent alveolar bone. It protects blood vessels and nerves by absorbing mechanical force such as mastication force and provides proprioception. Moreover, it mediates orthodontic tooth movement under compressional or tensional force and plays a critical role in recovery from periodontal disease [1]. Mesenchymal stem cells were isolated from PDL [2] and have shown potential for periodontal regeneration [3]. PDL regeneration can also be enhanced under mechanical stimulus [4][5][6]. Low-magnitude, high frequency mechanical vibration promotes human PDL stem cells differentiation [7], and mechanical shear stress promotes the osteogenic differentiation of dental stem cells, including those derived from the pulp and the PDL [8]. An in vivo study using rats showed

Patient Selection
Subjects free from underlying disease were selected from orthodontic patients who had visited Yonsei University Dental Hospital between June 2017 and May 2018 for orthodontic extraction of at least two premolars. The inclusion criteria were as follows: a sound premolar in one quadrant without caries or restorations, and no history of taking any steroidal or non-steroidal anti-inflammatory drugs during force application. The exclusion criteria were root malformation such as dilaceration, which might result in difficult extraction. The minimum sample size of 10 subjects was estimated to overcome the individual variation in gene expression in PDL cells on the basis of a previous study using RNA-seq to analyze differences in gene expression between periodontitis-affected and healthy sides [19]. The study protocol was approved by the institutional review board of Yonsei University Dental Hospital (IRB No. 2-2017-0028) and performed after obtaining written informed consent from every participant.
Based on the inclusion/exclusion criteria, 14 patients were initially enrolled. Samples from three patients were excluded at initial quality check, and the samples from 11 subjects (5 men and 6 women; mean age, 22.4 years; age range, 17.5-31.0 years) were used in this study (Table A1). Two premolars were extracted in seven subjects; three premolars (two premolars in the experimental group and one premolar in the control group) were extracted in one subject (subject #4); and four premolars (two premolars in the experimental group and two premolars in the control group) were extracted in 3 of 16 three subjects (subjects #7, #8, and #9). In cases of 3-or 4-premolar extraction, the two teeth on the same side were assigned to the same group.

Mechanical Stimulus (Applying Orthodontic Force)
To stimulate the PDL, orthodontic force was applied to the first or second premolar on one side (the experimental group; n = 11), whereas the tooth on the other side was left untreated (the control group; n = 11). The study was a split-mouth, randomized, controlled trial with a 1:1 allocation ratio. A 0.016-inch nickel-titanium wire (Tomy International, Yokohama, Japan), which is known to deliver light continuous force (100× g) with a deflection of 0.5-1.8 mm, was used to apply mechanical stimulus in the experimental group ( Figure 1). The orthodontic force was applied for three weeks [20], and thereafter, the teeth in both groups were carefully extracted. PDL cells were washed with cold PBS after gentle scrapping from the middle third of the root and prepared for RNA-seq.

Mechanical Stimulus (Applying Orthodontic Force)
To stimulate the PDL, orthodontic force was applied to the first or second premolar on one side (the experimental group; n = 11), whereas the tooth on the other side was left untreated (the control group; n = 11). The study was a split-mouth, randomized, controlled trial with a 1:1 allocation ratio. A 0.016-inch nickel-titanium wire (Tomy International, Yokohama, Japan), which is known to deliver light continuous force (100× g) with a deflection of 0.5-1.8 mm, was used to apply mechanical stimulus in the experimental group ( Figure 1). The orthodontic force was applied for three weeks [20], and thereafter, the teeth in both groups were carefully extracted. PDL cells were washed with cold PBS after gentle scrapping from the middle third of the root and prepared for RNA-seq.

RNA Sequencing
To evaluate changes in gene expression after mechanical stimulus, total RNA was extracted from PDL cells using TRIzol Reagent (Invitrogen, Waltham, MA, USA). The RNA was stored at −70 °C and measured at an optical density of 260 nm. The mixtures of total RNA were incubated with Oligo dT (Gibco BRL, Rockville, NY, USA).
The library was constructed, sequenced using an Illumina HiSeq2500 sequencer (Illumina, CA, USA), and the data obtained through RNA-seq. A gene set underlying mechanical stimulus was analyzed using gene set enrichment analysis (GSEA). Differentially expressed gene (DEG) analysis was also performed to determine differences in gene expression between the experimental and control groups. For GSEA, libraries from 11 subjects in each experimental or control group were merged into one library for each group, and then, the experimental and control groups in total, were compared. Based on the GSEA results, an enrichment map was obtained. DEG analyses were performed to compare the libraries between the experimental and control groups in each sample as well as in the merged samples. Based on the DEG results, we selected target genes based on gene expression fold-changes: more than 1.5 fold-change in gene expression in the experimental group compared to the control group was defined as upregulated, and less than −1.5 fold-change was defined as downregulated [21][22][23].

Real-Time PCR
Real-time PCR was performed for the selected upregulated and downregulated genes to verify RNA-seq results. A cDNA synthesis reaction was performed using a mixture of AccuPower PCR PreMix (Bioneer, Daejeon, Korea) and water. The cDNA was synthesized from total RNA obtained from both groups using the SuperScript First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA) and was amplified with an ABI-7300 (Applied Biosystems, Mortlake, Waltham, Massachusetts, USA). The amplified cDNA was detected with SYBR Green PCR Master Mix Reagent Kit (Takara, Seoul, Korea). PCR conditions were as follows: incubation for 10 min at 95 °C, followed by 40 cycles of 10 s denaturation at 95 °C, and annealing for 60 s at 60 °C. The reaction mixture lacking cDNA was used as a negative control in each run. Primer sequences are summarized in Table 1. Ratios of the intensities of the target genes and GAPDH signals were used as a relative measure of the expression Control Experimental

RNA Sequencing
To evaluate changes in gene expression after mechanical stimulus, total RNA was extracted from PDL cells using TRIzol Reagent (Invitrogen, Waltham, MA, USA). The RNA was stored at −70 • C and measured at an optical density of 260 nm. The mixtures of total RNA were incubated with Oligo dT (Gibco BRL, Rockville, NY, USA).
The library was constructed, sequenced using an Illumina HiSeq2500 sequencer (Illumina, CA, USA), and the data obtained through RNA-seq. A gene set underlying mechanical stimulus was analyzed using gene set enrichment analysis (GSEA). Differentially expressed gene (DEG) analysis was also performed to determine differences in gene expression between the experimental and control groups. For GSEA, libraries from 11 subjects in each experimental or control group were merged into one library for each group, and then, the experimental and control groups in total, were compared. Based on the GSEA results, an enrichment map was obtained. DEG analyses were performed to compare the libraries between the experimental and control groups in each sample as well as in the merged samples. Based on the DEG results, we selected target genes based on gene expression fold-changes: more than 1.5 fold-change in gene expression in the experimental group compared to the control group was defined as upregulated, and less than −1.5 fold-change was defined as downregulated [21][22][23].

Real-Time PCR
Real-time PCR was performed for the selected upregulated and downregulated genes to verify RNA-seq results. A cDNA synthesis reaction was performed using a mixture of AccuPower PCR PreMix (Bioneer, Daejeon, Korea) and water. The cDNA was synthesized from total RNA obtained from both groups using the SuperScript First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA) and was amplified with an ABI-7300 (Applied Biosystems, Mortlake, Waltham, Massachusetts, USA).
The amplified cDNA was detected with SYBR Green PCR Master Mix Reagent Kit (Takara, Seoul, Korea). PCR conditions were as follows: incubation for 10 min at 95 • C, followed by 40 cycles of 10 s denaturation at 95 • C, and annealing for 60 s at 60 • C. The reaction mixture lacking cDNA was used as a negative control in each run. Primer sequences are summarized in Table 1. Ratios of the intensities of the target genes and GAPDH signals were used as a relative measure of the expression level of the target genes. To ensure accuracy of the experiments, primer specificity was confirmed by the dissociation curve after PCR, and real-time PCR assays were performed in triplicate for each sample. The mean fold-change in expression in the experimental group compared with the control group was calculated from the Ct values, and the range of the fold-changes was represented by standard deviations of the values [21]. Table 1. Primer sequences used for real-time quantitative PCR.

Genes
Transcript ID Forward Reverse

Statistical Analysis
For differential expression analysis, gene level count data were generated using HTSeq-count v0.5.4p3 tool [22] with the option "-m intersection-nonempty" and -r option considering paired-end sequence. Based on the calculated read count data, DEGs were identified using the R package called TCC [23]. TCC package applies robust normalization strategies to compare tag count data. Normalization factors were calculated using the iterative DEGES/edgeR method. Q-value was calculated based on the p-value using the p.adjust function of R package with default parameter settings. Differentially expressed genes were identified based on the q value threshold less than 0.05.

Extracted Premolar Samples
Two premolars were extracted in seven subjects; three premolars (two premolars in the experimental group and one premolar in the control group) were extracted in one subject (subject #4); and four premolars (two premolars in the experimental group and two premolars in the control group) were extracted in three subjects (subjects #7, #8, and #9). In cases of 3-or 4-premolar extraction, the two teeth on the same side were assigned to the same group.

Gene Set Enrichment Analysis of Mechanically Stimulated PDL Cells
GSEA enrichment maps show enriched gene sets, most of which are related to the cell cycle pathway ( Figure 2). Additionally, pathways for DNA replication, immune system, and metabolism were represented among upregulated genes. The node size, which indicates the number of enriched genes, was large in the following gene sets: cell cycle, cell cycle mitotic, DNA replication, antigen processing ubiquitination proteasome degradation, and class I MHC-mediated antigen processing presentation.
GSEA enrichment maps show enriched gene sets, most of which are related to the cell cycle pathway (Figure 2). Additionally, pathways for DNA replication, immune system, and metabolism were represented among upregulated genes. The node size, which indicates the number of enriched genes, was large in the following gene sets: cell cycle, cell cycle mitotic, DNA replication, antigen processing ubiquitination proteasome degradation, and class I MHC-mediated antigen processing presentation.

DEG Analysis of Mechanically Stimulated PDL Cells
A heatmap was obtained by comparing the merged libraries of 11 subjects between the experimental and control groups ( Figure 3). The 11 heatmaps obtained from each subject are presented in the Appendix Materials ( Figures A1-A11). Fifty-nine genes were selected by q value (q < 0.05), then 13 up-regulated and 20 down-regulated genes were finally selected based on a fold-change of 1.5 (Tables 2 and 3).

Discussion
This prospective study investigated the transcriptional expression in human PDL cells stimulated with an orthodontic force with RNA-seq. Because we intended to identify gene transcripts potentially related to the remodeling of PDL, this study focused on the differential expression of PDL genes using GSEA, DEG analysis, and real-time PCR. From these results, we found 11 significant DEGs when orthodontic force was applied on the PDL.
In the early phase of orthodontic tooth movement, mechanically stimulated PDL cells recruit local factors related to homeostasis and remodeling [16,24] of the surrounding periodontal tissue, which may result from the upregulated cell cycle pathways. Under tension force, PDL cells participate in DNA synthesis of osteoblasts; are differentiated to osteoblasts; or are released from G2 block, which contributes to the alveolar bone formation. Moreover, morphological deformation of the PDL induced by mechanical force plays a key role in recovery of its original shape via the remodeling process [25]. Therefore, the upregulated cell cycle of the PDLs would make orthodontic

Discussion
This prospective study investigated the transcriptional expression in human PDL cells stimulated with an orthodontic force with RNA-seq. Because we intended to identify gene transcripts potentially related to the remodeling of PDL, this study focused on the differential expression of PDL genes using GSEA, DEG analysis, and real-time PCR. From these results, we found 11 significant DEGs when orthodontic force was applied on the PDL.
In the early phase of orthodontic tooth movement, mechanically stimulated PDL cells recruit local factors related to homeostasis and remodeling [16,24] of the surrounding periodontal tissue, which may result from the upregulated cell cycle pathways. Under tension force, PDL cells participate in DNA synthesis of osteoblasts; are differentiated to osteoblasts; or are released from G 2 block, which contributes to the alveolar bone formation. Moreover, morphological deformation of the PDL induced by mechanical force plays a key role in recovery of its original shape via the remodeling process [25]. Therefore, the upregulated cell cycle of the PDLs would make orthodontic tooth movement possible by facilitating the remodeling process.
Among the 11 genes, CPNE3, OPHN1, and PPM1F are specially interesting, which were likely related to PDL remodeling based on the function and the signaling pathway of each gene. CPNE3 is a calcium-dependent membrane-binding protein that plays a major role in the synthesis of phospholipids and the immune system. CPNE3 belongs to the copine family, which can bind both calcium ions and phospholipids simultaneously. In particular, CPNE7 is involved in PDL regeneration [26]. The gene has a functional role in attachment of PDL cells to cementum and promotes the physiological arrangement of PDL fibers. Therefore, CPNE3, in the same gene family as CPNE7, may be involved in the regeneration of PDL. Further studies are needed to investigate whether and how CPNE3 affects PDL regeneration.
OPHN1 encodes the protein oligophrenin-1 that has a Rho-GAP domain shown to negatively regulate RhoA, Rac, and Cdc42 in vitro and in nonneuronal cells [27]. Rho GTPases participate in important cell biological processes, including cell growth control, cell motility, and development [28]. As expression of OPHN1 increases, the rhoA, Rac, and Cdc42 proteins, which are known to activate caspase 3 and finally activate apoptosis, decrease, which leads to reduced apoptosis of PDL cells. Therefore, upregulation of OPHN1 may indirectly contribute to the maintenance of PDL cells.
PPM1F is a member of the serine/threonine protein phosphatase family. This gene is known to be involved in caspase-dependent apoptosis [29]. When orthodontic force is applied, apoptosis in PDL cells may be inhibited by downregulation of PPM1F. It can be assumed that downregulation of PPM1F has a positive effect on PDL regeneration.
The upregulated and downregulated genes identified in the DEG analysis were confirmed by qPCR, although the qPCR ratio of each sample was different. Furthermore, there was a great variety of gene expression between subjects. It was assumed that limitations of the clinical study using human samples might cause the variety. As seen in the demographic features of the subjects (Table A1), different age, sex, tooth allocation, and occlusion of each subject would result in a wide variety of gene expressions unlikely in the animal models. In addition, the orthodontic force in the experimental group and the occlusal force in the control group would not be exactly the same in each tooth, which may result in differences in gene expression [4]. Moreover, the nonlinear properties of the PDL under different loading directions [30] would affect the outcome. The unstable nature of RNA and different elapsed times from tooth extraction to RNA extraction from the PDL may be another reason for the difference because proliferation of PDL fibers changes depending on post-extraction time [31].
The small numbers of DEG, 13 upregulated and 20 downregulated genes, were selected by comparing the merged libraries of 11 subjects between the experimental and control groups in the present study. Although the GSEA analysis demonstrated not only the cell cycle pathways but also the pathways related to metabolism, immune system, disease, and programmed cell death, the present study focused only on the genes related to the cell cycle. There are other factors that limited the number of DEG: high cutoff value, short duration of mechanical stimulation, and non-differentiation between compression and tension sides of the root surface, which might prevent the genes involved in PDL regeneration from expressing differentially. Moreover, the 33 genes could not be fully verified by qPCR in this study. Although 22 genes were discarded because they were not related to the cell cycle, some of the genes might play a role in PDL regeneration indirectly.
Mechanical stimulus contributes to PDL regeneration in vitro [4,7] and prevention of ankylosis after transplantation or replantation in in vivo animal models [9]. To the best of our knowledge, this is the first report of gene expression changes in the PDL following mechanical stimulus, in humans. We confirmed the upregulation of the cell cycle pathways after mechanical stimulus and found novel genes related to the process. Future studies are needed to identify the biological role of these genes and to obtain clinical relevancies such as saving the tooth using mechanical stimulus in trauma or transplantation cases.