Sub-clustering in skeletal class III malocclusion phenotypes via principal component analysis in a southern European population

The main aim of this study was to generate an adequate sub-phenotypic clustering model of class III skeletal malocclusion in an adult population of southern European origin. The study design was conducted in two phases, a preliminary cross-sectional study and a subsequent discriminatory evaluation by main component and cluster analysis to identify differentiated skeletal sub-groups with differentiated phenotypic characteristics. Radiometric data from 699 adult patients of southern European origin were analyzed in 212 selected subjects affected by class III skeletal malocclusion. The varimax rotation was used with Kaiser normalization, to prevent variables with more explanatory capacity from affecting the rotation. A total of 21,624 radiographic measurements were obtained as part of the cluster model generation, using a total set of 55 skeletal variables for the subsequent analysis of the major component and cluster analyses. Ten main axes were generated representing 92.7% of the total variation. Three main components represented 58.5%, with particular sagittal and vertical variables acting as major descriptors. Post hoc phenotypic clustering retrieved six clusters: C1:9.9%, C2:18.9%, C3:33%, C4:3.77%, C5:16%, and C6:16%. In conclusion, phenotypic variation was found in the southern European skeletal class III population, demonstrating the existence of phenotypic variations between identified clusters in different ethnic groups.

Class III malocclusion has been commonly described as having a retruded or hypoplastic upper maxilla, a prognathic or hyperplastic mandible, or a combination of both 1 . These shared characteristics in skeletal class III malocclusion vaguely define the different sub-phenotypes clearly identified in recent studies [2][3][4][5][6][7][8] .
The current diagnosis of skeletal radiographic analysis includes multivariate techniques that allow the phenotypic characterization of different groups of class III skeletal malocclusions. Principal component analysis (PCA) enables the large number of variables used in cephalometric analysis to be reduced to fewer components, thus grouping the variables of greatest interaction on each axis 9 . This has made it possible to create phenotypic groups (clusters) for this specific malocclusion using cluster analysis to ascertain the least dispersion in each group and concomitantly achieve the greatest difference between each of the groups found 10 .
The number of distinguishable skeletal class III sub-phenotypes varies substantially between different studies and populations described in the literature: from seven different sub-phenotypes in Korean adults 2 to three 3,6,7 clinically distinguishable clusters in growing Caucasian patients; and four 8 , five 5 , and 14 4 skeletal class III clusters (C) have also been obtained via cluster analysis (CA) and described in the literature. These studies clearly show great variations within the same skeletal malocclusion and can be classified in different phenotypic groups with specific characteristics that differ from one ethnic group to another, and even within the same ethnic group [2][3][4][5][6][7][8]11 . The aforementioned studies have extensively described different numbers of class III malocclusion sub-phenotypes by taking into consideration skeletal and dental, dentoalveolar, and soft tissue measurements in different ethnicities and age groups 11 . The use of skeletal cephalometric measurements alone to perform these analyses may ensure that only skeletal clusters are identified, thus avoiding the interaction of dental components and soft tissue variations that could mask an accurate skeletal perception 12 . The main purpose of the present Study population, eligibility criteria, and subject inclusion process. The study population comprised 699 patients of southern European ethnic origin from Spain seeking orthodontic treatment. Eligible subjects were recruited from patients affected by skeletal class III malocclusion attending the Orthodontic Postgraduate Program at the Complutense University of Madrid and other private dental practices nearby. All the participants had been previously diagnosed with class III malocclusion by their doctors at each of these centers. All subjects were of southern European origin, had practically completed their growth as assessed by cervical vertebral maturation stage (CVMS) 14 and given informed consent to participate in the research.
Final selection for inclusion was made after compliance with all other specific criteria (Supplementary Table 1). Of the initial 699 participants, a final selection of 212 adult skeletal class III participants was made (97 men and 115 women; mean age for men was 28 ± 10.66 and 29 ± 10.12 for women) for cluster generation.
Sample size estimation was based on information available from previous studies 5 . The sample size was calculated with a confidence level (1 − α) of 95%, statistical power of 90%, precision (d) of 0.30, and variance (S2) of the quantitative variable of reference group of 0.69. This calculation determined a total sample of 131 individuals. The sample size was then adjusted for losses; the expected proportion of losses (R) was fixed at 15%, and a sample size of 154 subjects was finally estimated.
Clinical and diagnostic records. Pre-treatment lateral cephalometric radiographs of every patient were obtained prior to orthodontic treatment (Orthoceph by Siemens, Sirona Orthophos plus DS Ceph, Gendex Orthoralix). All patients were placed 10 cm away from the radiographic plate. All of the lateral radiographs were imported into Dolphin Imaging software (11.0, Dolphin Imaging & Management Solutions, Chatsworth, CA, USA) and calibrated and standardized with a 10 mm digital ruler. Chronological age was recorded and CVMS was determined 14 in all patients to verify that they had no remaining significant craniofacial growth. In addition, patient gender and ethnic origin were noted as part of the required criteria for inclusion in the present research.

