Giorgia Girotto

Nicola Pirastu

Rossella Sorice

Ginevra Biino

Harry Campbell

Adamo P D'adamo

Nicholas D Hastie

Teresa Nutile

Ozren Polasek

Laura Portas

Igor Rudan

Sheila Ulivi

Tatijana Zemunik

Alan F Wright

Marina Ciullo

Caroline Hayward

Mario Pirastu

MD Paolo Gasparini

Hearing function and thresholds: a genome-wide association study in European isolated populations identifies new loci and pathways

Keywords: Molecular basis of hearing system, GWAS, meta-analysis, isolated populations, "in silico" pathways

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The hearing system is characterized by three structures: a) the outer part, b) the middle ear, and c) the inner ear, that all play a role in hearing function. The hearing system is difficult to study through biochemical routes, due to the small amounts of tissue available for analysis and by key molecules that may be present in only a few tens of copies per cell, thus compounding the difficulty [START_REF] Steel | A genetic approach to understanding auditory function[END_REF] . Hair cells of the inner ear are constituted by stereocilia arranged in bundles at their upper surface, packed with actin filaments, and deflected by the vibrations of sound [START_REF] Steel | A genetic approach to understanding auditory function[END_REF] . This activity opens ion channels modulating potential within the cell. Cell activation releases neurotransmitters to synaptic junctions between hair cells and neural fibers of the auditory nerve. The neural spike subsequently propagates in the auditory nerve fiber. Nerve impulses are finally perceived by the brain, primarily in the temporal lobe, where they can be processed and assigned meaning.

Little is known about the molecular basis of variation of normal hearing function.

Several molecules have been identified as having a role in auditory function and hair cell transduction because they are specifically expressed in or around the stereocilia and mutations in their genes lead to hearing impairments in either humans or mice models [START_REF] Quint | Use of mouse genetics for studying inner ear development[END_REF] . These dysfunctional proteins are involved in impaired molecular-physiologic processes of potassium and calcium homeostasis, apoptotic signaling, stereocilia linkage, mechano-electric transduction, electromotility, and many other processes [START_REF] Van Laer | Nonsyndromic hearing loss[END_REF] . Briefly, these molecules include myosins which represent one of the largest group of deafness associated molecules, adhesion protein such as cadherins, member of the ferlin family, components of the tectorial membrane, genes involved in ion homeostasis such as connexins [START_REF] Martínez | Gap-junction channels dysfunction in deafness and hearing loss[END_REF] and many others [START_REF] Hilgert | Function and expression pattern of nonsyndromic deafness genes[END_REF] . According to the Hereditary Hearing Loss homepage (see HHL), more than 140 loci for nonsyndromic HHL have been mapped, and approximately 80 genes identified in humans.

In animal models ~70 loci for nonsyndromic HHL and ~60 genes have been so far described.

Hearing loss can also be multifactorial or complex in causality such as age-related hearing impairment (ARHI) [START_REF] Reiss | Presbyacusis: pathogenesis and treatment[END_REF] and noise-induced hearing loss (NIHL) [START_REF] Konings | Genetic studies on noise-induced hearing loss: a review[END_REF] , reflecting the interaction of a number of genetic and environmental factors. Despite some relevant efforts done to identify the molecular bases of these conditions, up to now, only a few genes have been associated with both ARHI [START_REF] Reiss | Presbyacusis: pathogenesis and treatment[END_REF] and NIHL [START_REF] Konings | Genetic studies on noise-induced hearing loss: a review[END_REF] . Almost all genes so far identified in all forms of hearing loss are those directly related to hearing as qualitative traits (i.e. disease genes mainly involved in monogenic inherited forms), while almost nothing is known about genes implicated in defining hearing as quantitative trait (i.e. thresholds) in normal hearing.

The use of isolated populations to reduce heterogeneity of complex and/or quantitative traits has already proven very useful in identifying DNA polymorphisms associated with these traits [START_REF] Varilo | Isolates and their potential use in complex gene mapping efforts[END_REF] .

Nevertheless, the argument about the advantage of the use of such population is still an open issue [START_REF] Shifman | The value of isolated populations[END_REF] . In principle, the inbreeding, typical of small communities, reduces genetic heterogeneity, increases homozygosity, and reduce environmental factors providing greater power for detection of susceptibility genes [START_REF] Varilo | Isolates and their potential use in complex gene mapping efforts[END_REF][START_REF] Lowe | Genome-wide association studies in an isolated founder population from the Pacific Island of Kosrae[END_REF] . Moreover, such populations could be extremely useful to detect rare variants.

Here, we combine the power of studying isolated populations with a) a general screening of hearing function, b) high-throughput SNP analysis, c) Genome Wide Association Studies (GWAS), and d) up-to-date in silico pathways construction, to analyze hearing as a continuous trait and detect underlying genes and networks. To reach this goal, we analyzed different hearing thresholds as well as Pure Tone Averages (PTA) and principal components from a Principal Component (PC) analysis of hearing traits performing a large meta-analysis on isolated populations coming from different European regions. Results have been then used to construct a highly significant biological pathway.

Materials and Methods

Subjects

Within the International consortium called G-EAR, we recruited 3417 subjects coming from several isolated villages: Carlantino located in South Eastern Italy (267) [START_REF] Bedin | Age-related hearing loss in four Italian genetic isolates: an epidemiological study[END_REF] , Friuli Venezia Giulia Genetic Park [START_REF] Nelis | Genetic structure of Europeans: a view from the North-East[END_REF] , characterized by 6 villages located in North Eastern Italy (968), Korcula [START_REF] Polasek | Genome-wide Association Study of Anthropometric Traits in Korčula Island[END_REF] , an island in the Adriatic sea (Croatia) (795), Campora and Cardile two isolated villages located in the Cilento National's park , characterized by 2 different isolated villages in South Western Italy (421) [START_REF] Bedin | Age-related hearing loss in four Italian genetic isolates: an epidemiological study[END_REF] , and, Talana (469), an isolated village from Ogliastra Genetic Park in central part of Sardinia, Italy [START_REF] Fraumene | High resolution analysis and phylogenetic network construction using complete mtDNA sequences in sardinian genetic isolates[END_REF] , plus Split (497) an outbred not isolated population located on the Dalmatian coast (Croatia). All tests were performed using standard audiometers. Subjects underwent pure-tone audiometry, tympanogram, and acoustic reflex testing in both ears. Measurements were all obtained after any acoustically obstructing wax had been removed. The analysis of hearing function was done calculating the pure-tone average of airconduction (PTA at the lower 0.25, 0.5 and 1 kHz, medium 0.5,1 and 2 kHz, and high frequencies 4,8 kHz).

A questionnaire to obtain socio-demographic information, as well as data on physical activity (i.e. job, sport, etc.), lifestyle (e.g smoking, alcohol consumption, coffee intake, diet including taste and food preferences, etc.), clinical examinations (psicological, neurological, cardiological, etc.) clinical chemistry including blood count and more than 20 parameters, drugs, diseases and other information regarding the health status (BMI, bone density, blood pressure, etc.) have been collected for each subject. Only subjects aged 18 or older were included in the analysis. Clear familial forms of severe hearing loss have been excluded from the study.

Phenotypes

For the analyses we used only the threshold from the best ear, defined as the ear with lowest value of hearing loss for each individual. Since age explains a big part of the hearing variance, for this reason each trait was first linearly regressed against age. Residuals from this regression were cleaned from outliers (mean +/-6 SD) and, since their distributions were skewed, normalized using rank normal transformation and used for the association analysis using sex as a covariate. The following quantitative traits have been tested:

• Pure Tone Average at low, medium and high frequencies (PTAL, PTAM, PTAH)

• Seven different thresholds (250 Hz, 500 Hz, 1 KHz, 2 KHz, 4 KHz, and 8 KHz)

• The first 3 components from PCs (PC1, PC2, PC3) estimated on all of the frequencies as previously described [START_REF] Van Laer | A genome-wide association study for age-related hearing impairment in the Saami[END_REF] . We did not regress out sex or age before running PCA analysis and we checked that none of the components used reflected either one. Each component describes a different pattern in the data [START_REF] Van Laer | A genome-wide association study for age-related hearing impairment in the Saami[END_REF] : PC1 is a "size variable" that represents an overall measure of a subject's hearing ability, PC2 and PC3 are instead "shape variables". In particular PC2 shows the ratio of hearing between the high and the low 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.

