Next Generation Sequencing and Animal Models Reveal SLC9A3R1 as a New Gene Involved in Human Age-Related Hearing Loss

Age-related hearing loss (ARHL) is the most common sensory impairment in the elderly affecting millions of people worldwide. To shed light on the genetics of ARHL, a large cohort of 464 Italian patients has been deeply characterized at clinical and molecular level. In particular, 46 candidate genes, selected on the basis of genome-wide association studies (GWAS), animal models and literature updates, were analyzed by targeted re-sequencing. After filtering and prioritization steps, SLC9A3R1 has been identified as a strong candidate and then validated by “in vitro” and “in vivo” studies. Briefly, a rare (MAF: 2.886e-5) missense variant c.539G > A, p.(R180Q) was detected in two unrelated male patients affected by ARHL characterized by a severe to profound high-frequency hearing loss. The variant, predicted as damaging, was not present in healthy matched controls. Protein modeling confirmed the pathogenic effect of p.(R180Q) variant on protein’s structure leading to a change in the total number of hydrogen bonds. In situ hybridization showed slc9a3r1 expression in zebrafish inner ear. A zebrafish knock-in model, generated by CRISPR-Cas9 technology, revealed a reduced auditory response at all frequencies in slc9a3r1R180Q/R180Q mutants compared to slc9a3r1+/+ and slc9a3r1+/R180Q animals. Moreover, a significant reduction (5.8%) in the total volume of the saccular otolith (which is responsible for sound detection) was observed in slc9a3r1R180Q/R180Q compared to slc9a3r1+/+ (P = 0.0014), while the utricular otolith, necessary for balance, was not affected in agreement with the human phenotype. Overall, these data strongly support the role of SLC9A3R1 gene in the pathogenesis of ARHL opening new perspectives in terms of diagnosis, prevention and treatment.


