Potential interactions among single nucleotide polymorphisms in bone‐ and cartilage‐related genes in skeletal malocclusions

© 2020 The Authors. Orthodontics & Craniofacial Research published by John Wiley & Sons Ltd 1Department of Orthodontics, University Medical Centre of Regensburg, Regensburg, Germany 2Department of Clinic and Surgery, School of Dentistry, Federal University of Alfenas, Alfenas, Brazil 3Department of Dentistry, School of Dentistry, Univille (Joinville University), Joinville, Brazil 4Department of Stomatology, Federal University of Paraná, Curitiba, Brazil 5Department of Pediatric Dentistry, School of Dentistry of Ribeirão Preto, USP Universidade de São Paulo, Ribeirão Preto, Brazil 6Department of Oral Diagnosis, School of Dentistry of Piracicaba, University of Campinas (UNICAMP), Campinas, Brazil


| INTRODUC TI ON
Skeletal malocclusions are complex craniofacial growth and developmental problems. 1 They are a set of human craniofacial morphologic characteristics that either exceed or exhibit deficiency of maxillary and mandibular dimensions, resulting in an improper relationship of the jaws that distorts the balance of the face. 2 Evidence gained especially from family and twin studies has demonstrated that genetic factors are strongly involved in the aetiology of skeletal malocclusions 3,4 Genes encoding proteins involved in bone and cartilage biology and skeletogenesis are candidate for skeletal malocclusions. The Bone Morphogenetic Protein (BMP) family is the largest subfamily of the structurally conserved transforming growth factor-beta (TGF-β) superfamily. BMPs are multi-functional growth factors that regulate the development, proliferation and differentiation 5,6 of mature osteoprogenitor cells into osteoblasts. 6 A review including several studies has shown the involvement of BMP2 in bone formation, 5 and BMP4 is involved in cell differentiation during skeletogenesis. 7 BMPs are also among the key pathways regulating craniofacial development and facial patterning. They regulate postnatal craniofacial growth and are associated with dental structures. 5 BMPs are involved during cartilaginous development, as are SMADs. 7,8 SMADs are important signalling pathway proteins that regulate the transcription of TGF-β superfamily genes. SMAD6 inhibits BMP signalling by interacting with transcription repressors. 9 RUNX2 (Runt-related transcription factor) is a key transcription factor associated with osteoblast differentiation and is considered a master regulator of skeletogenesis. RUNX2 is essential for the differentiation of pluripotent mesenchymal cells into osteoblasts, but also acts in mature osteoblasts maintaining the expression of bone matrix protein genes. 10 RUNX2 is expressed in different craniofacial tissues such as cartilage during the proliferation and maturation of chondrocytes. 10,11 Wnt signalling, one of the key cascades regulating development, is involved in important aspects of organogenesis such as control of cell polarity, fate and migration, 12 during cartilage and bone development, repair and regeneration. 13 The expression of Wnt signalling pathway genes during craniofacial development has been extensively investigated. 14 Genetic factors are involved in the aetiology of skeletal malocclusions and in vertical, sagittal, and transverse interrelationships of the dental arches. 2 Single-nucleotide polymorphisms (SNPs) are the most frequent variations in the human genome. SNPs in many genes were associated with different skeletal malocclusion phenotypes in different populations. [15][16][17][18][19][20][21][22] Therefore, in this study we investigated SNPs in bone-and cartilage-related genes in the aetiology of sagittal and vertical skeletal malocclusions.

| Sample
The Human Ethics Committee of the (identifying information) approved this study. Informed consent was obtained from all patients/children and/or their parents/legal guardians (in the case of minors).

Following the Strengthening the Reporting of Genetic
Association study (STREGA) statement checklist, 23 we evaluated genomic DNA (gDNA) extracted from saliva samples and pre-treatment lateral cephalograms from self-reported Caucasians as previously described. 15 All the included patients were enrolled in the orthodontic treatment at the (identifying information).
Biologically unrelated patients with no underlying syndromes nor congenital alterations and those without previous orthodontic and/or orthopaedic treatments were consecutively included in this study from 2015 to 2017. All patients that met the inclusion criteria were invited to participate in the study. Among 146 assessed individuals, one oral cleft patient and two patients, whose siblings were already included in the study, were excluded, yielding a total of 143 included patients.

