Pendred Syndrome, or Not Pendred Syndrome? That Is the Question

Pendred syndrome (PDS) is the most common form of syndromic Hearing Loss (HL), characterized by sensorineural HL, inner ear malformations, and goiter, with or without hypothyroidism. SLC26A4 is the major gene involved, even though ~50% of the patients carry only one pathogenic mutation. This study aims to define the molecular diagnosis for a cohort of 24 suspected-PDS patients characterized by a deep radiological and audiological evaluation. Whole-Exome Sequencing (WES), the analysis of twelve variants upstream of SLC26A4, constituting the “CEVA haplotype” and Multiplex Ligation Probe Amplification (MLPA) searching for deletions/duplications in SLC26A4 gene have been carried out. In five patients (20.8%) homozygous/compound heterozygous SLC26A4 mutations, or pathogenic mutation in trans with the CEVA haplotype have been identified, while five subjects (20.8%) resulted heterozygous for a single variant. In silico protein modeling supported the pathogenicity of the detected variants, suggesting an effect on the protein stabilization/function. Interestingly, we identified a genotype-phenotype correlation among those patients carrying SLC26A4 mutations, whose audiograms presented a characteristic slope at the medium and high frequencies, providing new insights into PDS. Finally, an interesting homozygous variant in MYO5C has been identified in one patient negative to SLC26A4 gene, suggesting the identification of a new HL candidate gene.


Introduction
Hearing loss (HL) is the most common sensory disorder affecting about 1-3 children in over 1000 [1][2][3]. In the developed countries, the most prevalent causes of HL are genetic defects that impair the production of protein with essential roles in the hearing pathway. HL is generally classified in non-syndromic HL (NSHL), affecting at least 70% of the patients, and syndromic HL (SHL), when other organs abnormalities are present alongside the HL [2]. To date, over 400 syndromes associated with HL have been described [4,5], suggesting the high degree of phenotypic and genetic heterogeneity that characterized the affected patients.
Pendred syndrome (PDS) is considered the most common form of SHL accounting for~4-10% of all cases of hereditary deafness [6,7]. PDS is inherited with an autosomal recessive pattern and it is associated with sensorineural HL, inner ear malformations, consent and their next of kin provided it in case of minors. All research was conducted according to the ethical standard defined by the Helsinki Declaration.
The individuals underwent a careful clinical evaluation. In detail, patients' past medical records were investigated by clinical geneticists to exclude HL cases that could be ascribed to non-genetic causes (e.g., infections or trauma). The patients underwent a careful clinical evaluation specifically focused on dysmorphic features and other congenital abnormalities not specifically related with the syndrome itself. Moreover, for all the patients, pure tone audiometric testing (PTA) or auditory brainstem response (ABR) was performed to define the degree of HL. To define the presence of EVA the Cincinnati criteria were followed (e.g., the vestibular aqueduct is defined enlarged if it has a midpoint width (1.0 mm) or opercular width (2.0 mm) greater than the 95th percentile) [16,17]. Thus, to perform the statistical analysis Magnetic Resonance Imaging (MRI) and Computerized Tomography (CT) scan were performed as well. Finally, thyroid function assessments were carried out. The levels of serum thyroid stimulating hormone (TSH), serum free thyroxine (T4), and serum free triiodothyronine (T3) were evaluated. Ultrasound was also performed whenever altered levels were found on the blood test or if thyroid clinical examination suggested the presence of goiter or nodules. The thyroid was considered enlarged if the measurements reveled a mean volume greater than 8 to 10 mL (range 3 to 20 mL) [18,19].

DNA Extraction and Quantification
Genomic DNA was extracted from peripheral whole blood samples using the QIAsym-phony®SP instrument with QIAsymphony®DNA Midi kit (Qiagen, Venlo, The Netherlands) and DNA concentration measured using Nanodrop ND 1000 spectrophotometer (Nan-oDrop Technologies Inc., Wilmington, DE, USA).

GJB2 and GJB6 Analyses
Sanger method was employed to sequence the entire coding region of GJB2 gene (primers available upon request). In particular, DNA was analyzed on a 3500 Dx Genetic Analyzer (ThermoFisher, Waltham, MA, USA), using ABI PRISM 3.1 Big Dye terminator chemistry (ThermoFisher, Waltham, MA, USA) according to the manufacturer's instructions. Subsequently, GJB6 deletions (i.e., D13S1830-D13S1854) were screened by multiplex PCR, as previously described [20].

