The effect of genetic polymorphisms on treatment duration following premolar extraction

To elucidate genetic factors affecting orthodontic treatment duration, we employed targeted next-generation sequencing on DNA from the saliva of 117 patients undergoing orthodontic treatment after premolar extraction. The clinical characteristics of patients are summarized, and the association of clinical variables with treatment duration was assessed. Patients whose treatment duration deviated from the average were classified into an extreme long group or an extreme short group. We identified nine single nucleotide polymorphisms (SNPs) of six genes that significantly differed in the two groups via targeted sequencing. The frequency of the CC genotypes of WNT3A, SPP1 (rs4754, rs9138), and TNFSF11, TT genotype of SPP1 (rs1126616), and GG genotype of SFRP2 was significantly higher in the extreme long group than in the short group. In the extreme short group, the TC genotype of SPP1, AA genotype of P2RX7, CT genotype of TNFSF11, and AG genotype of TNFRSF11A tended to exhibit higher frequency than in the long group. Taken together, we identified genetic polymorphisms related to treatment duration in Korean orthodontic patients undergoing premolar extraction. Our findings could lead to further studies predicting the prolongation of the orthodontic treatment duration, and will be of great aid to patients as well as orthodontists.


Results
Correlation between clinical parameters and orthodontic treatment duration. Among the clinical parameters of patients, horizontal anterior retraction was positively correlated (Beta = 0.1026, P-value = 0.0014) with treatment duration, while crowding and displacement exhibited a negative correlation (Beta = − 0.1241, P-value = 0.0029; Beta = − 0.1104, P-value = 0.0041, respectively). Age, sex, Angle classification, overjet, overbite, A point-nasion-B point (ANB), peer assessment rating (PAR) score, Frankfurt mandibular plane angle (FMA), and U1 to SN did not show any significant association. Multivariate linear regression was performed after adding all clinical variables as covariates and the analysis indicated horizontal anterior retraction as the variable most strongly correlated with orthodontic treatment duration ( Table 1).
Classification of extreme phenotype groups. The subjects were classified based on treatment duration. After creating a correlation plot between horizontal anterior retraction and orthodontic treatment duration, values outside the standard residual ± 1 were set as the extreme groups. Values above the standard residual (+ 1) were classified into the 'extreme long group' with a treatment duration longer than the average, and values below the standard residual (− 1) were classified into 'extreme short group' with a treatment duration shorter than the average. The number of patients in each group was 18. Clinical characteristics of the subjects. The clinical characteristics of the subjects are shown in Table 2.
There was no significant difference between groups with regard to age, gender, Angle classification, crowding, displacement, overjet, overbite, ANB, FMA, U1 to SN, as well as PAR score. As expected, treatment duration was significantly different between groups, as they had been classified based on orthodontic treatment duration ( Table 2).
Association between SNPs and orthodontic treatment duration. Targeted capture sequencing identified 142 variants with a minor allele frequency of > 3% in the candidate regions. No common variants were significantly associated with orthodontic treatment duration after Bonferroni correction for multiple testing (cut-off P-value = 0.05/142 = 3.52 × 10 -4 ). Only nine SNPs showed nominal significance (P < 0.05) in frequency between the extreme short and long groups (Table 3).
In the extreme long group, the frequency of the CC genotypes of WNT3A, SPP1 (rs4754, rs9138), and TNFSF11 (rs12585229, rs931273), the TT genotype of SPP1 (rs1126616), as well as the GG genotype of SFRP2, was significantly higher than that in the short group. In contrast, the frequency of the TC genotype of SPP1 (rs4754 and rs1126616), the AC genotype of SPP1 (rs9138), the AA genotype of P2RX7, and the CT genotype of TNFSF11 was higher in the extreme short group when compared with that in the extreme long group. www.nature.com/scientificreports/ Furthermore, the heterozygote of TNFRSF11A (AG genotype) showed overdominance in the extreme short group. Three SNPs of SPP1 (rs4754, rs1126616, rs9138) were linked together and rs9138, as an effective SNP, remained statistically significant after haplotype analysis was performed. The SNPs of TNFSF11 (rs12585229, rs931273) were also linked together. Rs752107 of WNT3A and rs9138 of SPP1 were identified as 3′-UTR SNPs. www.nature.com/scientificreports/ The rs3751143 SNP of P2RX7 was a missense variant, while the SNPs of TNFSF11 (rs12585229, rs931273) and TNFSF11A (rs4524034) were intron variants. The rs3810765 SNP of SFRP2 was located in a splicing region involved in transcription. Remaining SNPs were synonymous variants encoding the same amino acid sequences. In all subjects, genetic association with orthodontic treatment duration for significant SNPs of extreme phenotype sampling study was showed on Supplementary Table 1. The SNPs significantly related to the treatment duration are represented using scatter plots. The CC type of SPP1 (rs4754, rs9138) and the TT type of SPP1 (rs1126616) were associated with longer orthodontic treatment duration when compared to other types (Fig. 1). The AA type of TNFRSF11A (rs4524034) and the CC type of TNFSF11 (rs931273) were associated with an increase in treatment duration (Fig. 2). Other SNPs did not exhibit any significant difference in orthodontic treatment duration based on genotypes in all subjects (P > 0.05).
The SNPs discovered above and the clinical variables related to the orthodontic treatment duration in the univariate analysis were added as parameters for multivariate analysis to confirm the effect of clinical variables and genetic polymorphisms on the treatment duration. As a result, it was found that rs9138 of SPP1 and rs4524034 of TNFRSF11A affect the orthodontic treatment duration (Supplementary Table 2). For more details, genetic information of significant SNPs in extreme phenotype sampling study and allele frequency in all subjects were described in Supplementary Table 3.