| Phenotypes definition
Pre-orthodontic lateral cephalograms with the mandible in centric relationship were used, and digital cephalometric tracings performed by a calibrated orthodontist using the software Dolphin Imaging version 8.0 (Dolphin Imaging, Chatsworth, CA, USA). 16 Steiner's ANB, SNA, SNB angles and Ricketts' NBa-PtGn angle were measured to determine the sagittal skeletal jaw relationship  Supplementary Table S1. The genotyping was blindly performed using the Taqman ™ method for real-time PCR (ABI PRISM ® 7900HT Sequence Detection System, Foster City, CA, USA). Additionally, 10% of the sample was genotyped twice and an agreement of 100% was observed. The reaction was previously described. 15

| Statistical analysis
Chi-squared test was used to estimate the Hardy-Weinberg equilibrium and to compare genotype distribution among groups. Binary logistic regression analysis adjusted by gender and age was also performed. The established alpha for the exploratory results was P < .05. We also used the formal threshold for statistical significance after Bonferroni correction for multiple testing P < .0055 (0.05/9 SNPs).
Multifactor dimensionality reduction (MDR), a model-free and non-parametric method, was used to identify SNP-SNP interactions 24 using gender and age (before and after the growth spurt) as co-variables. MDR analysed all possible SNP combinations. A 10-fold cross-validation (CV) was performed, which calculated the ratio for each combination, separating 'high' or 'low' risk genotypes for each phenotype. The 1000 permutation test adjusted and determined statistical significance of the analysis. Models with 9/10 or 10/10 CV consistency, the testing balancing accuracy (TBA) >0.55 and P ≤ .05 were considered best models. Entropy values were obtained by the Jakulin and Bratko (2003) 25 formula, and MDR created dendrograms and interaction graphs using these values.

| RE SULTS
The mean age was 15.2 years (standard deviation = 7.3), 69 males and 74 females. The sample distribution according to the phenotypes is presented in Table 1.
All genotype distributions are demonstrated in Table 2 For the SNPs and phenotypes with suggestive association (P < .05), a logistic regression analysis was also performed using age and gender as co-variates (Table 3).  To explore high-order SNP-SNP interactions for each phenotype, we performed MDR analyses (Supplementary Table S2). Table 4 summarizes the MDR analysis and demonstrates the best MDR-predicted interaction models for the phenotypes that present significant models. Entropy measures among SNPs were calculated to obtain epistatic effects. Figure 1 shows the interactions between SNPs (dendrogram and interaction map) for phenotypes with interaction models. with the cranial base. 3,4 In the past decades, many genetic studies have been evaluating the association between different genes and maxillary/mandibular discrepancies as well as face morphology [15][16][17][18][19][20][21][22] These molecular genetic studies have mainly focused on skeletal class III and prognathism phenotypes. 17 More recently, some studies expanded the craniofacial phenotypes evaluated, including other sagittal and vertical patterns 16,[18][19][20][21][22] Therefore, in the present study, we decided to evaluate sagittal and vertical craniofacial phenotypes and each dental arch separately, in order to evaluate whether any gene/SNP acts in the maxilla, mandible or both jaws discrepancies.

