GWAS of Post-Orthodontic Aggressive External Apical Root Resorption Identified Multiple Putative Loci at X-Y Chromosomes

Personalized dental medicine requires from precise and customized genomic diagnostic. To conduct an association analysis over multiple putative loci and genes located at chromosomes 2, 4, 8, 12, 18, X, and Y, potentially implicated in an extreme type of external apical root resorption secondary to orthodontic forces (aEARR). A genome-wide association study of aEARR was conducted with 480 patients [ratio~1:3 case/control]. Genomic DNA was extracted and analyzed using the high-throughput Axiom platform with the GeneTitan® MC Instrument. Up to 14,377 single nucleotide polymorphisms (SNPs) were selected at candidate regions and clinical/diagnostic data were recorded. A descriptive analysis of the data along with a backward conditional binary logistic regression was used to calculate odds ratios, with 95% confidence intervals [p < 0.05]. To select the best SNP candidates, a logistic regression model was fitted assuming a log-additive genetic model using R software [p < 0.0001]. In this sample the top lead genetic variants associated with aEARR were two novel putative genes located in the X chromosome, specifically, STAG 2 gene, rs151184635 and RP1-30E17.2 gene, rs55839915. These variants were found to be associated with an increased risk of aEARR, particularly restricted to men [OR: 6.09; 95%CI: 2.6–14.23 and OR: 6.86; 95%CI: 2.65–17.81, respectively]. Marginal associations were found at previously studied variants such as SSP1: rs11730582 [OR: 0.54; 95%CI: 0.34–0.86; p = 0.008], P2RX7: rs1718119 [OR: 0.6; 95%CI: 0.36–1.01; p = 0.047], and TNFRSF11A: rs8086340 [OR: 0.6; 95%CI: 0.38–0.95; p = 0.024]), found solely in females. Multiple putative genetic variants located at chromosomes X and Y are potentially implicated in an extreme phenotype of aEARR. A gender-linked association was noted.