Discussion
In the current work, we explored factors that influence orthodontic treatment duration, with a focus on genetic traits. Univariate linear regression was employed to analyze the clinical parameters of subjects, and horizontal anterior retraction exhibited a positive correlation with treatment duration. In other words, the greater the Table 3. SNPs significantly associated with orthodontic treatment duration (leveling and retraction period) in targeted sequencing study. CI Confidence interval, Beta Beta coefficient. P value was adjusted for gender and age at the beginning of treatment. Bold values denote statistical significance at the P < 0.05 level. *Nucleotide location number was assigned according to the WNT3A (Transcript ID:   www.nature.com/scientificreports/ horizontal retraction of anterior teeth, the longer the treatment duration. In contrast, crowding was negatively associated with orthodontic treatment duration. Crowding was previously reported to influence the increase in leveling duration 9 . However, after the completion of leveling, crowding is resolved, and the extraction space decreases, leaving less available space for anterior tooth retraction. According to simple linear regression data on the degree of crowding related to leveling and retraction period (data not shown), the decrease in anterior retraction duration was greater than the increase in leveling duration. Thus, simple linear regression suggested that crowding was associated with a shorter orthodontic treatment duration. The multivariate linear regression analysis performed thereafter indicated that horizontal anterior retraction was the only clinical variable significantly associated with orthodontic treatment duration. These results are consistent with the lack of consensus between previous studies on clinical factors affecting orthodontic treatment duration. Conflicting results have been reported regarding the effects of age, sex, and Angle classification on treatment duration 2,4,10,20 . There have also been discrepancies on the association of pretreatment ANB, FMA, PAR score, overjet, and overbite with treatment duration 7,11,21,22 . Although clinical parameters related to orthodontic treatment duration exhibited a low correlation in the current study, no conclusion can be drawn regarding the relationship of analyzed clinical parameters and treatment duration, as further clinical research is necessary.
Targeted next-generation sequencing revealed nine SNPs in six genes with a significantly different frequency between the extreme groups classified based on horizontal anterior retraction. The frequency of the C allele SNP of WNT3A (rs752107) was higher in the extreme long group compared with that in the extreme short group. Wnt signaling is crucial for the growth and maintenance of various tissues, including bone 23 . The WNT3A protein Rs752107, as a 3′-UTR SNP, may contribute to an increase in bone formation through its effect on microRNA (miRNA) binding affinity. That is, it alters the binding site of miRNA, which normally inhibits WNT3A translation by binding the WNT3A mRNA transcript harboring the C allele. In contrast, the T allele may compromise miRNA binding and enhance WNT3A protein translation. Higher WNT3A protein levels would be expected to negatively regulate osteoblast proliferation and differentiation, in turn decreasing bone formation 25,26 . Thus, patients with the CC genotype would be expected to exhibit increased bone formation and higher bone density, consistent with the association observed in the present study 27 . Increased bone density may subsequently slow orthodontic tooth movement, prolonging treatment duration 12 .
There are few studies on the association between WNT3A genetic polymorphisms and orthodontic treatment duration. The current findings suggest the possibility of rs752107 being related to bone remodeling, tooth movement, and thus orthodontic treatment duration. However, further studies are required to elucidate the molecular mechanism through which this SNP regulates bone formation as well as to determine whether rs752107 is in linkage disequilibrium with other functional polymorphisms.
SPP1 encodes OPN, which mediates alveolar bone resorption at the compression side when orthodontic force is loaded. In a recent study, OPN upregulation via ERK-and p38 MAPK-mediated signaling was reported at the tension side during orthodontic tooth movement 28 . OPN is associated with the shortening of orthodontic treatment duration via the acceleration of orthodontic tooth movement through the regulation of alveolar bone remodeling at both the compression and tension sides in response to orthodontic stress 19 .
Studies have previously reported the association between genetic polymorphisms of SPP1 and bone density. The C allele of SPP1 (rs4754) was suggested to be related to high bone density and low risk of fracture 29 , while the T allele of rs4754 was associated with low bone mineral density and a high risk of osteoporosis in women 30 . In the current study, the frequency of the CC genotype of SPP1 SNPs (rs4754, rs9138) and the TT genotype of rs1126616 was high in the extreme long group. The C allele of rs4754 is associated with reduced SPP1 function, resulting in decreased bone resorption and increased bone density. High bone density may slow the speed of orthodontic tooth movement and thus lengthen orthodontic treatment duration 12 .
The effects of other SNPs (rs1126616, rs9138) on bone biology and orthodontic tooth movement have not yet been studied, with the exception of a report on the significant relationship of rs9138 with external apical root resorption (EARR) in orthodontic treatment 31 . As a 3ʹ-UTR SNP, rs9138 can significantly affect miRNAmediated SPP1 degradation and its protein levels. Further research is required to assess the association of SPP1 polymorphisms and orthodontic treatment duration.
We observed a high frequency of the GG genotype of the SFRP2 SNP rs3810765 in the extreme long group. SFRP2 encodes secreted frizzled-related protein-2 (sFRP-2), one of the sFRP family proteins that bind directly to Wnt proteins and appear to antagonize their effects by sequestering them from their receptor. However, the role of sFRP-2 in teeth, bone, and supporting tissues remains elusive 32 . sFRP-2 was reported to enhance the inflammatory response and inhibit bone formation 33 ; however, another study showed that SFRP2 activates the osteogenic differentiation of apical papilla stem cells, increasing bone formation 34 . Although there is no research on the effect of rs3810765 in tooth and bone tissues, as a splice region SNP, rs381076 is expected to affect mRNA www.nature.com/scientificreports/ splicing, thereby altering RNA stability, and influencing SFRP2 gene expression. Further molecular studies are necessary to determine how the G allele would contribute to longer orthodontic treatment duration. The purinergic P2X7 receptor (P2RX7) is a member of the P2X subfamily of purinergic receptors activated by extracellular ATP. P2RX7 expression has been reported in both osteoclasts and osteoblasts where it regulates bone formation and resorption. Further, P2RX7 was demonstrated to be involved in mechanically-induced signaling between osteoblasts and osteoclasts 35 . Loss of P2RX7 function leads to increased osteoclast numbers and decreased trabecular bone mass 36 . In P2RX7 KO mice, alveolar bone formation was not significantly affected, but the susceptibility to external apical root resorption was increased with the application of orthodontic force 37 .
As a missense variant, P2RX7 SNP (rs3751143) is related to various diseases of bones, teeth, and periodontal tissues. The C allele of rs3751143 might contribute to low bone density and a high risk of osteoporosis in postmenopausal women 38 . In the current study, the A allele frequency of rs3751143 was high in the extreme short group, whereas that of the C allele was not significantly different in the extreme long group. Nonetheless, no conclusion can be drawn regarding the association between rs3751143 and orthodontic treatment duration, as reliable results from a greater number of subjects as well as functional studies are required.
In this study, the frequency of the C allele of TNFSF11 SNPs (rs12585229, rs931273) was significantly increased in the extreme long group, which might be associated with lower expression of the TNFSF11 gene. TNFSF11 encodes RANKL (receptor activator of NF-κB ligand), and TNFRSF11A encodes RANK (receptor activator of NF-κB). RANKL binds to its receptor RANK to induce the differentiation of osteoclasts through RANKL/RANK signaling 39 . During orthodontic tooth movement, the expression of TNFSF11 is increased in osteocytes of the periodontal ligament where the orthodontic force is applied, and bone resorption occurs due to osteoclast activation via enhanced RANKL levels 40 . Thus, RANKL promotes orthodontic tooth movement by inducing bone resorption in the direction of tooth movement. When a drug formulation that slowly releases RANKL was injected into the mesial and distal areas of rat premolars, followed by the application of orthodontic force, faster tooth movement was observed 41 .
As the SNPs of TNFSF11 (rs12585229, rs931273) and of TNFRSF11A (rs4524034) are intron variants, they are generally expected to have little effect on gene function. However, even an intron variant can affect gene expression at the transcription or translation level, as indicated by previous studies, wherein these SNPs were shown to be related to external apical root resorption 42 . While there have been some studies on genetic polymorphisms of TNFSF11 that affect bone metabolism, no research on TNFSF11 SNPs rs12585229 and rs931273 has been published. There is only a report relating the rs4524034 SNP of TNFRSF11A to menarche in women 43 . Nevertheless, the genetic polymorphisms of TNFSF11 and TNFRSF11A are worth investigating in depth, considering their important role in orthodontic tooth movement.
This study has several limitations, including a small sample size, setting of the orthodontic treatment duration that excluded finishing period, and retrospective nature of the study. There were also practical difficulties encountered in the process of selecting patients suitable for the set criteria. In addition, this study did not elucidate the functional mechanism of how the SNPs are related to the treatment duration. Future studies are required to elucidate these mechanisms, validate the effect of the SNPs on the orthodontic treatment duration, and determine their influence in other cohort groups.

Materials and methods
Definition of orthodontic treatment duration. In this study, the duration of orthodontic treatment was defined as the leveling and anterior retraction period for space closing, as these processes are the most relevant to the duration of orthodontic treatment in clinical practice.
Subjects. This retrospective study was approved by the institutional review board of Yonsei University Dental Hospital (number: 2-2016-0023). All subjects were Korean participants and unrelated. All clinical examinations were conducted in accordance with the Declaration of Helsinki, and written consent was obtained on paper with a description of the study prior to enrolling all subjects. Initially, 122 patients who received orthodontic treatment with maxillary premolar extractions at the Department of Orthodontics at Yonsei University Dental Hospital in Seoul, Korea from February 2008 to May 2020 were enrolled. Among them, 117 patients (24 males, 93 females, with an average age of 19.8 years) were finally selected; five patients were excluded owing to complications in the genetic analysis due to contamination, or because the patients could not proceed with the treatment due to personal reasons such as long-term absence, poor collaboration, and frequent breakage of brackets or appliances.
The selection criteria were as follows: (1) treatment with fixed orthodontic appliances in the permanent dentition, (2) maxillary premolar extracted on both sides and anterior teeth retracted for space closing, (3) panoramic and lateral cephalometric radiographs taken before and after orthodontic treatment, (4) no orthognathic surgery, (5) no systemic diseases affecting tooth movement, bone metabolism, or maxillofacial malformation.
All 117 patients were treated using the straight-wire appliance (SWA) technique, from a 0.016-in nickeltitanium to a 0.019 × 0.025-in stainless steel (G & H Orthodontics, Franklin, Ind.) with conventional brackets of 0.022 slots.

Measurement of clinical parameters related to treatment duration.
To determine the clinical parameters related to orthodontic treatment duration, 18 clinical variables of the subjects such as age, sex, Angle classification, horizontal anterior retraction, crowding, displacement, overjet, overbite, PAR score, ANB (T1: before treatment, T2: after treatment), FMA (T1, T2), and U1 to SN (T1, T2) were measured using pre-treatment dental models, lateral cephalometric radiographs, and chart review. Lateral cephalometric radiographs were acquired before and after orthodontic treatment by a trained radiographer using Cranex3 + (Soredex, Helsinki, www.nature.com/scientificreports/ Finland) at the Department of Oral and Maxillofacial Radiology, Yonsei University Dental Hospital. These lateral cephalometric radiographs were traced using V-Ceph™ 5.5 (Osstem, Seoul, South Korea). Displacement, measured in mm, was defined as the sum of the absolute values of crowding and spacing to quantify abnormal tooth position. PAR score was obtained from the PAR index to quantify the degree of initial malocclusion 44 . Horizontal anterior retraction was defined as the difference between the tip of maxillary anterior teeth parallel to the HRP (horizontal reference plane) before (T1) and after (T2) orthodontic treatment in the lateral cephalometric radiograph, as measured in mm. HRP, proposed by Burstone 45 , was the line minus 7° from the Sella-Nasion line, and it was taken as the horizontal reference of orthodontic tooth movement in the study (Fig. 3).
Reliability of the measurement method. Ten patients were randomly selected to check for errors in radiographic measurement, and twenty lateral cephalometric radiographs were acquired from them. Each radiograph was examined by the same operator after a 2-week interval to determine reliability of the measurement. The intraclass correlation coefficient (ICC) between the two examinations was over 0.9, indicating considerable consistency. The difference between the first and second assessments was not significant. Genetic analysis. Genomic DNA was isolated from saliva using the Oragene DNA self-collection kit (Genotek, Ottawa, Ontario, Canada) and extracted with the prepIT-L2P DNA extraction kit (DNA Genotek, Ottawa, Ontario, Canada) according to the manufacturer's recommendations. DNA quantification was performed using Qubit.
For targeted next-generation sequencing analysis, previously reported genes and their SNPs related to tooth movement or bone metabolism were screened after a literature review (Supplementary Table 4). Cytokine genes www.nature.com/scientificreports/ that induce the inflammatory response around periodontal ligament tissues were selected as genes related to tooth movement. Genes involved in Wnt signaling and the RANK/RANKL pathway were also selected owing to their influence on orthodontic tooth movement and regulation of alveolar bone formation or resorption. Genes included only the coding sequence (CDS) region. For SNPs outside the CDS region, additional target probes were designed. Hybridization capture-based next-generation sequencing was employed, and sequencing data were obtained using a genome analysis tool kit. Sequencing reads obtained from Illumina NextSeq 500 platforms were further analyzed using BWA-MEM, Picard (v1.115), and Samtools (v1.1), and GATK (v4.0.4.0) was used to call single nucleotide variants.
Statistical analysis. Statistical analyses were performed using R 3.5.2 (R Foundation) and GraphPad Prism software (GraphPad Inc., La Jolla, CA, USA). Univariate and multivariate linear regression were performed to analyze the correlation between parameters and orthodontic treatment duration. Differences in clinical parameters between extreme long group (showing long treatment duration) and extreme short group (showing short treatment duration) were analyzed using the χ 2 test or t-test, as appropriate. The dominant, codominant, and recessive model was used in genetic association analysis. Genotype frequencies of each SNP were tested using the Hardy-Weinberg equilibrium. All P-values were based on two-sided comparisons, and P < 0.05 was considered statistically significant. The degree of linkage disequilibrium (LD) of SPP1 and TNFSF11 SNPs was examined using Haploview 4.2 software.