| D ISCUSS I ON
Here, we also decided to perform MDR analysis to evaluate SNP-SNP interactions. This approach has proven to be a powerful tool for a variety of medical genetic studies. [29][30][31] In our study, MDR analysis allowed us to identify some interactions potentially involved in different craniofacial phenotypes.
Although non-syndromic mandibular retrognathism is a relatively common type of malocclusion, which refers to an abnormal posterior position of the mandible as a result of a developmental alteration, few studies explored the genetic aetiology of this condition. 32,33 Muscles are known to have extensive mutual effects on bones, and associated genes are candidate genes for skeletal mal- which is a condition that includes mandibular retrognathism (micrognathism). 35 Thus, it is plausible to assume that both SNPs in BMP2 studied here, the intronic and the missense (Arginine > Serine) variants, are involved in non-syndromic mandibular retrognathism.
The MDR analysis also suggested an interaction between rs235768 (BMP2) and rs1200425 (RUNX2). Interestingly, Runx2 was identified to be an important mediator of Bmp2 expression during cranial bone development 36. Furthermore, an association between rs59983488 in RUNX2 and maxillary protrusion was suggested in the genotypic distribution. Mandibular retrognathism and maxillary protrusion are traits of the skeletal class II phenotype.
Skeletal class II and skeletal class III are both anteroposterior discrepancies between the maxilla and mandible. In our study, the rs708111 in WNT3A was suggested as a protective factor for skeletal class III phenotype, which was previously associated with the palatal rugae pattern. 37 The rs708111 in WNT3A was also involved with skeletal class II, when SNP-SNP interaction was analysed via MDR analysis. The Wnt signalling pathways interact in an extensive network during bone formation regulated by a variety of molecules. 38 Our MDR results of skeletal class II phenotype reflect this complex interaction. A recent study identified SNPs associated with skeletal class II and skeletal class III phenotypes, two of these contained binding sites of RUNX2; however, the association was observed for skeletal class III. 18 Although we were not able to observe SNP-SNP interactions for skeletal class III, it is important to mention that the sample size of this group could be a limitation to identify such interactions.
In fact, the sample size is an important limitation of our study, which used a convenience sample to explore the genetic background of maxillary and mandibular discrepancies. Although this convenience sample allowed us to perform an exploratory study, the association of some SNPs with uncommon phenotypes may not be observed due to sample size limitations. This is particularly true in low penetrance SNPs. After performing a Bonferroni correction, many SNP associations became statistically insignificant. Although a correction for multiple variables reduces the chance of a type I error, it also increases type II error in a small sample and SNPs with small effects. For these reasons, our results should be interpreted with caution, but warrant and should prompt future investigations.
However, the 1.000 permutation test performed in MDR analysis is also an approach to adjust multiple tests, like Bonferroni correction, estimating type error I and power at 0.05 significance level. 39 Another limitation that should be highlighted is the fact that the population stratification correction was not performed to analyse the genetic association of our self-reported Caucasian population. In admixed populations, this could lead to associations with SNPs unlinked to the condition. Therefore, independent replication studies in different populations should be performed.
The twin-method study performed by Šidlauskas et al 2016, 3 suggested that the shape and sagittal position of the dental arches are under stronger genetic control. However, heritability is also involved in vertical morphology of the face3 , . 4 In our study, the vertical morphology of the face was also evaluated here and some interesting interactions were suggested for both dolichofacial and brachyfacial phenotypes. MDR analysis is a data approach that aims to identify multi-locus combinations of genotypes that are associated with either high-risk or low-risk combinations. Therefore, it is also possible that the same SNPs/ genes may be involved in both dolichofacial and brachyfacial phenotypes, however, with different risk genotypes. This was observed in the MDR analysis, which elected the same SNPs in WNT11 and WNT3A and also in SMAD6 for both dolichofacial and brachyfacial phenotypes.
SMAD6 is essential to regulate BMPs during cartilage development. 8 BMP signalling is complex, and there are multiple potential cross-talks, including Smad signalling 8 and Wnt signalling, 38 different BMPs either enhance or antagonize Wnt-induced osteogenic differentiation. 40 The RUNX2 rs1200425 was also included in the MDR model for the dolichofacial phenotype. Haploinsufficiency of RUNX2 causes cleidocranial dysplasia in humans, which is characterized by vertical morphology alteration and heterozygous Runx2-deficient mice present a similar phenotype, suggesting that the expression level of Runx2 influences skeletal facial phenotype. 41 A more detailed understanding of the cross-talk between the signalling important for postnatal craniofacial growth and development will help to elucidate the SNP-SNP interactions involved in the facial phenotypes, and further studies are necessary in order to investigate whether these SNPs are involved in the aetiology of sagittal and vertical skeletal malocclusions in humans.

| CON CLUS ION
Our results suggest that SNPs in BMP2, BMP4, SMAD6, RUNX2, WNT3A and WNT11 genes and their interaction could be involved in the aetiology of both sagittal and vertical skeletal malocclusions.

ACK N OWLED G EM ENTS
Open access funding enabled and organized by Projekt DEAL.
[Correction added on 10 November 2020, after first online publication: Projekt Deal funding statement has been added.]

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

AUTH O R CO NTR I B UTI O N S
ECK, PNF, PP and CK involved in conceptualization and designed the study. ECK, PP, RDC, RS and CK involved in funding support. MANM collected the sample and determined the orthodontic phenotypes. F I G U R E 1 Entropy Analysis (EA). Dendrograms and interaction maps of the EA for each outcome statistically significant in MDR analysis. The keys farther to the right in the dendrogram represent a strong interaction between the SNPs, and the keys farther to the left represent a weak interaction. In the interaction map, the values inside the box are the percentage of entropy of each SNP individually, and the values connecting the boxes indicate percentage entropy resulting from the combination. When the SNP-SNP entropy effect is bigger than the sum of values of each SNP, this indicates synergy, represented by the red and orange connections. Green and blue connections on the other hand indicate a redundancy effect, when the SNP-SNP entropy effect is smaller than the sum of values of each SNP. A, EA for skeletal class II. B, EA for mandibular retrusion. C, EA for dolichofacial phenotype. D, EA for brachyfacial phenotype MANM and PNF performed the cephalometric analysis. ECK, AOP, JC and RDC performed the laboratorial analysis. CLBR and RS performed the statistical analysis. ECK, CLBR and CK wrote the manuscript. All authors corrected and approved the final version of the manuscript.