BMP2 rs1005464 is associated with mandibular condyle size variation

This study aimed to evaluate the association between single nucleotide polymorphisms (SNPs) in endochondral development-related genes and mandibular condyle shape, size, volume, and symmetry traits. Cone-beam Computed Tomographies and genomic DNA from 118 individuals were evaluated (age range: 15–66 years). Data from twelve 3D landmarks on mandibular condyles were submitted to morphometric analyses including Procrustes fit, principal component analysis, and estimation of centroid sizes and fluctuating asymmetry scores. Condylar volumes were additionally measured. Seven SNPs across BMP2, BMP4, RUNX2 and SMAD6 were genotyped. Linear models were fit to evaluate the effect of the SNPs on the mandibular condyles’ quantitative traits. Only the association between BMP2 rs1005464 and centroid size remained significant after adjusting to account for the false discovery rate due to multiple testing. Individuals carrying at least one A allele for this SNP showed larger condylar size than common homozygotes GG (β = 0.043; 95% CI: 0.014—0.071; P value = 0.028). The model including BMP2 rs1005464, age and sex of the participants explained 17% of the variation in condylar size. Shape, volume, and symmetry were not associated with the evaluated SNPs. These results suggest that BMP2 rs1005464 might be associated with variation in the mandibular condyles size.


Genotyping quality control
The overall genotyping success rate was 100%.Allele frequencies for all the SNPs were in Hardy-Weinberg equilibrium, and the minor alleles showed a frequency above 15% (Table 1).There was 100% agreement between the duplicate genotyping and the original calls.

Main findings
The sample consisted of 118 participants (36 males and 82 females) with a mean age of 32.1 ± 14. 4  Six principal components (PC) were identified for each aspect of shape variation, symmetric and asymmetric, explaining a cumulative variance of 84.2% and 75.5%, respectively (Fig. 1).PC wireframes are shown in Figs. 2  and 3. Supplementary Table S3 describes the different shape configurations of the mandibular condyles for each PC, ranging from those that represented more negative and positive scores.Data on PC scores, centroid sizes and Mahalanobis shape fluctuating asymmetry scores showed an adequate distribution.The analyses showed an individual effect on the condylar size, shape, and volume (P < 0.001).Also, a side effect (right/left) on shape (P < 0.001) but not on condylar size (P = 0.230) and volume (P = 0.461) was observed.
The complete results of the implemented statistical models for genetic evaluations are presented in Supplementary Table S4.The following phenotype-genotype relationships were significant at the nominal level (P < 0.05; F test P value < 0.05; Table 2): BMP2 rs1005464-symmetric PC2, symmetric centroid size, shape fluctuating asymmetry, and condylar volume; BMP4 rs17563-shape fluctuating asymmetry; and SMAD6 Table 1.Alleles and genotypes frequencies in the current sample.SNP -single nucleotide polymorphism, MAF-minor allele frequency, H-W-Hardy Weinberg.* (1 = minor allele / 2 = major allele).rs2119261-asymmetric PC4.The models for the effects of BMP2 rs1005464 on symmetric centroid size and condylar volume showed the most significant explanatory power (adjusted R 2 : ~ 15-20%).The studied SNPs did not significantly influence the symmetric PC1, PC3-PC6; asymmetric PC1-PC3, PC5 and PC6; asymmetric centroid size; and right-left condylar volume difference.The association between BMP2 rs1005464 and symmetric centroid size was the only one that remained significant even after Benjamini-Hochberg adjustment (P ≤ 0.05; F test P value < 0.001; power = 0.99; Table 2).Individuals carrying at least one A allele for this variant showed significantly larger condylar size than GG common homozygotes (β = 0.043; 95% CI 0.014-0.071;adjusted P value = 0.028; F test P value < 0.001).Linear models, including BMP2 rs1005464, age, and sex of participants as predictor variables, explained 17% of the variation in condylar size.