DNA sampling and genotyping

All studies had appropriate ethical consent and consent forms for clinical and genetic studies have been signed by each participant in the study. Blood samples were collected and used to extract DNA using standard protocols. After measuring quantity and quality of DNA, samples were genotyped with Illumina 370k platform (Carlantino, FVG Genetic Park, Cilento, Korcula and Split) or Affymetrix 500K (Talana). Genotype quality control and data cleaning were conducted independently by each study group and are summarized in supplementary Table 2.

Genotypes were then imputed to the 2.5M HapMap CEU SNP set v22 (summary of imputation can be found in Supplementary Material). The same map was also used to look at regions of Linkage Disequilibrium.

Statistical analysis

Association analysis was carried out through a mixed model linear regression where the variance/covariance matrix is the genomic kinship. The analysis was implemented in GenABEL [START_REF] Aulchenko | GenABEL: an R library for genome-wide association analysis[END_REF] package for genotyped SNPs and ProbABEL 17 for imputed data. Meta-analysis was conducted using the inverse variance model as implemented in the MetABEL 16 R library. For PCs traits fixed effects meta-analysis was conducted where Z-scores were estimated from p-values and weighted on the sample-size as implemented in METAL (see METAL). SNPs with imputation quality (Rsq in MACH) less than 0.3 or with less than 30 copies of an allele in each population, were excluded. P-values were estimated using the Wald test. After quality control, 3417 subjects and approximately 2.2 million SNPs were used for meta-analysis.

Pathway analysis

From each GWA result all SNPs with P-values lower than 1x10 -4 have been selected. After removing duplicate loci, tagged by different SNPs, and keeping only the most significant ones, all genes that could be in linkage disequilibrium (LD) with these SNPs based on the HapMap CEU population have been taken into account, generating a list of 1276 genes. Identification of molecular network interactions and pathway analysis of most significant GWAS loci was completed using the Ingenuity Pathway Analysis (IPA) [START_REF] Sajan | Toward a Systems Biology of Mouse Inner Ear Organogenesis: Gene Expression Pathways, Patterns and Network Analysis[END_REF] tools from Ingenuity Systems® (Redwood City, California, USA; http://www.ingenuity.com). Briefly, the top 48 genes (cutoff of 5x10-5) arising from GWAS data were functionally characterized using IPA. Networks with a maximum of 70 genes or proteins were constructed, and scores were computed based on the likelihood of the genes being connected together due to random chance. A score of 2 indicates that there is a 1/100 chance that these genes are connected in a network due to random chance. Therefore, any networks with a score of 2 or above are considered statistically significant (with >99% confidence). IPA then map these genes to a global molecular network developed from information contained in the Ingenuity knowledge base (a manually curated database of experimentally proven molecular interactions from published literature). IPA determines the most significantly enriched biological functions and/or related diseases by calculating the P-value using Fisher's exact test. Using similar methods, significantly represented canonical pathways in a set of focus genes were also determined using IPA. To run the analysis only direct interactions were taken into account. After the definition of the most significant pathways, expression data within the ear have been obtained by searching the following databases: NCBI Gene for human and mouse, Eurexpress and Jackson Laboratory for mouse.

Results

The meta-analysis has succeeded in identifying some genome-wide suggestive loci associated with hearing traits, plus several additional ones strongly suggestive. All these trait-locus associations represent novel findings.

For all quantitative traits the top SNPs resulting from this analysis are listed in Table 1 (PTA), Table 2 (PCs), Table 3 andTable 4 (thresholds). Eight loci were strongly associated with the analyzed traits (max p=2.1x10 -7 ), while many others were in a group showing a maximum p-value of 1x10 -6 . Table 1 reports the most suggestive significant GWAS data obtained on PTA. For PTAL the strongest association signal was given by rs248626 (P-value 3.1x10 -6 ) located on chromosome 5, in a region containing DIAPH1 and a cluster of protochaderin genes [START_REF] Mansour | Hearing loss in a mouse model of Muenke syndrome[END_REF][START_REF] El-Amraoui | Cadherins as Targets for Genetic Diseases[END_REF] . Suggestive significant association with this locus has been also obtained at 500Hz threshold. Another interesting SNP is rs4603971 (chromosome 3) that is in LD with KCNMB2 a potassium large conductance calcium-activated channel [START_REF] Gabashvili | Ion Channel Gene Expression in the Inner Ear[END_REF] . It is known that several potassium channels are essential for hearing pathways [START_REF] Martínez | Gap-junction channels dysfunction in deafness and hearing loss[END_REF] . Moreover, KCNMB2 interacts with KCNMA1 that plays a key role in controlling the tuning of hair cells in the cochlea, regulation of transmitter release and innate immunity [START_REF] Wallner | Molecular basis of fast inactivation in voltage and Ca2+-activated K+ channels: a transmembrane beta-subunit homolog[END_REF] .

For PTAM, rs898967 (chromosome 16) is located within CMIP gene, which is expressed in the ear (see NCBI). Suggestive significant associations with this gene have been also found at 250 Hz, 1 KHz and at 2 KHz. Another positive SNP is rs641113 located on chromosome 10 is in LD with OPTN, an important gene required for myosin VI localization at the Golgi complex [START_REF] Sahlender | Optineurin links myosin VI to the Golgi complex and is involved in Golgi organization and exocytosis[END_REF] .

Myosin VI is essential for auditory and vestibular function in mammals and genetic mutations lead to hearing impairment and vestibular dysfunction in both humans and mice [START_REF] Sahlender | Optineurin links myosin VI to the Golgi complex and is involved in Golgi organization and exocytosis[END_REF] . Finally, for PTAH, there are two interesting SNPs: rs6673959 (chromosome 1) that is in LD with the DFFB gene, a pro-apoptotic related gene expressed in the cochlear cells [START_REF] Hu | Differential expression of apoptosis-related genes in the cochlea of noise-exposed rats[END_REF] and rs10936160 located in the MFSD1 gene whose expression in the ear is already known [START_REF] Carninci | The Transcriptional Landscape of the Mammalian Genome[END_REF] . The full list of SNPs with a P-value lower than 1x10 -6 for the three PTAs, is reported in Table1.

Concerning PCs a significant association was detected at PC1 (rs2687481, P-value 3.2 x10 -7 ) with GRM8 gene, a glutamate receptor that inhibits adenylyl cyclase decreasing the formation of cAMP [START_REF] Scherer | GRM8 The human metabotropic glutamate receptor 8 (GRM8) gene: a disproportionately large gene located at 7q31.3-q32.1[END_REF] . Additional associations were detected with several SNPs located within the FGF14 gene, a member of the fibroblast growth factors family and rs7182802 (P-value 9.9 x10 -6 ) located inside GABRG3, a ion channel gene differentially expressed in the aging ear [START_REF] Krejci | Molecular pathology of the fibroblast growth factor family[END_REF][START_REF] Souza | Novel approach to select genes from RMA normalized microarray data using functional hearing tests in aging mice[END_REF] .

Concerning PC2, rs6669265 (P-value 6 x10 -7 ), a SNP in LD with NR0B2 gene, whose protein is an unusual orphan receptor that contains a putative ligand-binding site, has been identified [START_REF] Cellanetti | Insights into the binding mode and mechanism of action of some atypical retinoids as ligands of the small heterodimer partner (SHP) LPPR5[END_REF] . The protein encoded by NR0B2, has been shown to interact with retinoid and thyroid hormone receptors, inhibiting their ligand-dependent transcriptional activation [START_REF] Bryant | Sensory organ development in the inner ear: molecular and cellular mechanisms[END_REF] . Two additional loci at P-value<10 -7 have been detected: one on chromosome 11 and the other one on chromosome 17.

Additional interesting SNPs are in LD with OTX2 (P-value 5 x10 -6 ) that acts as a transcription factor and may play a role in brain and sensory organ development [START_REF] Sajan | Toward a Systems Biology of Mouse Inner Ear Organogenesis: Gene Expression Pathways, Patterns and Network Analysis[END_REF] . Another relevant SNP is rs3783041 on chromosome 13 (P-value 5 x10 -6 ) that is in LD with DIAPH3. This gene belongs to the diaphanous family and promotes actin polymerization. It is required for cytokinesis, stress fiber formation, and transcriptional activation of the serum response factor [START_REF] Mansour | Hearing loss in a mouse model of Muenke syndrome[END_REF] .

Finally, regarding the different thresholds (Table 3), rs10815873, a SNP within PTPRD gene was found to be strongly associated at 125 Hz threshold. Data available from the literature show that some members of this gene family play an important role in the hearing system [START_REF] Goodyear | A receptor-like inositol lipid phosphatase is required for the maturation of developing cochlear hair bundles[END_REF] . The SNP rs248626, already detected for PTAL (see PTA results section and Table 1), was also found associated with 500 Hz. As already mentioned, it is in LD with a cluster of protochaderins genes.

Similarly, rs898967 located within CMIP gene (see PTA results and Table 1) was found to be significantly associated with 1KHz. Another interesting SNP at the 250 Hz is rs1849287 on chromosome 15 located in NMB gene; a member of this family (NMB1) blocks mechanotransducing ion channels in cochlear hair cells [START_REF] Drew | High-threshold mechanosensitive ion channels blocked by a novel conopeptide mediate pressure-evoked pain[END_REF] .

Others very interesting data have been obtained at 4KHz. The SNP rs9574464 (chromosome 13) shows an association (P-value 3.2x10 -7 ). This SNP is in LD with DCLK1 (see NCBI), a member of the protein kinase superfamily and the doublecortin family which is expressed in the inner ear. Additional suggestive SNPs are: a) rs2660178, located on chromosome 10, which is in LD with PCDH15, a gene already involved in causing hearing loss [START_REF] Hilgert | Function and expression pattern of nonsyndromic deafness genes[END_REF] ; b) rs1719101, in LD with GRM8, a glutamate receptor already detected at PC1 and c) rs16939415 in LD with PXMP3 a gene related with GATA3 transcription factor, an early regulator of auditory system development [START_REF] Milo | Genomic analysis of the function of the transcription factor gata3 during development of the mammalian inner ear[END_REF] . At 8 KHz we identified rs9899183, on chromosome 17, that is in LD with TNFS; some members of the TNF family are important players in the initiation of acute cochlear apoptosis [START_REF] Hu | Differential expression of apoptosis-related genes in the cochlea of noise-exposed rats[END_REF] . Some of the genes/loci here described map to already known HLL loci such as DFNA30, DFNA37, DFNA42 and DFNA54 for dominant forms and DFNB17, DFNB48, DFNA57, DFNB60 for recessive forms (see HHL) whose genes are not yet characterized.

