The Enigmatic Etiology of Oculo-Auriculo-Vertebral Spectrum (OAVS): An Exploratory Gene Variant Interaction Approach in Candidate Genes

The clinical diagnosis of oculo-auriculo-vertebral spectrum (OAVS) is established when microtia is present in association with hemifacial hypoplasia (HH) and/or ocular, vertebral, and/or renal malformations. Genetic and non-genetic factors have been associated with microtia/OAVS. Although the etiology remains unknown in most patients, some cases may have an autosomal dominant, autosomal recessive, or multifactorial inheritance. Among the possible genetic factors, gene–gene interactions may play important roles in the etiology of complex diseases, but the literature lacks related reports in OAVS patients. Therefore, we performed a gene–variant interaction analysis within five microtia/OAVS candidate genes (HOXA2, TCOF1, SALL1, EYA1 and TBX1) in 49 unrelated OAVS Mexican patients (25 familial and 24 sporadic cases). A statistically significant intergenic interaction (p-value < 0.001) was identified between variants p.(Pro1099Arg) TCOF1 (rs1136103) and p.(Leu858=) SALL1 (rs1965024). This intergenic interaction may suggest that the products of these genes could participate in pathways related to craniofacial alterations, such as the retinoic acid (RA) pathway. The absence of clearly pathogenic variants in any of the analyzed genes does not support a monogenic etiology for microtia/OAVS involving these genes in our patients. Our findings could suggest that in addition to high-throughput genomic approaches, future gene–gene interaction analyses could contribute to improving our understanding of the etiology of microtia/OAVS.

Microtia is a minimal expression of the OAVS clinical spectrum, since both entities have variable phenotypic expression, asymmetric involvement of facial structures, rightside predominance, and male predilection. Moreover, familial occurrence with incomplete penetrance may be seen for microtia and/or related anomalies, such as preauricular tags and pits [9][10][11].

Study Population
Forty-nine unrelated Mexican patients (31 males, 18 females; ages 0 to 18 years old) who were evaluated between 2015 and 2019 at the National Institute of Pediatrics (Mexico) and given a clinical diagnosis of microtia/OAVS (microtia or anotia with or without HH, structural alterations in the spine and/or kidneys) were included. Patients who met the clinical criteria for microtia/OAVS but had congenital malformations distinct from those involving the spine and/or kidneys and/or an alteration in somatic growth (low or high height) or intellectual disability were excluded. Those with reported prenatal exposure to specific teratogens associated with microtia/OAVS (e.g., alcohol, retinoids, or maternal diabetes) were also excluded. Each patient underwent a systematic physical examination by clinical geneticists and a detailed family history was obtained (25 familial cases and 24 sporadic cases). Imaging studies (computed tomography of the inner ear, orthopantomography, complete spinal radiography, and renal ultrasound) were performed on the patients and their parents. Detailed clinical information related to all included patients was previously published [41]. This research protocol was approved and registered by the Ethics, Research and Biosafety Committees of the National Institute of Pediatrics (Mexico City, Mexico, registry number 004/2017).

NGS of the Five Microtia/OAVS Candidate Genes
Genomic DNA isolated from peripheral blood or buccal cell samples of all patients was analyzed with a targeted five-gene NGS panel consisting of HOXA2 (NM_006735.3), TCOF1 (NM_001135243.1), SALL1 (NM_002968.2), EYA1 (NM_000503.5), and TBX1 (NM_080647.1). The studied sequences covered the promoter region, all coding exons, and the intron-exon boundaries (200 bp). NGS libraries were prepared with an IDT xGen lockdown probe customized panel kit according to the manufacturer's protocol. Libraries were sequenced on an Illumina MiSeq 2 × 150 platform (San Diego, CA, USA) through the Admera Health Company (South Plainfield, NJ, USA). Our in-house bioinformatics analysis pipeline included a quality evaluation and trimming of low-quality raw reads, alignment against the GRCh38 human reference sequence, and calling of single nucleotide variants (SNVs) and small insertion-deletion variants. The GATK Toolkit (version 4.2.6.1) [42] and the Alamut Batch (version 1.4), Focus (version 1.0), and Visual (version 2.7.2, Interactive Biosoftware, Rouen, France) software packages were used for variant annotation and filtering. All clinically relevant variants were confirmed by Sanger automated sequencing in the index cases and their available relatives. The identified gene variants were classified according to the criteria of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (ACMG/AMP) [43].

