Rare single-nucleotide variants in oculo-auriculo-vertebral spectrum (OAVS).

BACKGROUND: Oculo-auriculo-vertebral spectrum (OAVS) is a craniofacial developmental disorder that affects structures derived from the first and second pharyngeal arches. The clinically heterogeneous phenotype involves mandibular, oral, and ear development anomalies. Etiology is complex and poorly understood. Genetic factors have been associated, evidenced by chromosomal abnormalities affecting different genomic regions and genes. However, known pathogenic single-nucleotide variants (SNVs) have only been identified in MYT1 in a restricted number of patients. Therefore, investigations of SNVs on candidate genes may reveal other pathogenic mechanisms.
METHODS: In a cohort of 73 patients, coding and untranslated regions (UTR) of 10 candidate genes (CRKL, YPEL1, MAPK1, NKX3-2, HMX1, MYT1, OTX2, GSC, PUF60, HOXA2) were sequenced. Rare SNVs were selected and in silico predictions were performed to ascertain pathogenicity. Likely pathogenic variants were validated by Sanger sequencing and heritability was assessed when possible.
RESULTS: Four likely pathogenic variants in heterozygous state were identified in different patients. Two SNVs were located in the 5'UTR of YPEL1; one in the 3'UTR of CRKL and one in the 3'UTR of OTX2.
CONCLUSION: Our work described variants in candidate genes for OAVS and supported the genetic heterogeneity of the spectrum.

INTRODUCTION

Oculo‐auriculo‐vertebral spectrum (OAVS, [OMIM :164210]) is a developmental disorder characterized by alterations in craniofacial morphogenesis mainly involving structures derived from the first and second pharyngeal arches (Beleza‐Meireles, Clayton‐Smith, Saraiva, & Tassabehji, 2014). The clinically heterogeneous and complex OAVS phenotype includes hemifacial microsomia, asymmetric ear anomalies, such as preauricular tags and microtia, ocular defects (e.g., epibulbar dermoid and microphthalmia) as well as vertebral anomalies (Beleza‐Meireles et al., 2014). Other important features may be associated, for example, cardiac, genitourinary malformations, and central nervous system defects (Barisic et al., 2014).

OAVS' etiology is complex and still not completely understood. Several environmental factors have been associated with the spectrum, such as multiple gestations, gestational diabetes, and the use of vasoactive drugs in gestation (Barisic et al., 2014; Werler, Sheehan, Hayes, Mitchell, & Mulliken, 2004). Genetic factors have also been revealed due to the existence of familial cases with autosomal dominant (Ballesta‐Martínez et al., 2013; Beleza‐Meireles et al., 2015; Tasse et al., 2005, 2007; Vendramini‐Pittoli & Kokitsu‐Nakata, 2009) and recessive inheritance (Beleza‐Meireles et al., 2015; Farra‐Awwad, 2011). Furthermore, genetic heterogeneity in OAVS is supported by the occurrence of several chromosomal abnormalities involving different genomic regions and genes (Beleza‐Meireles et al., 2015; Bragagnolo et al., 2018; Rooryck et al., 2010). However, to date, MYT1 (OMIM: 600379) has been the only gene implicated in this condition with two variants found in a cohort of 169 patients with the spectrum (Lopez et al., 2016) and other identified in one of 57 patients part of the Brazilian cohort studied here (Berenguer et al., 2017). Furthermore, functional studies by transient knockdown of MYT1's homolog in zebrafish revealed craniofacial cartilage alterations (Lopez et al., 2016).

Thus, in the present study, we sequenced exonic and untranslated regions (UTR) of the 10 selected candidate genes in 73 OAVS' patients allowing the identification of four likely pathogenic rare heterozygous variants in three of these genes (YPEL1, CRKL, and OTX2).

MATERIALS AND METHODS

Editorial Policies and Ethical Considerations

Written informed consent was obtained from patients or their parents as approved by the local ethics committee.

Enrollment

We studied a cohort composed by 73 patients, 65 index patients, and eight affected relatives. The patients were included when they fulfilled the minimal OAVS diagnostic criteria modified from Tasse et al. (2005) and Gougoutas, Singh, Low, and Bartlett (2007): hemifacial microsomia and unilateral or asymmetric microtia and/or preauricular tags or pits. Patients with skeletal anomalies, affecting bones other than ribs and spine, were excluded.

Candidate genes selection

