Introduction

Age-related hearing impairment (ARHI), or presbycusis, is the most frequent sensory disability among the elderly.1, 2 The involvement of genetic factors in the development of this multifactorial disorder was demonstrated by several heritability studies.3, 4, 5 Only recently have the first studies aiming at the identification of these genetic factors been conducted. Hitherto, three linkage studies have been performed,6, 7, 8 but only one region on chromosome 8q24.13–q24.22 reached genome-wide significance.8 In addition, several association studies based on a candidate gene approach and two genome-wide association studies (GWAS)8, 9 have been conducted. Early candidate gene-based association studies are reviewed by Van Eyken et al.10 Recently, based on a candidate gene approach, we reported associations between ARHI and variants in GRHL2,11 encoding a transcription factor that is also involved in monogenic hearing loss. In addition, based on a pooling GWAS, we reported association with variants in the glutamate receptor-7 gene (GRM7),9 encoding a metabotropic glutamate receptor assumed to modulate hair cell excitability and synaptic efficacy. Acquired mitochondrial mutations have also been proposed as a cause of ARHI in humans. These include both the ‘common’ 4977-bp mitochondrial DNA deletion associated with ageing12 and accumulation of many different mitochondrial mutations.13

Despite these first series of studies, the genetic aetiology of ARHI remains elusive. In this paper, we present the results of a GWAS for ARHI in the Finnish Saami population. The Saami, previously known as Lapps, live across Northern Scandinavia and the Kola Peninsula. Traditionally, Saami have been fishers, hunters and reindeer herders. A minority of the Saami continue to be engaged in these traditional economic activities, whereas the majority has adopted other occupations and a more westernised lifestyle. Among European populations, the Saami are considered a population genetic ‘outlier’. Despite the fact that both the Saami and Finnish languages belong to the Uralic language group, there is strong genetic differentiation between the Saami and Finns.14 Several genetic studies have been conducted to shed light on the genetic origin of the Saami.15, 16, 17, 18 These studies suggested that the Saami are descendants of Europeans who, several thousand years ago, reached northern Fennoscandia from both a western and an eastern route and subsequently got admixed. They also found a small East Asian contribution to the gene pool.18, 19 Several genetic studies have further suggested that the extent of linkage disequilibrium (LD) in the Saami is dramatically higher compared to that of outbred populations.20, 21, 22, 23 Higher levels of LD are observed in other population isolates as well.24, 25 Most of the genetic studies of genetically isolated populations were conducted involving young founder populations, in which the high extent of LD can be attributed to a founder effect. However, in ancient, small, genetically isolated populations without evidence of demographic expansion, high levels of LD are believed to have been generated by genetic drift.26, 27 The Saami represent such a population. Therefore, the choice of this study population was based on its anticipated statistical power advantage for association mapping, offering better genomic coverage for a given set of markers, as compared to outbred populations.

Owing to the expected higher levels of LD in the Saami population and based on what was known about the structure of the human genome at the time this study was designed, we anticipated that the Affymetrix 100 K array pair would entail sufficient genomic coverage for identification of putative ARHI-susceptibility genes.

Materials and methods

Subjects

Saami subjects, aged between 50 and 75 years, were enrolled through the public population register from four municipalities in Northern Finland (Ivalo, Utsjoki, Vuotso and Enontekiö) in three stages. In the first stage, geographical criteria were applied; only subjects from an area with a high probability of finding Saami individuals were selected. In the second stage, putative study subjects were invited based on an expert evaluation of Saami communities. Finally, the Saami identity of the subject in question was confirmed through a discussion involving the subject and the researcher. Except for the age range that was used, the ARHI selection criteria were as described previously.28 Briefly, all volunteers completed an extended questionnaire detailing medical history and exposure to environmental factors putatively influencing hearing. Subjects suffering from one or more medical pathologies from a standardised list of conditions that could potentially affect hearing ability were excluded. All subjects who passed this initial screening underwent an otoscopic investigation before pure-tone audiometry. Air conduction thresholds were measured at 0.125, 0.25, 0.5, 1, 2, 3, 4, 6 and 8 kHz, and bone conduction at 0.5, 1, 2 and 4 kHz for all participating volunteers. Audiological exclusion criteria were an air-bone gap of more than 15 dB averaged over 0.5, 1 and 2 kHz in at least one ear or asymmetrical hearing impairment with a difference in air conduction thresholds exceeding 20 dB in at least two frequencies between 0.5, 1 and 2 kHz. Written informed consent was obtained from all volunteers. The study was approved by the Finnish National Advisory Board on Health Care Ethics and by the ethical committees or the appropriate local institutional review boards at all participating institutions.