Statistical and Gene Interaction Analyses
The allele frequencies of the variants identified in the five studied genes were obtained using the allelic counting method. We first determined if the variant was in Hardy-Weinberg equilibrium using the two-tailed Fisher's exact test, which assessed whether the distribution of genotypes was as expected in our population. When information was available, an association analysis was performed to compare the genotype frequencies in our study population with those reported in Mexican individuals from Los Angeles (1000 Genomes Project, phase 3; https://www.ncbi.nlm.nih.gov/genome/gdv/, accessed on 17 May 2022) using Armitage's trend test as applied by the deFinetti online software (https://ihg.helmholtz-muenchen.de/cgi-bin/hw/hwa1.pl, accessed on 17 May 2022).
To identify possible intergenic and/or intragenic interactions between the variants in the analyzed loci, we used the Multifactor Dimensionality Reduction (MDR) ver. 3.0.2 software (Vanderbilt University Medical School, Nashville, TN, USA) [36].

Results
The median read depth for the five-gene panel was 639X (range 86X-1940X), and the coverage was 99.9%. Thirty-nine gene coding variants were identified in TCOF1, SALL1, EYA1, and TBX1; of them, 17 were predicted to alter the amino acid sequence (non-synonymous). These included 15 missense mutations, 1 in-frame microdeletion, and 1 in-frame microduplication that were classified as benign (B, n = 12) and likely benign (LB, n = 5) according to ACMG/AMP criteria (see Table 2). Furthermore, 22 synonymous or non-coding variants were identified in our patients (see Table 3). No variant was identified in the coding region of the HOXA2 gene.

Association Analysis
When we compared the allelic frequencies (AFs) of benign (B) or likely benign (LB) variants observed in our population with those reported in Mexican individuals from Los Angeles (1000 Genomes Project), a statistically significant difference (p < 0.05) was observed only for the p.(Ser159del) SALL1 variant, which was found specifically in microtia/OAVS patients (n = 5, 3 familial and 2 sporadic cases). Since the missense p.(Val1275Ile) SALL1 variant was present in a homozygous state in all cases and all individuals of the reference group, we did not perform an association analysis for this variant (see Table 2). Comparing the AFs of synonymous variants among our patients with those of the reference group, a statistically significant difference (p < 0.05) was observed only for the p.(Gly587=) TCOF1 variant (see Table 3).

MDR Interaction Analysis
We analyzed possible interactions among the 17 genotypes of the variants found in TCOF1, SALL1, EYA1, and TBX1 for our 49 microtia/OAVS cases and those present in the 63 individuals belonging to the reference group of a Mexican population from Los Angeles (1000 Genomes Project). Five synonymous variants (TCOF1: p.(Thr210=); SALL1: p.(Pro558=) and p.(Val190=); EYA1: p.(Ile195=); and TBX1: p.(Pro45=)) could not be subjected to MDR interaction analysis due to the lack of information in the reference group.
One significant gene-gene interaction related to the presentation of microtia/OAVS was identified: that between the non-synonymous p.(Pro1099Arg) TCOF1 variant and the synonymous p.(Leu858=) SALL1 variant (Figure 1). In this MDR interaction analysis, the balance accurate cross-validation testing result was 0.8175 and the cross-validation consistency value was 10/10. In a permutation analysis performed with 10,000 repetitions, the obtained p-value was statistically significant (p < 0.001).
Furthermore, using the network created by STRING, we added the intergenic interaction that we identified in this study between the variants p.(Pro1099Arg) TCOF1 (rs1136103) and p.(Leu858=) SALL1 (rs1965024) and incorporated in this network the following based on reported evidence for other gene-gene interactions: (a) The co-expression of SIX1 and Eya1 synergistically regulates the expression of SALL1 during kidney development [44] and also play keys roles in ear determination [45,46]. (b) In the mouse limb bud, Hoxa13 and Hoxd13 inhibit the expression of Sall1 [47]. In mouse embryonic stem cells, Sall1 appears to inhibit various Hox genes, including Hoxd13 and Gsc [48]. (c) Although it is not known whether Sall1 has a regulatory relationship with Hoxa2 in the branchial arches, the former is expressed early in head mesenchyme and then becomes restricted around the first branchial cleft in close proximity of Hoxa2, which encodes an important transcription factor in the external ear morphogenesis of mice [40]. (d) Hoxa2 coordinates the downregulation of Gsc, which acts as a transcriptional repressor in wild-type cartilage during mouse embryogenesis [49]. (e) In zebrafish, gsc downregulates the expression of bapx1 in the second pharyngeal arch [50,51]. (f) Both increased and decreased RA signaling could induce craniofacial abnor-malities, such as those found in OAVS [13]. (g) MYT1 overexpression reportedly induces the downregulation of RA receptor β (RARB), whereas mutated MYT1 does not [23] (see Figure 1).
This interaction network was built by the STRING V.11.0 software (https://stringdb.org/, accessed on 20 September 2022) and includes the proteins encoded by the five genes studied herein (in dark circles) plus those related to RA (underlined in blue) and craniofacial disorders (underlined in pink). The thickness of a gray line indicates the strength of the data compatibility based on the STRING evidence. The solid red line indicates the intergenic interaction that we herein identified between p.(Pro1099Arg) TCOF1 (rs1136103) and p.(Leu858=) SALL1 (rs1965024) variants, and the colored dotted lines represent the interactions documented in the literature.

Discussion
Since microtia/OAVS shows a heterogeneous etiology, incomplete penetrance, and variable expressivity, its clinical and molecular diagnoses have proven challenging. In the literature, the inclusion criteria for patients with this spectrum are diverse and not always well described or defined, which could be considered a limitation when interpreting and discussing results. We consider that the clinical inclusion criteria utilized herein allowed us to study a more homogeneous population and thereby avoid biases.
Some groups have previously analyzed genetic factors related to microtia/OAVS [22-24,26-28,52-56], but no previous gene-gene interaction analysis has been performed in these patients. Given the genomic differences found across populations and the high This interaction network was built by the STRING V.11.0 software (https://string-db. org/, accessed on 20 September 2022) and includes the proteins encoded by the five genes studied herein (in dark circles) plus those related to RA (underlined in blue) and craniofacial disorders (underlined in pink). The thickness of a gray line indicates the strength of the data compatibility based on the STRING evidence. The solid red line indicates the intergenic interaction that we herein identified between p.(Pro1099Arg) TCOF1 (rs1136103) and p.(Leu858=) SALL1 (rs1965024) variants, and the colored dotted lines represent the interactions documented in the literature.

Discussion
Since microtia/OAVS shows a heterogeneous etiology, incomplete penetrance, and variable expressivity, its clinical and molecular diagnoses have proven challenging. In the literature, the inclusion criteria for patients with this spectrum are diverse and not always well described or defined, which could be considered a limitation when interpreting and discussing results. We consider that the clinical inclusion criteria utilized herein allowed us to study a more homogeneous population and thereby avoid biases.
Some groups have previously analyzed genetic factors related to microtia/ OAVS [22][23][24][26][27][28][52][53][54][55][56], but no previous gene-gene interaction analysis has been per-formed in these patients. Given the genomic differences found across populations and the high prevalence of this disorder in Latin-American populations worldwide, it is important to obtain a more precise knowledge in patients with microtia/OAVS of this ethnic origin. Gaining a better understanding of the genetic etiology of this malformation could enable clinicians to provide more accurate genetic counseling to families [17][18][19][20][21]37].

Association Analysis
When we performed the association study between the variants identified in our patients and the reference group, we identified that the in-frame p.(Ser159del) microdeletion of SALL1 appears to be a risk factor for microtia, as evidenced by the statistically significant difference in the AF between our patients and the reference group and the presence of this variant in only our microtia/OAVS cases. Similar associations have been identified between variants of certain genes and the risk of developing malformations, such as congenital heart defects, biliary atresia, pyloric stenosis, hypospadias, and microtia [53,54,57].
The genome-wide association study (GWAS) approach has been successful in identifying new susceptibility loci for common structural congenital defects, such as oral clefts, congenital heart defects, biliary atresia, pyloric stenosis, hypospadias, craniosynostosis, and clubfoot [57]. However, congenital ear abnormalities, including anotia/microtia, have not previously been addressed by GWAS. The sole exception to this was a genome-wide linkage analysis performed on two families with OAVS [58]. In one family, the authors identified a highly suggestive linkage to a region harboring the GSC (Goosecoid homeobox) gene, which was considered to be a good candidate gene for this entity. However, coding-region changes and gross rearrangements were excluded in these two OAVS familial cases and in 120 additional sporadic cases [58].
Synonymous variants may influence the development of various human diseases, including birth defects [59,60]. A statistically significant difference was identified for the p.(Gly587=) TCOF1 variant. The AF for this variant was greater in the reference group, suggesting that it confers protection against or a decreased risk for microtia/OAVS in our Mexican population. There is no single mechanism by which a synonymous change could exert a biological effect. Accumulating evidence shows that biological systems take advantage of the degeneracy of the genetic code to control gene expression, protein folding efficiency, and coordinated expression across several gene families. The most obvious and well-characterized mechanism by which synonymous changes can exert a deleterious biological effect is by perturbing pre-mRNA splicing [59]. However, a synonymous variant could also be in linkage disequilibrium with deleterious functional variants located nearby. To our knowledge, the p.(Gly587=) of TCOF1 could be the first described synonymous variant associated with the microtia/OAVS trait.

MDR Interaction Analysis
In a very recent review of genetic and non-genetic factors involved in the development of microtia/OAVS, there was no mention of data related to gene-gene interactions [32]. As the interaction analysis for benign and/or synonymous variants could suggest a genetic protection or susceptibility factor for complex traits [34], such as ear malformations, we decided to apply an MDR analysis, which is a nonparametric model-free method for identifying epistasis using the identified variants [61]. This strategy identified a single statistically significant intergenic interaction between the non-synonymous p.(Pro1099Arg) TCOF1 variant and the synonymous p.(Leu858=) SALL1 variant (Figure 1). At the statistical level, combining the TCOF1 and SALL1 genotypes allowed us to discriminate between cases and controls, indicating that there is an interaction between these two genes. Although this interaction has not previously been reported, the TCOF1 and SALL1 proteins are related to craniofacial disorders [23] (Figure 1). The available evidence generally supports the involvement of the analyzed genes/proteins and various other genes/proteins in the development of microtia/OAVS. Notably, SALL1 interacts with most of these genes/proteins and thus appears to play a central role. In addition, it was the gene in which the most non-synonymous variants were identified in our group of patients (Figure 1), suggesting that SALL1 warrants future study in this regard.
Moving forward, the inclusion of a larger number of microtia/OAVS-related genes or the use of whole-exome sequencing should identify new variants that can be considered in future gene-gene interaction studies [62]. Epigenetic inheritance also has been suggested as a possible pathogenic mechanism [5]. For example, a histone acetylation-dependent imbalance in the allelic expression of the key craniofacial development gene, BAPX1 (also called NKX3-2, MIM*602183), was observed in five patients with OAVS [51]. Thus, the contribution of epigenetic mechanisms to the etiology of microtia/OAVS deserves attention in future genetic studies.
The published evidence and our present findings collectively support the complexity of ear embryogenesis and microtia/OAVS development, which involves the temporal and spatial expression of different proteins and signaling by multiple pathways. We did not identify any pathogenic variant in the five studied genes, but we found a gene-gene interaction between TCOF1 and SALL1. This highlights the need to identify other genes and genotype interactions that contribute to the etiology of craniofacial disorders, including microtia/OAVS [33]. Future research is also needed to assess the involvement of the RA pathway in the genetic etiology of microtia/OAVS.

Conclusions
Although gene-gene interactions are known to play an important role in the etiology of many complex diseases, no previous study has addressed gene interactions in patients with microtia/OAVS. Our finding of a gene interaction between TCOF1 and SALL1 in a group of Mexican patients with this entity supports the complex nature of ear embryogenesis and the development of microtia/OAVS. Further research is warranted, such as the inclusion of more candidate loci, which should lead to the identification of new gene-gene interactions underlying microtia/OAVS.