Introduction
============

Age-related Hearing loss (ARHL) is the predominant sensory impairment in the elderly ([@B24]; [@B27]), affecting millions of people worldwide ([@B37]). Moreover, projections suggest that, in the next decade, the number of patients will double, largely due to an increased lifespan ([@B10]).

Age-related Hearing loss is characterized by bilateral and progressive hearing loss that usually starts at the high frequencies ([@B24]) causing communication difficulties associated with cognitive decline, social isolation and depression ([@B10]). It is a complex disease in which genetic and environmental risk factors (i.e., noise, smoking, alcohol, etc.) interplay ([@B9]; [@B40]; [@B54], [@B55]). 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) ([@B17]; [@B52]; [@B18], [@B19]; [@B58]; [@B56]), 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 ([@B57]). 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 ([@B20]; [@B39]; [@B1]; [@B3]). Moreover, it has been shown that approximately 70% of human genes have, at least, one zebrafish ortholog ([@B23]) 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 ([@B12]).

Starting from these considerations, we developed a combined strategy, based on: (1) Targeted Re-Sequencing (TRS) of a panel of 46 ARHL candidate genes ([@B36]) 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 ([@B26]), allow us to propose *SLC9A3R1* as an extremely promising new ARHL candidate gene.

Materials and Methods {#s1}
=====================

Ethics Statement
----------------

### 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 ([@B18], [@B19]; [@B58]; [@B56]); (b) literature updates and (c) animal models ([@B36]).

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) ([@B38]; [@B11]; [@B46]; [@B2]) and conservation across species (PhyloP) ([@B42]). 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](#fn01){ref-type="fn"}^ as well as in 1000 Genomes Project^[2](#fn02){ref-type="fn"}^, NHLBI Exome Sequencing Project (ESP) Exome Variant Server^[3](#fn03){ref-type="fn"}^, ExAC Browser^[4](#fn04){ref-type="fn"}^ and gnomAD browser^[5](#fn05){ref-type="fn"}^. 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, ([@B7])\]. 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 ([@B29]).

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 ([@B51]; [@B50]; [@B22]). To solvate the protein, the SPC3 water model was used within a cubic box sized 1.5 nm ([@B5]). 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 ([@B21]) algorithm was used to constrain all the bond lengths and the SETTLE ([@B35]) 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 ([@B4]). The structures were collected at regular interval for every 100 ps to trace the trajectories and the tools in PyMOL ^[6](#fn06){ref-type="fn"}^, 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 ([@B28]). 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](#fn07){ref-type="fn"}^. 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 ([@B31]).

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 ([@B59]), thus specifically labeling a subset of retinal ganglion cells, hair cells of the inner ear and the lateral line neuromasts ([@B15]).

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 ([@B41]).

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 ([@B49]).

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.

Once dissected, embryos were fixed in 4% paraformaldehyde (PFA) over night (O/N), and then dehydrate with increasing concentrations of Metanol 0,1% Tween20 in PBS (PBT) for long-term storage (25, 50, and 75%). Embryos were then rehydrated and treated with proteinase K (Sigma-Aldrich). Afterwards embryos were incubated with a hybridization mix \[i.e., 50% Formamide, 5× saline-sodium citrate (SSC), 0.1% Tween 20, Citric acid to adjust HM to pH 6.0 (460 μl of 1 M citric acid for 50 ml of mix), 50 μg/ml heparin, 500 μg/ml tRNA\] containing the riboprobe (dilution: 1:100) and then with antibody against digoxigenin \[Peroxidase-conjugated anti-DIG Fab (Roche), dilution 1:1000\]. Nitroblue Tetrazolium (NBT: 3,5 μl of 50 mg/ml NBT solution in 1 ml of alkaline Tris buffer) was used with the alkaline phosphatase substrate 5-Bromo- 4- Chloro-3-Indolyl Phosphate (BCIP: 4,5 μl of 50 mg/ml BCIP solution in 1 ml of alkaline Tris buffer).

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 ([@B13]). HR has been enhanced through the co-injection of Rad51 protein ([@B48]) 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 gene-specific oligo (60 bps). After annealing, DNA was fill-in by T4 polymerase (Gene-specific oligo: TAATACGACTCACTATA-N20-GTTTTAGAGCTAGAAATAGCAAG; Constant oligo: AAAAGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACGGACTAGCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC).

CRISPR Injection and Screening
------------------------------

Tg(brn3c:mGFP)or Tg(isl3:GFP) embryos were injected with a mix of sgRNA, Cas9 (mRNA), ssOligo and Rad51 human Protein. 20--25 embryos of 48 hpf were selected and pooled for genomic DNA extraction. Genotyping was performed using specific diagnostic primers, a pair of primers to amplify the surrounding area (SLC9A3R1_Fw/SLC9A3R1_Rv = 379 bps; SLC9A3R1_Fw: CCCTTTGTAGGACTCGAAGAACGAG; SLC 9A3R1_Rv: TGTTCCAAACTAAGCCAGAGCAGAAC) and one primer to discriminate the mutated region (SLC9A 3R1_KI_Fw/SLC9A3R1_Rv = 278 bps; SLC9A 3R1_KI_Fw: CTTCTTGCGCAAATGGTTGTTG) ([Figure 1C](#F1){ref-type="fig"}).

![CRISPR-Cas9 strategy for the generation of the R180Q-*slc9a3r1* K/I model. **(A)** Comparison of human and zebrafish SLC9A3R1 amino acid sequence. 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.](fgene-10-00142-g001){#F1}

The remaining injected animals were grown until sexual maturity. F1 animals were Sanger sequenced to identify positive targeted mutants generated through CRISPR KI. These animals were incrossed to obtain F2 homozygous larvae to perform their phenotypical analysis.

Genotyping
----------

Genomic DNA of single larvae from the different experiment was extracted using Extract-N-AmpTM Tissue PCR Kit (Sigma). For the *Slc9a3r1* KI a first PCR reaction was performed using the primers SLC9A3R1_Fw and SLC9A3R1_Rv and the following conditions: 94°C 2 min, 25× (94°C 15 s, 58°C 30 s, 72°C 1:10 min), 72°C 7 min, 4°C hold. Then, a nested PCR reaction was performed using the primers SLC9A3R1_Fw, SLC9A3R1_Rv, SLC9A3R1_INT_MUT_Fw (5′-CTCATCCACCGCCCTGATG-3′), SLC9A3R1_INT_WT_Rv (5′-CCGGCCAGTATATTCAAGCTGTC- 3′) and the same thermal protocol.

*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 ([@B6]). 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.

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 Characterization
--------------------------

### 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 ([@B45]).

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 ([@B25]) 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 ([@B45]).

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.

Statistical Analysis
--------------------

Behavioral and otolith differences between *slc9a3r1*^+/+^, *slc9a3r1*^+/R180Q^, and *slc9a3r1*^R180Q/R180Q^ were evaluated with an unpaired *T*-test with Welch's correction.

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 ([@B36]). 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](#F2){ref-type="fig"}) and all *in silico* predictor tools classified it as damaging ([Supplementary Figure S1B](#SM1){ref-type="supplementary-material"}). 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](#SM1){ref-type="supplementary-material"}) without vestibular signs or symptoms. Additional clinical data include: (1) patient 582130: normal BMI, high pressure treated with drugs. During his life he was affected by Tuberculosis (TBC), lung cancer and dislipidemia. Habits: no smoking and limited alcohol consumption; (2) patient 5934486: overweight, high pressure treated with drugs, asthma treated with drugs, sarcoidosis. Habits: no smoking and no alcohol consumption. No other clinical features were detected in the patients nor exposure to noise.

![Molecular dynamics (MD) simulations and conformational analysis of NHERF1. **(A)** Schematic representation of NHERF1 protein domains. **(B)** Secondary structural details of the PDZ-II domain. **(C,D)** Secondary structural elements of the representative structures of the WT and mutated proteins. **(E,F)** Modification of the entire domain's structure. **(G)** Representation of the total number of hydrogen bonds of the WT residue compared to the mutated one \[from 67 (WT) to 62 (p.(R180Q)\].](fgene-10-00142-g002){#F2}

The c.539G \> A, p.(R180Q) variant has a frequency of 0.002 in our cohort of 464 ARHL patients and it is described in gnomAD browser with a minor allele frequency (MAF) of 2.886e-5 (rs146832150) (date of access 21/08/2018).

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](#F2){ref-type="fig"}) show that residue R180 is located in the surface of a well-structured β sheet ([Figure 2C](#F2){ref-type="fig"}) underlying its importance for the protein structure and binding activity. MD simulations revealed that secondary structural elements of both the wild-type (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: I179 and N167) in the nearby region ([Figure 2C,D](#F2){ref-type="fig"}). The outcome is a series of local network changes leading to rearrangements of neighboring loops of the mutational spot ([Figure 2C,D](#F2){ref-type="fig"}, see the triangles for the impact) and following modifications on the surface of the domain itself whose entire conformation is significantly modified ([Figure 2E,F](#F2){ref-type="fig"}). This finding correlates with the change in the total number of hydrogen bonds from 67 (WT) to 62 (p.(R180Q)) ([Figure 2G](#F2){ref-type="fig"}). 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](#SM2){ref-type="supplementary-material"}).

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](#SM2){ref-type="supplementary-material"}).

*In vivo* Experiments: Zebrafish KI Model
-----------------------------------------

### *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](#F3){ref-type="fig"}). *slc9a3r1* mRNA was detected broadly across the larvae, but its expression was enriched in the hematopoietic lineage, liver primordium ([Figure 3B](#F3){ref-type="fig"}) and in the inner ear ([Figure 3C--E](#F3){ref-type="fig"}).

![*slc9a3r1* expression analysis by *in situ* hybridization (ISH). **(A)** Schematic of 5 dpf inner ear cellular organization. The position of transversal views **(C--E)** is outlined by dotted lines. **(B)** Lateral view of the anterior region of the 5 dpf larvae. Yellow circle delimits the inner ear location. **(C,C',C")** Transversal view of the anterior region of the inner ear. The inner ear is outlined by a white dotted line. Arrows point to hair cell patches location. **(C)** *slc9a3r1* ISH. **(C')** brain3c:GFP. **(C")** *slc9a3r1* and brain3c:GFP merged image. **(D,D',D")** Transversal view of the medial region of the inner ear. The inner ear is outlined by a white dotted line. Arrows point to hair cell patches location. **(D)** *slc9a3r1* ISH. **(D')** brain3c:GFP. **(D")** *slc9a3r1* and brain3c:GFP merged image. **(E,E',E")** Transversal view of the posterior region of the inner ear. The inner ear is outlined by a white dotted line. Arrows point to hair cell patches location. **(E)** *slc9a3r1* ISH. **(E')** brain3c:GFP. **(E")** *slc9a3r1* and brain3c:GFP merged image.](fgene-10-00142-g003){#F3}

### 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 ([Figure 1A--D](#F1){ref-type="fig"}). Overall 135 six days post 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 *slc9a3r*1^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](#F4){ref-type="fig"}) while no variation in the ability to respond to mechanical stimuli was observed in comparison to the wild type siblings.

![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).](fgene-10-00142-g004){#F4}

### 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 ([@B14]). No differences have been observed for the *slc9a3r1*^+/R180Q^ and *slc9a3r1*^R180Q/R180Q^ models (see [Figure 4B](#F4){ref-type="fig"}). 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 ([@B43]). The alteration of the size, shape or the absence of this structure could have striking effects on the hearing ability ([@B25]). 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) ([Figure 5A,B](#F5){ref-type="fig"}). These results thus explain the hearing impairment observed in the *slc9a3r1*^R180Q/R180Q^ KI larvae.

![Measurement of the otolith size. **(A)** Measurement of the posterior otolith in *slc9a3r1*^+/+^ and *slc9a3r1*^R180Q/R180Q^ larvae. **(B)** Measurement of posterior otolith area of *slc9ar1* larvae (number of tested animals: *slc9a3r1*^+/+^ *N* = 57, *slc9a3r1*^+/-^ *N* = 70, *slc9a3r1*^R180Q/R180Q^ *N* = 32). ^∗∗^*p*-value \< 0.005 and ^∗∗∗^*p*-value \< 0.0005.](fgene-10-00142-g005){#F5}

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 ([Figure 6A,B](#F6){ref-type="fig"}). 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](#F6){ref-type="fig"}) 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](#F6){ref-type="fig"}), 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.

![Hair cells number, morphology and orientation in Zebrafish larvae. **(A)** 3D reconstruction of the saccular macula. Hair cells are marked with GFP. FM1-43FX dye labels only mature and fully functional hair cells. **(B)** Saccular hair cells count in *slc9a3r1* larvae (number of tested animals: *slc9a3r1*^+/+^ *N* = 9, *slc9a3r1*^+/-^ *N* = 10, *slc9a3r1*^R180Q/R180Q^ *N* = 6). **(C)** 3D reconstruction of the saccular hair cell layer. brn3c:mGFP signal labels the membrane of the hair cells while red signal labels the kinocilia. **(D)** Hair cell polarization within the saccular macula of *slc9a3r1* KI and their wt siblings.](fgene-10-00142-g006){#F6}

### 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.

In 6 dpf larvae (3 *slc9a3r1^+/+^*, 6 *slc9a3r1*^R180Q/R180Q^, 6 *slc9a3r1*^+/R180Q^ animals used) few differences (although not statistically significant) in the volume of the total SAG (42.439 μm^3^ for *slc9a3r1^+/+^*, 35.454 μm^3^ for *slc9a3r1*^R180Q/R180Q^, 40.237 μm^3^ for *slc9a3r1*^+/R180Q^), and no alterations in the volume of the posterior SAG (13.676 μm^3^ for *slc9a3r1^+/+^*, 11.404 μm^3^ for *slc9a3r1*^R180Q/R180Q^, 12.709 μm^3^ for *slc9a3r1*^+/R180Q^) were detected ([Figure 7A,B](#F7){ref-type="fig"}), thus suggesting that the hearing impairment is not due to a reduction of neurons in this structure.

![Stato-acoustic ganglion (SAG) characterization. **(A)** 3D reconstruction of the SAG in *slc9a3r1* KI and their WT siblings (number of tested animals: *slc9a3r1*^+/+^ *N* = 3, *slc9a3r1*^+/-^ *N* = 6, *slc9a3r1*^R180Q/R180Q^ *N* = 6). The anterior part of the SAG is depicted in pink, while in green is indicated the posterior part which innervates the saccular macula and it is related to the auditory function. **(B)** The graph shows the total SAG volume (volume of the pink and green area) and the volume of the posterior SAG of *slc9a3r1* animals. Although small differences (not statistically significant) no alterations were observed in the volume of the posterior stato-acoustic ganglion.](fgene-10-00142-g007){#F7}

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 ([@B16]) 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 function ([@B3]; [@B30]). 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 ([@B12]).

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 ([@B26]). 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 [@B34]. 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 ([@B47]). 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 ([@B44]) 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 ([@B32]).

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 ([@B8]). 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 ([@B33]). 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 ([@B53]; [@B8]) 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 ([@B26]). 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.

Conflict of Interest Statement
==============================

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**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.

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

<https://www.ncbi.nlm.nih.gov/projects/SNP/>

<http://www.1000genomes.org/>

<http://evs.gs.washington.edu>

<http://exac.broadinstitute.org/about>

<http://gnomad.broadinstitute.org/>

[www.pymol.org](http://www.pymol.org)

<http://bioinfo.ut.ee/primer3/>

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>

###### 

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).

###### 

Click here for additional data file.

###### 

*In vitro* experiments: mRNAs quantification and proteins analyses. After cells transfection with expression vectors containing either the WT or the mutant *SLC9A3R1* cDNA, mRNA and protein analyses were performed. **(A)** qRT-PCR on mRNA showed no difference in the expression levels of the mutant (^∗^) compared to the WT. **(B)** Western blot analysis demonstrated that WT and mutated protein are equally expressed. **(C)** Immunolocaliation showed that both the WT and the mutated proteins have a cytoplasmic expression.

###### 

Click here for additional data file.

[^1]: Edited by: Mike Mikailov, U.S. Food and Drug Administration, United States

[^2]: Reviewed by: Theodora Katsila, University of Patras, Greece; GuangJun Zhang, Purdue University, United States

[^3]: This article was submitted to Genetic Disorders, a section of the journal Frontiers in Genetics