Discussion
Variations in the morphometric characteristics of the mandibular condyle might be involved in the development of skeletal malocclusions, occlusal disorders and joint dysfunction [23][24][25][26][27] .Little is known about the genetic basis of common craniofacial morphological variations.The complex background involved in the morphological determination of craniofacial phenotypes must be revealed.Understanding the role of genes in normal condyle development is the basis for exploring the molecular mechanisms involved in pathological conditions that affect the temporomandibular joint complex and craniofacial morphology.This is the first genetic association study that widely explored the phenotypic variability of this structure and evaluated its association with SNPs in endochondral development-related genes.Our findings suggest that BMP2 rs1005464 would be associated with the size of the mandibular condyles.www.nature.com/scientificreports/Morphometric analyses revealed wide variation in the 3D configuration of the mandibular condyles.These structures were observed to present diverse shapes with variation in the width, length, and height of the mandibular condyle (neck and head, or only head) and different spatial orientations of the condylar head (i.e., pitch, yaw, roll).Although complex, these data appear more informative than commonly used analyses on radiographs or bidimensional projections from cone beam computed tomographies (CBCT) 28 .No significant correlation was observed between the present study's shape PC scores, size, and volume measures.These findings suggest that each condylar phenotype evaluated would provide relevant, different, and complementary information.
The mandibular condyle is subject to significant shape and size changes during active growth 29 , and even during adulthood 30 .This is due to its adaptive capacity that allows anatomical modification to ensure morphological, functional, and occlusal stability 31 .It has been reported that a rounded shape is more common in young adults, while the flat or angled shape is more observed in older people 32,33 .Although with several other complex changes added, our findings also showed that with the increase in age, the condylar head would adopt a more flattened configuration (i.e., characteristics compatible with negative β for the PC5 symmetric component; Fig. 2).
In another aspect, the evidence is controversial regarding the influence of sex on condylar shape 30,[32][33][34][35] ; however, concerning to its size and volume, studies consistently show that men have larger condyles 21,[36][37][38][39] .Our findings were in line with this previous evidence.Although the differences in shape between men and women were not evident, greater size and condylar volume were observed in males.
Regarding side effects, previous studies have shown that shape, size, and volume present differences between the right and left condyles 36,37,39,40 .In the present sample, only a side effect on the mandibular condyle shape was detected.We hypothesize that the right-left shape differences could be due to masticatory function occurring mainly on a preferred side 41 .The absence of right-left differences for volume and size could be because there were no individuals with evident facial skeletal asymmetries (chin deviation greater than 4 mm) in the present sample.Since some condylar traits could bring symmetry information not detected by other measures and considering that it has been reported that diverse condylar phenotypes are important determinants of mandibular configuration and facial morphology 22,42,43 , our data supports the need to carry out a 3D shape, size, and volume analysis for a better characterization of these structures.
Genetic analyses unveiled a significant association between BMP2 rs1005464 and mandibular condyle size.BMPs are important growth factors that belong to the TGFβ superfamily with critical roles in skeletal development and chondrogenesis 44 .Deletion of BMP receptors in the embryonic or early postnatal periods results in non-development of secondary mandibular cartilages or attenuation of condylar cartilage extension, respectively [45][46][47] ; demonstrating the importance that BMP signaling would have in the formation and growth of this structure.Specifically, BMP2 is known to participate in the formation of condensations, chondrocytes proliferation and differentiation, and extracellular matrix synthesis on articular tissues 44,48,49 .This molecule is highly expressed in the condylar cartilage anlagen 50,51 .It has been shown that adding exogenous BMP2 in condylar explants of Runx2-deficient animals induces chondrogenic differentiation, suggesting that this molecule would be an important factor for secondary cartilage development 50 .Moreover, BMP2 is required for postnatal maintenance of condylar cartilage integrity.Deletion of Bmp2 in chondrocytes leads to early degeneration and decreased mineralization of the mandibular condyle, decreased cell proliferation, and increased expression of degenerative markers 10 .Based on this information, we assume that the association between rs1005464 and the mandibular condyle size is due to the influence of this variant on BMP2-mediated processes of formation, growth, and postnatal maintenance of condyle integrity.
rs1005464 has already shown a previous association with mandibular retrognathism and, in interaction with other SNPs in endochondral development genes, also contributed to presenting the dolichofacial pattern 19 ; both features related to a small mandibular condyle.The study showed that GG homozygotes were more likely to have mandibular retrognathism 19 .Our findings aligned with this result, showing that GG individuals had smaller mandibular condyle centroid sizes.Considering that BMP2 participates not only in the endochondral ossification process but also in the intramembranous one 52 , it cannot be answered, based on our findings, if this variant influences only the size of the condyle or also contributes to the size of the mandibular body and rami.Unfortunately, several CBCTs analyzed did not fully include these structures (i.e., mandibular body and rami), preventing this evaluation.It is important to mention that BMP2 rs1005464 has additionally been associated with phenotypes like tooth size 53 and dentoalveolar size discrepancies 54,55 , suggesting that it is likely that this variant has a pleiotropic effect influencing different size-related craniofacial features.
The present study has some limitations.Although the analyses were adjusted for the age and sex of the participants, and even though malocclusions did not influence the evaluated traits for the present study, the presence of selection bias due to other factors related to the study's participants cannot be ruled out, since a convenience strategy for sampling and recruitment was implemented.In addition, since we worked with patients' chart data, obtaining additional information to control for confounding factors was impossible.Considering that factors such as mandibular functionality, muscular activity, and static and dynamic occlusal contacts strongly influence the mandibular condyle 3,4,56 , these may have affected our results.Similarly, the participants' body height could also have affected the results of condyle size and volume [13][14][15][16] .
Regarding another aspect, the number of SNPs evaluated was limited.SNPs in other genes crucially involved in mandibular condyle development (e.g., SOX9, TGFβ, DLX5, SHOX2, FGFs, TWIST1, PTHrP and IHH) 5,6 need to be investigated in future studies.In addition, it should be emphasized that the absence of association of the other SNPs in BMP4, RUNX2 and SMAD6 does not mean that these genes are not involved in the configuration of the mandibular condyle; different variants in these genes could likely have a significant effect.On the other hand, some strengths should also be mentioned.Although the sample size is relatively small, the implemented models reached a power greater than 80% in most analyses.The reported BMP2 rs1005464 association remains significant after applying the Benjamini-Hochberg adjustment, decreasing the probability that this result is a type 1 error.Another strength is the method implemented to assess the mandibular condyles, extracting the greatest possible information about its phenotypic variability.

