DNA-based prediction of external ear morphology in the Chinese population: an exploratory study

Samples collection

A total of 675 Chinese volunteers, aged 18–93 years (429 males and 246 females), participated in this study. DNA was extracted from oral swabs using the phenol–chloroform isoamyl alcohol (Butler, 2012) and quantified by NanodropTM 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Five digital photographs of the head were taken at the Frankfort horizontal plane using a Canon 80D camera: the left side (90°), the left angle (45°), a frontal view (0°), the right angle (45°), and the right side (90°) (Adhikari et al., 2015). This study was approved by the Ethics Committee at Sichuan University (k2020032). This work was carried out in accordance with the principles of the Declaration of Helsinki, and all donors gave written informed consent.

Human ear phenotypes

Right side, right angle, and frontal photographs were used to assess 13 ear traits on a three-point categorical scale (level_0, level_1, and level_2). The ear traits evaluated included ear protrusion, lobe attachment, lobe size, tragus size, antitragus size, intertragic incisure, superior helix rolling, posterior helix rolling, folding of the antihelix, antihelix curvature, crus helix expression, superior crus of antihelix expression, and Darwin’s tubercle (Fig. 1). These traits and the scoring standards were in accordance with previous studies (Rani, Krishan & Kanchan, 2022; Adhikari et al., 2015; Verma, Bhawana & Vikas, 2014). Intraclass correlation coefficients (ICCs) were calculated to assess inter- and intra-observer reliability (Shrout & Fleiss, 1979). The detailed methods were described in Table S1. All photographs of the volunteers were scored by the same rater (X.W.). The frequency distribution of the three categories for each ear trait was presented in Fig. S1. The Spearman’s rank correlation coefficients were calculated to assess the correlations among the various ear traits.

Candidate marker selection

