Multiethnic GWAS Reveals Polygenic Architecture of Earlobe Attachment.

The genetic basis of earlobe attachment has been a matter of debate since the early 20th century, such that geneticists argue both for and against polygenic inheritance. Recent genetic studies have identified a few loci associated with the trait, but large-scale analyses are still lacking. Here, we performed a genome-wide association study of lobe attachment in a multiethnic sample of 74,660 individuals from four cohorts (three with the trait scored by an expert rater and one with the trait self-reported). Meta-analysis of the three expert-rater-scored cohorts revealed six associated loci harboring numerous candidate genes, including EDAR, SP5, MRPS22, ADGRG6 (GPR126), KIAA1217, and PAX9. The large self-reported 23andMe cohort recapitulated each of these six loci. Moreover, meta-analysis across all four cohorts revealed a total of 49 significant (p < 5 × 10-8) loci. Annotation and enrichment analyses of these 49 loci showed strong evidence of genes involved in ear development and syndromes with auricular phenotypes. RNA sequencing data from both human fetal ear and mouse second branchial arch tissue confirmed that genes located among associated loci showed evidence of expression. These results provide strong evidence for the polygenic nature of earlobe attachment and offer insights into the biological basis of normal and abnormal ear development.

Introduction

Earlobe attachment (MIM: 128900) is often presented as an example of a readily observable Mendelian phenotype in educational materials and continues to be studied as a Mendelian phenotype in the contemporary primary literature (for example, see Ordu et al.). Yet, as early as 1937, Wiener pointed out that earlobe attachment is likely to be a polygenic trait exhibiting a continuous phenotypic distribution. Although earlobe attachment is a neutral morphological trait, understanding its genetic etiology is valuable in that it offers a glimpse into the biological basis of ear development, improving our understanding of genes potentially involved in developmental defects. Moreover, it serves as an instructive example of simple versus polygenic inheritance in an accessible trait.

Recent genome-wide association studies (GWASs) have investigated earlobe attachment4, 5 and reported significant associations with variants in EDAR (MIM: 604095) and SP5 (MIM: 609391). Although promising, these and other suggestive associations have yet to be replicated in independent samples. Of note, ethnic differences in the frequency of lobe attachment are well documented, suggesting that genetic heterogeneity might underlie the trait and that deciphering its genetic architecture might require trans-ethnic studies. This notion is supported by the fact that one of the two previously reported associations was with a missense EDAR variant that is common in Asian and American populations but absent or infrequent in European and African populations.

We propose that large-scale genetic studies of normal human morphological traits can provide insights into the genes and pathways involved in developmental malformations. For example, the external human ear exhibits a highly complex morphology that develops from the first and second branchial arches and requires precise spatial and temporal coordination of tissue proliferation, fusion, and apoptosis. Disruption of these processes can cause birth defects, such as nonsyndromic microtia (congenital under-development of the external ear [MIM: 600674]), which is a fairly common developmental defect that differs in frequency across populations. Moreover, because the jaw and associated masticatory musculature are also derived from the branchial arches, a number of craniofacial syndromes involving arch deficiencies, such as hemifacial macrosomia (MIM: 164210) and Treacher Collins syndrome (MIM: 154500), are characterized by external ear abnormalities. Understanding the genetic factors that contribute to normal structural variation in human ears could provide critical insights into ear morphogenesis, as well as morphogenetic processes in general. In this report, we focus on one aspect of external ear morphology: the lobe.

Specifically, we scanned the genome for variants associated with lobe attachment in 74,660 individuals from four independent and ethnically distinct samples. These included (1) European American (n = 1,791), Latin American (n = 5,062), and Chinese (n = 2,857) cohorts in which lobe attachment was scored as a tripartite (free, partially attached, or fully attached) phenotype by an objective expert rater and (2) a European-ancestry cohort (n = 64,950) comprising research participants who were customers of 23andMe, a personal genomics company, and self-reported lobe attachment as a binary (free or attached) phenotype. We performed two nested meta-analyses to combine results across the three expert-rater-scored cohorts and subsequently across all four cohorts. All participants were genotyped on an Illumina genome-wide array and imputed to the 1000 Genomes Project reference panel, and a GWAS was performed after adjustment for necessary covariates and principal components (PCs) of ancestry (see Table S1). Manhattan and quantile-quantile plots for the GWAS results in individual cohorts are shown in Figure S1.

Subjects and Methods

Study Design