Description of the phenotype

For each subject, only hearing thresholds from the best ear, defined as the ear with the lowest average threshold over all measured frequencies, were used in the analyses. Supplementary Figure 1 shows that hearing thresholds increase with age and that, for the higher frequencies, thresholds for men are higher and display greater variability (i.e., there is heteroscedasticity) than those for females. We therefore adjusted hearing thresholds for age and sex. The effect of age was removed by regressing log-transformed threshold values (log10(20+threshold)) for each of the measured frequencies on age, age squared and age cubed. To resolve the heteroscedasticity problem and to correct for sex, this adjustment for age was carried out for males and females separately, after which the residuals of each regression were scaled and subsequently combined. Note that adding the age and sex covariates in the model for association testing does not deal with this problem satisfactorily.

To deal with the multivariate audiometric phenotype, we applied principal component (PC) analysis as a data-reduction technique. Performing association tests for each frequency separately would have resulted in a large number of tests, worsening the multiple-testing problem. In addition, such analysis would have completely ignored the correlation structure of the audiometric data, that is, the shape of the audiogram. The latter would also apply when simple pure-tone averages (PTAs), which are more widely used in hearing research, would have been analysed. Our first PC (PC1) was highly correlated with the PTA over all frequencies, adjusted for sex, age, age square and age cubed (Spearman correlation=0.96). Spearman correlations between PC1 and sex- and age-adjusted PTAs for the high (2–8 kHz) and low frequencies (0.125–1 kHz) were 0.88 and 0.61, respectively. The association testing results for PC1, therefore, would not differ much from those that would be obtained when PTAs would have been analysed.

Table 1 lists the eigenvalues and eigenvectors for the first three PCs, PC1 to PC3. The results of separate calculations for males and females yielded similar values, justifying data pooling (results not shown). Each of the PCs provides independent bits of information. The first three PCs were retained for the analysis as they captured approximately 80% of the total variation. The coefficients of the eigenvectors (loadings) are directly proportional to the correlations between the transformed hearing thresholds corrected for age and gender, and each of the PCs (see, for example, Johnson and Wichern29). PC1 is clearly a ‘size’ variable, providing an overall measure of a subject's hearing ability. PC2 and PC3 are ‘shape’ variables. PC2 contrasts low with high frequencies and is a measure of the slope of the audiogram. PC3 contrasts the middle frequencies with the lower and higher frequencies and can be considered a measure of the concavity of an audiogram (Table 1 and Figure 1). These phenotypes are biologically meaningful30, 31 and we have previously shown that they have significant heritabilities.8

Table 1 Eigenvectors and eigenvalues for the first three PCs calculated from age- and sex-corrected hearing thresholds
Full size table
Figure 1
figure 1

Comparison of typical audiograms between subjects with high and low values for PC1, PC2 and PC3. For each of the three traits, three subjects from the lower and upper 10% extremes of the distribution of the trait are shown. The audiograms for both ears are given, to illustrate the extent of intra-subject variability. The left and right ears correspond to the filled and open circles, respectively. For PC1, on average, hearing loss is much more pronounced in individuals from the higher extremes (top row) as compared with those from the lower extremes (second row from the top). Subjects with low values for PC2 (fourth row from the top) tend to have a much more sloping audiogram, compared with those with high values for PC2 (third row from the top). Subjects with high values for PC3 (second row from the bottom) tend to have a more concave audiogram, compared with those in the lower extreme for PC3 (bottom row). Note that the displayed audiometric data have not been adjusted for sex and age.