Radiological landmarks, measurements, and axis model generation.
A single experienced operator (L.F.V.) measured the craniofacial landmarks and cephalometric measurements via the radiometric software used in this study to avoid additional potential errors. The examiner underwent a calibration period prior to initiation of this research using a subset of duplicated radiological landmark positioning and measurement assessment.
Overall, 102 cephalometric measurements and 67 radiological landmarks were located in each patient (see Supplementary Table 2 for full description). Of these landmarks, 29 were on pure skeletal structures, 12 were on dental positions, 15 were in soft tissue structures, and the remaining 11 were related to airway structures. A total of 21,624 radiographic measurements were obtained as part of the cluster model generation. Each patient was examined and described via 102 radiographic variables (Table 1): 55 were skeletal measurements, 25 were dental parameters, and 22 were soft tissue (15) and airway (7) parameters. Among the skeletal parameters, 23 were angular measurements, 24 were distances, and 8 were proportional measurements ( Table 1).
The skeletal variables underwent PCA to generate the best suited axis model to find the most significant components that might explain most of the variance found in our data set regarding skeletal craniofacial structures, reducing the dimensionality of the data and ranking them by weighted importance. PCA was conducted to reduce the number of initial variables used in the cephalometric analysis (n = 102) to a smaller set of axes correlated to each other. To assess the components, the axes were rotated so that each axis had large correlations with few variables. The varimax rotation was used and Kaiser normalization was applied to prevent variables with more explanatory capacity from affecting the rotation (IBM SPSS Statistics version 25.0).

Statistics
Reliability, accuracy, and method error. The same operator obtained 102 radiographic measurements of each patient. The reliability of the measurements was tested using a replica of the measurements of 15 patients randomly selected and separated in time by a 3-week interval. Specifically, the method error (ME) was calculated once all of the cephalometric traces were completed using the intraclass correlation coefficient, two-way mixed effects model for absolute agreement ICC 15 . The paired Student's t-test with a p value > 0.05 considered nonstatistically significant differences between the original and repeated measurements. The accuracy of the method was calculated using the Dahlberg formula in the replicated sample 16

Sub-phenotypes characterization and clusters configuration.
A priori data processing explained 92.7% of the variance using 10 main axes (Fig. 1). These enabled the application of cluster analysis and other multivariate techniques. The cluster analysis was performed (Coheris Analytics SPAD version 9.1) to identify and group subjects with similar phenotypes to make comparisons between them. A best suited model of 6 clusters was defined based on Ward's criterion, defining the distances between groups and establishing the least possible dispersion within groups to ensure the greatest homogeneity for each cluster (Fig. 2). To obtain a graphical representation, the generated skeletal cluster model was represented using the mean value of the main explanatory cephalometric variables of each component.  www.nature.com/scientificreports/ The potential interaction of CVMS 14 and sex on each of the generated clusters was also evaluated using the chi-square test. The supplementary cephalometric variables measured on the lateral cephalometric radiographs (dental, soft tissues, and airways) were recorded to observe how those variables produced each cluster.