Whole Exome Sequencing (WES)
WES was carried out on an Illumina NextSeq 550 instrument (Illumina Inc., San Diego, CA, USA). According to the manufacturer's instructions. Briefly, 50 ng of genomic DNA were enzymatically fragmented, and, after an end repair and dA-tailing reaction, each fragment was ligated to a universal adapter and then amplified using the Unique Dual Index primers (Twist Bioscience). Afterwards, genomic libraries were prepared using the Twist Human Core Exome + Human RefSeq Panel kit (Twist Bioscience, South San Francisco, CA, USA) that allows to cover 99% of the protein-coding genes. In conclusion, the hybridized fragments have been captures, amplified and sequenced.
This workflow, designed for Illumina paired-end sequencing data, allows the generation of a final VCF file containing information regarding germline variants, such as Single Nucleotide Variants (SNVs), short insertion/deletions (INDELs) and exon-level copy number variations (CNVs). Finally, VCF files were analyzed on EnGenome Expert Variant Interpreter (eVai) software (evai.engenome.com) that allows variant annotation and interpretation. In particular, eVai combines Artificial Intelligence with the American College of Medical Genetics (ACMG) guidelines [21] to classify and prioritize every genomic variant, suggesting a list of possible related genetic diagnoses.
SNVs and INDELs were excluded if they lead to synonymous amino acid substitutions that were not predicted as damaging or did not affect splicing or highly conserved residues. Furthermore, variants with a quality score (QUAL) < 20 or called in off-target regions were excluded as well.
Variants previously reported as polymorphism were removed thanks to a comparison between the identified genetic variants and data reported in NCBI dbSNP build153.
(  [25] and dbscSNV score [26] were used to evaluate the pathogenicity of novel variants.
Finally, on a patient by patient basis, variants were discussed in the context of phenotypic data at interdisciplinary meetings. As the last step, and the most likely disease-causing variants were confirmed by direct Sanger sequencing.

CEVA Haplotype Analysis
Considering that the subjects involved in the study are Caucasian, we employed Sanger sequencing to test the presence of the CEVA haplotype in M1 patients. We followed the procedure as already described in 2.3. (primers available upon request). In particular, we sequenced fragment of 200-400 bp containing the twelve SNP of the haplotype (i.e., rs17424561, rs79579403, rs17425867, rs117113959, rs17349280, rs117386523, rs80149210, rs199667576, rs9649298, rs117714350, rs199915614, rs150942317).

Multiplex Ligation Probe Amplification (MLPA)
MLPA analysis was performed searching for deletion/duplication in the SLC26A4 gene. The SALSA ® MLPA ® probe mixes P280-100R SLC26A4 (MRC-Holland, Amsterdam, the Netherlands) was employed, according to the manufacturer's protocol. For the data analysis, the software Coffalyser.Net was used in combination with the lot-specific MLPA Coffalyser sheet. The following cutoff values for the dosage quotient (DQ) of the probes were used to interpret MLPA results: 0.80 < DQ < 1.20 (no deletion/duplication), DQ = 0 (deletion), and 1.75 < DQ < 2.15 (duplication).

Prediction of Membrane Topology and 3D Molecular Model of Pendrin
TOPCONS [27] was used to predict the consensus prediction of membrane protein topology of pendrin (encoded by SLC26A4). No homologous transmembrane protein was detected in the Protein Data Bank. Overall, the software predicted 13 transmembrane (TM) regions, 7 intracellular regions and 6 extracellular regions ( Figure 1) with different confidence.
A BLASTP run of the sequence of human pendrin using the Protein Data Bank as a reference database highlighted some hits, the best one covering 91% of the protein sequence (human SLC26A9 chain A; PDB: 7CH1; s.i. 35.13%). We also used the neural network based Alphafold algorithm [28] for the prediction of human Pendrin three-dimensional structure (https://alphafold.ebi.ac.uk/entry/O43511 (accessed on 7 September 2021). The performance was satisfactory, especially for the TM regions, for which very high (predicted Local Distance Difference Test, pLDDT < 90) to confident (70 < pLDDT < 90) levels characterized the prediction of the TM regions. Superposition of the experimental (CryoEM) structure of chain A of the dimeric SLC26A9 assembly (PDB: 7CH1) [29] to the structure Genes 2021, 12, 1569 5 of 14 originated by Alphafold led to relatively low (2.73 Å) root-mean squared deviation (RMSD); therefore, we used the Alphafold model for the analysis of the putative effects of amino acid substitutions in pendrin. TOPCONS [27] was used to predict the consensus prediction of membrane protein topology of pendrin (encoded by SLC26A4). No homologous transmembrane protein was detected in the Protein Data Bank. Overall, the software predicted 13 transmembrane (TM) regions, 7 intracellular regions and 6 extracellular regions ( Figure 1) with different confidence. Figure 1. Membrane topology prediction for the protein expressed by SLC26A4 (pendrin) according to TOPCONS [27]. Boxes refer to the transmembrane (TM) region (in vs. out, grey; out vs. in, white); red lines refer to intracellular regions; blue lines refer to extracellular regions. X-axis refers to the position within protein sequence. Reliability estimate of the predictions is reported in the Yaxis.
A BLASTP run of the sequence of human pendrin using the Protein Data Bank as a reference database highlighted some hits, the best one covering 91% of the protein sequence (human SLC26A9 chain A; PDB: 7CH1; s.i. 35.13%). We also used the neural network based Alphafold algorithm [28] for the prediction of human Pendrin three-dimensional structure (https://alphafold.ebi.ac.uk/entry/O43511 (accessed on 7th of September 2021). The performance was satisfactory, especially for the TM regions, for which very high (predicted Local Distance Difference Test, pLDDT < 90) to confident (70 < pLDDT < 90) levels characterized the prediction of the TM regions. Superposition of the experimental (CryoEM) structure of chain A of the dimeric SLC26A9 assembly (PDB: 7CH1) [29] to the structure originated by Alphafold led to relatively low (2.73 Å ) root-mean squared deviation (RMSD); therefore, we used the Alphafold model for the analysis of the putative effects of amino acid substitutions in pendrin.

Results
Twenty-four patients presenting with HL and temporal bone abnormalities suggesting PDS have been considered in the study.
All the clinical features of the enrolled patients are reported in Table 1. CT and MRI were performed in 17 out of 24 patients, while CT alone was carried out for 23 patients out of 24. Briefly, the in-depth morphological evaluation highlighted EVA, ESE, Mondini malformation, or IP-II, and semicircular canal and vestibular abnormalities in the vast majority of them. In particular, EVA, and eventually ESE, alone were present in 21% of our patients (ID6, ID8, ID12, ID18, ID23). On the other hand, both EVA and Mondini malformation were detected in 79% of the patients (ID1, ID2, ID3, ID4, ID5, ID7, ID9, ID10,  ID11, ID13, ID14, ID15, ID16, ID17, ID19, ID20, ID21, ID22, ID24). As regards thyroid function tests (serum TSH, serum-free T4, and serum-free T3 levels), data were available for 23 out of 24 patients and one (ID6) showed subclinical hypothyroidism (TSH 5.07 mU/L with a normal value ranging from 0.27 and 4.20) at the age of 5 years. Ultrasound examination of the thyroid was performed in two out of 24 patients.  [27]. Boxes refer to the transmembrane (TM) region (in vs. out, grey; out vs. in, white); red lines refer to intracellular regions; blue lines refer to extracellular regions. X-axis refers to the position within protein sequence. Reliability estimate of the predictions is reported in the Y-axis.

Results
Twenty-four patients presenting with HL and temporal bone abnormalities suggesting PDS have been considered in the study.
All the clinical features of the enrolled patients are reported in Table 1. CT and MRI were performed in 17 out of 24 patients, while CT alone was carried out for 23 patients out of 24. Briefly, the in-depth morphological evaluation highlighted EVA, ESE, Mondini malformation, or IP-II, and semicircular canal and vestibular abnormalities in the vast majority of them. In particular, EVA, and eventually ESE, alone were present in 21% of our patients (ID6, ID8, ID12, ID18, ID23). On the other hand, both EVA and Mondini malformation were detected in 79% of the patients (ID1, ID2, ID3, ID4, ID5, ID7, ID9, ID10,  ID11, ID13, ID14, ID15, ID16, ID17, ID19, ID20, ID21, ID22, ID24). As regards thyroid function tests (serum TSH, serum-free T4, and serum-free T3 levels), data were available for 23 out of 24 patients and one (ID6) showed subclinical hypothyroidism (TSH 5.07 mU/L with a normal value ranging from 0.27 and 4.20) at the age of 5 years. Ultrasound examination of the thyroid was performed in two out of 24 patients. All the patients underwent a first round of genetic testing, including the analysis of GJB2 and GJB6 genes, which resulted negative.
WES allowed to identify the molecular cause of PDS in four patients.
In particular, WES highlighted the presence of homozygous and compound heterozygosis mutations in the SLC26A4 (NM_000441.1) gene in four patients: In detail, ID1, a 29-year-old female, is affected by profound HL, with a notable worst performance at medium and high frequencies (Figure 2), together with EVA, Mondini, ESE, semicircular canal/vestibular abnormalities, and no thyroid function dysfunction.
WES revealed the presence of a homozygous mutation (c.1225C>T, p.(R409C)) predicted as pathogenic by several in silico tools and already described as pathogenic in association with EVA [30] (Table 2).   WES revealed the presence of a homozygous mutation (c.1225C>T, p.(R409C)) predicted as pathogenic by several in silico tools and already described as pathogenic in association with EVA [30] (Table 2).   Table 2) that has already been classified as pathogenic and associated with bilateral EVA [32].
Patient ID22, a 29-year-old female, shows an audiometric pattern characterized by a profound slope at the medium and high frequencies, and mild/moderate HL at the low ones ( Figure 2) in addition to bilateral EVA, ESE, IP-II and semicircular canal/vestibular abnormalities (Table 1). At the age of 28, her last thyroid ultrasound showed mild enlargement of the thyroid gland with a heterogeneous echogenicity without nodules. TSH, T3 and T4 levels have been normal throughout her follow-up. ID22 carries a homozygous mutation (c.1489G>A, NM_000441.1, p.(G497S)) classified as pathogenic and detected in patients displaying HL and EVA [34] (Table 2).
A molecular diagnosis was also provided for patient ID12, a 12-year-old female who carries two SLC26A4 mutations at the compound heterozygous state. The patient displays moderate and asymmetric HL at the low frequencies and severe/profound HL at the medium and high frequencies ( Figure 2) together with bilateral EVA (Table 1) without thyroid dysfunction. The two known mutations (c.1001+1G>A and c.1149+3A>G), already classified as pathogenic, alter the splicing process [35,36] (Table 2).
Moreover, six additional patients resulted carriers of only one heterozygous mutation and were initially classified as M1: All the variants were classified as pathogenic by different in silico prediction tools and have been previously associated with PDS or EVA [30,[33][34][35][36][37][38], with the only exception of the one carried by ID19 ( Table 2). All the patients, aged between six and 16, are affected by moderate to profound HL and displayed both bilateral EVA and IP-II. The only exception is ID10, who displays asymmetrical HL and monolateral EVA and IP-II, all affecting mainly the left ear (Table 1, Figure 2, Table S1). For all six patients, we were able to measure TSH, T3, and T4 levels and so far, no one showed hypothyroidism or subclinical hypothyroidism. However, one patient, ID17, had a thyroid ultrasound showing a normal thyroid volume associated with a heterogeneous echogenicity of the thyroid gland without focal lesions.
Interestingly, for patient ID19, a 16-year-old male, who displays profound bilateral HL, EVA, IP-II and semicircular canal/vestibular abnormalities, the analysis of the CEVA haplotype revealed its presence in trans with the pathogenic variant ((c.600G>A, NM_000441.1, p.(Q200Q))) ( Table 2). The predicted effects of this allele on the encoded protein are both synonymous aminoacidic substitution and splicing alteration; in fact, the nucleotide change involved the last nucleotide of exon 5, likely altering the normal splicing site and thus impacting on the final structure of the pendrin protein.
Further, we performed MLPA analysis targeting the SLC26A4 gene in the M1 and M0 patients, and all the individuals resulted negative for deletion or duplication in the SLC26A4 gene.
In order to strengthen the pathogenic effect of the identified variants, we performed an in silico protein modeling. To date, no experimental structure has been solved for pendrin. The performance of the novel neuronal network-based algorithm Alphafold [28] is extremely high in predicting the three-dimensional fold of proteins even when the sequence similarity is too low to allow for the building of reliable homology models. The Alphafold-predicted three-dimensional structure of pendrin is reported in Figure 3.
The overall topology is in line with predictions based on the widely used TOPCONS methodology (see methods) and shows the typical fold of a transmembrane transporter. It is worth noting that the protein sequence of pendrin shares 32% sequence identity with the human SLC26A9 Cl − and Na + transporter, for which the dimeric cryo-EM structure has been recently solved [29]. The structure represented in Figure 3 is therefore a prediction limited to the monomeric subunit of a putatively dimeric assembly.
All mutated sidechains in the missense variants found in this study except V577A, located in the cytoplasmic domain, are in the transmembrane region of pendrin and their predicted orientation is towards the interior of the transporter rather than towards the lipid milieu ( Figure 3A, B, viewed from the top). It is therefore plausible that the mutations do not affect the supramolecular assembly of pendrin into dimers, but rather may lead to destabilization of the protein or prevent its function of sodium-independent transporter of chloride and iodide. SLC26A4 gene.
In order to strengthen the pathogenic effect of the identified variants, we performed an in silico protein modeling. To date, no experimental structure has been solved for pendrin. The performance of the novel neuronal network-based algorithm Alphafold [28] is extremely high in predicting the three-dimensional fold of proteins even when the sequence similarity is too low to allow for the building of reliable homology models. The Alphafold-predicted three-dimensional structure of pendrin is reported in Figure 3. The overall topology is in line with predictions based on the widely used TOPCONS methodology (see methods) and shows the typical fold of a transmembrane transporter. It is worth noting that the protein sequence of pendrin shares 32% sequence identity with the human SLC26A9 Cland Na + transporter, for which the dimeric cryo-EM structure has been recently solved [29]. The structure represented in Figure 3 is therefore a prediction limited to the monomeric subunit of a putatively dimeric assembly.
All mutated sidechains in the missense variants found in this study except V577A, located in the cytoplasmic domain, are in the transmembrane region of pendrin and their predicted orientation is towards the interior of the transporter rather than towards the lipid milieu ( Figure 3A, B, viewed from the top). It is therefore plausible that the mutations do not affect the supramolecular assembly of pendrin into dimers, but rather may lead to destabilization of the protein or prevent its function of sodium-independent transporter of chloride and iodide.
Interestingly, some pendrin mutations identified here correspond to residues playing important functional roles in SLC26A9. In particular, Arginine in position 185, substituted Interestingly, some pendrin mutations identified here correspond to residues playing important functional roles in SLC26A9. In particular, Arginine in position 185, substituted by a Threonine in p.(R185T) is the analog of R169 in SLC26A9, a basic residue involved in an electrostatic interaction with a Cl − ion. Moreover, the side chain of R409 protrudes to a region that, in SLC26A9 is occupied by a Na + ion and the R409C substitution would lose any potentially conserved electrostatic intersection in pendrin. It is not clear what could be the effect of the V577A substitution in the cytoplasmic milieu of pendrin, but the fold of that region is predicted with less confidence by Alphafold and structural differences with SLC26A9 are more prominent that in other regions. As for the other variants, the G497S substitution might induce some steric clash between adjacent helices, while the Q235R substitution might induce electrostatic repulsion with K447, thus destabilizing the protein core. Finally, replacing the aliphatic L445 side chain with a bulky, aromatic Trp could also lead to significant structural rearrangement in a packed hydrophobic region. As for the p.(Q200Q) variant, the synonymous mutation is obviously not predicted to affect either protein structure or function.
Finally, 14 individuals did not display any causative mutation in SLC26A4 (M0), nor in other PDS related genes. The comparison between their audiometric profiles with those of M2/M1 individuals revealed a milder phenotype in M0 patients, suggesting the possible involvement of other causative genes. (Table S1; Figure S1). In this light, we searched for variants in genes already known for being causative of SHL/NSHL (Hereditary Hearing Loss Homepage; http://hereditaryhearingloss.org/, date of access: 7 September 2021), or in genes likely involved in the auditory function, without finding any interesting result, with the only exception of patient ID24. The subject carries an interesting homozygous variant in MYO5C (NM_018728.4, c.3592C>T p.(R1198C)) ( Table 2), a novel gene that has not been associated with HL yet. The variant is predicted as pathogenic by the in silico tools used during data analysis and is described with an ultra-rare minor allele frequency in the gnomAD database (rs200717531), with no homozygous individuals reported. ID24 displays bilateral EVA, without ESE, bilateral IP-II, semicircular canal and vestibular abnormalities ( Table 1). The newborn hearing screening revealed a monolateral hearing impairment (right ear), that progressed rapidly to bilateral and profound HL (Figure 4) by the age of three. The patient did not manifest hypothyroidism or subclinical hypothyroidism. 2021), or in genes likely involved in the auditory function, without finding any interesting result, with the only exception of patient ID24. The subject carries an interesting homozygous variant in MYO5C (NM_018728.4, c.3592C>T p.(R1198C)) ( Table 2), a novel gene that has not been associated with HL yet. The variant is predicted as pathogenic by the in silico tools used during data analysis and is described with an ultra-rare minor allele frequency in the gnomAD database (rs200717531), with no homozygous individuals reported. ID24 displays bilateral EVA, without ESE, bilateral IP-II, semicircular canal and vestibular abnormalities (Table 1). The newborn hearing screening revealed a monolateral hearing impairment (right ear), that progressed rapidly to bilateral and profound HL (Figure 4) by the age of three. The patient did not manifest hypothyroidism or subclinical hypothyroidism.

Discussion
When PDS was first described, the diagnosis was made only based on clinical signs (e.g. deafness, goitre, thyroid dysfunctions, temporal bone malformations). However, since the discovery of the main disease-causing gene, SLC26A4, the diagnostic procedures have improved, and today it is possible to define the molecular basis of the patients' disease [9]. Nevertheless, PDS patients are characterized by high clinical variability and the detection rates of PDS causing mutations remain limited, highlighting the importance of a careful clinical examination.
According to literature data, the detection rate of at least one mutation ranges between 25-50% [12,13,15]. Among our cohort, we identified five M2 patients (20.8%) (including one patient carrier of a pathogenic mutation in trans with the CEVA haplotype) and five M1 subjects (20.8%), for an overall SLC26A4-mutations carriers rate of 41.7%. All the mutations here identified are predicted as damaging and our in silico modeling suggested their possible effect on protein stabilization and function.
The lack of a second pathogenic allele explaining the clinical phenotype of M1 patients might be due to several reasons, some of them tested in our work.
In a limited series of patients, the hypothesis of a digenic inheritance involving other genes, such as FOXI1, KCNJ10, and EPHA2, has been shown [13,14]; however, none of the subjects of our cohort carried a likely pathogenic variant in any of the above-mentioned genes.
We also checked for the presence of any CNV in the SLC26A4 gene using MLPA, with negative results. In this light, after excluding the presence of the CEVA haplotype, the most reasonable explanation implies the presence of variants in regions not explored with WES, or complex structural variants not detectable with the methodologies used in this work. In addition, a new pathophysiological model for PDS has also been recently hypothesized [39]. Hosoya et al. proposed that a single mutation on the SLC26A4 gene might cause the formation of misfolded pendrin aggregates and, if defects of the ubiquitinproteasome system and/or autophagy pathways are also present, this could ultimately result in inner ear cells death. Briefly, they suggest that some M1 patients' phenotype might be due to a decreased clearance of the mutant protein causing degeneration, rather than a channelopathy. Thus, specific genetic backgrounds, but also environmental factors, might exacerbate the stress induced by monoallelic SLC26A4 mutation leading to inner ear damage, which could ultimately cause the PDS phenotype in a specific subset of patients.
Based on these data, additional efforts for understanding the pathological mechanism responsible for the clinical phenotype of these patients are needed, including, for instance, the application of whole genome sequencing, RNA studies to detect quantitative or qualitative defects in SLC26A4 transcripts or further studies involving human-induced pluripotent stems cells from PDS patients.
Given the availability of an extreme accurate clinical evaluation, we tried to outline the main features of our PDS patients.
As regards the radiological observations, a clear distinction between M2/M1 and M0 patients could not be made. Thus, the identification of bilateral EVA or Mondini should not be considered as a proof of the presence of SLC26A4 gene variants.
Data regarding the thyroid hormone levels were available for almost all the patients (23 out of 24) and except for one patient (ID6, 13-year-old) who shows subclinical hypothyroidism, all the other subjects seem to be euthyroid. As mentioned before, goitre generally becomes apparent in adolescence (more or less after the age of ten) and in our cohort ten out of 23 patients are aged <10 years. Moreover, environmental factors, such as iodine intake can play an important role since literature data report that the prevalence of goitre is inversely proportional to daily iodine intake [40]. Considering that all the recruited patients live in Trieste, a seaside city, it is reasonable to hypothesize that they are exposed daily to high iodine concentrations.
Concerning the audiological data, interestingly, we highlighted a new genotypephenotype correlation that, to our knowledge, has been previously described but never deeply explored due to a lack of accurate audiological data of PDS patients [3,[40][41][42][43][44]. In particular, M1 and M2 subjects display highly similar audiometric patterns. The majority of these patients are characterized by mild-moderate HL at the low frequency, which gradually worsens, becoming moderate-severe at the medium frequencies and profound at the high frequencies. Only few M1/M2 patients display profound HL with flat audiograms. Furthermore, audiometric data of M1/M2 patients seem to differ from the M0 cases who show a milder phenotype, confirming the hypothesis that different genetic causes should be searched in these patients. Furthermore, we identified an unexpected phenotype in a M2 patient carrier of the CEVA haplotype. In fact, according to literature data the concomitant presence of a heterozygous mutation in trans with CEVA haplotype is associated with milder forms of PDS [15]. However, patient ID19 displays a profound hearing impairment at all frequencies, suggesting a higher phenotypic variability associated with the CEVA haplotype.
The remaining 14 patients of our cohort do not carry any mutation in PDS genes, nor in other known HL genes. Moreover, eight out of 14 M0 patients displayed fluctuation and/or progression of the HL, which are signs typically associated with PDS. This consideration might suggest that the HL of these eight patients could be ascribed to mutations located in genomic regions not explored with WES, such as SLC26A4 introns. Alternatively, the causative mutations might be found in novel disease genes. Interestingly, among the M0 individuals, we identified a patient carrier of a homozygous missense variant in MYO5C gene. MYO5C encodes the myosin-Vc protein that belongs to the class V myosins, one of the most ancient and distributed group of myosins, hypothesized to act as motors for actin-dependent organelle transport [45]. Myosin-Vc seems to be mostly expressed in epithelial and glandular tissues and shares~50% identity with the two other classes of myosins V in vertebrates, myosin-Va (Myo5a) and myosin-Vb (Myo5b) [46]. Even though no mutations in MYO5C have been linked to any known disease yet, other myosins have been associated with HL before, such as MYO15A, MYO3A, MYO6, MYO7A, etc [47]. In this light, the identification of an ultra-rare homozygous missense variant, predicted as pathogenic by several in silico tools, represents an interesting finding worth mentioning. Additional studies clarifying the role of this gene will be crucial to define if it indeed could be involved the etiopathogenesis of HL. Nevertheless, the evaluation of MYO5C in other PDS-like subjects may lead to the identification of other patients carriers of pathogenic variant, reinforcing the hypothesis of a novel candidate gene for PDS.
Overall, these findings highlight the efficacy of WES and the analysis of the CEVA haplotype together with an accurate clinical evaluation, both from the audiological, and radiological point of view. In fact, a deep phenotypical evaluation might allow a precise definition of the essential symptoms common to all PDS patients leading to the definition of consistent phenotype-genotype correlations that will ultimately increase the detection rates.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes12101569/s1, Table S1: Patients' hearing thresholds. Figure S1: Audiometric pattern of carriers and non-carriers of pathogenic alleles.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study and written informed consent has been obtained from the patient to publish this paper.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.