We used a nested genome-wide meta-analysis approach to identify and replicate variants associated with ear attachment. A GWAS was performed separately in four independent cohorts representing populations with distinct ancestries. For three of these cohorts (European American, Latin American, and Chinese), lobe attachment was scored by objective raters as a tripartite phenotype. For the fourth cohort (European-ancestry individuals in the 23andMe sample), lobe attachment was self-reported as a dichotomous phenotype. Because of differences in phenotype definition (tripartite versus dichotomous) and collection method (rater-scored versus self-reported), we performed two genome-wide meta-analyses of (1) the three rater-scored cohorts and (2) all four cohorts.

Data Access

Summary statistics for the 10,000 most significant SNPs from the meta-analyses are provided in Table S2. The individual-level genetic data for the expert-rated European American cohort are available from dbGaP: phs000949.v1.p1. Full summary statistics for all SNPs in the European American, Latin American, and Chinese cohorts are available upon request. Summary statistics for the 23andMe cohort can be requested directly from 23andMe and will be made available to qualified researchers under the terms of a data-transfer agreement with 23andMe to protect the privacy of the participants. Please contact David Hinds for more information and to apply to access the data.

Recruitment and Phenotyping

The European American cohort consisted of 1,791 unrelated European-ancestry individuals aged 3–49 years and recruited from Pittsburgh, Seattle, Houston, and Iowa City as part of a 3D Facial Norms Project.11, 12 Participants were screened for conditions affecting craniofacial morphology, including a history of congenital malformations, trauma, and surgery. The Latin American cohort comprised 5,062 participants from the Consortium for the Analysis of the Diversity and Evolution of Latin America (CANDELA) recruited from Brazil, Chile, Colombia, Mexico, and Peru. The Chinese cohort comprised 2,857 ethnic Han Chinese participants recruited from Taizhou in the Jiangsu Province of China as part of the Taizhou Longitudinal Study. All participants provided informed consent, and all study protocols were approved by the institutional review boards of the pertinent research institutions.

In the European American, Latin American, and Chinese cohorts, earlobes were classified as free, partially attached, or attached. An individual was considered to possess attached earlobes if at least one ear was rated as attached. For the European American cohort, two independent observers examined the ears of participants from 3D craniofacial surface images captured by digital stereophotogrammetry. For the Latin American cohort, the same rater scored lobe attachment according to digital photography (Nikon) of the right side (45° angle) and front of the face. For the Chinese cohort, two independent observers scored lobe attachment from 2D digital photography (Canon EOS 600D) of both sides at a 45°angle.

An additional cohort, composed of research participants from the consumer base of 23andMe as previously described,5, 15 was included. For this study, the 23andMe sample comprised 64,950 unrelated individuals of European ancestry, and data on ear attachment were collected via self-reporting in online surveys using a dichotomous (attached or free) phenotype definition. Example imagery of attached and free earlobes was provided to participants for reference. 23andMe research participants provided informed consent and answered surveys online according to a human subjects protocol approved by Ethical and Independent Review Services, an external institutional review board.

Genotyping, Quality Control, and Imputation

Genotyping was performed separately in the four cohorts. For the European American cohort, DNA was extracted from saliva samples and genotyped along with 72 HapMap control samples for 964,193 SNPs on the Illumina HumanOmniExpress+Exome v.1.2 array by the Center for Inherited Disease Research. Genetic data cleaning and quality control have been described in detail previously. In brief, samples were interrogated for sex, chromosomal aberrations, relatedness, genotype call rate, and batch effects. SNPs were interrogated for call rate, discordance among 70 duplicate samples, Mendelian errors among HapMap controls (parent-offspring trios), deviations from Hardy-Weinberg equilibrium, and sex differences in allele frequencies and heterozygosity. For the Latin American cohort, DNA was extracted from blood samples obtained by a certified phlebotomist and genotyped for 730,525 SNPs on the Illumina HumanOmniExpress array. Quality filters included genotyping call rates per participant and per SNP and minor allele frequency (MAF). Because of admixture within the sample, filters for Hardy-Weinberg equilibrium were not implemented. For the Chinese cohort, DNA was extracted from peripheral-blood samples and genotyped for 887,270 SNPs on the Illumina HumanOmniZhonghua-8 array. SNP-level quality filters were applied for missing call rate, MAF, deviation of genotype frequencies from Hardy-Weinberg equilibrium, and technical filters (shown in Table S3).