Introduction
External apical root resorption (EARR) resulting from orthodontic forces represents one of the most undesirable iatrogenic effects secondary to mechanical strain during orthodontic movement, provoking an irreversible loss of root structure and tooth attachment in the apical third of the tooth root [1,2]. EARR is mostly manifested in its mild to moderate forms [3][4][5]; however, the most aggressive phenotype, with a frequency <1-5% and >5 mm apical loss, might critically compromise tooth viability [6,7]. EARR of any degree represents a complex pathological trait with multilevel etiological-risk factors that have not been completely elucidated to date [8]. Several diagnostic and clinical factors have been associated with EARR [9]. Several factors, such as treatment time, apical displacement, and gender-specific risk have been found to be associated with EARR, but these findings systematically show some degree of inconsistency and controversy in literature [10][11][12]. In fact, EARR occurrence and severity remains unpredictable and is not fully explained by clinical evidence alone. In this context, a genetic component and its contribution to this pathological feature, has been a critical issue that has recently garnered considerable attention [13,14]. Particularly, risk loci in somatic chromosomes have been targeted extensively, whereas very few studies are available regarding regions potentially associated with EARR at sexual chromosomes [15]. To date, a limited number of candidate gene association studies [15][16][17][18][19][20][21][22][23][24] have provided evidence suggesting that some genetic variants might exert a positive or negative influence over EARR of a mild to moderate degree at the level of the IL1 gene cluster [15,[20][21][22][23][24], TNFRSF11B [25], P2RX7 [17,18], SSP1 [19,23], or TNFRSF11A [15,26] at an autosomic level but also at IRAK1 gene [27] on sexual chromosome X. Despite the above mentioned studies, there is no available scientific evidence regarding how genetic factors might be specifically associated with the most severe phenotype of EARR, i.e., aggressive EARR (aEARR). Moreover, whether aEARR is positively/negatively associated with previous genetic variants and the studied clinical/diagnostic factors remains to be elucidated.
Therefore, the present study aimed to perform the first genome-wide association study conducting an association analysis over multiple putative loci and genes located at somatic chromosomes 2,4,8,12,18, and at sexual chromosomes X and Y, potentially implicated in aEARR.

Study Design
We performed a genome-wide association study of aEARR, a derived extreme and well-delimited root resorption phenotype, in UCM 3D g consortium participants ( Figure 1). The UCM 3D g consortium database is based on a general-population cohort of roughly 0.01 million patients aged 9-67 years old, recruited across Spain between 2005 and 2019.
J. Pers. Med. 2020, 10, x FOR PEER REVIEW 2 of 17 aggressive phenotype, with a frequency <1-5% and >5 mm apical loss, might critically compromise tooth viability [6,7]. EARR of any degree represents a complex pathological trait with multilevel etiological-risk factors that have not been completely elucidated to date [8]. Several diagnostic and clinical factors have been associated with EARR [9]. Several factors, such as treatment time, apical displacement, and gender-specific risk have been found to be associated with EARR, but these findings systematically show some degree of inconsistency and controversy in literature [10][11][12]. In fact, EARR occurrence and severity remains unpredictable and is not fully explained by clinical evidence alone. In this context, a genetic component and its contribution to this pathological feature, has been a critical issue that has recently garnered considerable attention [13,14]. Particularly, risk loci in somatic chromosomes have been targeted extensively, whereas very few studies are available regarding regions potentially associated with EARR at sexual chromosomes [15]. To date, a limited number of candidate gene association studies [15][16][17][18][19][20][21][22][23][24] have provided evidence suggesting that some genetic variants might exert a positive or negative influence over EARR of a mild to moderate degree at the level of the IL1 gene cluster [15,[20][21][22][23][24], TNFRSF11B [25], P2RX7 [17,18], SSP1 [19,23], or TNFRSF11A [15,26] at an autosomic level but also at IRAK1 gene [27] on sexual chromosome X. Despite the above mentioned studies, there is no available scientific evidence regarding how genetic factors might be specifically associated with the most severe phenotype of EARR, i.e., aggressive EARR (aEARR). Moreover, whether aEARR is positively/negatively associated with previous genetic variants and the studied clinical/diagnostic factors remains to be elucidated.
Therefore, the present study aimed to perform the first genome-wide association study conducting an association analysis over multiple putative loci and genes located at somatic chromosomes 2,4,8,12,18, and at sexual chromosomes X and Y, potentially implicated in aEARR.

Study Design
We performed a genome-wide association study of aEARR, a derived extreme and welldelimited root resorption phenotype, in UCM3D g consortium participants ( Figure 1). The UCM3D g consortium database is based on a general-population cohort of roughly 0.01 million patients aged 9-67 years old, recruited across Spain between 2005 and 2019.

Sample Size Calculation
Sample size estimation was based on the minimum sample size required for a genetic association study based on: (a) the prevalence of the disease: 0.01; (b) ratio of cases to controls (~1:3); (c) alpha error 1 × 10 −4 /beta error: <0.20; (d) base relative risk: 2.9; and (f) frequency of risk allele: 0.25. Sample size was calculated using the freely available Genetic Power Calculation software [28]. It was determined that 450 patients would be needed to establish an association between the presence of a single nucleotide polymorphism (SNP) of interest under these conditions and the appearance of an advanced state of EARR, 100 cases for the aggressively affected cohort and 350 controls, with 3% overestimation to include expected dropouts were required. A pre-piloted protocol was followed for data filtering and extraction from the UCM 3D g and SALUD R databases. As detailed in Supporting Information File 2, the diagnostic variables, clinical setting variables, and radiological values were collected and subjected to QC (accession registered URI+i ref#:5-201119). Manifest-codifying data errors in the database were removed by selecting and identifying implausible values in each category of diagnostic, clinical, or radiological variables (e.g., Discrepancy index (DI) = 400) or others that were not adequately recorded in terms of the type of unit or quantitative/categorical format (e.g., age = a). Extreme values of duration of mechanical loading (treatment time > 72 months) justified elimination as they were more likely to be data errors or non-physiological extremes rather than feasible variable inputs. The quality-checked data were processed in the next steps ( Figure 1).

Phenotyping and Radiographic Measurements
All subjects selected for final inclusion in the study were assigned to the affected or control cohorts of patients according to radiological screening, using radiographic measurements performed in duplicate and in a double-blinded (P.I. and R.S) manner on orthopantomographic and teleradiographic projections that were already available and used at the clinic for routine diagnosis and treatment. Subjects were classified and divided into these two groups, based on the presence or absence of the phenotype of aggressive post-orthodontic EARR of more than 5 mm in blinded radiographic measurements. The affected cohort included patients with severe EARR > 5 mm and the control group was composed of patients with EARR < 5 mm [7].
The following methods have been detailed previously; an aEARR phenotype was assessed after the roots were measured from before and after treatment on panoramic radiographs focusing on the maxillary central and lateral incisors [4,5,15]. All pre-and post-treatment images were calibrated beforehand and a correction factor for magnification was applied in all cases. Measurements were performed on digital radiographs using diagnostic software (Adobe Photoshop CS8, Adobe Systems Incorporated, San Jose, CA, USA) enabling the image filters to provide maximum precision when localizing the terminal points of the roots. Accordingly, the tooth with the highest EARR value was selected as the dependent variable of interest for that subject using the method described by Linge and Linge [30], modified by Brezniak et al., (2004) [31]. Pre-and post-treatment radiographs through the initial and final root (r1 and r2, respectively) and crown (c1 and c2, respectively) lengths were used to determine the changes in dental and root length. If the root became shorter during treatment, the EARR value resulted from r1-r2 [c1/c2]).
Apical displacement and variation in tooth inclination were quantified using the superimposition of radiographic measurements on a lateral radiograph, using a modified version of the method described by Baccetti et al., (1998) [32].
DNA samples were genotyped using the high-throughput Axiom platform with the GeneTitan ® MC Instrument (Axiom Genome-Wide Human Assay technology, CeGen). This method has been extensively validated in literature [33,34].
Data were analyzed, using the Axiom Analysis Suite 4.0 software for genotype clustering and calling. We then implemented a quality control for the data, where we checked possible sample stratification, missing SNP genotype, SNP monomorphic status, and SNP minor allele frequency (see Supporting Information File 3) from 687,133 markers. Next, to achieve our objectives, we selected 14,377 SNPs located in the X and Y-chromosomes along with other candidate genes at chromosomes 2, 4, 8, 12, and 18 for the analysis (see Supporting Information File 4).

Overall Statistical Analysis of Clinical and Radiological Variables
A descriptive analysis of the data for quantitative and categorical variables based on diagnostic or treatment factors was performed (mean, SD, ranges, frequencies, and distributions). Backward conditional binary logistic regression was used to assess the extent to which specific diagnostic and treatment parameters interfere within the observed aEARR process; odds ratios (OR) with a 95% confidence interval were also calculated. SPSS was used for data analysis (version 22.0; LEAD Technologies, Chicago, IL, USA), with statistical significance set at a value of p < 0.05 [35].

Genetic Association Tests
To select our best SNP candidates, we fitted a logistic regression model adjusted by the type of treatment, total length of treatment (load duration), and gender, and assumed a log-additive genetic model. These statistical analyses were performed using R software (R Core Team, Vienna, Austria. version 3.6.1) [35].

Reliability and Accuracy of the Measurement Method
To prevent inter-observer variation, the same experienced operator (R.S) carried out all the measurements defined previously. However, a second experienced examiner (P.I) replicated the measurements on a subset of 50 patients to calculate the inter-observer error. The kappa coefficient was analyzed to determine concordance between EARR-affected and non-affected group assignment, based on radiological screening, and the result was a value of one. The method error was also calculated for measurements acquired from radiographic records, comparing double measurements of 20 randomly chosen subjects at an interval of 20 days. A paired Student's t-test was used for calculations, with an absence of significance being regarded as indicative of concordance between repeated measurements; the intraclass correlation coefficient (ICC) for absolute agreement was also calculated for both intra-and inter-observer errors. Accuracy of measurement was obtained from the equation: where d is the difference between the double measurements and n is the number of paired double measurements [36].

Phenotyping: Reliability of the Radiographic Measurement Methods and the Associated Errors
aEARR was phenotyped according to radiographic measurements with a threshold for root resorption value greater than 5 mm. Reliability of the measurements retrieved no statistically significant differences between replicated assessments (p > 0.05) and intraclass correlation coefficient (ICC: 0.930), and the concordance index resulted in optimal values (k = 1.00) for both intra-and inter-examiner measurements, respectively [37]. Method error for measurements obtained from panoramic radiographs following the described method was calculated to be below <0.04 mm.

Sample Characteristics, Description and Analysis of aEARR Risk Associated with Clinical Features
The flow chart diagram (Figure 1) describes the sampling filtering steps and the following research strategy. The present study sample included a cohort of patients with relatively homogenous diagnostic characteristics, who had undergone mechanical load during a full orthodontic treatment, as detailed in Table 1. The mean ages of patients treated in the affected and control cohorts were~21 ± 12 and~23 ± 12 years old, respectively, with a fair balance found for sex distribution. The American Board of Orthodontics (ABO) discrepancy index was found to be quite homogenous in both groups [~16 ± 9 and~15 ± 8] with a mean treatment time of~27 ± 9 and~25 ± 8 months, respectively (Table 1). Results from the associative testing by regression analysis are detailed in Table 1. The results showed that recorded clinical factors do not confer an additional risk of aEARR when compared to the control cohort. Specifically, differences in treatment time [OR: 0.974; 95%CI: 0.944-1.005; p = 0.095] or treatment type, extraction vs. non-extraction [OR: 1.668; 95%CI: 0.839-3.316; p = 0.145], did not confer an additional risk for aEARR process when results from both the study cohorts in the regression analysis were compared.

Genotype Distributions and Analysis of aEARR Risk Associated with Genetic Variants at Multiple
Putative Loci on Chromosomes 2, 4, 8, 12, 18, X, and Y.
In addition to diagnostic and treatment-related co-variables, patients included in the present study were genotyped for specific novel target genetic variants and other SNPs explored in previous studies with regard to a risk of any degree of EARR found along the genome in chromosomes 2, 4, 8, 12, 18, X, and Y (see Supporting Information File 3).
When the whole cohort sample was explored within the associated risk for aEARR, a gender-dependent effect was detected to influence the results. The top lead genetic variants associated with aEARR [p value < 1 × 10 −4 ], after gender stratification, focused on two novel putative genes located in the X chromosome, specifically, STAG 2 gene, stromal antigen 2 gene, rs151184635 (prior to and after adjustment for confounding factors) and RP1-30E17.2, clone-based (Vega) gene rs55839915 (after adjustment for confounding factors). These two target variants were found to be associated with an increased risk of aEARR; this effect was particularly restricted to men [OR: 6.09; 95% CI: 2.6-14.23 and OR: 6.86; 95% CI: 2.65-17.81, respectively] (Tables 2 and 3).   In addition to these two top genetic variants (as compiled in Tables 2 and 3), a total number of 27 novel genetic variants out of 14,717 [p value < 0.001] were identified with marginal association values, specifically in the sexual chromosomes; some of them were associated particularly with an independent susceptibility risk of experiencing aEARR in men or women secondary to mechanical orthodontic load. None of the previously studied genetic variants located at chromosomes 2, 4, 8, 12, and 18, also included in the present study, were found to be marginally associated with a positive or negative risk of aEARR [p value > 0.001]. In this respect, when a p value threshold was set to a conventional value of <0.05, the number of potential associations increased substantially as provided in Supporting

Discussion
The present study offers, for the first time, significant valuable data on genetic variants located on chromosomes X and Y that are potentially implicated in aEARR facilitated by orthodontic forces. aEARR has been defined as a permanent loss of apical dental root structure of more than 5 mm. Severe forms of EARR secondary to orthodontic forces are the least frequent type of EARR, which is more often detected as of mild or moderate degrees [7]. Thus, aEARR describes a clearly radiographically identifiable phenotype of EARR with an apical third loss of more than 5 mm that is produced within a limited period, as this is the case of a mean orthodontic treatment length of~20 months [38]. Despite being relatively rare, aggressive phenotypes are prone to exacerbate patient morbidity and are more likely to provoke mechanical and functional disabling consequences for the tooth, potentially inducing irreversible pulp or periodontal damage and triggering an inflammatory and/or infectious process that might result in eventual tooth loss [6,39].
Therefore, from a clinical perspective, it is clearly urgent to suitably define the triggering factors that might contribute to the etiology of this complex aggressive pathology [7]. In line with this, several diagnostic and clinical factors have been previously associated with moderate degrees of EARR secondary to orthodontic forces [2,40]. Demographic factors such as age and gender, morphological factors such as root shape type or clinical factors such as treatment duration, magnitude of orthodontic forces, previous dental trauma, maxillary expansion degree, direction of tooth movement, extraction treatment, use of intermaxillary elastics, or even appliance type have been suggested [9,[41][42][43][44]. Nevertheless, a significant number of controversial related results are also found in literature. Recent meta-analysis and systematic reviews have suggested a greater relevance to some of them, but these are usually based on scientific evidence with low certainty levels [9,45]. Hence, several differing clinical management protocols for aEARR have been suggested [1].
Regression analysis performed in the present study suggests that both cohorts are unlikely to be marginally influenced (p > 0.05) by some of the collected covariables. In this regard, the treatment time has been described as a risk factor influencing EARR severity within a degree-time dependent effect [44]; however, this suggestion is not supported for the case of aEARR based on the present research data. A robust risk association was not directly imputed to treatment length by itself according to comparisons between both cohorts in the current study [44][45][46]. It is thus plausible that some confounding factors associated with treatment length and not treatment time itself, might explain some of the differences in some cases [38,47]. In this respect, some of the lengthiest treatments are very often associated with non-conventional tooth-movement rates or uncharacteristic treatment mechanics, which might be linked with infra-diagnosed conditions or factors, i.e., this may be the case for basal differences in bone remodeling metabolism and rates, and also for differences in craniofacial muscular responses which might affect occlusal loading and even tooth micro-trauma during orthodontic treatment and may even interfere with medications such as AINES and bisphosphonates [48][49][50]. Although a meticulous medical history was recorded; no unquestionable certainty should be presumed from the patients' responses. Further non-compliance, round tripping movements, and movements near the bone cortex might also be associated with lower orthodontic tooth movement rates and/or extended treatment times [45,51]. Even considering these reflections, heterogeneity between study designs, ethnic differences, and the differences within phenotypic characterization might contribute to some of the observed results [17,23,52,53].
Apart from clinical-related factors, accumulating evidence supports the influence that genetic factors exert over the occurrence of post-orthodontic EARR with moderate severity. Previously published studies have provided useful evidence regarding the suggestive role of some specific genetic variants within the EARR process following candidate gene approaches [15][16][17][18][19][20][21][22][23]. In this respect, none of the previously studied genetic variants, also examined in the current study, showed robust statistically significant associations with aEARR [p > 0.0001], while a few of them (SSP1: rs11730582, P2RX7: rs1718119, TNFRSF11A: rs8086340) showed marginal associations [p < 0.05]. Interestingly, these marginal associations were found just in the female group, not showing the same trend in males, prior to and after adjustment with the clinical confounders. This might suggest a gender-specific effect. In this respect, this is the first study to perform a deep analysis regarding the influence of more than 14,000 specific genetic variants located on sexual chromosomes within an aggressive phenotype of EARR in the context of mechanical orthodontic loading. Moreover, this paper offers very valuable data regarding multiple novel putative loci and the genes potentially implicated in this type of extreme phenotype located not chromosomes X and Y, as well as at the level of other candidate genes with less power imputation located at autosomes 2, 4, 8, 12, and 18. Specifically, just two variants located at chromosome X, STAG 2 rs151184635 and RP1-30E17.2 rs55839915, where identified as the best associated SNPs [p value < 0.0001]. Interestingly, these associations showed a gender-dependent association limited to men [54,55]. We detected several statistically significant sex-specific interactions for many SNPs in the targeted genes. This sexual dimorphism is present in a vast majority of human pathologies based on genetic and hormonal differences that might modulate gene expression increase or minimize disease risk and progression [56][57][58][59]. In this regard, differences in bone remodeling rate have been described to be highly influenced by sex-specific features along with potential differences in bone mineral density and metabolism, or even hormone balance status [48,51,60]. Moreover, other gene expression factors directly linked to cytotoxic protection mechanism or related cytokine secretion might be underlying in the sex-specific differences as in other pathological entities [61][62][63]. Future studies should thus clarify if these results imply a potential interaction of these candidate genes with specific pathways related to gonadal steroids or additional gender-dependent mediators [64].
The top identified associations in the current study targeted rs151184635 and rs55839915 variants. RP1-30E17.2 [Ensembl: ENSG00000225689.1] contains 216,374 genetic variants with at least four splice variants. The SNP rs55839915 associated in the present study overlaps in three different transcripts being a long intergenic non-coding RNA type (lincRNA) that lacks coding potential [65]; however, lncRNAs have been described as regulators of gene expression, scaffold formation, and epigenetic control mediating within different pathological complex entities and physiological processes [66,67]. No specific underlying mechanisms have been described associated with this variant, and so, further studies should provide some potential explanatory hypothesis for aEARR. Meanwhile, STAG 2 encodes stromal antigen 2 protein, a subunit of the cohesin complex [Ensembl: ENSG00000101972 MIM: 300826]. This gene is linked to separation of chromatids during cell division, and its inactivation is associated with several types of human cancer [68]. The rs151184635 variant codes a non-coding transcript variant occurring within an intron overlapping three different transcripts at STAG2, TEX13D, and SH2D1A, all long non-coding RNA genes [69][70][71][72]. Although no direct robust functional consequence has been described in literature, SH2D1A, which is a mediator in cytolytic pathways as key activator of T-and NK cell-cytotoxicity, showed differential gene expression associated with specific immunopathologies as is the case of systemic juvenile idiopathic arthritis, associated with the potential onset of macrophage activation syndrome [73]. In this regard, extremely high overproduction of pro-inflammatory cytokines, IFN-γ, IL-1, IL-18, TNF, IL-2, IL-6, and macrophage colony-stimulating factor (M-CSF), as well as repressors such as circulating TNF receptors and IL-1ra are correlated with this macrophage-based pathology [74]. In connection with this, it has been described that prolonged mechanical strain within the periodontal ligament (PDL) around the tooth root is associated with an increase in CD68+, iNOS+ M1-like macrophages, an imbalance in the M1 > M2 ratio polarization and root resorption pathology linked to IFN-γ oversecretion by T-cells and PDL stem cells [75,76]. Whether any potential functional consequence in terms of alternative splicing or gene expression modulation might have a direct/indirect effect on the onset of aggressive phenotype of EARR, remains to be a remote hypothesis that should be explored deeply based on robust future molecular studies. Interestingly, X-chromosome activation of the vast majority of genes linked to the X chromosomes occurs in one of the two X chromosomes of any cell in women. Random or quasi-random selection of which X chromosome will remain inactivated in females follows a different process in males, which might lead to differential expression or repression of some inherited X-linked genes. This raises the possibility that expression of specific long non-coding RNAs, as should be the case for STAG2, TEX13D, or SH2D1A, might exert an influence in modulating gene expression in specific domains that escape X silencing, as occurs in other species, in a gender-dependent way. This should partially explain whether the genetic variants (rs151184635 and rs55839915) in the present study are at least marginally associated in males but not in females [77][78][79].

Limitations
The present study recruited a representative cohort of severely EARR affected patients that were radiographically screened by means of panoramic projections measurements >5 mm as assessed by two experienced examiners [7]. The radiographic diagnostic method used for patient assignment to affected and control groups is not exempt from limitations in terms of accuracy compared to other x-ray methods such periapical radiographs, Cone Beam Computerized Tomography (CBCT), Computerized Tomography (CT), or ex vivo methods (histological and/or histochemical) [80]. Nevertheless, panoramic radiograph assessments might produce sufficient reliability to perform absolute linear measurements when the head position is controlled in a range of 10 • of the inclination in regard to the horizontal plane. Moreover, panoramic radiographs have been associated with underestimation of the tooth root lesion in mild lesions; however, this type of screening for assessing such types of aggressive phenotypes (>5 mm) is less prone to misclassification as shown by optimal reproducibility and error method results retrieved by the two experienced operators in the present study.
Secondly, although the UCM3Dg consortium database provides a unique opportunity to investigate the underlying influence of genetics over aEARR, the sample size is still relatively moderate, once the obtained results are re-analyzed, and the present design did not use an independent external cohort as a replication study due to the particularity of this aggressive phenotype. Thus, spurious associations remain a challenge and the external validity of the present findings need to be confirmed in further investigations by using different external cohorts to test the model.
Lastly, differences observed between study results in terms of the diagnostic, clinical, and genetic factors associated with EARR of different degrees are a major concern that should be addressed. With respect to this, it seems plausible that the heterogeneity found between study designs, ethnic differences [52], absence of a standardized and unique phenotypic characterization of EARR in terms of diagnostic methods and diagnostic criteria values/thresholds between studies might have contributed to some of this controversy. Nevertheless, this is a critical issue that should be clearly revisited and global consensus-based standards should be achieved in this regard to grant a next generation improvement, not only in internal validity of the studies but equally critical, in external validity of the research in the field. This should support bench-to-clinic research findings with increased certainty levels [1,81,82].
Despite all the aforementioned relevant concerns, this study offers novel and extremely valuable data regarding multiple putative loci and genes located at chromosomes X and Y potentially implicated in this type of extreme phenotype, along with other candidate genes with less imputation power at chromosomes 2, 4, 8, 12, 18, X, and Y. To the best of our knowledge, this is the first XY-chromosome-wide association study (XYWAS) to investigate sex-specific genetic effects on XY chromosomes within one of the most severe phenotypes of EARR in the context of orthodontic forces.

Conclusions
Multiple putative genetic variants located at chromosomes X and Y are potentially implicated in an extreme phenotype of aEARR. Particularly, STAG2 genetic variants rs151184635 and RP1-30E17.2 genetic variant rs55839915 were found to be associated with an increased risk of being afflicted with aEARR, only in men.