INTRODUCTION
Age-related Hearing loss (ARHL) is the predominant sensory impairment in the elderly (Huang and Tang, 2010;Kidd Iii and Bao, 2012), affecting millions of people worldwide (National Health and Nutrition Examination Survey [NHANES], 2015). Moreover, projections suggest that, in the next decade, the number of patients will double, largely due to an increased lifespan (Bowl and Dawson, 2015).
Age-related Hearing loss is characterized by bilateral and progressive hearing loss that usually starts at the high frequencies (Huang and Tang, 2010) causing communication difficulties associated with cognitive decline, social isolation and depression (Bowl and Dawson, 2015). It is a complex disease in which genetic and environmental risk factors (i.e., noise, smoking, alcohol, etc.) interplay (Bovo et al., 2011;Ohgami et al., 2013;Vuckovic et al., 2013Vuckovic et al., , 2014. To date, it is still unclear whether ARHL is caused by (1) rare Mendelian gene variants with large effect size or (2) multiple variants each contributing to the disease. Despite significant research efforts carried out during the last decade, the genetics risk factors involved in ARHL are still mainly unknown. To date, only few ARHL susceptibility genes have been detected by genome-wide association studies (GWAS) (Friedman et al., 2009;Van Laer et al., 2010;Girotto et al., 2011Girotto et al., , 2014Wolber et al., 2014;Vuckovic et al., 2015), and thus there is a strong need to plan for new research activities aimed at understanding the molecular mechanisms underlying this diseases and to define possible targets for therapeutic and preventive plans.
The use of animal models is a powerful tool for the discovery and/or the validation of human disease-genes. In particular, zebrafish has become an attractive model for the study of the development and function of the vertebrate inner ear (Whitfield et al., 2002). Although the zebrafish ear does not contain an equivalent of the mammalian cochlea, many features (e.g., the organization and morphology of the supporting cells and hair cells, the mechanical stimulation of the ear, etc.) are conserved and analogous with other species (Haddon and Lewis, 1996;Nicolson, 2005;Abbas and Whitfield, 2009;Baxendale and Whitfield, 2016). Moreover, it has been shown that approximately 70% of human genes have, at least, one zebrafish ortholog (Howe et al., 2013) and a number of genes required for the hearing function in zebrafish have been also associated with auditory defects in mammals (i.e., mice and humans) further supporting their role across different species (Coimbra et al., 2002).
Starting from these considerations, we developed a combined strategy, based on: (1) Targeted Re-Sequencing (TRS) of a panel of 46 ARHL candidate genes (Morgan et al., 2018) in 464 patients, followed by (2) zebrafish models of the most interesting variants identified.
Here, we report the results obtained using this approach demonstrating the presence of a SLC9A3R1 (also named Na + /H + Exchange Regulatory Cofactor, NHERF1) pathogenic variant in two unrelated ARHL patients and finally that this variant has deleterious effects in a zebrafish knock-in (KI) model. These findings, together with a previous study demonstrating the involvement of Slc9a3r1 in hearing loss in mouse (Kamiya et al., 2014), allow us to propose SLC9A3R1 as an extremely promising new ARHL candidate gene.

Human
The study was reviewed and approved by the Ethics Committee of the Burlo Garofolo children's hospital in Trieste (Italy) (2007 242/07). Written informed consent was obtained from each participant to the study and all the research was carried out according to the ethical standards defined by the Helsinki declaration.

Zebrafish
Zebrafish procedures were performed in accordance with Spanish and European laws, guidelines and policies for animal (Real Decreto 1201/05 (BOE 252, October 21, 2005 and European Directive 2010/63/EU on the protection of animals used for scientific purposes of October 20, 2010). This project was approved by the institutional Ethical Committee for Animal Experimentation of the PRBB, where ZeClinics conducts all experimental work (approval number CEA-OH/9421/2).

Patients Recruitment
A large cohort of 464 ARHL patients coming from inbred (Friuli Venezia Giulia-FVG-Cohort, Carlantino Cohort) and outbred (Milan, Naples, Trieste) Italian populations was analyzed by TRS. The cohort consists of 258 males and 206 females, all aged over 50 and with high frequencies bilateral HL developed around the 5th decade of life [high-frequency pure tone average (PTAH) > 40dB]. No disorders in the external and middle ear, as well as any vestibular problem or syndromic features were present.
Furthermore, 350 healthy individuals (all aged ≥ 50 y.o with PTAH ≤ 25dB) coming from five out of six villages of the FVG cohort (Resia, Sauris, Clauzetto, San Martino del Carso, Illegio) were used as an internal control.

Targeted Re-sequencing
A total number of 46 ARHL-candidate genes, including SLC9A3R1, were sequenced using Ion Torrent PGM TM (Life-Technologies). Genes were selected according to data from (a) GWAS meta-analyses on isolated and outbred populations of European and Asian ancestry (Girotto et al., 2011Wolber et al., 2014;Vuckovic et al., 2015); (b) literature updates and (c) animal models (Morgan et al., 2018).
Briefly, 10 ng of genomic DNA were used to construct DNA libraries using Ion AmpliSeq Library Kit 2.0 (Life Technologies). Template Ion Sphere Particles were prepared using Ion PGM Template OT2 200 kit and a single end 200 base-read sequencing run was carried out using Ion PGM sequencing 200 kit v2 (Life Technologies), on Ion Torrent PGM TM (Life Technologies). Ten indexed patients' libraries were sequenced simultaneously on each Ion 318 Chip. Sequencing data were then analyzed according to the Ion Torrent SuiteTM v3.6. The annotated SNVs\INDELS were evaluated according to several in silico predictor tools (SIFT, Polyphen2, MutationTaster, LRT) (Ng and Henikoff, 2003;Chun and Fay, 2009;Schwarz et al., 2010;Adzhubei et al., 2013) and conservation across species (PhyloP) (Pollard et al., 2010). Variants were classified as ultra-rare (MAF < 0.001), rare (MAF < 0.01) or common (MAF > 0.01) based on the frequencies reported NCBI dbSNP build142 1 as well as in 1000 Genomes Project 2 , NHLBI Exome Sequencing Project (ESP) Exome Variant Server 3 , ExAC Browser 4 and gnomAD browser 5 . Finally, those variants most likely to be disease causing (i.e., rare and ultra-rare variants predicted as damaging by all in silico predictor tools) were analyzed by Sanger sequencing and tested in controls.

SNP Genotyping
SLC9A3R1 variant was checked in 350 healthy controls using a TaqMan SNP genotyping assay (Assay ID: C_164690872_10, Thermo Fisher Scientific). Reactions were performed according to the manufacturer's instructions. Data were analyzed with the Taqman Genotyper Software (Thermo Fisher Scientific).

Mutational Protein Modeling
The 3D structure of PDZ2 domain of the NHERF1 protein was retrieved from protein data bank [ID: 2KRG, (Bhattacharya et al., 2010)]. This structure was considered as wild type (WT) and used in discovery studio [DS, (Accelrys Inc., San Diego, CA, United States)] to produce a mutant model (R180Q) as previously described (Krishnamoorthy et al., 2011).

Molecular Dynamics Simulations
The 3D structures of WT and mutant were used in the Groningen machine for chemical simulations (GROMACS) to perform molecular dynamics (MD) simulations. The atoms of the systems were prepared with GROMOS96 force field (van Gunsteren, 1996;Van Der Spoel et al., 2005;Hess et al., 2008). To solvate the protein, the SPC3 water model was used within a cubic box sized 1.5 nm (Berendsen et al., 1981). The counter ions were added to neutralize the systems and the periodic boundary conditions were applied in all directions. The prepared systems consist of 30202 and 30204 atoms in total. The LINCS (Hess et al., 1997) algorithm was used to constrain all the bond lengths and the SETTLE (Miyamoto and Kollman, 1992) algorithm was applied to constrain the geometry of the water molecules. A twin range cut-off was used: 0.8 nm for Van der Walls and 1.4 nm for electrostatic interactions, for managing long range interactions. The steepest descent algorithm was applied to energy minimize the systems with the tolerance of 2000 Kj/mol/nm. Subsequently, these structures were pre-equilibrated with 100 ps simulation before performing the production MD simulation for 25 ns with a time-step of 2 fs at constant temperature (300 K), pressure (1 atm) and number of particles, without any position restraints (Berendsen et al., 1984). The structures were collected at regular interval for every 100 ps to trace the trajectories and the tools in PyMOL 6 , DS and GROMACS were utilized for analyzing the structures and their molecular interactions.

Cluster Analysis and Surface Map
To represent the structure from the MD simulations, each trajectory with 2500 structures were classified into clusters based on their structural deviations. A structure from the top ranked cluster was chosen for representation as it was frequently occurring conformation. The representative structures were further used in DS for mapping the surface and to analyze the modifications on the binding surface (potential sites for binding partners).

In vitro Molecular Cloning
The impact of the identified mutation on mRNA and protein levels was tested by transient transfection in HEK 293 cells using expression clones containing either the wild-type (WT) or the mutant cDNA. cDNAs were cloned into a pCMV6-Entry vector (Origene, Rockville, MD, United States), Myc-tagged.
The calcium phosphate transfection method was used (Kingston et al., 2003). Forty-eight hours after transfection total cell proteins and RNAs were prepared and analyzed by Western Blot (WB) and quantitative Real Time PCR (qRT-PCR), respectively.

Western Blot Analysis
For protein analysis, HEK 293 cells were lysed in IPLS buffer (50 mM Tris-HCL pH7.5, 120 mM NaCl, 0.5 mM EDTA and 0.5% Nonidet P-40) supplemented with proteases inhibitors (Roche). After sonication and pre-clearing, protein lysate concentration was determined by Bradford Assay (Bio-Rad). An 8% polyacrylamide gel was used for protein electrophoresis. After blotting, membranes were blocked with 5% skim milk in Tris-buffered saline, 0.1% Tween 20 (TBST) and then incubated with primary c-Myc Antibody 9E10 monoclonal (Santa Cruz) overnight. Secondary antibodies [anti-mouse antibody (Santa Cruz)] were diluted in blocking buffer and incubated with the membranes for 45 min at room temperature. Proteins were detected with the ECL detection kit (GE Health Care Bio-Sciences).
Housekeeping proteins (e.g., β-actin or Hsp90) were used as an internal control for protein loading as well as for reference in the WB analysis.

Quantitative Real-Time PCR (qRT-PCR)
RNA was extracted from cell pellets using High Pure RNA isolation Kit (Roche). Total RNA (1 µg) was reverse transcribed to cDNA using Transcriptor First Strand cDNA Synthesis kit (Roche). qRT-PCR was performed using standard PCR conditions in a 7900HT Fast Real Time PCR System (Applied Biosystems) (i.e., 95 • C for 10 , 40 cycles of 95 • C for 15 and 60 • C for 1 , followed by a dissociation stage of 95 • C for 15 , 60 • C for 15 and 95 • C for 15 ) with Power SYBR Green PCR Master Mix (Thermo Fisher Scientific). Gene-specific primers were designed by using Primer3Web software 7 . All experiments were performed in biological triplicate. Expression levels have been standardized to Neo gene expression and all data have been analyzed using the 2 − CT Livak Method (Livak and Schmittgen, 2001).

Zebrafish Husbandry
Zebrafish KI lines used in this study were generated, analyzed phenotypically and maintained at ZeClinics (Barcelona, Spain). Two different zebrafish transgenic lines have been used as a background for the KI lines generation: (1) Tg(brn3c:mGFP), which is a stable reporter line where a GFP tagged to the membrane is expressed under the control of the brn3c promoter (Xiao et al., 2005), thus specifically labeling a subset of retinal ganglion cells, hair cells of the inner ear and the lateral line neuromasts (Di Donato et al., 2013).
(2) Tg(isl3:GFP), which is a stable reporter line where a GFP is expressed in the cytoplasm of sensory neurons under the control of the isl3 promoter (Pittman et al., 2008).
Embryos were maintained in petri dishes at 28.5 Celsius. Developmental stages were evaluated as hours and days post-fertilization (hpf and dpf). All zebrafish experimental protocols were approved by the Generalitat de Cataluña.

Gene Expression in Zebrafish Larvae
Gene expression in Zebrafish larvae (5 dpf) was performed by whole mount in situ hybridization (ISH). At this stage, inner ear is developed to allow functional hearing, which implies hair cells involved in that function are mature. Synthesis of antisense RNA and whole-mount in situ hybridization were performed as previously described (Thisse et al., 2004).
Specific riboprobes intended to recognize slc9a3r1 mRNA were designed (Slc9a3r1_Fw:ATGTCCAGCGACCTCAGGCC; Slc9a3r1_Rv_T7: TAATACGACTCACTATAGGGCTCCAGTCC ATCTGCGGAGCTC). cDNAs were amplified by PCR, using Expand High Fidelity PCR System (Roche), from a custom Zebrafish cDNAs library obtained by RT-PCR from (SuperScript III Rev Transcript kit, Invitrogen) a mRNA pool coming from 5 dpf Zebrafish larvae [Trizol mRNA extraction protocol, Trizol (Sigma-Aldrich)]. A T7 sequence linker in reverse primers was included to directly use the synthesized PCR products as templates to amplify the reverse digoxigenin-labeled riboprobe to be used for ISH.
In order to check gene expression in hair cells, marked with GFP through its activation by the Brn3c promoter, secondary antibody against GFP was used (rabbit anti-GFP, Torrey Pinnes; 1:400).
Stained embryos were processed for imaging through 2 different methods: (1) Whole embryos were imaged with a bright field stereoscope to determine the overall expression pattern.
(2) Transversal sections were acquired though a cryostat (20 µm width) to determine precisely the mRNA location within the inner ear across three different anterioposterior positions. Gene expression pattern was compared with the position of hair cells, marked with GFP through its activation by the Brn3c promoter. Sections were imaged with Leica DM5 light/fluorescent microscope at 40× magnification.

Generation of slc9a3r1 Knock-In (KI) Zebrafish Line
To generate the KI, we designed a strategy where double-strand breaks (DSBs), generated by CRISPR technology, are repaired through homologous recombination (HR) DNA repair mechanism (Cong et al., 2013). HR has been enhanced through the co-injection of Rad51 protein (Takayama et al., 2017) and providing a single-stranded oligonucleotide (ssOligo) that contains two homologous arms flanking the desired mutation.

Synthesis of slc9a3r1 sgRNA
After assessing Slc9a3r1 protein conservation with its human orthologous, single guide RNAs (sgRNAs) were synthesized in vitro from double strand DNA that contained a T7 promoter region, a specific CRISPR recognition site (20N spacer) and a guide sequence that allows recognition by CAS9. This double strand DNA was synthesized by the partial annealing of a constant oligo (80 bps) and a genespecific oligo (60 bps). After annealing, DNA was fill-in by T4 polymerase (Gene-specific oligo: TAATACGACTCACTATA-N20-GTTTTAGAGCTAGAAATAGCAAG; Constant oligo: AA AAGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACGG ACTAGCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC).

slc9a3r1 KI Sound Response
Six day-old larvae obtained by pairwise mating of adult tg(brn3c:mGFP;R180Q-slc9a3r1 +/R180Q ) were tested to determine their ability to respond to sound stimuli at different frequencies, measuring the acoustic startle response (Bhandiwad et al., 2013). The EthoVision XT 12 software and the DanioVision device from Noldus Information Technologies (Wageningen, Netherlands) were used. This closed system consists of a camera placed above a chamber with circulating water and a temperature sensor that is set at 28 • C. A 96-wells plate is placed in the chamber, which can provide different stimuli (light/dark environment, tapping, sound) controlled by the software. For each sound frequency tested, larvae were left for 10 min in dark for acclimation, 2 min with the light on (natural locomotor behavior of zebrafish is active in dark and immobile in light) and then the sound stimulus was provided. Finally, a tapping was given to the 96-wells plate in order to verify the absence of general locomotor defects. The system provides tapping stimuli that range from 1 to 8 arbitrary units (AU). An intensity of 3 AU was chosen for the experiments (1 and 2 were too weak to induce any response in normal larvae and higher intensities could over stimulate larvae and mask possible subtle effects).
The final readout to assess the capacity to respond to sound was the percentage of stimuli responding 6 day-old larvae.
Four different sound stimuli (300, 325, 350, and 375 Hz) that showed the highest percentage of responding larvae were tested.
Then genomic DNA was extracted from single larvae to identify their specific genotype. Purple rectangle outlines R180 conservation among human and zebrafish. (B) Scheme displaying the nucleotide and amino acid translated sequence present at the target genome (below) and the ssOligo (above) used for modifying the target sequence. Below, green arrow in genome sequence outlines the sgRNA site. Above, red squares in ssOligo outline homologous recombination right and left arms. (C) Scheme displaying the PCR strategy used to identify potential mutants. (D) Example of R180Q-Slc9a3r1 mutant identified through the comparison of Sanger sequences obtained from a wild type (above) and a heterozygous individual (below). Both sequences display nucleotide and translated amino acid sequence. Above, targeted amino acid R180 is displayed in green. Below, modified amin oacid Q180 is displayed in red and modified nucleotides are displayed in red and lower case.
To further verify the normal locomotion activity of slc9a3r1 R180Q/R180Q larvae, we analyze their movement during the entire trial and compare it with the wild-type.

Hair Cell Number
In order to test the number of functional hair cells in the saccular macula (the sensory patch responsible for hearing) tg(brn3c:mGFP;R180Q-slc9a3r1 +/R180Q ) adult carrier were incrossed and 6 dpf progeny were intra-ear injected with FM1-43FX, a vital dye that penetrates in only mature and fully functional hair cells through calcium channels. In particular, 6 dpf larvae were anesthetized using Tricaine (20 µM) and 1 nl of FM1-43FX (stock solution 300 µM) was microinjected into the inner ear lumen (intra-ear). Once injected, larvae were fixed with 4% PFA 2 h at room temperature and then incubated with PBS-Triton 2% O/N, in order to dissolve the otoliths. The larvae were then embedded on their lateral sides in 1% low melting point agarose and the hair cells expressing membrane GFP (mGFP) and that internalized the FM1-43X were imaged using a SP8 Leica confocal microscope and Z-stacks spanning the entire saccular macula were taken (one z-plane imaged every 1 µm). Raw data was analyzed with FIJI software (Schindelin et al., 2012).
Confocal imaging of the saccular macula of slc9a3r1 +/+ , slc9a3r1 +/R180Q and slc9a3r1 R180Q/R180Q was done and images where analyzed. An in-house FIJI-based macro allows us to 3D reconstruct the entire saccular macula, measure the GFP (brn3c) and DsRed (FM1-43FX) signals and estimate the total number and the number of mature hair cells.

Hair Cells Morphology and Polarity
To test sensory hair cells morphology and orientation, 6 dpf larvae from tg(brn3c:mGFP;R180Q-slc9a3r1 +/R180Q ) adult pairwise mating were cryosectioned and an immunohistochemistry against the acetylated-tubulin, which labels the kinocilium, was performed. This staining, together with the membrane GFP signal help to detect alteration in hair cell morphology or in the orientation of the hair cell bundles.
Six dpf larvae were fix with 4% PFA 2 h at room temperature and then incubated with PBS-Triton 2% O/N in order to dissolve the otoliths. After several washes with 0.1%-Tween-20 PBS, larvae heads were incubated 1 h in 15% sucrose (in PBS) then in 15% sucrose/7.5% gelatin and placed in cryomold in the desired orientation (transversally for hair cell morphology, laterally for hair cell bundle orientation), while tails were used for genotyping. Blocks were frozen in 2-Methylbutane for tissue preservation and cryosectioned at 30 µm on a Leica CM 1510-1 cryostat. Sections were collected on Superfrost slides. Then the inner ear containing slides were blocked with blocking solution [0.1%Tween-20 in PBS (PBT), 2% bovine serum albumin (BSA) and 10% goat serum] for 1.5 h at RT. Mouse anti-acetylated-tubulin (Sigma; 1:1000) and rabbit anti-GFP (Torrey Pines; 1:400) were incubated overnight at 4 • C in blocking solution. After washing with PBT for the whole day, anti-rabbit Alexa488 and anti-mouse Alexa648 (Invitrogen; 1:400) were incubated overnight at 4 • C in blocking solution. Slides were then mounted with mowiol and imaged using a SP8 Leica confocal microscope and Z-stacks spanning the entire saccular macula were taken (one z-plane imaged every 1 µm). Raw data were analyzed with FIJI software in order to generate 3D reconstruction of the saccular macula or analyze the hair cell bundle orientation. To quantify and generate the graphs showing the ciliary orientation a Python-based in-house program was used.
For every single cell, we determine the orientation based on the position of the kinocilia and plot every cell on an Angle's plot (Inoue et al., 2013) based on their polarization and finally generate a Polar plot to represent the percentage of the cells oriented in the dorsal, ventral, anterior, and posterior direction.

Stato-Acoustic Ganglion Characterization
The number of neurons in the posterior stato-acoustic ganglion (SAG) (which innervate the saccular macula) was tested.
Adults Tg(isl3:mGFP;R180Q-slc9a3r1 +/R180Q ) were incrossed to obtain progeny. Six dpf larvae were fixed with 4% PFA 2 h at room temperature and then incubated with PBS-Triton 2% O/N in order to dissolve the otoliths. The larvae were then embedded on their sides in 1% low melting point agarose and the sensory neurons expressing GFP were imaged using a SP8 Leica confocal microscope and Z-stacks spanning the entire posterior SAG were taken (one z-plane imaged every 1 µm). Raw data were analyzed with FIJI software (Schindelin et al., 2012).
Using a FIJI-based in-house macro the entire posterior SAG was 3D reconstructed the number of the sensory neurons was calculated.

Otolith Analysis
To test the size of the otolith, a crystalline structure contacted by the hair cell ciliary array and essential to drive the movement of the hair cell bundle, we incrossed tg(brn3c:mGFP;R180Q-slc9a3r1 +/R180Q ) adult carriers and 6 dpf progeny was anesthetized using Tricaine (20 µM), placed on their side, and imaged with Leica M165 Stereoscope before genotyping.
The images were analyzed with FIJI software and the area of the posterior otolith was measured.
The significance threshold was set at 0.05.

RESULTS
A total of 464 ARHL Italian patients were screened with a TRS panel of 46 ARHL candidate genes (Morgan et al., 2018). A mean of 33.5 megabases of raw sequence data were available for each subject. The coverage, on 95% of the targeted region, was at least 20-folds, with a 270-fold mean-depth total coverage. On average 333 single nucleotide variants (SNVs) and small insertions/deletions (INDELs) were called for each patient. After applying the filtering procedure described in the Section "Materials and Methods, " a missense variant at the heterozygous state, c.539G > A, p.(R180Q) in SLC9A3R1 (NM_004252.4) gene was detected in two unrelated ARHL male patients (582130, aged 79, and 593486, aged 63). It affects the PDZ2 domain of NHERF1 protein (Figure 2A) and all in silico predictor tools classified it as damaging (Supplementary Figure S1B). Both subjects come from the same isolated village located in North Eastern Italy (Friuli Venezia Giulia-FVG cohort) and they were recruited thanks to a huge project on isolated community in Italy (INGI consortium) based on a volunteer study. Both patients show severe to profound high-frequency hearing loss resembling a sensory ARHL phenotype and the age of onset was 48 y.o.
(582130) and 51 y.o. (593486) (Supplementary Figures S1A,B We confirmed the absence of the variant in all the available relatives of the patients (2 healthy daughters of 582130) as well as any other ARHL patient and control from the same small community (15 ARHL cases and 37 healthy controls). Moreover, 350 additional controls coming from the other five villages of the FVG cohort have been tested, confirming the absence of the SLC9A3R1 variant.

Protein Simulations and Conformational Analysis
Protein modeling and molecular dynamics (MD) simulations were used to test the impact of the mutation at molecular level. The secondary structural details of the PDZ2 domain ( Figure 2B) show that residue R180 is located in the surface of a well-structured β sheet ( Figure 2C) underlying its importance for the protein structure and binding activity. MD simulations revealed that secondary structural elements of both the wildtype (WT) and the mutant are stable during simulations. However, a major change was observed in the side chain of Q180 compared to R180, which further affects their hydrogen bonding partners (R180: H169, Q196 and N167 and Q180: Frontiers in Genetics | www.frontiersin.org I179 and N167) in the nearby region (Figures 2C,D). The outcome is a series of local network changes leading to rearrangements of neighboring loops of the mutational spot (Figures 2C,D, see the triangles for the impact) and following modifications on the surface of the domain itself whose entire conformation is significantly modified (Figures 2E,F). This finding correlates with the change in the total number of hydrogen bonds from 67 (WT) to 62 (p.(R180Q)) ( Figure 2G). Overall, the modeling suggests a significant conformational change that can alter the normal function of the domain and the protein itself.

In vitro Experiments: qRT-PCR, Western Blot, Immunostaining
In order to understand the effect of the mutated allele on RNA and protein stability, HEK 293 cells have been transfected with expression vectors containing either the WT or the mutant cDNA. qRT-PCR and WB analysis 48 h after the transient transfection did not reveal any difference in the levels of both mRNA and protein between WT and mutant (Supplementary Figure S2).
Moreover, immunostaining assay on transfected HeLa cells revealed no differences in the cellular localization of the mutated protein compared to the WT (Supplementary Figure S2).

In situ Hybridization
Prior to the generation of the KI zebrafish model, slc9a3r1 expression was tested in 5 dpf zebrafish larvae, by whole mount in situ hybridization (Figure 3). slc9a3r1 mRNA was detected broadly across the larvae, but its expression was enriched in the hematopoietic lineage, liver primordium ( Figure 3B) and in the inner ear (Figures 3C-E).

KI Hearing Phenotype
CRISPR/Cas9 technology has been used to generate the KI model (R180Q-slc9a3r1), in which has been introduced a precise nucleotide modification mimicking the human mutation (Figures 1A-D). Overall 135 six days post Frontiers in Genetics | www.frontiersin.org fertilization (dpf) larvae obtained by pairwise mating of adult tg (brn3c:mGFP;R180Q-slc9a3r1 +/R180Q ), were tested to determine their capacity to respond to sound stimuli at different frequencies (300, 325, 350, and 375). In general, for all stimuli, the response of slc9a3r1 R180Q/R180Q KI larvae (n = 31) was lower as compared to that of heterozygous slc9a3r1 +/R180Q (n = 47) and WT animals (n = 57). Interestingly, a statistical significant difference in sound perception at 325 Hz was definitely evident between the WT and the homozygous mutants (P = 0.0132), and a worsening hearing perception was also noticed in the heterozygous animals ( Figure 4A) while no variation in the ability to respond to mechanical stimuli was observed in comparison to the wild type siblings.

Behavioural Test
Locomotion was studied after light/darkness stimuli to verify that response differences to sound stimuli among WT and mutant siblings were not promoted by general defects in locomotive behavior.
Wild type Zebrafish larvae show a relatively high swimming activity in darkness when compared to the activity in light (de Esch et al., 2012). No differences have been observed for the slc9a3r1 +/R180Q and slc9a3r1 R180Q/R180Q models (see Figure 4B). The graph shows normal behavior and response to visual stimuli of homozygous slc9a3r1 mutants compared to heterozygous and WT, demonstrating that the impairment observed in sound response is specific to hearing function and not to defects in central nervous system.
To further investigate the role of this gene/mutation, several functional tests in the KI model have been performed and, among them, a significant reduced otolith size has been detected.

Measurement of the Otolith Size
Despite fish vestibular-acoustic system presents several functional elements that are very similar to the human ones (e.g., hair cells, supporting cells, etc.), they do not have a counterpart of the cochlea. Moreover, in fish inner ear the otoliths are essential structures for balance and hearing. Otoliths are crystalline structure contacted by the hair cell ciliary array and are essential to drive the movement of the hair cells bundle and start the electrical signaling in case of sound (Riley and Phillips, 2003). The alteration of the size, shape or the absence of this structure could have striking effects on the hearing ability (Inoue et al., 2013). For this reason, we looked for any difference in the otoliths of WT compared to that of the mutated animals.
Results shows a 5,8% reduction of the otolith size in slc9a3r1 +/R180Q compared to WT (P = 0.0014) and an even more increased 10.2% reduction in the size of the otolith of slc9a3r1 R180Q/R180Q compare to WT (P < 0.0001) (Figures 5A,B). These results thus explain the hearing impairment observed in the slc9a3r1 R180Q/R180Q KI larvae.
Other functional tests demonstrated no alterations in the hair cells number and morphology as well as in the count of neurons, as reported below.

Hair Cells Number and Morphology
The reduced capacity of slc9a3r1 mutant larvae to respond to acoustic stimuli could have several explanations such as: (a) a reduced number or a defective maturation of the sensory hair cells leading to not functional hair cells in the saccular macula, (b) an abnormal morphology or an altered orientation of the sensory hair cells. As regards to (a), 9 WT animals compared to 10 heterozygous mutants and 6 homozygous mutants were tested. No differences in the total number (35.79 ± 6.11 cells for WT, 38.32 ± 5.11 cells for R180Q/R180Q, 43.18 ± 6.53 for ± total cells) or fully mature hair cells (31.16 ± 5.22 cells for WT, 32.26 ± 3.25 cells for R180Q/R180Q, 35.56 ± 3.27 for +/R180Q total cells) among slc9a3r1 R180Q/R180Q, slc9a3r1 +/R180Q mutants and their Wt siblings were observed (Figures 6A,B). As regards to (b) 6 dpf larvae (slc9a3r1 +/+ n = 9, slc9a3r1 +/R180Q n = 10, slc9a3r1 R180Q/R180Q n = 6) were cryosectioned and an immunohistochemistry against the acetylated-tubulin, that labels the kinocilium, was performed. Data analysis revealed no differences in hair cell body morphology or size, hair cell bundle organization or kinocilia/stereocilia length suggesting that the KI do not alter the hair cell morphology (Figure 6C) but induce a reduced response to sound stimuli due to other events. Moreover, no statistically significant differences were observed in the orientation of the hair cells within the sensory patch (see Figure 6D), suggesting that the hearing impairment in KI is not due to hair cell altered polarization. All together these findings indicate the KI hearing impairment is not caused by alterations in the hair cell population.

Count of Neurons
Since the hair cells of slc9a3r1 mutants do not show any defect, the number of neurons in the posterior stato-acoustic ganglion (SAG), which innervate the saccular macula and could be related to a hearing loss, was checked.

DISCUSSION
Despite significant research efforts carried out, the genetics risk factors involved in ARHL are still mainly unknown and only few ARHL susceptibility genes have been detected so far.
Zebrafish mutants can show morphological and functional defects similar to those of other species (Ernest et al., 2000) and their vestibular-acoustic system is very similar to that of human beings, thus representing a cost and time-saving alternative to mouse model allowing a rapid assessment of inner ear FIGURE 4 | Auditory/sensory response and light/dark behavioral response in Zebrafish larvae. (A) Response of slc9a3r1 KI animals compared to the wild-type ones to all the different sound stimuli tested (number of tested animals: slc9a3r1 +/+ N = 57, slc9a3r1 +/-N = 47, slc9a3r1 R180Q/R180Q N = 31). Statistically significant difference indicated with "*" ( * p-value < 0.05). (B) Behavior and response activity to visual stimuli of slc9a3r1 R180Q/R180Q compared to slc9a3r1 +/+ and slc9a3r1 +/-. Tapping intensity is indicated as arbitrary unit (AU).  function (Baxendale and Whitfield, 2016;Leventea et al., 2016). Furthermore, a number of genes required for hair-cell function in the zebrafish have been already associated with auditory defects in mice and humans, thus revealing their conservation across species (Coimbra et al., 2002). Using a combination of genetics/genomics studies followed by functional studies we deeply investigated the role of SLC9A3R1 in ARHL. In particular, the identification of an ultra-rare human mutation [c.539G > A, p.(R180Q)] in two ARHL patients from an isolated community of the Italian Alps led to a series of functional studies including the generation of the first KI zebrafish model for hearing phenotype developed so far demonstrating the pathogenic effect of this variant and thus a possible relevant role of SLC9A3R1 in ARHL. SLC9A3R1 encodes the Na + /H + exchange regulatory factor 1 (NHERF1) protein (also called Ezrin-radixin-moesin-binding protein of 50 kDa, Ebp50), which belongs to NHERF family of scaffolding proteins. These proteins play key roles in mediating a number of transports and signaling phenomena, interacting with different molecular partners. Interestingly, several protein-protein interactions involving NHERFs take place in the cochlea, and accordingly the Nherf1 KO mouse shows hearing defects, due to hair cells anomalies (Kamiya et al., 2014). NHERF1 interacts with many different molecular partners through its PDZ domains. The missense variant here identified affects the second PDZ domain of the protein (PDZ2) and involves the R180 residue, which is essential for protein-protein interaction, as demonstrated by Mamonova et al. (2012). Accordingly, our protein modeling confirmed this data, suggesting that the p.(R180Q) mutation induces changes in the network of interactions by altering the total number of hydrogen bonds and thus compromising the binding activity of the protein itself. Moreover, the KI zebrafish model further support the pathogenic effect of the p.(R180Q) variant, showing that it causes a reduction in size of the saccular otolith, a structure corresponding to human otoconia whose abnormalities might cause vertigo, dizziness and isolated hearing loss. Zebrafish otoliths are biomineralised structures required for balance and hearing. While the utricular otolith is mainly involved in balancing, the saccular one is implicated in the sense of hearing (Stooke-Vaughan et al., 2015). Interestingly, genetic screenings analyzing several zebrafish mutants with hearing phenotype showed abnormalities in the otolith volume and formation, and revealed that otolith defects are not necessarily associated with other morphological changes of the inner ear (Schibler and Malicki, 2007) as it happens to our slc9a3r1 KI model. This finding clearly justifies the absence of any balance problem in our KI model and this is in agreement with the ARHL phenotype of the two ARHL patients here described (i.e., high-frequencies hearing loss without any vestibular symptoms). Furthermore, data obtained for SLC9A3R1 gene are also in agreement with those reported for other genes involved in the otoconia development such as Claudin and Tectorin, in which vestibular signs are frequently absent suggesting a functional compensation by other proteins (Lundberg et al., 2015).
Molecular mechanism underlying KI morphological changes still need to be understood, however, it has been shown that different Na + /K + -ATPase genes are essential for otolith formation in zebrafish (Blasiole et al., 2006). In this light, a relevant role might be played by the Na + /K + -ATPase α1 subunit that binds the PDZ2 domain of NHERF1 in opossum kidney cells (Mahon et al., 2003). Since the strong homology of this subunit between zebrafish and opossum sequences (∼71%) it is reasonable to assume that this binding could also happen in zebrafish. Taking into account that our KI model carries a mutation in the PDZ2 domain in a residue critical for protein-protein interaction, it is highly probable that this function could be severely compromised (Voltz et al., 2001;Blasiole et al., 2006) thus leading to an altered otolith formation and a subsequent hearing loss. Furthermore, this mechanism of action could also explain the phenotype's difference between the KO mouse (that shows anomalies in the shape of the hair cells bundles) and our KI zebrafish model, which specifically reproduces the variant identified in our patients.
Despite carrying the same variant, the KI zebrafish model displays early-onset hearing loss, while the two patients here described are affected by ARHL. This difference might be explained by several hypotheses.
First of all, even though animal models are extremely useful for recapitulating many human diseases, it is important to keep in mind that in some case the phenotypes' spectra may be different between the two species. Having said that, the slc9a3r1 +/R180Q model does not show a statistically significant reduction of the hearing thresholds, despite presenting an altered otolith's volume compared to the WT. In this light, considering that the human patients carry the SLC9A3R1 variant at the heterozygous state, it is reasonable to think that also in this case, despite the presence of an anatomical defect, the hearing phenotype does not occur at an early age.

CONCLUSION
Our results confirmed the role of SLC9A3R1 in the hearing system function and development, strengthening the previously reported data of the mouse KO model (Kamiya et al., 2014). Moreover, since the KI zebrafish model generated in the present study carries the specific variant identified in two ARHL human patients, this work provides evidence that SLC9A3R1 gene might play a role in human ARHL. Overall these findings suggest that SLC9A3R1 could act as a likely Mendelian gene, with a large effect size, in the etiopathogenesis of the late-onset hearing loss form detected in our patients. In this light, larger patients cohorts should be screened at molecular level for the presence of SLC9A3R1 gene variants/mutations to understand the overall contribution of this gene to the disease itself. Moreover, the identification of other ARHL patients carrying SLC9A3R1 variants would be also essential for the development of new plans for disease's prevention and treatment.

DATA AVAILABILITY
Sequencing data are available at the European Genome-phenome Archive (EGA) at the following link https://ega-archive.org/ studies/EGAS00001003072.

AUTHOR CONTRIBUTIONS
GG: study design assessment, data production and analysis, and writing the manuscript. AM: data production and analysis and writing the manuscript. NK: protein modeling and molecular dynamics simulation. MC: raw data analysis and quality control. BM: raw data analysis and quality control. SB: in vitro studies. MLB: data production. MDS: support of functional studies and conceiving the experiments. PG: conceiving the experiments and writing the manuscript.

FUNDING
This research was supported by RBSI14AG8P-SIR2014 to GG (http://sir.miur.it/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
We gratefully acknowledge ZeClinics company for technical contribution in Zebrafish models.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2019.00142/full#supplementary-material FIGURE S1 | Clinical features and DNA sequence chromatograms of the two ARHL patients. (A) Audiograms of patients 582130 and 5934856; the downward slope indicates that high frequencies are severely affected. (B) The figure displays DNA sequence chromatograms showing the nucleotide variant identified in the two ARHL patients and all the details of the mutation. rs ID, Reference SNP ID from dbSNP database. Frequency, frequency of gnomAD Exome all is reported. Predictor tools, GERP++ (higher number is more conserved, >0 is generally conserved), PhyloP (Pathogenicity score: conserved > 0.95, not conserved < 0.95), Polyphen-2 (Pathogenicity score: probably damaging: D, possibly damaging: P, benign: B), SIFT (Pathogenicity score: D: disease causing, N: polymorphism, P: polymorphism automatic), MutationTaster (Pathogenicity score: D:disease causing, N:polymorphism, P: polymorphism automatic), CADD Phred (Pathogenicity score: >10 predicted to be deleterious).