For the 23andMe cohort, DNA extraction and genotyping were performed on saliva samples by Laboratory Corporation of America clinical laboratories certified by the Clinical Laboratory Improvement Amendments and accredited by the College of American Pathologists. Samples from this cohort were genotyped on one of four Illumina platforms: two versions of the HumanHap550 chip plus 25,000 custom SNPs, the HumanOmniExpress plus custom content to increase overlap with the HumanHap550 platforms, or a fully custom-designed array. Participants with samples that failed to reach 98.5% call rates were re-contacted for a replacement sample and were re-analyzed. Quality filters were applied for genotype call rate, MAF, Hardy-Weinberg equilibrium, and artifact effects by date.

For all studies, unobserved variants were imputed with haplotypes from the 1000 Genomes Project as the reference (phase 1 for the Latin American and 23andMe cohorts and phase 3 for the European American and Chinese cohorts). Pre-phasing (using SHAPEIT2 for the rater-scored cohorts and the company’s own tool according to the Beagle algorithm for the 23andMe cohort) was performed before imputation. Imputation was performed with IMPUTE2 for rater-scored cohorts and with Minimac2 for the 23andMe cohort. For the European and Latin American cohorts, masked variant analysis, in which genotyped SNPs were imputed for assessment of imputation quality, indicated high accuracy of imputation. Table S4 shows imputation quality filters.

Population Structure

To assess population structure, we performed principal-component analysis (PCA) within each cohort by using subsets of uncorrelated SNPs. Plots of the top PCs of ancestry for the European American, Chinese, and 23andMe cohorts are shown in Figure S5. The complex population structure in the Latin American sample was the focus of a previous investigation. On the basis of scatterplots of the PCs and scree plots of the eigenvalues, we determined that adjustment for 4, 5, 0, and 5 PCs was necessary for the European American, Latin American, Chinese, and 23andMe cohorts, respectively. In the 23andMe cohort, we compared phased genomic segments with reference data across 31 populations to assign the mostly likely ancestry source of each segment. We aggregated local-ancestry assignments to determine the overall proportions of ancestry of each individual. Of the 68,965 consenting 23andMe participants with available phenotypes, 64,950 individuals were determined to have >97.5% European ancestry and were included in this study. In general, genetically determined European ancestry closely matched the self-reported ancestry of the 23andMe participants.

Association Analyses

Earlobe attachment was analyzed as a semiquantitative phenotype (coded 0, 1, or 2 for free, partially attached, or attached earlobes, respectively) separately in the European American, Latin American, and Chinese cohorts. We tested genetic association while adjusting for necessary covariates (such as age, sex, height, and body mass index; see Table S1) and PCs of ancestry by using linear regression under the additive genetic model. For the analysis of the X chromosome, we coded genotypes as 0, 1, or 2 per the additive genetic model for females and as 0 or 2 for males in order to maintain the same scale between sexes. GWAS results for each study were combined via inverse-variance-weighted meta-analysis.

For the 23andMe cohort, earlobe attachment was analyzed as a binary phenotype (coded 0 for free or 1 for attached). We used logistic regression including adjustment for age, sex, genotyping platform, and the top five PCs to test for genetic association under the additive genetic model. Results across all four cohorts were combined via Stouffer’s method of inverse-weighted meta-analysis (based on p values, direction of effect, and sample size). This method of meta-analysis was chosen because it is robust to differences in the scale of the effect estimates among the expert-rater-scored cohorts and between the expert-rater-scored and self-reported 23andMe cohort as a result of differences in phenotype assessments. Compared with meta-analysis methods that use effect sizes and standard errors, Stouffer’s method results in only a small loss of efficiency, which is outweighed by its robustness to known and unknown phenotype differences across cohorts, and it does not require an assumption that effects are the same across cohorts. We used the binomial test (i.e., sign test) to model the consistency of direction of expert-rater-scored effects with the 23andMe cohort.

Functional Annotation

We used HaploReg to query the lead SNP (i.e., the SNP with the smallest p value) at each associated locus in order to extract evidence of functional variation (promotor and enhancer histone marks, DNase hypersensitivity, expression quantitative trait locus [eQTL] information) for all SNPs in linkage disequilibrium (LD; r2 > 0.8) with the lead SNP. 351 genes of interest were defined on the basis of a physical proximity of 500 kb to the lead SNP at each locus. These genes were queried in a number of online databases. We used Mouse Genome Informatics (MGI) to annotate expression in relevant tissues and phenotypic consequences and used the VISTA Enhancer Browser to annotate active enhancer elements in relevant tissues. We used OMIM, PubMed, DECIPHER, and ClinVar to annotate human phenotypic information. We performed genomic enrichment analyses by using the Genomic Regions Enrichment of Annotations Tool (GREAT).

Tissue Collection, RNA Isolation, and Sequencing