Full size image

Genotyping and quality control of data

Genomic DNA was extracted using standard protocols and diluted to 50 ng/μl. The quality of the genomic DNA was assessed by spectrophotometric analysis and gel electrophoresis. Each sample was genotyped with the Affymetrix GeneChip Human Mapping 100 K array pair (Santa Clara, CA, USA), as described in the GeneChip Human Mapping 100 K Assay Manual. Genotype calling was performed using the BRLMM algorithm of the Affymetrix GeneChip Genotyping Analysis Software (GTYPE) version 4.1. Quality control of data was performed based on the following criteria: a sample call rate of >94%, a marker call rate of >90%, a Hardy–Weinberg equilibrium P-value of >0.001 and a minor allele frequency of >5%. In total, 352 samples from the Saami subjects were genotyped with the Affymetrix 100 K array. One subject was excluded because of low sample call rate. The total genotyping rate in the remaining samples was 0.9931. One additional sample was removed owing to an unintentional duplication event, and three further subjects were excluded owing to missing phenotypic information. After data quality control, 347 subjects and 83 381 SNPs were used for statistical analysis. Genotyping for SNP rs11928865 was performed using a Taqman assay (Applied Biosystems, Foster City, CA, USA) with a Lightcycler480 real-time PCR machine (Roche, Basel, Switzerland).

Association testing

Estimates of genome-wide pairwise identity-by-descent sharing, provided by PLINK,32 revealed a substantial degree of undocumented relatedness among the subjects. This was also reflected in an excess of SNPs that were not in Hardy–Weinberg equilibrium. The results of a linear regression for each of the three quantitative traits would, therefore, not be trustworthy as dependent data lead to underestimation of the standard errors for the regression coefficients, yielding P-values that are too small. In addition, population stratification may confound the results of the analysis. To account for the dependency of the data, as well as for population stratification, we analysed the data using a mixed model.33 With this approach, the covariance matrix of the phenotype is modelled as:

where y is the column vector with the phenotypes of the 347 subjects, K is a matrix with kinship coefficients estimated based on the genome-wide data, I is an identity matrix and σg2 and σe2 are the genetic and the residual variances, respectively. This was implemented using the MIXED procedure of SAS 9.1.3 within SAS macro, which was based on that used by Yu et al.33 The 2K matrix with the proportion of alleles shared identical-by-descent was estimated using PLINK.32 The model also corrected for origin. As samples originated from four municipalities, three dummy variables were included as covariates in the model. An additive genetic model was assumed for each SNP. As it was computationally infeasible to run this analysis for the entire set of SNPs, for each trait only the top 4000 SNPs (approximately 5% of the total number of SNPs passing quality control) resulting from an initial screening based on linear regression analysis (naive analysis) were considered for the mixed model analysis. To verify that no true association signals were lost by using this analytical strategy, we first examined how well the P-values of the naive analysis correlated with those of the mixed model analysis. For this purpose, we analysed 8000 randomly selected SNPs meeting the quality control standards, representing approximately 10% of the SNPs that passed the quality control, with the mixed model. Comparison of P-values between naive analysis and mixed-model results showed that no SNPs with highly significant P-values in the mixed-model analysis had a non-significant P-value in the naive analysis. A reciprocal relationship would be attributed to inflated test statistics. Therefore, for the data at hand, by analysing the top 5% (4000 SNPs) of a naive linear regression analysis, the chance of missing any SNP that would not be significant in the naive analysis, but highly significant in the mixed-model analysis seems very low (Supplementary Figure 2). Quantile–quantile plots in Supplementary Figure 2 show that residual population stratification or cryptic relatedness is of minor concern. Genomic control inflation factors varied from 1.02 to 1.08. The graph on the right of each panel in Supplementary Figure 2 further shows that P-values obtained from the naive analysis tend to be too low. For example, for PC1, the most extreme −log10(P-value) from the naive analysis is 4, whereas it is 3.5 for the mixed model.