Methods
The present study followed an analytical observational cross-sectional design involving evaluating patient clinical records and subsequent analysis of phenotype-genotype relationships in eligible individuals.The Research Ethics Committee of the School of Dentistry of Ribeirao Preto, University of São Paulo, Brazil, approved and supervised the proper conduct of this study (n.3.150.551).Research was conducted after approval of the Institutional Ethics Committee and all experiments were performed in accordance with regulations of the latest version of Declaration of Helsinki guidelines and its amendments.All participants and / or legal guardians gave written informed consent to participate in the study.Recommendations for Strengthening the Reporting of Genetic Association Studies were followed for the report 57 .

Participants
Chart data of 403 patients from graduate and private dental clinics in two cities (Ribeirão Preto and Sorocaba), state of São Paulo, Brazil, with an indication for CBCT between 2008 and 2019, were screened to determine their eligibility.Ribeirão Preto and Sorocaba are in the southeastern region of Brazil.In this region, European ancestry predominates (60.7%), followed by African (32.0%) and Amerindian (7.3%) 58 .
CBCT from all individuals were taken for clinical purposes.Biologically unrelated individuals who had likely passed the peak of pubertal growth spurt of the jaws (≥ 15 years old) 59 , whose CBCT included both mandibular rami and condyles, were selected by a non-probabilistic convenience sampling.The exclusion criteria were diagnosed syndromes, dentofacial anomalies or metabolic diseases; history of facial trauma; degenerative temporomandibular joint disorders; parafunctional or chewing habits; previous orthopedic and / or ortho-surgical treatment; tooth loss affecting vertical dimension; and low-quality images.
After participant recruitment, 118 were considered eligible for the present study.Due to the relatively small number of participants, it was decided to analyze the total sample and retrospectively calculate the obtained power.