The pathway analysis was carried out using Ingenuity software and 48 seed molecules arising from GWAS data. The analysis was able to generate only one significant pathway (Figure 1) characterized by a very high P-value (1x10 -79 ) and 49 different molecules. This score reflects the extremely low probability that genes were represented in this particular network by chance alone. Thirty out of them are present within the output of our GWAS further strengthening the impact of this work, while the remaining 19 (added by the software) are other proteins that interact with them at different levels to create this predicted pathway. We then searched for different human and mouse expression databases demonstrating that 34 genes (69%) contained within this pathway are expressed in the ear (in color in Figure 1). These findings suggest that most of the genes here identified might have an interesting role in studying normal variation of hearing function. Finally, this hypothetical pathway is also characterized by a series of IPA canonical pathways, i.e. standardized pathways very well known in the literature having a molecular role in a biological system. The most represented in terms of members are "colorectal cancer metastasis signaling", "glucocorticoid receptor signaling" and "G-protein coupled receptor".

Discussion

GWAS became the tool of choice for the identification of genes for complex and quantitative traits, since they are able to analyze large amounts of data [START_REF] Friedman | GRM7 variants confer susceptibility to age-related hearing impairment[END_REF] . Despite recent progress, almost nothing is known about hearing thresholds and the molecular bases of variation of normal hearing, apart from genes identified as being directly involved in HHL (see HHL).

Here, we present the first GWAS performed on hearing traits.

Our scan resulted in the identification of genes that have a realistic biological role in hearing function and that might be considered as good candidates for further research activities. Among them, DCLK1 (4KHz), doublecortin-like kinase 1, a member of the protein kinase superfamily and the doublecortin family which is expressed in the inner ear. The encoded protein is involved in several different cellular processes, including neuronal migration in the developing brain and in maturation of nervous system. A possible role of DCLK-1 in the maturation of the nervous system could be interesting on hypothesizing an important role of this gene in the development of sensitive neurons that we know essential for the hearing function. Another interesting gene is PTPRD (125 Hz), a member of the protein tyrosine phosphatase (PTP) family. PTPs are known to be signaling molecules that regulate a variety of cellular processes including cell growth, differentiation, mitotic cycle, and oncogenic transformation. Some members of this family play an important role in the hearing system. For example, studies of homologous genes in others species suggest the role of PTPs family in promoting neurite growth and regulating neurons axon guidance, both mechanisms important for neuronal development. A third gene is GRM8 (PC1), a member, together with GRM7, of the group III of metabotropic glutamate receptors family which has been divided into three groups on the basis of sequence homology, putative signal transduction mechanisms, and pharmacologic properties [START_REF] Scherer | GRM8 The human metabotropic glutamate receptor 8 (GRM8) gene: a disproportionately large gene located at 7q31.3-q32.1[END_REF], 34 . GRM8 and GRM7 proteins show 87% of homology and 76% of identity using BLAST analysis. Very interestingly, both associations with GRM8 and GRM7 have been detected at PC1 in two independent studies, one case controls focused on ARHI 6 (qualitative trait) and the present one on genetic bases of variation of normal hearing function (quantitative trait). A fourth significant candidate is NR0B2 (PC2), whose product interacts with retinoid and thyroid hormone receptors, inhibiting their ligand-dependent transcriptional activation [START_REF] Cellanetti | Insights into the binding mode and mechanism of action of some atypical retinoids as ligands of the small heterodimer partner (SHP) LPPR5[END_REF] . Nuclear receptors for thyroid hormone and retinoic acid are expressed in the developing sensory epithelia of the inner ear and their ligands play roles in hair cell development [START_REF] Bryant | Sensory organ development in the inner ear: molecular and cellular mechanisms[END_REF] . The last candidate is CMIP at 1kHz (c-MAF inducing protein), a gene also expressed in the inner ear.

Additional strong biological candidates arose from results of SNPs within a group showing a maximum p-value of 1x10 -6 . In particular, we should mention FGF14 (PC1) [START_REF] Krejci | Molecular pathology of the fibroblast growth factor family[END_REF] . It has been demonstrated that Fibroblast growth factors have been implicated in a wide range of cellular processes. Although the role of FGF signaling in the maintenance of normal auditory function remains to be elucidated, some members of this family play a dosage-sensitive role in the differentiation of the auditory sensory epithelium and have a fundamental function in cochlea development. Moreover, the inhibition of FGFs signaling could cause a reduction in hair cells and support cells and an alteration of FGFRS, which signal are produced by auditory neurons, could modify the right development of the cochlea [START_REF] Krejci | Molecular pathology of the fibroblast growth factor family[END_REF] ; ii) LPRP4 (PC3) and a quite long list of other genes such as PCHD15, KCNMB2, PIBF1, OPTN, DFFB, MFSD1, NMB, PXMP3, and TNFSF12 gene. They are good candidates for variation of normal hearing function, and could be eventually considered as candidates for both NIHL and ARHI.

Since previous studies on molecular basis of inherited hearing loss have detected loci/genes which show, as expected, a minimum overlap with those here reported, present data describe additional genes which might be involved in hearing function.

These results led to "in silico" building of an extremely interesting pathway characterized by a network of genes in which the vast majority of them are expressed in the inner ear. Despite it represents an hypothetical network, it is characterized by a series of significant functional relationships amongst proteins and molecules with 3 main canonical networks: Colorectal cancer metastasis signaling, Glucocorticoid receptor signaling and G-protein coupled receptor.

As regards to the first network, recently has been found that DFNA5 (a gene causing a non syndromic autosomal dominant type of hearing loss) is also involved in some type of cancer such as colorectal one [START_REF] Vrijens | Characterization of the murine Dfna5 promoter and regulatory regions[END_REF] . Concerning the second one, although the exact mechanism of glucocorticoid action on the inner ear are not known, the inner ear of both humans and experimental animals demonstrate an abundance of glucocorticoid receptors in both neuronal and non-neuronal tissues. Moreover, glucocorticoids are widely used to treat different hearing disorders [START_REF] Canlon | Glucocorticoid receptors modulate auditory sensitivity to acoustic trauma[END_REF] . Finally, as regards to the last network, some members of G-protein coupled receptor family are known to play a relevant role in hearing function.