Ten genes were selected for sequencing when they filled at least one of the three criteria: (a) genes previously described in the literature with possibly pathogenic variants and/or indicated as candidates for the main phenotypes of OAVS; (b) genes encompassed by potentially pathogenic copy‐number variations (CNVs) previously detected in our cohort (Bragagnolo et al., 2018); (c) genes associated with development of the first and second pharyngeal arches.

According to these criteria, the following genes were selected for sequencing studies: YPEL1 (OMIM: 608082, NM_013313.5), MAPK1 (OMIM: 176948, NG_023054.2), CRKL (OMIM: 602007, NG_016354.1), OTX2 (OMIM: 600037, NG_008204.1), GSC (OMIM: 138890, NG_034111.1), HMX1 (OMIM: 142992, NG_013062.2), NKX3‐2 (OMIM: 602183, NG_023192.1), MYT1 (OMIM: 600379, NM_004535.3), PUF60 (OMIM: 604819, NG_033879.1), and HOXA2 (OMIM: 604685, NG_012078.2).

Target gene sequencing

Sequencing panel was custom designed with the Ion AmpliSeq ™ Designer online tool (www.ampliseq.com) targeting exons and UTR from the 10 selected genes (padding: ± 5bp). This design allowed the analysis of 165 amplicons varying from 125 to 375 base pairs.

DNA samples were isolated from peripheral blood using Gentra Puregene Kit (Qiagen‐ Sciencesm). Libraries were constructed using Ion AmpliseqTM Library Kit 2.0 and normalized using Ion Library EqualyzerTM Kit according to the manufacturer's protocol. Amplified libraries were then submitted to emulsion polymerase chain reaction (PCR) on the Ion OneTouchTM 2 system using the Ion PGMTM Hi‐QTM OT2 Kit (Thermo Fisher) following manufacturer's protocol. Ion sphere particles were enriched on the Ion OneTouch ES (Thermo Fisher), loaded on an Ion 318 chip (Thermo Fisher) and sequenced on the Ion Torrent PGMTM (Thermo Fisher).

Exons and UTR from the 10 selected candidate genes were sequenced in the cohort of 73 patients in two separate runs.

Data analysis

Data were processed using Ion Torrent Suite (Thermo Fisher) software that applied quality controls and removed poor signal reads. Variant initial annotation and quality control were performed using Ion ReporterTM software (Thermo Fisher). All variants were annotated according to the Human reference genome GRCh37/hg19.

Filtering detected variants