The SNPs reported by large-scale GWASs that exhibited statistically significant associations with ear morphology were initially selected (Adhikari et al., 2015; Shaffer et al., 2017; Wang et al., 2022; Li et al., 2023). Duplicate variants identified across different studies were carefully eliminated. We then refined our selection by filtering out markers with a minor allele frequency (MAF) below 0.05, based on population-specific data from the Beijing Han Chinese (CHB) and Southern Han Chinese (CHS) of the 1000 Genomes Project Phase 3 (https://www.internationalgenome.org). LDlink online software (https://ldlink.nih.gov/?tab=home) was utilized to exclude SNPs that were in high linkage disequilibrium (LD) with each other in the CHB and CHS, based on the GRCh38.p14 reference genome (Machiela & Chanock, 2015). Specifically, an LD threshold of r² < 0.5 was set to ensure that our final list of SNPs was as independent as possible.

Multiplex SNaPshot assay

Primers for amplification and single-base extension (SBE) were designed following the operating instructions of the PyroMark Assay Design 2.0 software. To validate the specificity of the polymerase chain reaction (PCR) primers, in-silico analysis was performed using the UCSC genome browser (Hendling & Barišić, 2019). Meanwhile, both PCR and SBE primers were analyzed using Autodimer to detect potential hairpin and primer dimer formations (Vallone & Butler, 2004). The PCR primers were further examined according to the results of polyacrylamide gel electrophoresis (PAGE), and the SBE primers were examined according to the single-site SBE assay.

The multiplex PCR was performed in a five-µl reaction volume containing two µl PyroMark PCR Master Mix (Qiagen, Hilden, Germany), one µl of mixed PCR primers (Table 1), two ng of DNA, and RNase-free water. The PCR amplification was carried out under the following conditions: initial denaturing at 95 °C for 15 min, 20 cycles at 94 °C for 30s, 65 °C for 90s (drop down 0.5 °C each cycle), 72 °C for 30s, 19 cycles at 94 °C for 30s, 55 °C for 90s, 72 °C for 30s, followed by 72 °C for 10 min. The following PCR procedures were carried out as described in Jiang et al. (2023). Specifically, two µl shrimp alkaline phosphatase (SAP) and one µl of exonuclease I were added to the PCR products, incubating at 37 °C for 1 h to eliminate the single-strand primers and dNTP, followed by heat inactivation at 80 °C for 10 min to eliminate the enzyme. Then, the mini-sequencing was performed in five µl reactions containing one µl of the PCR products dealt with digestion, one µl of mixed extension primers (Table 1), 1.2 µl SNaPshot Multiplex Ready Reaction Mix, and 1.8 µl RNase-free water. It was conducted under 26 cycles of 96 °C for 10 s, 50 °C for 5 s, and 60 °C for 30 s. To remove the remnant ddNTP, one µl SAP was then added to the product by incubation at 37 °C for 1 h, followed by heat inactivation at 80 °C for 10 min.

Purified products (one µL) were combined with nine µL of HiDi™ Formamide sizing standard mixture. This sizing standard mixture was prepared by adding three µL GeneScan™ 120 LIZ^® Size Standard to one mL HiDi™ Formamide. Electrophoresis was performed on an ABI 3130xl Genetic Analyzer (Thermo Fisher Scientific, Waltham, MA, USA) with 36 cm capillary arrays and POP-7™ polymer. The specific reaction conditions for 3,130 were as follows. Samples were injected at 2.5 kV for 10s, and electrophoresed at 13 kV for 600 s at a run temperature of 60 °C. Raw data were analyzed in GeneMapper^® ID v3.2 software (Applied Biosystems, Foster City, CA, USA) with a peak detection threshold of 100 relative fluorescence units (RFU).

Association analysis

Allele frequencies for each SNP were calculated based on our dataset of the Chinese population. All genotype data were tested for Hardy-Weinberg (HW) equilibrium and LD analysis using the online software SHEsis (http://analysis.bio-x.cn/myAnalysis.php) (Shi & He, 2006).

DNA variants with an MAF higher than 0.05 were tested for their association with each ear phenotype on the entire dataset of 675 samples using multinomial LR (MLR) analysis. Allelic odds ratios (ORs) with 95% confidence intervals and respective p-values were determined for minor alleles categorized in an additive manner. The effect of SNPs on the respective phenotypes was measured by the ORs (Noreen et al., 2023). P-value < 0.05 was considered nominally significant. The association analyses were conducted in R version 4.4.1 using the ‘nnet’ package.

Prediction model

A comparative analysis was conducted to assess the efficacy of five predictive models for three-class classification: MLR, SVM, RF, AdaBoost, and k-nearest neighbors (KNN). The SVM model was implemented with a radial basis function kernel and a regularization parameter set to 1.0, with probability estimates enabled. The RF classifier was constructed with 100 decision trees. The KNN model employed three nearest neighbors, while AdaBoost model was configured with 100 estimators and utilized the SAMME algorithm. These parameter selections were intended to optimize each model’s performance in the multi-class classification setting. Ten-fold cross-validation was employed to guard against arbitrary partitions of the dataset. Specifically, the dataset was randomly split into ten equal subsets. During each fold of the cross-validation process, nine of these subsets were designated as the training set, while the remaining subset served as the testing set. For each category, positive predictive value (PPV), negative predictive value (NPV), sensitivity, specificity, and area under the curve (AUC) values were calculated based on the average performance across the ten testing sets obtained from the iterations.

Specific categories of ear traits, exhibiting good predictive performance in their respective three-class classification tasks, were selected. To simplify the classification process, binary prediction models were developed for each selected category. Additionally, interactions between genetic variants were examined using the MDR 3.0.2 software. The multifactor dimensionality reduction (MDR) method is a powerful strategy for detecting and interpreting statistical locus-locus epistasis (Moore et al., 2006; Moore, 2004). As a data mining technique, it is used to detect and characterize non-linear relationships between variables (Moore & Williams, 2009; Ritchie et al., 2001). The dendrogram graphs provided by MDR illustrate the presence, strength, and nature of epistatic effects (Moore et al., 2006). Furthermore, the genetic interactions were tested by including interactions into our binary prediction models. The predictive workflow was illustrated in Fig. S2. Alternative thresholds of probability for the phenotype prediction were tested, ranging from p > 0.5 to p > 0.85 with a 0.05 interval. Prediction modelling was carried out using Python version 3.7.

Qualitative ear phenotypes

The ICCs ranged from 0.47 to 0.97, indicating moderate to good intra-rater and inter-rater consistency (Table S1). The scores for ear traits examined in the Chinese population showed a weak correlation between them (Table S2), with Spearman correlation coefficients all below 0.3. Specifically, there was a correlation between antihelix fold and superior crus of antihelix expression (r = 0.218). In contrast, lobe size exhibited a significant negative correlation with the ear protrusion (r = −0.207).

Five ear traits exhibited weak but statistically significant correlations with sex: superior helix rolling (r = 0.343), lobe size (r = 0.219), ear protrusion (r = −0.209), lobe attachment (r = 0.127), and posterior helix rolling (r = −0.178). Separately, age was negatively correlated with the scores of four ear traits: antitragus size (r = −0.137), superior helix rolling (r = −0.136), antihelix curvature (r = −0.126), and lobe attachment (r = −0.124) (Table S3).

Candidate markers selection and multiplex SNaPshot assay

A list of 45 SNPs was initially selected following an exhaustive review of the literature on ear morphology. MAF screening and LD analysis were then applied to refine our selection, retaining 23 SNPs. However, designing feasible primers for six SNPs was challenging, including rs12695694, rs9496426, rs10211400, rs1602631, rs57788627, and rs62169502. Incorporating SNP rs2742261 into the multiplex system with other SNPs also posed difficulties. Furthermore, upon examination of our dataset, we found that the SNP rs10824309 was monomorphic and thus excluded from the multiplex assays.

Ultimately, after rigorous screening and examination, a total of 15 genetic markers (rs10198822, rs6802174, rs74030209, rs17023457, rs17034666, rs3827760, rs7812632, rs1948400, rs6699106, rs263156, rs3789101, rs7771119, rs62169501, rs1619249, rs1960918) were included in the development of the final multiplex SNaPshot assays (Fig. 2). All peaks were detected with a 100 RFU threshold and all blood samples were successfully genotyped, underscoring the successful establishment of the robust and reliable SNaPshot system.

Association analysis

The allele frequencies were shown in Table S4, and deviations from HWE were noted for eight SNPs (Table S5). Details of LD analysis are shown in Table S6, indicating that all the SNPs were in linkage equilibrium.

All SNP markers included in the association analysis, except for rs17034666, were found to be associated with ear traits in the Chinese cohort. Among these SNPs, rs6802174 exhibited significant associations with six ear traits: ear protrusion, lobe attachment, tragus size, antihelix fold, antihelix curvature, and Darwin’s tubercle. Subsequently, rs74030209 was significantly associated with five ear traits: ear protrusion, lobe attachment, lobe size, antihelix curvature, and crus helix expression. Similarly, rs1948400 demonstrated significant associations with five ear traits, including ear protrusion, lobe attachment, intertragic incisure, antihelix curvature, and Darwin’s tubercle (Fig. 3). Due to the large number of ear traits included in our study, we have summarized the SNP marker that showed the strongest association with each trait, while the detailed information has been provided in Tables S7 to S20.

Four genetic variants (rs6802174, rs74030209, rs1948400, and rs7771119) were significantly associated with ear protrusion, with rs6802174 showing the strongest association. Compared to those with the CC genotype in rs6802174, individuals with the CG genotype had 3.09 times higher odds of having a large protrusion (level_2) (p-value = 3.03E−05). Furthermore, individuals with the GG genotype had 2.36 times higher odds of having a large protrusion (p-value = 0.0468).

Three SNPs (rs6802174, rs74030209, and rs1948400) were significantly associated with lobe attachment. Compared to those with the CC genotype in rs6802174, individuals with the GG genotype had a 2.70-fold increased likelihood (p-value = 0.0122) of having average lobe attachment (level_1). Three SNPs (rs74030209, rs7812632, and rs263156) explained the variation in lobe size. Compared to individuals with the CC genotype in rs74030209, those with the CT genotype have a significantly lower likelihood of having large lobe size, with ORs of 0.48 (p-value = 0.0176) and 0.35 (p-value = 0.0019), respectively.

Five SNPs (rs6802174, rs3827760, rs263156, rs3789101, rs1960918) were significantly associated with tragus size. In rs3827760, individuals with the AA genotype were 0.02 times (p-value = 0.0125) less likely to have large tragus size (level_2) compared to those with the GG genotype. Two genetic variants, including rs17023457 and rs3827760, were observed to be associated with antitragus size. In rs17023457, individuals with the CC genotype had a 2.26-fold increased likelihood (p-value = 0.0155) of having an average antitragus size compared to those with the TT genotype.

Four genetic predictors (rs10198822, rs1948400, rs62169501, and rs1960918) were significantly associated with intertragic incisure. Individuals with the TC genotype in rs10198822 had a 0.44-fold decreased likelihood (p-value = 0.0011) of having horseshoe-shaped intertragic incisure (level_1) compared to those with the TT genotype. In contrast, only one SNP, rs10198822, was observed to be associated with superior helix rolling. Individuals with the TC genotype exhibited a notable increase in the likelihood of a higher degree of superior helix rolling compared to those with the TT genotype. The odds increased by 3.58 times (p-value = 0.0443) for partial folded superior helix rolling (level_1), and by 5.72 times (p-value = 0.0057) for over-folded superior helix rolling (level_2).

The statistical significance was obtained for four SNPs (rs6699106, rs7771119, rs62169501, and rs1619249), which explained the variation in posterior helix rolling. In rs7771119, the presence of allele A was associated with a decreased risk of having a higher degree of posterior helix rolling. In rs7771119, allele A was associated with decreased risk of high-degree posterior helix rolling. Compared to CC, CA had 0.43 times lower odds (p-value = 0.0083), and AA had odds 0.18 times those of CC (p-value = 0.0149). Three SNPs—rs6802174, rs7812632, and rs3789101—were associated with the antihelix fold. In rs7812632, individuals with the GG genotype were 0.27 times (p-value = 0.0372) less likely to have partial folded antihelix (level_1) compared to those with the CC genotype.

Four genetic predictors (rs6802174, rs74030209, rs1948400, and rs3789101) were significantly associated with antihelix curvature. In rs3789101, individuals with the CC genotype had a 0.27-fold decreased likelihood (p-value = 0.0303) of having a strong antihelix curvature (level_2) compared to those with the GG genotype. Two variants (rs74030209, rs6699106) were associated with crus helix expression. Compared to individuals with the CC genotype in rs6699106, those with the TT genotype had a significantly lower likelihood of having a higher degree of crus helix expression, with ORs of 0.38 (p-value = 0.0135) and 0.36 (p-value = 0.0359), respectively.

Three SNPs (rs17023457, rs3827760, and rs7812632) explained the variation in the superior crus of antihelix expression. Compared to those with the GG genotype in rs3827760, individuals with the GA genotype had an increased odds of having prominent antihelix superior crus (level_2) by 2.64 times (p-value = 0.0309). Five SNPs (rs6802174, rs17023457, rs7812632, rs1948400, and rs263156) were significantly associated with Darwin’s tubercle. In rs263156, individuals with the GT genotype had a 0.34-fold decreased likelihood (p-value = 0.0229) of having a prominent Darwin’s tubercle (level_2) compared to those with the GG genotype.

Prediction results

For the 13 ear traits, the micro-average AUC in the three-class prediction for ear phenotypes ranged between 0.50 and 0.60. The MLR and AdaBoost models exhibited better prediction accuracy, whereas the SVM and RF models showed reduced performance, and the KNN model performed the worst (Fig. 4). Notably, for absent tragus (level_0), the predicted AUC exceeded 0.70 in the MLR and AdaBoost models (Fig. 4). The indicators describing prediction accuracy for the five models were presented in Table S21.

For absent tragus (level_0), the AdaBoost model achieved a prediction accuracy with an AUC of 0.73, an NPV of 0.78, a PPV of 0.45, a sensitivity of 0.54, and a specificity of 0.71. In contrast, the MLR model demonstrated lower accuracy, with an AUC of 0.71, an NPV of 0.78, a PPV of 0.49, a sensitivity of 0.54, and a higher specificity of 0.75. For medium tragus size (level_1), the AdaBoost model demonstrated an AUC of 0.65, an NPV of 0.71, a PPV of 0.48, a sensitivity of 0.38, and a specificity of 0.76. Meanwhile, the MLR model showed a lower accuracy, with an AUC of 0.60, an NPV of 0.68, a PPV of 0.38, a specificity of 0.69, and a sensitivity of 0.35. When dealing with prominent tragus cases (level_2), the AdaBoost model demonstrated an AUC of 0.62, alongside an NPV of 0.73, a PPV of 0.47, a sensitivity of 0.47, and a specificity of 0.72. The MLR model showed comparable performance, with an AUC of 0.66, an NPV of 0.72, a PPV of 0.45, a sensitivity of 0.44, and a specificity of 0.72.

In order to reduce model complexity, the tragus size was further converted into binary classifications. Specifically, tragus size was categorized into present (level_1 and level_2) versus absent (level 0). Predictive models for the binary classification task were established using AdaBoost and binary LR (BLR) models. The AdaBoost model achieved an AUC of 0.73, with an NPV of 0.59, a PPV of 0.68, a sensitivity of 0.65, and a specificity of 0.62 (Table 2). Similarly, the BLR model exhibited an NPV of 0.63, a PPV of 0.70, a specificity of 0.67, a sensitivity of 0.67, and a lower AUC of 0.71 (Table 2).

MDR analysis revealed a redundant interaction for the absent tragus (level_0) between rs3827760 in EDAR and rs1960918 in LRBA, which was indicated by blue lines in Fig. S3. When incorporated into the AdaBoost model, this interaction led to an improvement in the AUC to 0.74 (ΔAUC = 0.0061) (Fig. 5A). Similarly, incorporating the interaction into the BLR model resulted in an AUC increase to 0.72 (ΔAUC = 0.0172) (Fig. 5B). The confusion matrices and precision–recall curves of BLR and AdaBoost were displayed in Figs. S4 and S5, respectively. For both models, prediction accuracy parameters improved with the increase of the probability threshold (Table S22). For the AdaBoost model, setting the threshold at p > 0.65 led to high sensitivity, specificity, and PPV, while keeping the overall prediction error relatively low. However, at higher thresholds, the AdaBoost model tended to adopt a more conservative approach, resulting in a greater number of samples being classified as uncertain. In contrast, the BLR model exhibited more stable performance across different thresholds. Specifically, when the threshold was set to p > 0.70, the BLR model not only achieved high specificity, sensitivity, and PPV but also struck a good balance between accuracy and the ability to handle uncertain cases.