Branchial arch two tissue was dissected from wild-type embryonic day 10.5 mice, as well as equivalently staged homozygous sbse and dmbo embryos, and snap frozen on dry ice. RNA was isolated with the QIAGEN RNeasy Mini Kit (74104, QIAGEN), and its quality and concentration were assessed with an Agilent 2200 TapeStation system at the Genomics Core of the Fred Hutchinson Cancer Research Center. All samples had an RNA integrity number greater than 8. For each genotype, RNA was pooled from three male pups, and a total of 1 μg RNA of each genotype was sent to the Genomic Services Lab at the Hudson Alpha Institute for Biotechnology for preparation of indexed directional libraries and ribosomal reduction RNA sequencing (RNA-seq). Samples were paired-end sequenced with 250 million reads on an Illumina HiSeq v.4 PE100. The Cufflinks software suite was used for transcriptome assembly and differential expression analysis. We used the DESeq2 package to weight expression with count data and the Integrative Genomics Viewer to visualize sequences.28, 29

For the analyses of human fetal RNA, ear tissue was obtained (after informed parental consent was provided) from material collected by the Birth Defects Research Laboratory (under approval by the institutional review board of the University of Washington). The gestational age of conceptuses, reported as the number of days after fertilization, was estimated from fetal foot length. The tissue was snap frozen, and RNA was processed as described above.

Results

Earlobe-Attachment Loci Observed in Trans-ethnic Meta-analysis

Rates of lobe attachment differed across the cohorts (Table S1), which was expected given the known differences across ethnic groups. Meta-analysis of the GWAS results from the three expert-rater-scored cohorts (see Figures 1A and 1B) yielded six loci that were significantly (i.e., p < 5 × 10−8) associated with earlobe attachment (Table 1; Figure 2): 2q13, 2q31.1, 3q23, 6q24.2, 10p12.2, and 14q13.1. These loci included the genes EDAR (2q13; lead SNP: rs3827760), SP5 (2q31.1; lead SNP: rs6756973), MRPS22 (MIM: 605810; 3q32; lead SNP: rs9866054), ADGRG6 (LOC153910 or GPR126 [MIM: 612243]; 6q24.2; lead SNP: rs58122955), KIAA1217 (MIM: 617367; 10p12.2; lead SNP: rs7096127), and PAX9 (MIM: 167416; 14q13.1; lead SNP: rs1950357). One missense variant, rs3827760 in EDAR, was observed among the six lead SNPs and other variants in high LD (r2 > 0.8) with the lead SNPs. The others were either intronic or intergenic, and as summarized by HaploReg, several showed evidence of DNase hypersensitivity and histone marks indicative of enhancer or promoter regulatory elements in skin and other cell types. None of the lead SNPs or variants in high LD with the lead SNPs were known eQTLs.

Figure 1

Genome-wide Scans

(Left) Manhattan plots showing the –log10-transformed p values (y axis) by physical genomic position (x axis) for each SNP in (A) the meta-analysis of the three rater-scored cohorts and (C) the meta-analysis of all four cohorts. The horizontal line represents the threshold for genome-wide significance (p < 5 × 10−8). (A) Six significant loci (green) were observed, and genes near the lead SNP in each locus are annotated. (C) 49 associated loci were observed: the same six loci in (A) are shown in green, and these reached genome-wide significance in more than one cohort; the 15 loci in blue showed genome-wide significance in one cohort and replication-level significance (p < 0.001) in at least one additional cohort or the meta-analysis of expert-rater-scored cohorts; and the 28 loci in red were observed via meta-analysis.

(Right) Quantile-quantile plots showing the observed distribution of –log10-transformed p values (y axis) against the expected distribution (x axis) under the null hypothesis of no association (diagonal line) for (B) the meta-analysis of the three rater-scored cohorts (genomic inflation factor = 1.066) and (D) the meta-analysis of all four cohorts (genomic inflation factor = 1.563). The presence of significantly associated loci is indicated by the deviation of observed p values from the tail of the null distribution, as shown by points above the diagonal in the upper right of the plots.

Table 1

Evidence of Association for the Lead SNP in Each Significant (p < 5 × 10−8) Locus Nominated in the Meta-analysis across the Rater-Scored Cohorts

	2q13	2q31.1	3q23	6q24.2	10p12.2	14q13.1