Some of the amplicons also contemplated a few intronic and downstream regions; however, in this study, only single‐nucleotide variants (SNVs) in exons and UTR were selected. Variants with low sequencing coverage (<30×) (Tarabeux et al., 2014) were filtered out. Rare variants were further selected and defined as those with a minor allele frequency (MAF) <1% in The Genome Aggregation Database (gnomAD, https://gnomad.broadinstitute.org/), using 15,708 genome samples (Lek et al., 2016) and in a Brazilian reference cohort database with 609 genome samples (http://abraom.ib.usp.br/) (Naslavsky et al., 2017).

Variants interpretation

Filtered SNVs were annotated using the WGSA annotation pipeline (v0.76) (Liu et al., 2016) that integrates several public databases and prediction tools, such as FATHMM‐XF (Rogers et al., 2017); CADD (v1.4) (Kircher et al., 2014); DANN (Quang, Chen, & Xie, 2015); ClinVar (accessed on 2018/09/30) (Landrum et al., 2014); GWAS catalog (accessed on r2018‐11–26) (Welter et al., 2014); GenoCanyon (v1.0.2) (Lu et al., 2015); gnomAD (r2.1) (Karczewski et al., 2019); TargetScan v7.2 (Agarwal, Bell, Nam, & Bartel, 2015); miRBase (v3) (Griffiths‐Jones, 2004); Enhancer Finder; and also dbNSFP (v4.0b1c) (Liu, Jian, & Boerwinkle, 2011) that includes SIFT; Polyphen‐2; and Mutation Taster. RegulomeDB (v1.1) (Boyle et al., 2012) was also consulted.

Pathogenicity of selected variants was ascertained based on prediction from four in silico tools: SIFT; Polyphen‐2; FATHMM‐XF; and Mutation Taster depending on the variant type and regions. For missense variants, SIFT, Polyphen‐2, FATHMM‐XF, and Mutation Taster were used. For synonymous variants and variants in UTR, only FATHMM‐XF and Mutation Taster were applied.

SNVs were considered likely pathogenic when predicted as pathogenic by all queried algorithms. ClinVar and GWAS catalog were consulted to investigate if these variants had been previously associated with diseases. Additionally, annotation highlighting regulation features (GenoCanyon, TargetScan, miRBase, Enhancer Finder, and RegulomeDB) were considered to improve pathogenic interpretation of the variants.

Sanger sequencing

Variants considered as likely pathogenic were validated by Sanger sequencing of PCR‐amplified products. Primers used for PCR amplification are available upon request. PCR products were purified and sequenced using the BigDye® Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems) and the Genetic Analyzer 3130xl. Chromatograms were visualized using Chromas (https://technelysium.com.au/wp/chromas/). Variants were assessed in patient's parents to verify heritability, when DNA samples were available. No affected relatives were referred. Patients 1 and 3 have unaffected half siblings, patient 2 has one unaffected brother and patient 4 does not have siblings. However, DNA samples for these individuals were not available.

Expression analysis of likely pathogenic variants harboring genes

Gene expression analysis was performed by real‐time (RT)‐qPCR for the genes with likely pathogenic variants. Whole blood RNA was isolated from patients with likely pathogenic variants and from 11 controls using the PAXgene BloodRNA MDx Kit (Qiagen‐Sciences). cDNA was synthesized using High Capacity cDNA Reverse Transcription kit (Thermo Fisher Scientific). TaqMan assays (Thermo Fisher Scientific) were selected for genes harboring likely pathogenic variants and presenting whole blood expression (YPEL1 and CRKL). GAPDH and ACTB expression were used as internal controls. Assays and whole blood expression are detailed in Table S1

Relative differences in expression were quantified using the 2− ∆∆Crt method with the patient's sample as calibrator. Gene expression was considered modified when controls' 2− ∆∆Crt mean presented a minimum fold regulation of ±2 when compared with the patient's 2− ∆∆Crt.

RESULTS

Clinical and cytogenetic evaluation

Normal G‐banding karyotypes at 550‐band level were found in all patients but one, who presented with a balanced inversion on chromosome 12–46,XX,inv(12)(q15q24.1). Array analysis revealed 14 CNVs (nine duplications and five deletions) with potential pathogenicity to OAVS in 12 of our patients. Clinical and cytogenomic evaluation for all but three patients had previously been reported elsewhere (Bragagnolo et al., 2018).

Detection of SNVs and their clinical impacts

Each of the two 318 chip used for the sequencing generated an average of 4.2 million reads which correspond to a mean coverage of 451 × per sequenced sample. Panel design enabled the sequencing of approximately 95% of exons and UTR, not contemplating only about 1,700 bp from these regions.

A total of 243 distinct variants were detected. About 123 distinct SNVs remained after filtering by coverage and selecting variants in exons and UTR. After filtering for variants' frequency in population cohorts (gnomAD and AbraOM), 57 rare SNVs were selected and proceeded to downstream functional interpretation.

According to the in silico pathogenicity prediction, four of 57 SNVs were considered as likely pathogenic, all of them being in heterozygous state in four different patients (Table 1). The four patients' carriers of the identified SNVs do not present any CNVs relevant to their phenotypes. (Bragagnolo et al., 2018) (Table 2).

Table 1

Likely pathogenic rare SNVs identified in OAVS patients from our sample

P	Gene	Variant Nomenclaturea	Variant IDb	Gene Element	Inheritance	MAF (gnomAD)	MAF (AbraOM)	FATHMM‐XF	Mutation Taster	CADD phred	DANN score	Geno Canyon
1	YPEL1	NC_000022.10:g.22065172C > G	rs1029869759	5’UTR	NI	—	—	0.976841	Disease causing	14.74	0.87	0.99
2	YPEL1	NC_000022.10:g.22065186G > A	rs186536735	5’UTR	Maternal	0.1%	0.2%	0.972	Disease causing	10.74	0.90	0.99
3	CRKL	NC_000022.10:g.21307535C > A	rs139209821	3’UTR	NI	0.03%	—	0.693715	Disease causing	19.42	0.86	0.99
4	OTX2	NC_000014.8:g.57268275T > G	rs916081360	3’UTR	NI	0.006%	0.08%	0.655947	Disease causing	20.6	0.85	0.99

Abbreviations: MAF, minor allele frequency; NI, No information, DNA from parents unavailable; P, patient.

According to the Human Genome Variation Society.

Variant ID from single nucleotide polymorphism database (b151).

Table 2

Clinical description, sex, age, and CNVs of the patients

P	Gene	Variant IDa	Sex, age	Clinical features	CNVsb	CNVs in other patientsc
1	YPEL1	rs1029869759	M, 4 y	Facial asymmetry, coloboma of upper left eyelid and left epibulbar dermoid, microtia I on the left ear and microtia II with the absence of lobe on the right, bilateral preauricular tags, and left nostril with slit in its anterior portion	No	Deletiond
2	YPEL1	rs186536735	F, 8 y	The patient presented hemifacial microsomia, microtia III on the right ear with atresia of external acoustic meatus and middle ear malformation, microtia II on the left ear and the absence of ossicular chain, moderate left sensorineural hearing impairment and retrognathia, right lung agenesis, single kidney, cervical anomalies, and situs inversus totalis.	No	Deletiond
3	CRKL	rs139209821	M, 14 y	Hemifacial microsomia, alopecia areata, right epibulbar dermoid, microtia III with moderate to severe mixed hearing loss, micrognathia, and scoliosis.	No	Deletiond
4	OTX2	rs916081360	F, 1 m	Facial asymmetry, microtia III, triangular face, cervical anomalies, horseshoe kidney, hydrocephaly, esophagus atresia, and dextrocardia	No	No CNVs

Abbreviations: CNVs, copy number variations; F, female; M, male; m, months; P, patient; y, years;.

Variant ID from single nucleotide polymorphism database (b151).

Pathogenic, likely pathogenic, or CNVs of uncertain significance identified in the patients.

CNVs encompassing these genes identified in Brazilian patients from Bragagnolo et al. (2018).

A deletion including YPEL1 and CRKL was identified in one patient from Bragagnolo et al. (2018).

Among these SNVs, two were in the 5’UTR of YPEL1; one was located in the 3’UTR of CRKL, and one was identified in the 3’UTR of OTX2 (Table 1). All four variants were validated by Sanger sequencing.

Three of four likely pathogenic SNVs presented CADD scores above 10 and the SNV located in the OTX2 presented score above 20, indicating that these variants are, respectively, among the top 10% and 1% more likely damaging variants in the genome. DANN scores for all variants were above 0.8 also indicating a high probability of being deleterious.

Considering the annotations of regulation features, GenoCanyon revealed a high score for all four variants (>0.9). Enhancer Finder indicated that none of the variants were located within predicted enhancer sites for general development, brain, heart, or limb. According to RegulomeDB, three of these variants do not appear to affect their DNA binding sites to regulatory proteins, but the variant located in the CRKL is likely to affect the binding of BACH1 protein as revealed by Chip‐Seq experiments in K562 cells (https://www.encodeproject.org/experiments/ENCSR000EGD/).

Gene expression analysis of YPEL1 and CRKL showed unmodified transcription levels in the patients when compared to controls (Figure S1–S3).

A summary of the SNVs considered as likely pathogenic and the patients' phenotypes are presented in Table 2.

DISCUSSION

Although little is known about the genetic basis of OAVS, recent studies leveraged evidence of a complex etiology, involving several genes and genetic factors (Lopez et al., 2016; Rooryck et al., 2010; Tasse et al., 2005; Vendramini‐Pittoli & Kokitsu‐Nakata, 2009). We aimed to identify rare variants with pathogenic potential in a cohort of 73 Brazilian patients affected by the spectrum. Given the Brazilian's historical mixed ancestry, the use of a reference population cohort composed of Brazilian individuals was crucial to filter variants based on their population frequencies.

We identified four likely pathogenic rare SNVs in three (YPEL1, CRKL, and OTX2) of 10 candidate genes in four unrelated patients, among the patients from our cohort. Both variants in YPEL1 were in its 5’UTR, while the other two variants were localized at 3’UTR. One variant was maternally inherited, suggesting incomplete penetrance, possibly with several genetic and environmental factors leading to phenotype. All variants were previously annotated by single nucleotide polymorphism database; however, no clinical features have been associated with any of them.

None of the four patients presented altered karyotypes or likely pathogenic or pathogenic CNVs at chromosomal microarray investigation. This reinforces that the SNVs detected may play an important role in the development of the disease in these patients.

YPEL1, located in the 22q11.21 region, encodes for a not well‐characterized protein called yippee‐like 1, that is believed to participate in cell division (Farlie et al., 2001). In mice, Ypel1 is expressed in the pharyngeal arches, indicating that this gene may play a role in the development of the face (Farlie et al., 2001). A zebrafish model with Ypel1 protein level knocked down presented several craniofacial defects, including underdeveloped jaw (Aerts et al., 2006). CNVs involving this gene have been previously associated with OAVS and the patients affected presented preauricular tags (Xu, Fan, & Siu, 2008), epibulbar dermoid (Tan et al., 2011), as observed in our patient P1, and hemifacial microsomia (Tan et al., 2011), as observed in both patients P1 and P2.

CRKL is also localized in the 22q11.21 region and encodes for a tyrosine kinase important for signaling pathways. The snoopy mice strain, carrier of a Crkl mutation that results in loss of protein expression, presents severe craniofacial anomalies, including mandibular and maxillary hypoplasia, microstomia, micrognathia, and ocular anomalies with variable expressiveness (Miller et al., 2014). Accordingly, our patient P3 presents hemifacial microsomia and micrognathia.

OTX2, localized in the 14q22.3, encodes for a homeobox protein that acts as a transcription factor. Mice with heterozygous mutation in this gene exhibit several craniofacial anomalies such as micrognathia, agnathia, and anophthalmia, showing that this gene plays a role in neural crest development in mice (Matsuo, Kuratani, Kimura, Takeda, & Aizawa, 1995). In a family with five individuals affected by hemifacial microsomia, a 1.3 Mb duplication of chromosome 14q22.3 has been identified (Zielinski et al., 2014). In this study, OTX2 was suggested as the candidate gene to the craniofacial phenotype. Furthermore, duplications involving this gene have been previously related to OAVS (Ballesta‐Martínez et al., 2013). Even though OTX2 is associated with ocular anomalies, such as microphthalmia and anophthalmia (Jones et al., 2016; Lonero et al., 2016; Shimada, Takagi, Nagashima, Miyai, & Hasegawa, 2016; Slavotinek et al., 2015), our patient (P4) does not present ocular alterations. A 120 kb deletion encompassing this gene was described in a patient with mandibular dysostosis, micrognathia, and external acoustic meatus agenesis. Similar to our patient, this case also did not present ocular alterations (Latypova et al., 2016).

Although the identified variants were predicted to be likely pathogenic, their functional consequences remain unknown. Our expression studies did not show a direct impact of the variants on whole blood transcriptional level. However, expression perturbations during the target tissues development and posttranslational modifications cannot be ruled out. 5’UTR and 3’UTR contain sequence elements that are important for gene expression and regulation (Mignone, Gissi, Liuni, & Pesole, 2002; Steri, Idda, Whalen, & Orrù, 2018). Synthesized mRNA is regulated by elements present in the UTR and by miRNAs interacting with their targets in the 3’UTR (Mignone et al., 2002; Steri et al., 2018).

To predict whether these variants could impact a regulatory mechanism, we considered the annotations and predictions from several tools. According to GenoCanyon, all variants presented a high score indicating a functional potential to these genomic locations (Lu et al., 2015). RegulomeDB predicted that the variant in the CRKL is likely to affect the binding of BACH1 protein to this genomic region. This binding site was evidenced by Chip‐Seq experiments performed in K562 cells, a human erythromyeloblastoid leukemia cell line (https://www.encodeproject.org/experiments/ENCSR000EGD/). However, this binding site is not replicated in human embryonic stem cell line H1 (https://www.encodeproject.org/experiments/ENCSR000EBQ/). The BACH1 (BTB and CNC homology) is mapped to the 21q21.3 region and encodes a transcription factor that contains leucine zipper and BTB‐zinc finger domains (Blouin et al., 1998). The zebrafish homolog Bach1b is ubiquitously expressed during early embryonic development (Jiang et al., 2017). In addition, the mice homologous gene was expressed in the pharyngeal arch (Reymond et al., 2002; Yokoyama et al., 2009) Therefore, we hypothesize that the variant identified in the CRKL could be interfering with the BACH1 protein binding during formation and development of the pharyngeal arches, culminating in CRKL expression dysregulation and the patient's phenotype.

It is worth mentioning that discordances between predicted and actual effect of variants have been reported (Miosge et al., 2015). These discordances can result in false‐positive clinical interpretations of variants or in the assumption of unknown significance to variants that have functional impacts. Additionally, the algorithms for the predictive tools differ from each other, which can cause different interpretations of the same variant. Thus, it is important that multiple tools predictions be considered together for the interpretation of genetic variants (Richards et al., 2015). In this work, when the predictions of the different tools were unanimous, we considered as strong evidence of the pathogenicity of the variant. This unanimity criterion made the classification more restricted and therefore more reliable.

In conclusion, our work described novel variants in candidate genes for OAVS. However, functional studies are essential to determine the potential role of these variants in OAVS' etiology. The fact that each variant was only identified in a single patient supports the genetic heterogeneity of the spectrum and strengthens the indications of a multifactorial inheritance with different factors involved in the disease, including environmental aspects.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

Click here for additional data file.

Click here for additional data file.