Moreover, among the molecules added by the Ingenuity software to build the final pathway we should mention BMP4, VDR, and PTHLH. Bone Morphogenetic Protein 4 (BMP4) is a member of the TGF-β superfamily and is known to be important for the normal development of many tissues and organs, including the inner ear. Recent studies have been demonstrated that Bmp4 heterozygous null (Bmp4+/-) mice are viable and some adults exhibit an inner ear defect [START_REF] Blauwkamp | The Role of Bone Morphogenetic Protein 4 in Inner Ear development and function[END_REF] . Regarding VDR (Vitamin D receptor), it seems to be important with a wide range of things and it has several important biological roles. VIT D deficiency, VDR malfunction, hypoparathyroidism and hypervitaminosis have been suggested to be potential caused of sensorineural hearing loss [START_REF] Zou | Progressive hearing loss in mice with a mutated vitamin D receptor gene[END_REF] . As a matter of fact, the VDR knock-out mouse shows a severe calcification in the thalamus causing an alteration of the connection between inferior colliculus and auditory cortex. As mentioned before, an alteration of expression of parathyroid hormone is related with hearing and, to be more precise, this alteration could contribute to the abnormal bone turnover in otoscleorosis; in this light PTHLH is another strong candidate [START_REF] Grayeli | Parathyroid hormone-parathyroid hormone-related peptide receptor expression and function in otosclerosis[END_REF] .

Of course, all the proteins present in the pathway can be taken into account as a basis for further in vitro functional experiments.

In conclusion, here we report a first series of data on hearing quantitative traits (i.e. normal hearing function). Candidate genes located in positive GWAS regions belong to several different gene families that show only a small overlap with those already identified as causing hearing impairments. However, present findings should be further confirmed and/or replicated in other populations/cohorts. Anyway, these results increase our knowledge of the molecular basis of normal hearing function and might open new perspectives for preventive and therapeutic strategies for hearing impairments.

Licence for Publication

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence on a worldwide basis to the BMJ