Gene candidate(s)	EDAR	SP5	MRPS22, FOXL2	ADGRG6 (GPR126)	KIAA1217, ARHGAP21	PAX9, NKX2-8
Lead SNP	rs3827760	rs6756973	rs9866054	rs58122955	rs7096127	rs1950357
Base position	109,513,601	171,542,573	138,997,688	142,921,276	24,506,439	37,209,698
Functional position	missense	intronic	intergenic	intronic	intronic	intronic
Data source	genotyped	imputed	genotyped	imputed	imputed	imputed
Minor/major alleles	G/A	C/G	A/G	A/G	C/T	C/A
European American
MAF	0.015	0.413	0.257	0.223	0.451	0.386
Beta	0.028	0.073	−0.041	0.119	−0.034	0.086
SE	0.081	0.020	0.022	0.023	0.019	0.020
p value	0.733	1.84 × 10−4	0.059	2.78 × 10−7	0.077	2.40 × 10−5
Latin American
MAF	0.404	0.661	0.550	0.185	0.490	0.359
Beta	0.062	0.088	−0.025	0.067	−0.043	0.027
SE	0.012	0.011	0.011	0.013	0.010	0.010
p value	7.67 × 10−8	7.29 × 10−16	0.019	2.34 × 10−7	2.06 × 10−5	9.41 × 10−3
Chinese
MAF	0.946	0.419	0.429	0.267	0.303	0.295
Beta	0.129	0.101	−0.128	0.051	−0.052	0.065
SE	0.032	0.015	0.015	0.017	0.016	0.016
p value	5.66 × 10−5	9.84 × 10−12	2.60 × 10−18	0.002	0.001	8.04 × 10−5
Meta-analysis for Rater-Scored Cohorts
p value	6.65 × 10−10	1.13 × 10−28	5.27 × 10−13	2.49 × 10−14	2.44 × 10−8	4.82 × 10−9
23andMe
MAF	0.007	0.410	0.249	0.238	0.451	0.389
OR	1.490	1.270	0.837	1.335	0.813	1.237
CI	(1.299, 1.709)	(1.238, 1.303)	(0.813, 0.862)	(1.296, 1.375)	(0.793, 0.834)	(1.206, 1.269)
p value	2.04 × 10−8	4.90 × 10−76	4.42 × 10−33	3.31 × 10−87	8.48 × 10−59	8.26 × 10−59
Meta-analysis for All Cohorts
p value	1.16 × 10−13	1.65 × 10−99	4.78 × 10−43	7.64 × 10−100	1.56 × 10−65	2.85 × 10−66

Abbreviations are as follows: MAF, minor allele frequency; SE, standard error of the beta-coefficient; OR, odds ratio; and CI, 95% confidence interval of the odds ratio.

Figure 2

Regional Association Plots Showing Significant Associations Observed in the Meta-analysis of the Three Rater-Scored Cohorts

Regional plots near (A) EDAR, (B) SP5, (C) MRPS22, (D) ADGRG6 (GPR126), (E) KIAA1217, and (F) PAX9 show –log10-transformed p values (left y axis) by physical position (x axis). Shading denotes the LD (r2) between each SNP and the lead SNP (purple). The blue overlay represents the recombination rate (right y axis). Gene positions are indicated under each plot.

Each of these six loci also showed significant evidence of association with earlobe attachment in the self-rated 23andMe cohort (Table 1). Moreover, in the meta-analysis across all four cohorts (Figures 1C and 1D), a total of 49 significant loci were observed (Figure S2; Table S5), which included the six loci observed in the meta-analysis of expert-rater-scored cohorts. Of the 49 loci, 15 showed significant associations (i.e., p < 5 × 10−8) in the 23andMe sample and replication-level p values of <0.001 in at least one additional cohort or in the meta-analysis across the three expert-rater-scored cohorts. Another 24 loci were driven primarily by significant associations observed in the 23andMe cohort; 14 of the 24 showed consistent directions of effects between 23andMe and all three expert-rater-scored cohorts (sign test p = 1.29 × 10−7), and 23 of 24 showed consistent directions of effects between 23andMe and at least two of the three expert-rater-scored cohorts (sign test p = 1.49 × 10−6). The remaining 4 of the 49 loci were significant only in the meta-analysis across all four cohorts but not in any individual cohort.

351 genes were located within 500 kb of the lead SNP across the 49 associated loci. Using public databases (see Subjects and Methods), we queried these 351 genes for documented expression and activity of enhancer elements in relevant tissues, as well as any known roles associated with ear phenotypes in human disorders or mouse models. Table S6 enumerates the genes at associated loci and the evidence substantiating their biological roles related to ear morphology. In total, 71 (20%) of the 351 genes were expressed in relevant tissues in mice (12 [3%] in the outer ear, 58 [17%] in the inner ear, and 39 [11%] in the branchial arches). Likewise, 21 (6%) of the genes were implicated in human syndromes manifesting with ear phenotypes, 22 (6%) were implicated in ear phenotypes in mouse models, and 16 (5%) were flanking active enhancer elements in relevant tissues. Overall, several plausible candidate genes were identified across the 49 associated loci identified via meta-analysis.