Results

After quality control, 347 subjects and 83 381 SNPs were used for association testing. The quantitative traits tested were the first three PCs calculated from the log-transformed hearing threshold values, adjusted for age and sex. To account for the multiple levels of relatedness present in our data and potential population stratification, the associations between SNPs and the traits were tested using a mixed-model approach, based upon genome-wide kinship estimates. Owing to computational constraints, only the top 4000 SNPs (approximately 5% of the total number of SNPs passing quality control) resulting from an initial screening based on a linear regression analysis were used in the mixed-model analysis. In the Materials and methods section, we have outlined that, using this strategy, it is unlikely that strong association signals will be lost. For each of the three PCs, the top 25 SNPs resulting from this analysis are listed in Tables 2, 3, 4.

Table 2 Top 25 SNPs resulting from the mixed-model analysis for PC1
Full size table
Table 3 Top 25 SNPs resulting from the mixed model analysis for PC2
Full size table
Table 4 Top 25 SNPs resulting from the mixed model analysis for PC3
Full size table

If, in agreement with other GWAS,34 we used a P-value cut-off of 5 × 10−7 to declare genome-wide significance, the top-ranked SNP for PC3 (rs457717; Table 4) survived this significance threshold (P-value=3.55 × 10−7). This SNP is located in an intron of the IQ motif-containing GTPase-activating-like protein (IQGAP2). The second-ranked SNP (rs1697845) is also positioned in IQGAP2, and is highly correlated with the first (D′=1; r2=0.965), indicating that this signal is not a genotyping artefact.

For the other two traits, PC1 and PC2, none of the SNPs reached this level of significance (Tables 2 and 3). Intriguingly, however, the eighth-ranked SNP for PC1 (rs161927; Table 2; P-value 0.000149) was positioned 40 kb downstream from the metabotropic GRM7 gene (Supplementary Figure 3). We have previously reported an association of a common variant in the GRM7 gene, rs11928865, with hearing ability based on a pooling-based GWAS of a European sample set.9 Therefore, we also genotyped this SNP in the Saami population and analysed these data using the mixed-model approach. This resulted in a P-value of 0.045. However, the A-allele was associated with an elevated hearing threshold in the Saami sample set, while in the European sample set, the associated allele was T.9 One of the two top-scoring GRM7 SNPs resulting from a GWAS of a Finnish sample set,9 rs779701, was present on the 100 K Affymetrix array. However, in the Saami sample set, this SNP did not result in a significant association (P-value of 0.735). The LD pattern within the GRM7 region deduced for the Saami (Supplementary Figure 3D) was similar to that observed in the HapMap for the European population (Supplementary Figure 3C). No obvious biological meaning related to the inner ear could be deduced for the gene regions of the seven higher-ranked SNPs for PC1.

For PC2, no strong biological candidates were present among the top 25 SNPs (Table 3). The presence of six SNPs from the same gene (KH domain-containing, RNA-binding, signal transduction-associated protein-2; KHDRBS2) in Table 3 is due to the fact that all of these SNPs reside in the same region of high LD.

Discussion

We conducted a genome-wide association scan for ARHI traits derived from audiometric data in the genetically isolated Finnish Saami population. We summarised the information contained in each audiogram by the first three components obtained through PC analysis. Previously, we have shown that these traits have significant heritabilities.8 Using these audiometric measures, we retain biologically important information on the shape of the audiogram that would be lost when simple PTAs or individual frequencies are analysed. Indeed, premortem-derived audiometric patterns have been correlated to pathophysiological findings in human temporal bones.30, 31 Nelson et al35 studied individuals with downward-sloping audiograms and reported that the severity of hearing loss, based on audiometric thresholds, was highly associated with a decrease in stria vascularis volume and degeneration of outer and inner hair cell and ganglion cell populations. The slope of the audiogram was significantly associated with the extent of ganglion cell degeneration.