Table 1 . Top SNPs resulting from the mixed-model analysis for PTAL (low), PTAM (medium) and PTAH (high

 1 ). Top SNPs obtained after meta-analysis of GWAS data are listed according different PTAs. Direction of effect in: Cilento, Split, Korcula, Carlantino, Friuli Venezia Giulia, Talana. Betas for single population have not been reported in value but in direction only. Abbreviations: Chr., chromosome, e.a., effect allele., SE ß., standard error ß.

	Publishing Group Ltd to permit this article (if accepted) to be published in Journal of
	Medical Genetics and any other BMJPGL products and sublicences such use and
	exploit	all	subsidiary	rights,	as	set	out	in
	(http://group.bmj.com/products/journals/instructions-for-authors/licence-		forms)."

Table 2 . Top SNPs resulting from the mixed-model analysis for PC1, PC2 and PC3

 2

. Top SNPs obtained after meta-analysis of GWAS data are listed according different PCs. Direction of effect in: Cilento, Split, Korcula, Carlantino, Friuli Venezia Giulia, Talana. Betas for single population have not been reported in value but in direction only. Abbreviations: Chr., chromosome, e.a., effect allele., SE ß., standard error ß.

Table 3 . Top SNPs resulting from the mixed-model analysis for 125Hz, 250Hz, 500Hz, 1kHz, 2kHz, 4kHz, 8kHz

 3 . Top SNPs obtained after meta-analysis of GWAS data are listed according different tresholds. Direction of effect in: Cilento, Split, Korcula, Carlantino, Friuli Venezia Giulia, Talana. Betas for single population have not been reported in value but in direction only. Abbreviations: Chr., chromosome, e.a., effect allele., SE ß., standard error ß.

	Supplemetary tables 1., 2., 3.

Acknowledgements

We thank Laura Esposito and Angela D'Eustacchio for technical support. The study was partially funded by Telethon Foundation (GGP09037), Fondo Trieste (2008), Regione FVG (L.26.2008), and Italian Ministry of Health (RC16/06) (to PG). The KORCULA-CROATIA study in the Croatian island of Vis was supported through the grants from the Medical Research Council UK to H.C., A.F.W. and I.R.; and Ministry of Science, Education and Sport of the Republic of Croatia to I.R. (number 108-1080315-0302). The SPLIT-CROATIA study was supported through the grants from the Medical Research Council UK; and Ministry of Science, Education and Sport of the Republic of Croatia. (number 108-1080315-0302). We would like to acknowledge the invaluable contributions of the recruitment team from Korcula, the Croatian Centre for Global Health, University of Split, the administrative teams in Croatia and Edinburgh and the people of Korcula and Split. The Cilento study was supported by grants from the Italian Ministry of Universities (FIRB -RBIN064YAT) and the Ente Parco Nazionale del Cilento e Vallo di Diano to MC. RS was supported by a fellowship of Regione Campania,

Italy. We thank the Ogliastra population and all the individuals who participated in this study. We are very grateful to the municipal administrators for their collaboration to the project and for economic and logistic support.

Web Resources

The URLs for data presented herein are as follows:

Database of human genes: Genecards, http://www.genecards.org/ Meta-analysis: METAL, http://www.sph.umich.edu/csg/abecasis/Metal/ U.S. National Institutes of Health (NIH) free digital archive of biomedical and life sciences: http://www.ncbi.nlm.nih.gov/ (NCBI) Hereditary Hearing loss Homepage (HHL), http://hereditaryhearingloss.org/ Ingenuity software for pathway analysis, http://www.ingenuity.com/ Eurexpress, http://www.eurexpress.org/ee/ Jackson Laboratory, http://www.informatics.jax.org/

Conflict of interest

The authors declare no conflict of interest

Supplemental Data

Supplemental data include three tables that can be found with this article online. rs9574476 13 35215647 C -0.16 0.03 ------


Bárbara J Henriques

Peter Bross

Cláudio M Gomes
email: gomes@itqb.unl.pt

Mutational hotspots in electron transfer flavoprotein underlie defective folding and function in multiple acyl-CoA dehydrogenase deficiency

Keywords: ETF, MADD, Molecular chaperones, Spectroscopy, Protein, Metabolic disorder MADD, multiple acyl-CoA dehydrogenase deficiency, WT, wild type, ETF, electron transfer flavoprotein, PDB, Protein Data Bank, CD, circular dichroism, SPR, surface Plasmon resonance, T m , midpoint unfolding temperature, C m , denaturant midpoint transition concentration

We have carried out an extensive in silico analysis on 18 disease associated missense mutations found in electron transfer flavoprotein (ETF), and found that mutations fall essentially in two groups, one in which mutations affect protein folding and assembly, and another one in which mutations impair catalytic activity and disrupt interactions with partner dehydrogenases. We have further experimentally analyzed three of these mutations, ETFβp.Cys42Arg, ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys, which have been found in homozygous form in patients and which typify different scenarios in respect to the clinical phenotypes. The ETFβ-p.Cys42Arg mutation, related to a severe form of multiple acyl-CoA dehydrogenase deficiency (MADD), affects directly the AMP binding site and intersubunit contacts and impairs correct protein folding. The two other variations, ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys, are both associated with mild MADD, but these mutations have a different impact on ETF. Although none affects the overall α/β fold topology as shown by far-UV CD, analysis of the purified proteins shows that both have substantially decreased enzymatic activity and conformational stability. Altogether, this study combines in silico analysis of mutations with experimental data and has allowed establishing structural hotspots within the ETF fold that are useful to provide a rationale for the prediction of effects of mutations in ETF.

A C C E P T E D M A N U S C R I P T

INTRODUCTION

Mutations in the ETFA/ETFB or ETFDH genes, which encode the two subunits of Electron Transfer Flavoprotein (ETF) and ETF-ubiquinone oxidoreductase (ETF-QO), respectively, cause multiple acyl-CoA dehydrogenase deficiency (MADD), a congenital metabolic disease with a broad clinical expression [START_REF] Frerman | Defects of Electron Transfer Flavoprotein and Electron Transfer Flavoprotein-Ubiquinone Oxidoreductase: Glutaric Acidemia Type II, The Metabolic & Molecular Basis of Inherited Disease McGrawHill[END_REF]. The systematic application of newborn screening programs has allowed for a wider coverage of MADD genotypes, that range from missense or nonsense mutations to frameshift. One of the proteins affected in this disorder is ETF, a key enzyme in the mitochondrial fatty acid beta oxidation and amino acid degradation pathways, which transfers electrons from at least 12 dehydrogenases via ETF-QO to the respiratory chain. As in many other metabolic disorders, clear genotype-phenotype relationships remain elusive in MADD: in particular for milder disease forms, it is expected that the mutant proteins are partially functional and that cellular and environmental factors such as temperature and cofactors can modulate disease expression [START_REF] Bross | Protein misfolding and degradation in genetic diseases[END_REF]. Progress in this issue has however also been constrained by the fact that only a very limited number of mutations has been characterised upon recombinant expression of the variant proteins. However, substantial steps forward can be achieved if the available data on clinical mutations are merged with a structural and functional analysis of the protein affected. In this paper we have investigated for mutational hotspots on the ETF structure, by analysing the protein for particularly susceptible regions in which MADD-related ETF variations detected in patients map.

Complementary experimental validation for the putative mutational hotspots identified was subsequently obtained by carrying out a biochemical and biophysical analysis of purified ETF variants with representative disease-causing mutations. The data gathered by this strategy contribute to a better understanding of the molecular factors underlying functional deficiency and may, in the long term, provide a framework for the design of novel therapies.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

MATERIALS AND METHODS

Chemicals All reagents were of the highest purity grade commercially available. Octanoyl-CoA, FAD, AMP and urea were purchased from Sigma. Isopropyl-β-D-thiogalactopyranoside (IPTG) was purchased from VWR International.

Structural Analysis. The crystallographic structures of human ETF (PDB code: 1efv, [START_REF] Roberts | Three-dimensional structure of human electron transfer flavoprotein to 2.1-A resolution[END_REF])

and ETF:MCAD complex (PDB code: 1T9G, [START_REF] Toogood | Extensive domain motion and electron transfer in the human electron transferring flavoprotein.medium chain Acyl-CoA dehydrogenase complex[END_REF]) were inspected using PyMOL (DeLano Scientific). Analysis of the molecular interactions, cofactor contacts, topological features and generation of models for the variants was carried out using the WhatIF web server [START_REF] Vriend | WHAT IF: a molecular modeling and drug design program[END_REF] and the PDBsum database [START_REF] Laskowski | PDBsum: a Web-based database of summaries and analyses of all PDB structures[END_REF]. The PolyPhen (Polymorphism Phenotyping) web server was used to predict the effect of amino acid variations on protein structure and function [START_REF] Ramensky | Human non-synonymous SNPs: server and survey[END_REF]. The input for this server was the UniProt accession number (ETFα: P13804 and ETFβ:P38117), sequence position and the two amino acid variants characterizing the desired mutation.

ETF expression in E. coli JM109 cells. The pK-lac-ETF-αβ3(C/T) plasmid encoding both

ETF subunits (here named pWt), described previously [START_REF] Bross | A polymorphic variant in the human electron transfer flavoprotein alpha-chain (alpha-T171) displays decreased thermal stability and is overrepresented in very-long-chain acyl-CoA dehydrogenase-deficient patients with mild childhood presentation[END_REF], was used as DNA template in the

ACCEPTED MANUSCRIPT

were the reverse complements of those shown. After PCR screening the sections containing the mutation were sequenced to exclude PCR errors and subcloned back into the original plasmid. The ETFβ-D128N mutant, also used in these studies, was already available [START_REF] Olsen | Clear relationship between ETF/ETFDH genotype and phenotype in patients with multiple acyl-CoA dehydrogenation deficiency[END_REF]. E. coli JM109 cells transformed with the respective plasmid, pWt, pC42R, pD128N, pR191C, were grown in LB (Luria-Bertani) supplemented with 10 µgml -1 kanamycin at 30ºC or 37ºC in a shaking incubator until OD 532 of 0.5-0.8 was reached. The cells were than induced with 1 mM IPTG for 4 hours. Cells were harvested by centrifugation, re-suspended in 10 mM Hepes, 10% ethylene glycol and 0.5 mM phenylmethylsulphonylfluoride (Roth) in presence of DNase (Applichem) and disrupted in a French press. The soluble extract was obtained by centrifugation at 12,000 g for 30 minutes and used for the ETF activity assays [START_REF] Rhead | Multiple acyl-coenzyme A dehydrogenation disorder responsive to riboflavin: substrate oxidation, flavin metabolism, and flavoenzyme activities in fibroblasts[END_REF] and for SDS-PAGE followed by western blotting analyses according to [START_REF] Bross | Co-overexpression of bacterial GroESL chaperonins partly overcomes non-productive folding and tetramer assembly of E. coli-expressed human medium-chain acyl-CoA dehydrogenase (MCAD) carrying the prevalent disease-causing K304E mutation[END_REF]. An ETF polyclonal antibody raised in rabbits against purified recombinant wild type ETF protein was used to identify ETF proteins and visualized with peroxidase conjugated secondary antibody reaction with 1-step TMB-Blotting (Pierce).

Co-expression of Molecular chaperones with ETFβ-p.Cys42Arg. Different conditions

were used to improve ETFβ-p.Cys42Arg mutant expression, always using LB as growth medium. First, E. coli cells transform with pC42R plasmid were expressed at 30ºC and 18ºC.

Induction was made with 1 mM IPTG. JM109 E. coli cells were co-transformed with the pC42R plasmid and a molecular chaperone plasmid, pOFX-BADSL2 and pOFX-BADKJE2, respectively for the GroES/EL complex and DnaK,J and GprE [START_REF] Castanie | A set of pBR322compatible plasmids allowing the testing of chaperone-assisted folding of proteins overexpressed in Escherichia coli[END_REF]. Cells transformed with the ETFβ-p.Cys42Arg and with pOFX-BAD2 empty plasmid, control vector, or with the plasmid encoding the respective molecular chaperone system, were grown in LB supplemented with 10 µgml -1 kanamycin and 10 µgml -1 tetracycline at 30ºC and 18ºC until OD 532 of 0.5-0.8 was reached. Cells were then induced with 0.2% arabinose and 30 minutes Flavin content was determined using the molar extinction coefficient ε 436nm =13400 M -1 .cm -1 reported for FAD bound to ETF [START_REF] Mckean | Subunit structure of electron transfer flavoprotein[END_REF]. ETF enzyme activity was measured following 2,6dichlorophenolindophenol (DCPIP) reduction at 600 nm in a coupled assay in which recombinant human MCAD and octanoyl-CoA were employed, as described in [START_REF] Rhead | Multiple acyl-coenzyme A dehydrogenation disorder responsive to riboflavin: substrate oxidation, flavin metabolism, and flavoenzyme activities in fibroblasts[END_REF]. One unit of catalytic activity is defined as nmol of DCPIP reduced per minute, in the conditions used in the assay. All specific activities reported are based on total flavin content. Pure ETF fractions with a 2.5 fold molar excess FAD were fast-frozen using liquid nitrogen and stored -80ºC.

Spectroscopic methods.

Before each experiment FAD excess added to buffers as a preservative was removed by extensive washing using ultra filtration/dilution, and all experiments were performed with pure proteins containing full occupancy of FAD site.

UV/visible spectra were recorded at room temperature in a Shimadzu UVPC-1601

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

spectrometer with cell stirring. Fluorescence spectroscopy was performed using a Cary Eclipse instrument. For tryptophan emission excitation wavelength was set at 280 nm, and FAD emission was followed setting excitation wavelength at 436 nm; slits were 5 and 10 nm for excitation and emission, respectively. Typically protein concentration was 1 µM. CD spectra were recorded on a Jasco J-815 spectropolarimeter with a cell holder thermostatically controlled with a Peltier. A quartz polarized 1 mm path length quartz cuvette (Hellma) was used, and protein concentrations were typically 0.1 mg.ml -1 . Thermal stability. Thermal unfolding with a linear temperature increase was followed using circular dichroism (ellipticity variation at 222 nm) and fluorescence spectroscopy, tryptophan emission (λ ex =280 nm; λ em =340 nm), FAD emission (λ ex =436nm; λ em =530 nm) and FRET from tryptophan emission to FAD cofactor (λ ex =280 nm; λ em =530 nm). In all experiments, a

Surface Plasmon

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

heating rate of 1ºC min -1 was used, and temperature was increased from 30 to 90ºC. Data were analysed according to a two-state model, and fits to the transition curves were made using OriginPro8.

Chemical Stability. The denaturation curves were measured diluting ETF (ETF-WT, ETFβ-p.Asp128Asn or ETFβ-p.Arg191Cys) in different urea concentrations and the tryptophan emission spectrum was recorded after 15 minutes of incubation at 25°C.

Transition curves were determined plotting the average emission wavelength against urea concentrations. Data were analysed according to a two-state model, and fits to the transition curves were made using OriginPro8.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

RESULTS AND DISCUSSION

Mutational hotspots in the ETF structure. A number of genetic defects in the ETFA/ETFB gene account for MADD. Most of the variations detected in patients are of the missense type, although frameshift and nonsense mutations are also described [START_REF] Stenson | The Human Gene Mutation Database: providing a comprehensive central mutation database for molecular diagnostics and personalized genomics[END_REF]. As the effect of modifying a single amino acid in ETF can be rather broad with respect to clinical expression, the possibility to predict whether a point mutation in a particular region of the protein would affect the folding, stability or activity, would be a potentially very useful parameter to establish genotype-phenotype relationships. In order to fill this gap we have set to define a structural mapping of mutations in ETF. For this purpose our strategy was twofold. First, we have compiled a list of known disease-associated missense mutations described for human ETF, whose predicted effects were analysed using a combination of in silico tools for mutagenesis predictions and structural analysis (Table 1). Secondly, we have analysed the ETF crystal structure in respect to specific structural factors that, if affected by mutations, are expected to perturb the structural integrity or the biological activity of the protein [START_REF] Roberts | Three-dimensional structure of human electron transfer flavoprotein to 2.1-A resolution[END_REF][START_REF] Toogood | Extensive domain motion and electron transfer in the human electron transferring flavoprotein.medium chain Acyl-CoA dehydrogenase complex[END_REF]. In this respect, we have analysed: i. positions that correspond to AMP-protein interactions; ii.

residues that are involved in H-bond connections at the ETFα/ETFβ contact surface; iii.

highly conserved segments throughout ETF ortologues; iv. amino acids involved in FAD binding or in the stabilization of the FAD moiety, and; v. positions involved in interactions with medium-chain acyl-CoA dehydrogenase (MCAD), both the anchor loop and the FAD domain. At this stage positions involved in the interaction with ETF-QO were not possible to analyze once there is no structural information available on the ETF:ETF-QO complex.

Merging the data obtained by these two complementary approaches shows that particular regions of the ETF structure are more susceptible to deleterious modifications than others, and the clustering of the known mutations in those structure segments is rather suggestive of structural hotspots for MADD mutations (Fig 1).

On the other hand, mutations in the second region have an impact on catalysis. This occurs either by disturbing directly the FAD protein interactions, or the second coordination sphere ligands, which are important modulators of the redox properties of the cofactor. In fact, the chemical and catalytic properties of the flavin are tightly controlled by the polypeptide, not only at the level of the residues involved in direct close-range interactions (i.e. within the first coordination sphere), but also by those which are further away from the organic cofactor but that are nevertheless connected to it through hydrogen-bond networks or other weak interactions (i.e. second coordination sphere). Also, this segment comprises positions of the FAD domain that are known to be involved in complex formation with MCAD [START_REF] Toogood | Stabilization of non-productive conformations underpins rapid electron transfer to electron-transferring flavoprotein[END_REF].

Therefore, mutations in this region are expected to affect mostly the biological activity of the protein, although concurrent destabilization may also occur.

In respect to the ETFβ subunit, the scenario is not so clear as a result of the fact that only 5 missense mutations are known in this subunit. Nevertheless, some highly susceptible regions can be defined: those involving two N-terminal segments of ETFβ which harbour simultaneously FAD and AMP binding residues (ETFβ Arg5-Ile20 and ETFβ Val34-Ala45), a region rich in AMP-protein interactions (ETFβ Gly123-Thr134), and segments involved in

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

interactions with partner proteins, namely the recognition peptide within the anchor domain that constitutes the primary interaction site of partner proteins with ETF (ETFβ Asp184-Lys200) [START_REF] Toogood | Extensive domain motion and electron transfer in the human electron transferring flavoprotein.medium chain Acyl-CoA dehydrogenase complex[END_REF], and also a short segment centred around ETFβ-Glu165, an essential residue that destabilises positions compatible with fast inter-protein electron transfer, thus ensuring high complex dissociation rates [START_REF] Toogood | Stabilization of non-productive conformations underpins rapid electron transfer to electron-transferring flavoprotein[END_REF], a strict requirement for subsequent interaction of reduced ETF with ETF:QO.

Analysis of prototypic mutations. In order to provide complementary experimental validation of the putative mutational hotspots identified in ETF, three mutations in ETFβ were selected for subsequent analysis: ETFβ-p.Cys42Arg, ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys. With the exception of ETFβ-p.Asp128Asn [START_REF] Olsen | Clear relationship between ETF/ETFDH genotype and phenotype in patients with multiple acyl-CoA dehydrogenation deficiency[END_REF][START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF], none of these mutations had been previously expressed and purified for a thorough in vitro functional and structural characterisation. All three mutations were detected in patients homozygous for the respective mutations allowing to relate the results to the phenotypes. Also, this selection fills a gap in respect to the characterisation of mutations in the ETFβ subunit, in which a lower number of MADD related mutations have been described compared to ETFα. Although an all-embracing division cannot be established, these mutations typify different functional and folding effects. To better predict the structural effect of these individual mutations we have built molecular models of these variants by homology modelling with respect to the wild type protein (Fig. 2). The ETFβ-p.Cys42Arg mutation, associated to a severe clinical expression, is expected to severely affect the interaction with the AMP cofactor, which plays a key role in the assembly of the ETF dimer [START_REF] Sato | In vitro assembly of FAD, AMP, and the two subunits of electron-transferring flavoprotein: an important role of AMP related with the conformational change of the apoprotein[END_REF][START_REF] Sato | In vitro refolding and unfolding of subunits of electrontransferring flavoprotein: characterization of the folding intermediates and the effects of FAD and AMP on the folding reaction[END_REF]. The ETFβ-p.Asp128Asn modification, which occurs at an outer layer of the FAD interacting moieties, has been shown to affect the plasticity of the tertiary structure and decrease directly the specific activity [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF], and its mutation likely affects

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

a stabilizing H-bond to ETFβ-Lys11 (Fig. 2). Finally, the ETFβ-p.Arg191Cys mutation is hypothesised to affect catalysis by impairing interaction with the partner dehydrogenase, as this residue is localized in the anchor domain. Also, analysis of the model of the mutant variation shows that removal of ETFβ-Arg191 may destabilize the fold by disrupting electrostatic interactions with the nearby ETFβ-Glu47, ETFβ-Arg51 and ETFβ-Glu54 (Fig. 2). The patients with the latter two mutations display mild forms of the disease, indicating that the modified conformation is nevertheless able to fold and assemble as a dimer, retaining some catalytic proficiency.

Mutations affect folding efficiency. The effect of each mutation on the efficiency of protein folding was investigated by heterologous expression of the different variants and differential analysis of the soluble versus insoluble protein produced, in respect to that of wild type ETF and comparing cell cultures with standardized induction protocols displaying superimposable cell growth characteristics. Therefore, any difference strictly results from distinct expression levels of the heterologously expressed protein.

Since it has been previously demonstrated that expression of ETFβ-p.Asp128Asn was improved upon cell growth at 30ºC [START_REF] Olsen | Clear relationship between ETF/ETFDH genotype and phenotype in patients with multiple acyl-CoA dehydrogenation deficiency[END_REF], over expression all of the analysed ETF variants was carried out both at 30 and 37°C (Fig. 3). The results obtained show that while wild type ETF and the mild mutant variants (ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys) were well expressed as soluble proteins, the severe phenotype mutation (ETFβ-p.Cys42Arg) was not (Fig. 3). The expression level of ETFβ-p.Arg191Cys was 87% of that of the soluble protein detected in the wild type whereas that of ETFβ-p.Asp128Asn was decreased to 48%. This suggests that while the expression of the variant ETFβ-p.Arg191Cys is almost unaffected, the ETFβ-p.Asp128Asn modification has a destabilizing effect on the folding process, making it rather susceptible to misfold and

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

aggregate, or be degraded. This observation during cell expression agrees with the data obtained in vitro using purified ETFβ-p.Asp128Asn for which a higher proteolytic susceptibility was determined [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF].

On the other hand, expression of the ETFβ-p.Cys42Arg variant yielded only 10% of ETF protein in the soluble fraction, indicating that this mutation affecting the AMP binding severely impairs ETF folding and assembly. Further expression assays lowering the E. coli growth temperature to 18ºC in an attempt to further improve conditions for in cell folding did not change the scenario, and even at this low temperature almost all the ETFβ-p.Cys42Arg

protein expressed went to the insoluble fraction (not shown). This result agrees with western blot analysis of skin fibroblasts from a patient carrying this mutation that showed no detectable ETF protein [START_REF] Curcoy | Late-onset form of beta-electron transfer flavoprotein deficiency[END_REF][START_REF] Schiff | Electron transfer flavoprotein deficiency: functional and molecular aspects[END_REF].

Molecular chaperones partly rescue ETFβ-p.Cys42Arg folding. An approach for the rescue of folding defects resulting from disease-causing mutations is that of stimulating the so called proteostasis network [START_REF] Powers | Biological and chemical approaches to diseases of proteostasis deficiency[END_REF]. This involves recruiting the protein quality control machinery and/or the action of small molecule substrates and cofactors as effectors of the folding process. In order to investigate if molecular chaperones could rescue the defective folding resulting from the ETFβ-p.Cys42Arg mutation, we have carried out co-expression experiments with the chaperonins GroEL and GroES (homologs of the human mitochondrial Hsp60 and Hsp10 chaperone system), and the chaperone DnaK (homolog of mammalian Hsp70s) and its cofactors DnaJ and GrpE. This strategy has been shown to result in the effective rescue of proteins with folding defects in other cases [START_REF] Bross | Effects of two mutations detected in medium chain acyl-CoA dehydrogenase (MCAD)-deficient patients on folding, oligomer assembly, and stability of MCAD enzyme[END_REF][START_REF] Salazar | Expression and characterization of two pathogenic mutations in human electron transfer flavoprotein[END_REF][START_REF] Andresen | The molecular basis of medium-chain acyl-CoA dehydrogenase (MCAD) deficiency in compound heterozygous patients: Is there correlation between genotype and phenotype?[END_REF][START_REF] Majtan | Rescue of cystathionine beta-synthase (CBS) mutants with chemical chaperones: purification and characterization of eight CBS mutant enzymes[END_REF]. The results obtained (Fig. 4) show that the molecular chaperones could in fact partially rescue the protein to the soluble fraction. However, enzymatic activity assays on the soluble fractions showed no

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

detectable ETF activity, indicating that although molecular chaperones successfully rescue some soluble conformations these are still be incorrectly folded and/or with defective cofactor insertion: presumably these soluble forms are still substantially destabilised and poor incorporation of FAD and/or deficient interaction with the partner dehydrogenases impair activity.

Structure and conformational properties of the ETFβ-p.Asp128Asn and

ETFβ-p.Arg191Cys variant proteins -The ETFβ-p.Arg191Cys variant was purified to homogeneity and spectroscopically characterised. The ETFβ-p.Asp128Asn was also purified and characterized as described in [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF] and the spectroscopic data from this variant were redrawn from published data [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF] in order to allow a direct comparison with the new variant described. The protein fold and secondary structure were analysed using far-UV circular dichroism (Fig. 5). The spectra obtained are typical of folded proteins with a α/β fold, and with the clear α-helix fingerprint with minima at 208 and 222nm. The mutant variants have superimposable spectra with respect to that of wild type ETF, indicating that the mutations do not disrupt the protein fold nor the secondary structure. Trp-emission spectra, which are indicative of tertiary structure interactions, were also recorded. ETF contains two Trp residues, and the fluorescence emission properties of this amino acid ranges typically from 315 to 350 nm and are related to solvent accessibility: a more solvent exposed Trp will emit at higher wavelengths. The results obtained show that both mutants have 10 nm red-shifted maxima of the emitting band, with maximum emission at 330 nm versus 320 nm of the wild type ETF (Fig. S2). This suggests that the structure of the mutant variants is not as constrained as that of the wild type ETF and that they displays a more dynamic tertiary structure. The spectroscopic fingerprint of the catalytic FAD cofactor was also analysed using visible absorption spectroscopy, which is very sensitive to changes in the H-bound interaction

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

network around the cofactor. From the results obtained we could conclude that the environment of the functional cofactor does not change significantly in any of the mutants, as the two characteristic absorption bands with typical emission maxima at around 376 nm and 436 nm, are retained. Overall, these results indicate that the ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys variants retain the fold and FAD interactions identical to those of the wild type protein but the differences in Trp fluorescence properties suggest decreased compactness and structural constraints in comparison to the wild type [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF]. The latter aspect may have an impact on the enzymatic activity of the mutants as it likely contributes to a less efficient interaction with the partner acyl-CoA dehydrogenases and to an increased protein instability.

Protein instability and MCAD interactions underpin functional deficiency -

The specific activity of the ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys variants was only 30% of that of the wild type protein (400 U/mg vs. 1300 U/mg), which shows that, in addition to the lower folding efficiency, these mutant variants also display a decreased specific activity [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF]. Since the mutant proteins retain structural features identical to those of the wild type protein under physiological conditions, we have investigated if differences in stability and interaction with the partner dehydrogenases could play a role. The interaction between the ETF variants and the medium chain acyl-CoA dehydrogenase (MCAD) was investigated in preliminary experiments by surface plasmon resonance, using MCAD immobilized in a NTA-sensor chip.

The results obtained showed that the ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys variants interact with MCAD, as an increase in the resonance signal is observed; however, the magnitude of the interaction at equilibrium was only 20-30% of that observed for wild type ETF, which is suggestive of a weaker interaction (Fig. S3).

A C C E P T E D M A N U S C R I P T

The effect of the β-p.Asp128Asn and ETFβ-p.Arg191Cys mutations on the conformational stability was also investigated spectroscopically using far-UV CD and fluorescence. Thermal denaturation of the proteins was monitored by following secondary and tertiary structure from far-UV CD and tryptophan and FAD fluorescence emissions (Fig. 6). All the techniques used yielded similar transition curves for each of the proteins analysed, showing that all methods are monitoring the same unfolding event with no noticeable intermediates formed. The analysis showed the same order of relative stabilities of the proteins involved: ETF-WT > ETFβ-p.Arg191Cys > ETFβ-p.Asp128Asn. The thermal stability studies showed that wild type ETF unfolds at an apparent midpoint unfolding temperature (T m ) of 60 ºC, and that ETFβ-p.Arg191Cys is just slightly destabilised (∆T m =-3°C, T m = 57 °C), whereas ETFβ-p.Asp128Asn exhibits a more significant destabilization (∆T m =-8°C, T m = 52 °C).

Chemical denaturation was followed using tryptophan emission, and the curves allowed the determination of apparent midpoint urea concentrations (C m ). The obtained values were 2.2, 2.9 and 3.9 M for ETFβ-p.Asp128Asn, ETFβ-p.Arg191Cys, and ETF-WT, respectively (Fig. 6B). Again, the order of stabilities was the same as in thermal denaturations indicating that the ETFβ-p.Asp128Asn variation has a stronger destabilizing effect than ETFβ-p.Arg191Cys. In fact it was already described that ETFβ-p.Asp128Asn exhibits some plasticity in tertiary structure and decreased specific activity as a function of time during incubation at 39°C, in an in vitro mimic of a fever episode [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF]. These differences among the studied variants can be rationalised from the molecular models made for the ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys variants (Fig. 2). The replacement of ETFβ-Asp128 propagates to the AMP coordination sphere by affecting an hydrogen bond between Asp128 and Lys11, therefore substantially destabilising the protein fold, as shown by thermal and chemical unfolding studies (Fig. 6). Likewise, the ETFβ-p.Arg191Cys variation also impacts on the conformational stability of the fold as removal of ETFβ-Arg191 disrupts

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

salt bridges with Glu47 and Glu57 at the nearby α1 helix. This case illustrates how a mutation in a functional region can also have an impact on overall protein stability.

CONCLUSIONS

In this work we have analysed the effects of missense mutations in ETF which correspond to different clinical presentations in MADD patients, combining in silico studies with experimental biophysical and biochemical data. Overall, this study showed that ETF mutations map to structural hotspots with respect to subunit interaction, cofactor binding and interaction with partner proteins, and this provides a framework that allows to predict more accurately the probability for a given mutation to affect protein function. Variations leading to severe phenotypes do not map to a particular structural region and a direct correlation with the clinical display can not be generalised, a general problem that has hampered the prediction of mutation effects for many disease-relevant proteins. However, it appears that amino acid changes that lead to drastic alterations of the chemical properties of residues involved in catalytic regions and subunit interactions have a higher probability to result in severe effects.

Among these is ETFβ-p.Cys42Arg which results in a protein that is unable to fold: such mutations will be clearly associated with severe clinical phenotypes. Other mutations will affect folding, cofactor binding, and assembly to different degrees, and most of them will be associated with mild phenotypes. The degree of impairment in respect to the clinical display is rather difficult to predict, as different effects can be observed: cofactor binding, interaction with partner proteins and conformational stability. A key aspect arising from our study is the fact that conformational instability of dynamic substructures may be transferred to other parts of the protein, affecting the protein as a whole [START_REF] Correia | Conformational stability of human frataxin and effect of Friedreich's ataxia-related mutations on protein folding[END_REF][START_REF] Correia | Dynamics, stability and iron-binding activity of frataxin clinical mutants[END_REF][START_REF] Kayatekin | Zinc binding modulates the entire folding free energy surface of human Cu,Zn superoxide dismutase[END_REF] . That is for example the case of ETFβp.Asp128Asn variant, in which the modification of a residue which has long range interactions with the cofactor, increases cofactor lability, decreases protein stability and

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

interaction with the partner dehydrogenase. Likewise, the ETFβ-p.Arg191Cys variation directly affecting the recognition loop also decreases protein stability, as removal of the positively charged Arg will disrupt electrostatic interactions with the acidic residues within helix α 1 (Fig. 2). This effect then propagates to the whole structure, affecting its conformational stability. The observations made for ETF represent a scenario, which is also relevant for many other proteins, as it underpins the necessity to perform and carefully evaluate experimental expression of clinical protein variants in order to collect more data that will contribute to establish general principles that can be used to more accurately predict the functional outcome of genomic variations in metabolic disease.

A C C E P T E D M A N U S C R I P T

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

ETFβ-p.Arg191Cys. Spectra were slightly offset for clarity of comparison. The data for the ETFβ-p.Asp128Asn variant was redrawn from [START_REF] Henriques | Role of flavinylation in a mild variant of multiple acyl-CoA dehydrogenation deficiency: a molecular rationale for the effects of riboflavin supplementation[END_REF]. See text and materials and methods for details. The PSIC score is the position-specific independent count parameter computed by PolyPhen [START_REF] Ramensky | Human non-synonymous SNPs: server and survey[END_REF].

  production of the two novel mutants, the ETFβ-p.Cys42Arg and ETFβ-p.Arg191Cys. The mutants were generated using the Quick change mutagenesis kit (for the ETFβ-p.Cys42Arg variant, and 5'-GAGGCTCAACGAGCCCTGCTACGCCACGCTG-3' for the ETFβ-p.Arg191Cys mutant (altered bases are underline and in bold type). The second mutagenic primers required A C C E P T E D M A N U S C R I P T

FiguresFigure 1 -

 1 Figures Figure 1 -Structural features mapped on ETF sequenceAnalysis of the ETF structure (PDB: 1efv) allowed the identification of the regions involved in FAD and AMP binding, MCAD interactions, and H-bonding between the ETFα and ETFβ subunits, which are represented as boxed segments in respect to the primary sequence positions. Stars denote MADD mutations[START_REF] Stenson | The Human Gene Mutation Database: providing a comprehensive central mutation database for molecular diagnostics and personalized genomics[END_REF] and are filled according to the clinical phenotype: severe forms (filled stars) and mild (open stars).

Figure 2 -

 2 Figure 2 -Structural details around the mutated positionsThe cartoons represent magnifications of the ETF structural regions in which the MADD variations are identified, overlaying the mutated and the original residues, as well nearby important . Straight dotted lines represent hydrogen bonds. Structures prepared using PyMOL.

Figure 3 -

 3 Figure 3 -Expression levels of ETF variants at 30 and 37°C. Expression of ETF variants at 30 and 37°C was quantified from Western Blot analysis (left) and expressed in a bar graph (right) in respect to that of the wild type protein. See text and materials and methods for details.

Figure 4 -

 4 Figure 4 -Effect of molecular chaperones on the expression of soluble ETFβp.Cys42Arg The amount of the ETFβ-p.Cys42Arg variant expressed in the soluble (s) and insoluble (i) form was quantified from western blots of cell extracts after inducing for 1 and 6h the coexpression of the molecular chaperones GroEL/ES (left) or DnaK/J/GprE (right). The bar graph shows the relative fraction of soluble (light grey) and insoluble (dark grey) protein. See text and materials and methods for details.

Figure 5 -

 5 Figure 5 -Spectroscopic analysis (A) far-UV CD (B) Visible absorptionThe different variants were analysed in respect to their far-UV CD (A) and visible absorption (B) spectra (the legend is common for both panels): a. ETFβ-p.Asp128Asn, b. ETFβ-WT, c.

Figure 6 -

 6 Figure 6 -Stability profiles of ETFβ-p.Asp128Asn and ETFβ-p.Arg191Cys variants The thermal (A) and chemical (B) denaturation profiles of the ETFβ-p.Asp128Asn (inverted triangles) and ETFβ-p.Arg191Cys (open circles) were compared in respect to that of wild type ETF (closed squares). The solid curves represent two-state sigmoid curves from which the apparent midpoint denaturations were determined: melting temperatures (T m ) of 52, 57 and 60°C and midpoint denaturant concentrations (C m ) of 2.2, 2.9 and 3.9M urea for ETFβ-p.Asp128Asn, ETFβ-p.Arg191Cys, and ETF-WT, respectively. See text and materials and methods for details.

Figure 5 (Figure 6 (

 56 Figure 1 (BW for Print)

Table 1 -Mutations in ETF related to MADD.

 1

		PSIC score	Structural rationale	Ref
	ETF α
	p.Leu95Val	0.085 Near inter subunit H-bonded region.	[27]
	p.Gly116Arg	2.667	Changes at buried site: hydrophobicity, overpacking and charge.	[22, 28]
			Change at buried site hydrophobicity. Closest
	p.Arg122Lys	1.507	contact is ETFβ-Gln146. Affects subunits	[19]
			interface.
	p.Phe144Ser	1.781	Changes at buried site: hydrophobicity and cavity creation.	[19]
			Closest contact is ETFβ-Ser223. Cavity
	p.Val157Gly	2.346	creation. Increase in flexibility. Affects subunits	[29]
			interface.
	p.Val165Ala	1.373 Near inter subunit H-bonded region.	[19]
	p.Leu212Pro	2.798 Change at buried site hydrophobicity.	[19]
	p.Arg249Cys	3.472	Disruption of FAD binding site. Stabilizes FAD semiquinone. Contacts with MCAD.	[19].
	p.Gly255Val	2.988 Maps at FAD interacting region.	[30]
			Disruption of FAD binding site.
	p.Thr266Met	2.834	Hydrophobicity change at buried site. Contacts	[7, 22, 28]
			with MCAD.
	p.Gly267Arg	2.982	Maps at FAD interacting region. Charged residue inserted in turn. Contacts with MCAD.	[30]
	ETF β
	p.Ala17Pro	1.612	Maps at FAD and AMP interacting regions. Maps nearby ETFα-Ile284.	[31]
			Disruption of AMP binding site. Changes at
	p.Cys42Arg	3.343	buried site: overpacking, hydrophobicity and	[18, 19]
			charge.
	p.Asp128Asn	2.438	Charge change at buried site. Maps at AMP interacting regions. Maps nearby ETFα-Tyr149.	[7, 15]
	p.Arg164Gln	2.492	Charge change at buried site. Near inter subunit H-bonded region. Maps nearby ETFα-Asn118.	[32]
	p.Arg191Cys	3.515	Removal of charged residue. Contacts with MCAD.	[19]

Acknowledgments

The work was supported by the Fundação para a Ciência e Tecnologia (FCT/MCTES, Portugal) through research grant PTDC/SAU-GMG/70033/2006 (to C.M.G.), and a PhD fellowship SFRH/BD/29200/2006 (to B.J.H.).

A C C E P T E D

M A N U S C R I P T