We performed genomic enrichment analysis by using GREAT to determine whether the gene set comprising the two nearest genes across the 49 associated loci (indicated by the position of the lead SNP) was enriched with relevant annotations across several ontologies. Of most relevance, we observed significant enrichment of over 16-fold for the human phenotype annotations “microtia,” “aplasia/hypoplasia of the external ear,” and “aplasia/hypoplasia of the ear” (p values < 1.5 × 10−5 for all), as well as 5-fold enrichment for “low-set ears” and “abnormal location of ears” (p values < 0.0001 for both). Significant enrichment was also observed for several embryonic-development- and morphogenesis-related Gene Ontology biological processes, many mouse morphology (including ear) terms, and mouse expression in many relevant tissues (notably branchial arch and ear). Detailed enrichment results are shown in Figure S3.

Expression Experiments

To confirm that genes located among associated loci are expressed in relevant tissues during development, we used RNA-seq to measure the expression of 174 genes located within 250 kb of one of the 49 lead SNPs. Gene expression was measured in two human fetal ears (at days 57 and 59 of development, when external ear structures are present but still developing) and in mouse embryonic day 10.5 branchial arch tissue isolated from two mutants, short body-short ear (sbse) and dumbo (dmbo), and sex- and background-matched wild-type (C57BL/6) controls. These mutants present with microtia (sbse) or low-set ears with “lobe duplication” (dmbo). For 4 of the 49 loci, no genes were located within 250 kb of the lead SNP. The majority of the remaining 45 loci had one or more genes that were robustly expressed in these relevant embryonic tissues (Figure S4). Human fetal ear tissue showed similar expression levels at days 57 and 59, and the greatest expression was observed on both days for PRRX1 (MIM: 167420; 1q24.2), a homeobox gene with relevant biology (see Discussion). Some of the genes located at the six loci (especially 2q31.1 and 2q13; Figure 3) identified in the GWAS of rater-scored cohorts and recapitulated in the 23andMe cohort were among the top ranking genes (of the 174 genes tested) in terms of expression in humans or differential expression between mutant and wild-type mice. For example, mouse orthologs of both SP5 (2q31.1) and EDAR (2q13) exhibited higher expression in dmbo mutant mice than in wild-type controls (log2 fold change > 0.5), and orthologs of GAD1 (MIM: 605363; 2q31.1), MYO3B (MIM: 610040; 2q31.1), and SH3RF3 (2q13) showed reduced expression in the mutant mice (log2 fold change < −0.5). Thus, separate genes at the 2q13 and 2q31.1 loci were both up- and downregulated (see Figures 2A and 2B for genetic association results for these loci). In a comparison of sbse mutants and wild-type mice, orthologs of ERICH2 (2q31.1) and NKX2-1 (2q13.1; MIM: 600635) were among the top overexpressed genes in sbse mutants among the 174 tested, and orthologs of GAD1 (2q31.1) and SH3RF3 (2q13) again exhibited reduced expression in the mutants (see Figures 2A, 2B, and 2F for genetic association results for these genes). Together, these results confirm the expression of genes located at associated loci in relevant human tissue and suggest that expression differs by ear phenotype for some genes at associated loci in a mouse model.

Figure 3

Heatmaps of Gene Expression

Heatmaps of gene expression in second branchial arch tissue from embryonic mice (wild-type mice and dmbo and sbse mutants) and fetal human pinna (at days 57 and 59) and for genes at the six loci observed in the meta-analysis of rater-scored cohorts; normalized total read counts are on the top, and fold changes between mutant mice and wild-type mice are on the bottom. The shading scale is shown for Z scores of expression and fold change. Genes for which expression data were not measured are shown in gray.

Discussion

We have performed the largest genome-wide study to date of earlobe attachment in 74,660 individuals from three ancestry groups. All six significant loci observed in the meta-analysis using the tripartite rater-scored phenotype were also significantly associated with earlobe attachment in the 23andMe cohort using the self-reported dichotomous phenotype. Furthermore, four of the six loci were the top ranking (by p value) associations observed in the 23andMe cohort (the other two loci ranked 7th and 42nd in the 23andMe GWAS). Thus, the GWAS of the 23andMe cohort recapitulated all of the major findings from the expert-rater-scored cohorts. Furthermore, four of the six loci contained genes showing notable expression (e.g., a high read count in humans and/or differential expression in mice) in relevant embryonic tissues.