Our association scan resulted in two SNPs in or near genes that are reasonable biological candidates and that are candidates for further research. The top-ranked SNP, rs457717 (P-value 3.55 × 10−7), that reached genome-wide significance, was associated with PC3 and was localised in an intron of the IQ motif-containing, GTPase-activating-like protein (IQGAP2). IQGAP2 is a member of the Ras superfamily of GTPases. This class of proteins regulates a wide variety of cellular signalling pathways. IQGAP2 is expressed inside the cochlea36 and has been implicated in cadherin-mediated cell adhesion.37 As cadherin-23 (CDH23) plays an important role in the sensory hair cells,38 and mutations in CDH23 are responsible for syndromic and non-syndromic autosomal recessive deafness in humans39, 40 and age-related hearing loss in inbred mouse strains,41 this warrants further research. Intriguingly, the SNP rs161927 (P-value 0.000149), that ranked eighth for PC1, was positioned immediately downstream from the metabotropic GRM7 gene. A previous GWAS of a European and a Finnish sample set, which was based on an age- and gender-corrected Z-score highly correlated with our PC1,9, 11, 28, 42 already suggested a role for GRM7 in ARHI.9 GRM7 encodes a metabotropic glutamate receptor that can be activated by L-glutamate, the primary excitatory neurotransmitter in the cochlear hair cell. GRM7 expression was observed in the inner and outer hair cells, and in the spiral ganglion nerve cell bodies of the inner ear.9 We have hypothesised that certain GRM7 variants may cause increased susceptibility to glutamate excitotoxicity, which may consequently lead to the development of ARHI. However, variants that were significant in our previous study did not replicate in our current study. The LD structure within the GRM7 region indicates that these are independent signals. These results, therefore, are consistent with allelic heterogeneity, which is a plausible scenario given the likely absence of any selective pressure for this disorder. Independent studies are needed to investigate the role of common variants in GRM7 in ARHI.

This study was conducted involving the genetically isolated Saami population. A number of studies provide evidence that the extent of LD is much higher in the Saami as compared to that in other European populations.20, 21, 22, 23 This would provide an important power advantage, as higher LD leads to better genomic coverage. It has been proposed that association studies based on a map of modest marker density within a population isolate such as the Saami would provide a first rough identification of the region associated with a particular complex disease. This first demarcation of the region should then preferably be refined by association studies of an outbred population using a denser map.22 We have published elsewhere the results of an elaborate genome-wide evaluation of the relative power for GWAS in the Saami.43

A possible limitation of this study is the limited sample size. Recent studies indicate that several thousand subjects need to be recruited in order to obtain sufficient power to detect the effect sizes typically expected for complex phenotypes (with each contributing variant explaining as little as 0.1–1.3% of the variation),44, 45 although successful GWAS have been conducted using only 96 cases and 50 controls.46 Obviously, large sample sizes can never be recruited in small isolated populations such as the Saami. A specific difficulty of working with isolated populations became apparent in the current study. Owing to the fact that isolated populations are often derived from a limited number of founders and that inbreeding is sometimes inevitable, an increased relatedness is observed within such populations. As a consequence, classical statistical methods that assume unrelatedness, cannot be applied to analyse an association study performed in a population isolate. In our exploratory phase, we investigated several approaches, including a naive linear regression analysis and a PC-based correction (EIGENSTRAT47). However, for the latter approach, the basic assumptions were violated and the inbreeding and relatedness present in our Saami sample created spurious axes of variation.48 Mixed models offered a solution. In this approach, we estimated the genetic relationships between the study subjects based on genome-wide data and subsequently accounted for them in the analysis.33

In conclusion, we present the results of a first GWAS performed in the genetically isolated Saami population. We identified IQGAP2 as a candidate ARHI gene. Independent replication is needed to confirm this lead. In addition, we found further evidence implicating GRM7 in ARHI.