Results
Reliability and accuracy of the method. Reliability tests of the cephalometric measurements determined a p value greater than 0.05 and ICC values higher than 0.8, with only one variable less than 0.85 (ANS-PNS/Me-Go (%)), indicating good reliability. The accuracy of the measurements retrieved an error value ranging from 0.02 mm for basal width and 0.2 mm for Wits (FOP) to a maximum value of 6.2 mm for anterior face height (N-Me) and 5.9 mm for length of the mandibular base (Go-Pg). For the angular radiographic measurements, we found error values of 0.02° for the facial taper and 0° for the inferior gonial angle (N-Go-Me) and maximum values of 1.5° for occlusal plane at SN and 1.2° for the articular angle. For the proportional measurements, we found errors of 0 to 3.1% for face height ratio (N-ANS/ANS-Me), ANS-PNS/Me-Go, and S-Ar/Ar-Go, respectively. Study sample distribution and types of axis models generated. A total of 699 lateral cephalometric radiographs were initially collected and revised. Of these, 487 had to be rejected due to lack of visible cephalometric points, poor quality, impacted teeth, retained teeth, dental agenesis, and patients with remaining growth or previous orthodontic treatment. Hence, a total of 212 subjects (97 men and 115 women) were included in the present research. Mean age was 28 ± 10.66 years for men and 29 ± 10.12 years for women. All the subjects were CVMS IV or V. Of the total sample, 45 subjects were CVMS IV and the remaining 167 were CVMS V. The total sample showed a mean ANB value of − 1.3 ± 1.5 and a Wits appraisal of − 5.2 ± 3.2. The MP-FH and y-axis variables, two of the most significant in the sample, had averages of 20° ± 5.9° and 66.25° ± 3.9°, respectively.
A 10-axis model was generated from the radiographic skeletal variables via PCA (Fig. 1). The results of PCA revealed that the first 10 principal components explained 92.7% of total variation. Specifically, the first three axes represented almost half of the variation. These axes were represented by sagittal and vertical measurements that  Table 2.
Skeletal class III sub-phenotypes: cluster characterization. The 10-axis model for radiometric variables led to a cluster configuration of a total of six skeletal class III malocclusion subphenotypes, or clearly distinguishable clusters (Fig. 3). The first skeletal cluster (C1) represented 9.9% (n = 21; mean age: 29.9 ± 10.39 years) of the sample. This skeletal class III malocclusion phenotype was characterized by an increased total anterior facial height, enlarged posterior cranial base and larger mandible size than the other groups (Fig. 3). The second skeletal cluster (C2) represented 18.9% of the sample (n = 40; mean age: 27.25 ± 8.47 years) and showed the vertical component. This group had an increased mandibular plane, as well as a lower proportion of posterior facial height relative to anterior facial height. Bi-maxillary retrusion was also found in this group. Both the upper and lower jaws were positioned in a more posterior location.
The third skeletal cluster (C3) was the most frequent phenotype in the sample, with 33% (n = 70; mean age: 28.49 ± 9.39 years) presenting characteristics intermediate between C2 and C4. C3 was described as a slight skeletal class III malocclusion. In this phenotype, reduced upper facial height was found and the smallest mandibular length of all groups. Midface distance and maxillary length were the second smallest of the six skeletal groups after C4. Slight mandibular projection was also shown. This cluster represented a boundary class III skeletal malocclusion (Fig. 3). www.nature.com/scientificreports/ The fourth skeletal cluster (C4) determined in this study represents the least frequent distribution of the skeletal class III cohort, 3.77% (n = 8; mean age: 28.64 ± 9.80 years) and exhibited a severe skeletal class III malocclusion with a reduced anterior and posterior cranial base. The variables for mandibular projection presented the most augmented values of all the skeletal clusters. In addition, upper maxilla hypoplasia was observed. This skeletal C4 had a more decreased posterior facial height than the other clusters and a severely reduced symphysis width (Fig. 3).
The fifth skeletal cluster (C5) (16%, n = 34), unlike C4 was characterized by increased maxilla size and increased length of the middle third of the face compared to the five other skeletal clusters. The mandibular structure had larger dimensions. And the variables characterizing the vertical dimension indicated a brachial phenotype, but of lesser degree than C6 (Fig. 3).
Finally, skeletal cluster 6 (C6) (16%, n = 34) represented a moderate skeletal class III malocclusion with a short face. Skeletal C6 had a considerably reduced anterior facial height (N-Me), a similarly decreased lower Table 2. Principal components generated based in a ten axes model. a variables with the highest contribution in each PC.   Table 3 shows the mean values of the skeletal radiographic variables of each of the groups. With respect to the distribution of gender and skeletal maturation stage in the generated clusters, the CVMS 14 was homogenous for all clusters (p > 0.05), while gender was asymmetrically distributed in four of the six clusters. C1, C3, C4, and C5 had differences in gender distribution (p < 0.05). C1 and C5 were mostly represented by male subjects, at 85.7% and 79.4%, respectively, while C3 and C4 were represented mainly by female subjects, at 74.3% and 100%, respectively. Similarly, in the male-dominated clusters, total linear measurements (Go-Gn, Go-Pg, NMe, SN, S-Ar, and ANS-PNS) were much higher than in the female-dominated clusters (C3 and C4). With