Among the six associated loci from the rater-scored meta-analysis, two (2q13 at EDAR and 2q31.1 at SP5) have been previously identified in the Latin American cohort, two (6q24.2 near ADGRG6 [GPR126] and 10p12.2 at KIAA1217) have previously shown suggestive evidence of association,4, 5 and one (3q23 near MRPS22) has been previously implicated in earlobe size. In contrast, the association at 14q13.1 near PAX9 did not show evidence of genetic association in previous studies.4, 5 These or nearby genes have known biological functions that indicate plausible roles in determining ear morphology. For example, EDAR encodes a cell-surface receptor important for the development of ectodermal tissues, including skin. The missense SNP rs3827760 (c.1109T>C [p.Val370Ala]) affects protein activity31, 32 and is associated with variation in tooth morphology, hair, sweat gland density, and facial morphology in Asians.33, 34, 35, 36, 37, 38, 39, 40, 41 This variant was the top SNP in the EDAR region in the meta-analyses, and its frequency differs dramatically across populations (e.g., G allele frequency of <1% in Europeans, 39% in Latin Americans, and >90% in Han Chinese in the 1000 Genomes Project). Adhikari et al. have shown that Edar is expressed in the mouse pinna and that, compared with wild-type mice, mouse mutants with loss of Edar function exhibit reduced ear protrusion and length, as well as a different shape.

SP5 encodes a transcription factor involved in the regulation of Wnt-mediated beta catenin signaling, which in turn is critical for multiple aspects of development, including that of the inner ear. ADGRG6 (GPR126) encodes a G-protein-coupled receptor whose disruption, via either mutation or morpholino, causes a swollen inner-ear phenotype in zebrafish. This locus was shown to be associated with earlobe size, a related phenotype, in the Latin American cohort. KIAA1217 is not known to be involved in ear development, but the associated variants are downstream of ARHGAP21 (MIM: 609870), variants in which have been associated with mandibular prognathism, a branchial arch defect.

The association at 3q23 occurs nearest MRPS22, which is implicated in a Cornelia de Lange-like phenotype including ear and skin dysmorphic features, and is upstream of FOXL2 (MIM: 605597), encoding a craniofacial transcription factor. This locus was previously implicated in earlobe size in the Latin American cohort. The 14q13.3 association is near PAX9, encoding a transcription factor involved in mouth and tooth development, as well as NKX2-8 (MIM: 603245), a homeobox candidate gene for microtia.

In addition to the six loci identified via meta-analysis of rater-scored cohorts and recapitulated in 23andMe, meta-analysis across all four cohorts yielded 43 additional significant associations, driven primarily by the large 23andMe cohort. Among the more promising candidates were genes implicated in both human and mouse ear phenotypes: TBX15 (MIM: 604127), PRRX1, and ZEB2 (MIM: 605802). TBX15 (the gene containing the lead SNP at this locus) is a transcription factor responsible for Cousin syndrome (MIM: 260660), in which ears are low-set and posteriorly rotated. In mouse, mutations in TBX15 cause abnormal ear position and the “droopy ear” phenotype, which resembles that seen in the dmbo mutants used in the expression analysis in this study. This locus was associated with ear phenotypes including antitragus size and folding of antihelix in the Latin American sample. PRRX1 (40 kb upstream of the lead SNP, which also showed the greatest expression in fetal human ear at day 57) is a homeobox gene implicated in Agnathia-otocephaly (characterized by severe malformations of the mouth, jaw, and ear; MIM: 202650), and copy number variants affecting PRRX1 have been observed in several patients with dysmorphic ear phenotypes. Moreover, PRRX1 is expressed in the inner, middle, and outer ear and first and second branchial arches in mouse, and mutations cause lower ear position and abnormal Meckel’s cartilage. ZEB2 (500kb downstream of lead SNP) is a homeobox gene implicated in Mowat-Wilson syndrome (MIM: 235730), in which the ears are cupped and the earlobes are upturned with central depression. In mouse, ZEB2 is expressed in the inner and middle ear, and first brachial arch, and knockouts lead to Mowat-Wilson-like features or lack of first branchial arch during embryogenesis. Other notable candidates include the growth factor BMP5 (MIM: 112265), and homeobox transcription factors DLX5 (MIM: 600028) and DLX6 (MIM: 600030), which are all expressed in ear and related tissues and are all implicated in ear phenotypes in mice.55, 56