Phenotyping
Participants' CBCTs were evaluated for obtaining 3D landmark coordinates and condylar volume measurements.Preliminarily, CBCT files generated in Digital Imaging and Communications in Medicine format (.dcm files) were converted to Guys Imaging Processing Laboratory format (.gipl files) in ITK-SNAP 3.6.0software (http:// www.itksn ap.org).Images were then resampled to 0.25 mm isotropic voxel size in 3D SLICER 5.0.3 software (http:// www.slicer.org).Subsequently, the 3D analyses were performed following the step-by-step below: 1. Head orientation matching the midsagittal plane with Nasion, Crista Galli and Basion and the axial plane with the Frankfurt horizontal plane (Transforms module, 3D SLICER).2. Construction of 3D volumetric label maps (segmentations) of approximately the upper two-thirds of the mandibular rami, including mandibular condyles and coronoid processes (Active Contour and Paintbrush modes, ITK-SNAP).3. Identification and pre-labelling of 3D landmarks (Paintbrush mode, ITK-SNAP).Fourteen landmarks on the right and left condyles and surrounding structures were initially evaluated in the entire sample (Table 3).
Anterior Condylion from both sides was excluded due to limitations during their identification.Hence, 12 landmarks were finally selected and included in further analyses (Fig. 4).www.nature.com/scientificreports/ 7. Cropping of the mandibular condyles for volume measurement.3D surface models of the mandibular condyles were cropped at the level of an oriented plane, passing through the Sigmoid notch and the tip of the Coronoid process on the right and left sides (Easy Clip module, 3D SLICER) (Fig. 5).8. Reconversion of cropped 3D surface models to cropped volumetric label maps (Mesh to Label Map module, 3D SLICER).9. Measurement of the volume (mm 3 ) of the mandibular condyles using the cropped segmentations (ITK-SNAP).
The aforementioned procedures were carried out by a single evaluator (GAMV) previously trained and calibrated by a specialized researcher (ACOR) with expertise in the described method 60 .Twenty percent of the sample was randomly selected to evaluate the method error.The same evaluator identified the 12 selected landmarks, extracted 3D coordinates, and measured the condylar volume again in this subsample after at least two weeks.Intra-observer repeatability was evaluated by the intraclass correlation coefficient.Random and systematic errors were assessed by Dahlberg's Formula and by estimating the proportion bias according to the Bland-Altman method, respectively.
The 3D coordinates of the landmarks were imported into MORPHOJ 1.07a software (https:// morph ometr ics.uk/ Morph oJ_ page.html) to conduct geometric morphometric analyses.These assessments were performed under two symmetry perspectives 61,62 : (a) object symmetry for shape evaluations (i.e., considering the relative disposition of the condyles to each other since the position of each condyle is an integral aspect of mandibular symmetry), and (b) matching symmetry for condylar size evaluations (i.e., disregarding the relative disposition of the condyles to each other, analyzing similarities and differences between both structures).
The following geometric morphometric analyses were conducted: 1. Full Procrustes fit to extract the shape variation of the data set.The analysis under object symmetry generated two data matrices for the symmetric and asymmetric components of variation.In contrast, the analysis under matching symmetry generated a data matrix containing the Procrustes coordinates.2. Generation of covariance matrices and subsequent PC Analysis to identify the most critical aspects of variation in the data sets (i.e., those explaining at least 5% of the variation).Individual scores were generated for the identified PCs.Wireframes were created to represent the mean configuration and deformation of the structures of interest (Fig. 6).3. Condylar size estimation via centroid size.As mentioned above, the analysis was performed under the matching symmetry perspective, where a data set was generated based on the individual means of the right and left condyles.4. Estimation of the condylar size asymmetry via asymmetric centroid size. 5. Estimating individual scores of shape fluctuating asymmetry (i.e., deviations from the mean asymmetry) in units of Mahalanobis distances (dimensioned concerning the variation of the sample asymmetry, independently of the directional asymmetry).6. Procrustes ANOVA to assess individual and side (right / left) effects on the condylar size and 3D shape variation.

Genotyping
Saliva samples containing squamous epithelial cells from the buccal mucosa were collected to obtain genomic DNA.The participants were asked to rinse their mouth with 5 ml of 5% saline solution for 60 s and expectorate the volume into 15 ml propylene tubes (CORNING Inc., Corning, NY, USA).Next, sterile disposable cytobrushes (Plus GT, Medscand, COOPERSURGICAL Inc., Trumbull, CT, USA) were swiped and rolled twice to three times at the tongue's base and inner surface of the cheek on the right and left sides to collect additional biological material.The tubes containing the samples were centrifuged at 550 rpm for 10 min to sediment the cell pellet.The supernatant was discarded, and the pellet was resuspended in 1 ml of extraction buffer (Tris-HCl 10 mM, pH 7.8; EDTA 5 mM; SDS 0.5%).The samples were stored at − 20 °C until processing.Genomic DNA was extracted from buccal cells following a previously reported method 63 .The concentration and purity of the extracted DNA were evaluated using a spectrophotometer (NanoDrop 1000, THERMO SCIENTIFIC Inc., Waltham, MA, USA).Seven single nucleotide polymorphisms (SNPs) across four genes (Bone Morphogenetic Protein 2 [BMP2], Bone Morphogenetic Protein 4 [BMP4], Runt-Related Transcription Factor 2 [RUNX2], and SMAD Family Member 6 [SMAD6]) were chosen based on a minor allele frequency > 5% by searching public databases, previously reported association with craniofacial phenotypes, and their location on genes with biological implications in the growth, development and maintenance of the mandibular condyles (Table 4).Genotyping was blindly performed by polymerase chain reactions (PCR) using endpoint analysis 64 and TaqMan technology on a real-time PCR system (Step One Plus Real-Time PCR System, APPLIED BIOSYSTEMS, Foster City, CA, USA).Primers, probes, and universal master mix were provided by APPLIED BIOSYSTEMS (Foster City, CA, USA).Ten percent of the sample was genotyped in duplicate to test the quality of the process.

Statistical analysis
Descriptive statistics were used to report participant characteristics, and alleles and genotypes distributions for each SNP.Hardy-Weinberg equilibrium was evaluated by applying the chi-square test (or chi-square with Yates's correction, when necessary).The distribution and normality of the phenotype data were assessed using histograms, the Shapiro-Wilk test, and skewness values.ANOVA evaluated individual and side (right/left) effects on condylar volume measures.Univariate linear regressions were used to evaluate the single effect of the variables sex, age, and malocclusion on the mandibular condylar traits (Omnibus ANOVA test).

Figure 1 .
Figure 1.Scree plots show the variance explained by symmetric (A) and asymmetric (B) PCs.

Figure 2 .
Figure 2. Mandibular condyles shape configurations of subjects with the most negative (− β) or positive (+ β) individual scores for symmetric PCs.S-superior, I-inferior, A-anterior, P-posterior, R-right, L-left.Light blue lines represent the average configuration, while dark blue lines represent the variation of interest.

Figure 3 .
Figure 3. Mandibular condyles shape configurations of subjects with the most negative (− β) or positive (+ β) individual scores for asymmetric PCs.S-superior, I-inferior, A-anterior, P-posterior, R-right, L-left.Light blue lines represent the average configuration, while dark blue lines represent the variation of interest.

Table 3 .
Definition of 3D landmarks.*Paired landmarks located on the right and left sides.† Landmarks used for delimitation of the mandibular condyles for volume measurement.most superior point from the line connecting the lateral and medial poles of the mandibular condyle Posterior Condilion (PCo)* The most posterior point from the line connecting the lateral and medial poles of the mandibular condyle Anterior Condilion (ACo)* The most anterior point from the line connecting the lateral and medial poles of the mandibular condyle Lateral pole (LP)* The most lateral point of the mandibular condyle Medial pole (MP)* The most medial point of the mandibular condyle Regions close to the mandibular condyle Sigmoid notch (SN)* † Deepest point of the mandibular Sigmoid notch Coronoid process (CP)* † The most superior point of the Coronoid process of the mandible

Figure 6 .
Figure 6.Wireframe of the mean configuration of 3D landmarks on the mandibular condyles.A (black letter)-X versus Y axis (side view); B-X versus Z axis (top view); C-Y versus Z axis (anterior view); Ssuperior; I-inferior; A (blue letter)-anterior; P-posterior; R-right; L-left; 1 and 3-right and left Sigmoid notch, respectively; 2 and 4-right and left Coronoid Process, respectively; 5 and 9-right and left superior Condylion, respectively; 6 and 10-right and left posterior Condylion, respectively; 7 and 11-right and left lateral pole, respectively; 8 and 12-right and left medial pole, respectively.