Discussion
In the present study, 6 different types of class III skeletal malocclusion phenotypes and 10 main components were obtained using grouping techniques. Both approaches, principal component generation and morphological class III cluster generation, form part of a sequential process in which the principal component axes are not real biological coordinates of the system but a means of grouping the variables for later construction of the class III clusters. Specifically, the analysis of the main components facilitated the generation of 10 axes, which were constructed exclusively from the 55 skeletal measurements of the subjects of the study. This allowed for extremely high coverage of variance (92.7% of the variance). We found differences in the size of the axis models generated compared with previously published studies 5,8,17 , which restricted the number of main components to five 17 and six 5,8 , respectively, and covered 67% 17 to 81.2% 5 of total variance in their sample. This may be because in the present study, only skeletal measurements were considered for the axis model generation.
Dental and soft tissue measurements were excluded to avoid potential interactions with the construction of the purely skeletal class III clusters, as these parameters could modify or mask the underlying craniofacial structure 12 . Hence, using a large number of critical measurements, 55 skeletal variables, we achieved a notable percentage of explained variation, which, although not excellent, shows minimal loss of information, which might be of critical importance for the configuration of skeletal clusters, although assuming, as with other diagnostic methods, that some uncovered features would be underrepresented.
As described in previously published studies 5,8,17 , the first three axes represented more than half [58.5%] of the total variation found in the sample, essentially referring to the description of sagittal and vertical parameters, while others 5,8,17 also represented the position of the lower incisor in the first half of the variation derived from the involvement of dental variables in the formation of PCs. Most of the parameters that represent sagittal and vertical measurements in our PC1 (26.6% of variation), PC2 (19.6% of variation), and PC3 (12.3% of variation) are equivalent to the parameters of PC1 (23.7% and 20.6%) and PC2 (17.3% and 19.34%) in other studies 5,8 , while PC1 in another study 17 is more like PC4 (9.3%) and PC5 (6.5%) (ANB, facial angle, and Pg-N perp) in the present research, despite the fact that the same sagittal and vertical parameters were included in these aforementioned components. The observed differences might be at least partially explained by the type of class III patients included and their distribution, essentially based on distinct ethnic origin and growing stage of the included subjects.
The 10 main models enabled the subsequent generation of skeletal class III malocclusion sub-phenotypes or clusters. According to the dendrogram analysis, potential construction with three or six skeletal clusters was Table 3. Craniometric mean values defined in each skeletal cluster. SD standard deviation.  www.nature.com/scientificreports/ possible (Fig. 2). Nevertheless, a six-cluster skeletal model was selected, since we observed substantial loss of representation of clearly detailed sub-phenotypes in the three-cluster model, such as C2 representing a vertical cluster with bimaxillary retrusion, or C4 representing maxillary hypoplasia with a prognathic mandible. The number of clusters found in previous studies that also used cluster analysis focusing on class III malocclusion classification varied between 4 8 and 7 clusters 2 , very close to the number of clusters established in the present research. Despite this, we also observed in the literature other studies with even more clusters than those found in our study, varying between 3 3,6,7 and 14 clusters 4 , which might provide high variation in coverage but completely lose efficiency in terms of practical utility in clinical settings. Other authors 3,4,6,7 used a diffuse cluster analysis when establishing a three-cluster model 6,7 , implying that the elements of analysis may belong to more than one cluster, while others used a hierarchical cluster analysis 3 instead of the mixed cluster analysis used in the present research. However, with respect to 14 different clusters 4 , the increased number of clusters could be derived, at least in part, from the inclusion of class III dental malocclusions as well as skeletal class III malocclusions, and the classification was also used to evaluate the effects of treatment 4 , which required large numbers in order to provide enough patients in each sub-group, while at the same time ensuring potentially higher homogeneity in intersubject characteristics in each of the constructed groups.
With respect to the skeletal clusters generated in the present research, the cluster representing the highest percentage of the sample was C3, at 33%. This cluster was characterized by a mesofacial pattern (y-axis: 66.7 ± 2.1, SN-GoGn: 29.0 ± 3.2, and NBa-PtGn: 90.1 ± 2.7) with slight maxillary retrusion (convexity: − 2.1 ± 1.4) and a slight anterior mandibular position (PgN perp: 5.4 ± 4.1) but with decreased total mandibular size (GoGn: 79.4 ± 3.9, GoPg: 69.8 ± 3.4, and CoGn: 114.4 ± 5.1) than the rest of the clusters. The characteristics of this cluster were similar to those found in C2 clusters in previous studies 5 , which was one of the most representative clusters in their sample. In the same research 5 and a previous study 17 , the group with the highest number of subjects was represented by a group characterized by a retrusive jaw (C5 and C3, respectively). In our study, the retrusive jaw model was the second most representative model in the sample (C2), representing 18.8% of the total sample. In another study, on the other hand, the model with maxillary deficiency alone, but without mandibular compromise, was not represented 8 . Interestingly, this must have derived from the ethnic origin of the sample population, since in studies that analyzed Caucasian subjects 5 and even a mixed sample with the highest percentage of subjects being of Caucasian origin 17 , this type of cluster was represented, as in the present research. In contrast, clusters generated from samples of Asiatic origin did not generate this type of sub-phenotype 8 , which clearly demonstrates the morphological differences between ethnic groups 11,18,19 and even between different populations of the same ethnic group when it comes to analyzing class III malocclusion sub-phenotypes.
In our study, the soft tissue, dental, and airway variables were excluded to avoid an inadequate interpretation of the generated groups, since it has been observed that soft tissues do not closely follow bone structures 12 . Despite this, several common characteristics related to the skeletal parts and soft tissues were observed.
In this respect, the soft tissue of the chin continued the underlying bone structure in the 6 sub-phenotypes found. Similarly, lower facial height in soft tissues (Sn'-Me') and hard tissues (ANS-Me) coincided in 4 (C1, C3, C4, and C6) of the 6 clusters. Conversely, skeletal facial height in C2 (69.89 ± 5.30) was 0.09 mm higher than in C5 (69.80 ± 3.45), but facial height in soft tissues (72.23 ± 5.06) was 0.73 mm lower than in C5 (72.96 ± 3.78). Despite this slight difference, our results agree with previous observations 12 , which concluded that the soft tissue of the lower facial third continues the skeletal profile, but not in the labial region 12 . In this study, as observed in previous ones [20][21][22] , a more significant relationship was found between the labial tissue and the position of the incisors.
The position of the lower lip with respect to the lower incisor was analyzed in several studies that found a close relationship between the two 23,24 , even observing a more significant relationship in the position of the lower incisor with respect to the lower lip compared to the relationship between the upper incisor and the upper lip 23 . This study found that in 5 of the 6 sub-groups, the position of the lower lip coincided with the position of the lower incisor. In the case of C6, the most retruded position of the lower incisor was observed with both the NB line (1.1 ± 2.3) and A-Pg (0.84 ± 2.62), coinciding with the most retruded position of the lower lip with respect to the S line (− 2.62 ± 2.22) and E line (− 5.29 ± 2.43). The relationship between the position of the lower incisor with the NB plane or A-Pg plane with respect to the position of the lower lip with the S line and/or E line was observed in C5, C4, C3, and C2, but not in C1. Individual characteristics may lead to different soft tissue behaviors with respect to the bones and/or teeth [23][24][25] , with different correlations between the two components. The method used for the analysis of soft tissues should also be taken into account, since recent studies support the need to conduct this analysis using facial photographs, since skeletal cephalometry and soft tissue cephalometry lead to different diagnoses of the patient's facial characteristics 23 .
Due to the considerable variations, caution must be taken when analyzing the soft profile, especially at the labial level. Parameters such as labial thickness or nasal size, which can vary widely in each individual 25 and do not follow a skeletal pattern 26 , cannot be ignored. These variations, not only individual, but also soft tissue variations between genders 25 , can lead to a distortion of the final labial position.
In previously published articles 3,[5][6][7][8]17 that conducted this type of classification of skeletal class III malocclusion, the samples contained similar gender percentages to our own sample: 40.8% were males and 58.6% were females. In other studies, the gender mix in the total sample ranged from 40.2% 8 to 48.7% 3 for males, and between 51.3% 3 and 59.7% 8 for females. Only one of the aforementioned studies 8 mentioned the gender percentage in each of the clusters found: two of the groups found in that study (C3 and C4) had substantially higher distributions of women than men. In our study, we observed another 3 clusters in which the difference in the female/male ratio was greater than 40% (C1: 14.3/85.7, C3: 74.3/25.7, and C5: 20.6/79.4). Although both clusters followed a similar trend in terms of gender distribution, a thorough analysis of the cluster characteristics found some differences in the configuration of skeletal structures derived from ethnicity of Asiatic and Caucasian origin 18,27,28 .