The myriad associations overwhelmingly demonstrate the polygenic nature of earlobe attachment, standing in contrast to previous notions regarding its Mendelian nature, which have been perpetuated through the primary literature and educational materials for nearly a century. In fact, the large number of confidently identified loci is on par with many continuous anthropometric traits such as height and body composition. Moreover, earlobe attachment is correlated with other aspects of lobe morphology, including earlobe size, so the overlap in associated loci with the previous study of lobe size is unsurprising. The effect sizes of variants observed in this study are fairly large for individual SNPs (resulting in a difference of up to 0.2 phenotype standard deviations per allele in the rater-scored cohorts and up to an odds ratio of 1.5 in the 23andMe cohort), although they are not large enough to cause Mendelian segregation. Moreover, differences in allele frequencies across ancestry groups were observed for some associated SNPs, which could explain part of the ethnic heterogeneity observed for earlobe attachment. Specifically, 16 of the lead SNPs of the 49 associated loci, notably the EDAR variant rs3827760 (which had a MAF difference of 0.93 between the European American and Chinese cohorts), showed MAF differences greater than 0.2 across the ancestry groups.

Consistent with a polygenic trait in a well-powered GWAS, we observed evidence of genomic inflation (e.g., genomic inflation factors greater than 1.0) separately in the meta-analyses and the 23andMe cohort. This occurs because the genomic inflation factor, although designed to be calculated from a set of independent null markers, is instead calculated in GWASs from the set of all of the SNPs tested, including truly associated SNPs in LD with causal variants. For polygenic traits (for which there are multiple truly associated loci) in well-powered studies (where small p values are obtained even for SNPs weakly correlated with true causal alleles), the lambda is expected to be greater than 1. We argue that this inflation is not due to population stratification because association models were adjusted for genetic ancestry estimated from the genetic data in each cohort. Moreover, inflation was not observed in the individual expert-rater-scored cohorts or in a subset of 2,000 participants from the 23andMe cohort, as would be expected if population stratification had caused epidemiological confounding. Instead, the inflation observed for earlobe attachment, which is similar to that observed for highly polygenic traits such as height, is expected given the contributions of numerous associated loci each tagged by many correlated SNPs.

Strengths of this work include the high-quality phenotyping based on digital imagery in the rater-scored cohorts, the inclusion of cohorts from different ancestry groups, the large sample size, and the method of meta-analysis, which was chosen to be robust to phenotype differences across the cohorts. Although phenotype data were collected via self-report in the 23andMe cohort, the fact that associations observed in the rater-scored cohorts were also identified in the 23andMe cohort suggests that the large sample size of the 23andMe cohort counterbalances noise (if any) as a result of the method of data collection. Despite these strengths, and because of the differences in phenotype assessment across cohorts, this study was limited by the fact that, within our testing framework, we were not able to directly test the heterogeneity of effects among cohorts. In addition, the large sample size of the 23andMe cohort, which benefited the statistical power of study, most likely had an outsized effect on the meta-analysis across all four cohorts. For this reason, we also reported results for meta-analysis across the three expert-rater-scored cohorts.

In conclusion, we have identified 49 associations with earlobe attachment, including 21 loci meeting the standard of genome-wide discovery (p < 5 × 10−8) plus independent replication (p < 0.001) and 28 loci showing evidence of discovery only (i.e., without independent replication), via meta-analysis. These genes provide insight into the complex biology of ear development. The fact that we observed several associated genes in which pathogenic variants are known to cause human syndromes with ear phenotypes is consistent with our hypothesis that whereas deleterious variants in genes can cause congenital defects and Mendelian conditions, regulatory variants in the same genes can influence normal phenotypic variation. Ultimately, understanding the genetics of normal human morphological traits can provide insights into the genes and pathways involved in developmental malformations.

Consortia

Members of the 23andMe Research Team include Michelle Agee, Babak Alipanahi, Adam Auton, Robert K. Bell, Katarzyna Bryc, Sarah L. Elson, Pierre Fontanillas, Nicholas A. Furlotte, David A. Hinds, Bethann S. Hromatka, Karen E. Huber, Aaron Kleinman, Nadia K. Litterman, Matthew H. McIntyre, Joanna L. Mountain, Elizabeth S. Noblin, Carrie A.M. Northover, Steven J. Pitts, J. Fah Sathirapongsasuti, Olga V. Sazonova, Janie F. Shelton, Suyash Shringarpure, Chao Tian, Joyce Y. Tung, Vladimir Vacic, and Catherine H. Wilson.

Conflicts of Interest

D.H. is employed by and owns stock or stock options in 23andMe.