Scientific Reports
| (2020) 10:17882 | https://doi.org/10.1038/s41598-020-74488-w www.nature.com/scientificreports/ Beyond differences in ethnicity and gender, skeletal class III malocclusion is strongly linked to genetic components and environmental factors 29 . Current phenotype and genotype studies of this malocclusion [30][31][32][33] have been reported in the literature to develop preventive and therapeutic advances in these patients 34,35 . Recent studies in the field of medicine focus on the search for preventive measures and individualized treatments based on the patient's morphological and genetic components 34,[36][37][38] . However, this will require future research, larger samples, and collaborative multicenter studies 5 to better understand the differences between regions and environmental factors involved. To reach this point, and bearing in mind that the distribution and frequencies of some specific genetic variations can differ substantially between different populations 39 , it is imperative to obtain an accurate and precise classification in the sub-phenotype diagnosis of each population and/or ethnic group. The present research offers significant data on skeletal class III malocclusion sub-phenotyping in subjects of southern European origin. To the best of our knowledge, this is the first scientific evidence in the field of skeletal clusters of subjects of European origin affected by class III malocclusion.

Conclusions
In the present study, 10 axes were responsible for 92.7% of variability derived from 55 pure skeletal radiographic measurements. Similar to the results of other articles, different clusters represented phenotypic variability in class III skeletal malocclusion. The present research identified six clusters that clearly represented different subphenotypes of skeletal class III malocclusion in a population of southern European origin.