Rare missense variant p.Ala505Ser in the ZAK protein observed in a patient with split-hand/foot malformation from a non-consanguineous pedigree.

Introduction

During embryonic development, an apical ectodermal ridge (AER) forms at the distal end of each limb bud and functions as a focus for regulatory molecules to guide limb development.[1] Dysregulation of the expression and activity of proteins in two zones within the AER, called the progression zone and the zone of proliferating activity, is thought to play a role in the etiology of five nonsyndromal forms of split-hand/foot malformation (SHFM) (Figure 1).[2-6] Disruption of these signaling processes results in ectrodactyly, syndactyly, polydactyly, or other alterations of the hands and feet, as previously reviewed.[1] Despite identification of multiple SHFM-causative loci, some patients present with a phenotype suggestive of SHFM but lack mutations in the associated genes. Both Mendelian and non-Mendelian inheritance have been reported for various etiologies of SHFM.

Figure 1.

Signaling within the developing limb bud relevant to etiologies of split-hand/foot malformation (SHFM) with possible associated hearing loss. The figure is oriented from proximal (left) to distal (right). Blue font denotes signaling thought to occur within the apical ectodermal ridge. (ΔNp63α is a transcriptional regulator protein.) Types of SHFM associated with hearing loss in this pathway are listed in red. See Kantaputra and Carlson (2018) for a thorough review.[8]

Pathogenic variant p.Phe368Cys in leucine zipper containing kinase AZK (ZAK) exhibited autosomal recessive inheritance in two separate consanguineous Pakistani cohorts.[7] Those unrelated consanguineous Pakistani families were studied by Spielmann et al. (2016) in a search for novel pathogenic loci. These pedigrees exhibited limb abnormalities consistent with SHFM and also had associated hearing loss.[7] The pedigrees suggested autosomal recessive inheritance, so autozygosity mapping and single nucleotide polymorphism (SNP) array genotyping were used to determine that chromosomal region 2q31 was likely to contain the disease locus.[7] Subsequent exome sequencing revealed missense variant p.Phe368Cys within a sterile alpha motif (SAM) domain of ZAK to be the causative mutation in these families.

The authors generated a ZAK knockout mouse model which caused all mice to die during embryonic development.[7] Deletion of the SAM domain (which contains p.Phe368Cys) by CRISPR/Cas caused a hindlimb defect that was similar to the defect observed in mice with mutant Trp63, a known SHFM gene.[7] These studies demonstrate the relevance of the SAM domain within ZAK to SHFM.

Given the recent demonstration of the role of ZAK in limb genesis and its association with hearing loss in two pedigrees, we hypothesized that a mutation in ZAK could be present in additional unrelated patients with ectrodactyly and associated hearing loss.

Methods

Patient samples

Biobanking of DNA samples for sequencing was performed with review and approval of the Institutional Review Committee of Self Regional Healthcare (Greenwood, SC, USA).

Phenotypic and historical data representing the entire cohort of SHFM patients at the Greenwood Genetic Center (GGC; Greenwood, SC, USA) were screened for possible autosomal recessive inheritance, greater expressivity in the feet than hands, and associated hearing loss (to best match the phenotype observed in the Pakistani cohort). Thirty-eight patients within the GGC’s cohort matched the selection criteria, had pedigrees that did not exclude autosomal recessive inheritance, and provided consent for DNA sequencing. All participants provided written informed consent.

DNA sequencing

ZAK-targeted PCR and subsequent Sanger DNA sequencing were performed on banked genomic DNA samples from the selected 38 individuals in an international cohort of unrelated patients with SHFM. PCR was used to amplify patient DNA as described previously.[9] The DNA primers targeting ZAK gene isoform 1 (as listed in Ensembl) are shown in Supplemental Table 1. Exon sequences within non-canonical isoforms (e.g., Ensembl ZAK transcript 201) were targeted with an additional set of primers (Supplemental Table 2). All RNA primers were conjugated to an appropriate M13 tag (forward or reverse).

ExoSAP-IT (Affymetrix, Santa Clara, CA, USA) was used to degrade contaminants, and the purified sample was labeled with BigDye (Thermo Fisher Scientific, Grand Island, NY, USA) as described.[10] BigDye Terminator v1.1 and v3.1 5× Sequencing Buffer (Thermo Fisher Scientific) were used to stabilize the sequencing reagents, which were subjected to purification using Qiagen DYE-EX 96-well processing plates (Qiagen, Venlo, Limburg, the Netherlands) to remove the BigDye terminators. Capillary electrophoresis using a 3730 DNA Analyzer (Thermo Fisher Scientific) was performed and the resulting files analyzed on the DNASTAR platform (DNASTAR, Madison, WI, USA).

Structural modeling for ZAK protein

The missense variant of interest, p.Ala505Ser, was not located in a conserved domain and no definitive template could be identified for homology modeling. Thus, protein models were constructed using Phyre2[11] and RaptorX[12] servers, which both use combination approaches to build protein models; Phyre2 uses combination ab initio threading and homology approaches, whereas RaptorX uses multi-template threading, template detection, and binding site prediction.

Our model for ZAK protein folding was further tested for predictions of folding free energy change by using several web servers, including mCSM, SDM, DUET, SAAFEC, and I-Mutant as previously described.[13-17]

Evolutionary conservation analysis

Multiple sequence alignment (MSA) among different species was performed to analyze the evolutional conservation of the wildtype (WT) residues involved in the mutation. The sequences of ZAK proteins from nine species (Figure 3e) were collected from UniProt. These sequences were then subjected to MSA using the T-Coffee server (https: //www.ebi.ac.uk/Tools/msa/tcoffee/).

Targeted genomic microarray analysis

DNA isolated from patient CMS15791 with missense variant p.Ala505Ser was analyzed using the Oxford Gene Technology Syndrome Plus Microarray system (Begbroke, UK). Analysis was performed using 105K DNA oligonucleotides distributed across the genome, with increased coverage at regions known to be relevant to SHFM, including 85 recognized microdeletion or microduplication syndromes and 41 unique sub-telomeric regions.

Results

DNA sequencing reveals two distinct missense mutations in ZAK

Two missense mutations in exon 18 of the ZAK gene were observed in 11 patients, and another 2 patients were found to have a silent mutation in exon 15, as listed in Table 1 and summarized in Figure 2, which portrays the canonical isoform of ZAK (Ensembl transcript 201). Given the low prevalence of SHFM, we predicted that only variants with rare allele frequencies would explain the phenotype. Thus, variant allele frequencies obtained from the University of California San Diego Human Genome Browser Database (UCSD HGB) are reported in Table 1. Of the 11 patients with missense mutations in exon 18, variant c.1515 G>T (p.Ala505Ser) in one patient was unreported in the USCD database (accessed April 2019) and was presumed to be rare. Ten other patients (26% of the cohort) contained variant c.1592 C>G, encoding p.Ser531Leu; however, 37% of the healthy population also contains this variant. Because all of the missense mutations observed in the cohort were present in a domain of ZAK with unknown function and not in the SAM domain associated with SHFM, in silico approaches were used to assess the potential effect upon protein stability.

Table 1.

Overview of all variants observed in the cohort of 38 patients within ZAK exons 15 and 18.

Exon	Alteration	Allele Freq.	refSNP (rs)#number	# of Patients: CMS1 number(s)	Predicted Protein Change
15	c.1263 T>C	7.488%	rs3739103	2 patients: 8763, 16898	synonymous (p.Ile421Ile)
18	c.1515 G>T	N/A	none reported	1 patient: 15791	p.Ala505Ser
18	c.1592 C>T	37.700%	rs3769148	10 patients: 3136, 3163A, 4699, 4798S, 7986S, 9710, 9924, 13206, 13795, 14980	p.Ser531Leu

Allele frequencies are reported in University of California San Diego Human Genome Browser Database and reading frame based on Ensembl Database. Certified Medical Service (CMS) numbers designate patients. N/A, not applicable.

Figure 2.

Overview of variants observed in ZAK canonical isoform 1. Exon 1 is an untranslated region (UTR) and exons 2 to 10 (numbers across top of figure) encode the protein kinase domain. The sterile alpha motif (SAM) spanning residues 339 to 407 (pink) is also highlighted in the ZAK structure because the SAM is implicated in the known causes of SHFM.

p.Ala505Ser results in mild destabilization of ZAK

In silico analysis focused on comparing the effects of the reported pathogenic substitution p.Phe368Cys with our observed rare substitution, p.Ala505Ser. Because Phyre2 uses a combination of ab initio threading and homology approaches, whereas RaptorX uses multi-template threading, template detection, and binding site prediction, the final models were compared and are shown in Figure 3. The WT sequences for ZAK are shown in Figure 3a and 3c, as generated by Phyre2 and RaptorX, respectively. Variant p.Ala505Ser is shown in Figure 3c and 3d, with a qualitative impact upon the predicted tertiary structure of the protein in the RaptorX models (Figure 3c and 3d). Next, to quantitatively compare the effects of these variants on protein stability, the models were used to calculate the predicted binding free energy change (ΔΔG), shown in Table 2. The known pathogenic variant p.Phe368Cys was confirmed to be destabilizing by most servers (median =−1.83 kcal/mol). Similarly, the observed rare substitution p.Ala505Ser was predicted by all servers to be destabilizing (median =−0.7 kcal/mol), albeit with lesser magnitude than p.Phe368Cys.

Figure 3.

ZAK protein structural changes caused by variant p.Ala505Ser. (a) Predicted wildtype ZAK protein model with highest score from Phyre2.[11] (b) Predicted Ala505Ser ZAK protein model with highest score from Phyre2.[11] (c) Predicted wildtype ZAK protein model with highest score from RaptorX.[12] (d) Predicted Ala505Ser ZAK protein model with highest score from RaptorX.[12] Red arrows denote qualitative differences in protein between WT and p.Ala505Ser. (e) Multiple sequence alignment of ZAK proteins among nine species (available from UniProt) in which variants were observed, with residues corresponding to the known pathogenic substitution p.Phe368Cys (top) and this variant p.Ala505Ser (bottom) shown in boxes.

Table 2.

Predicted binding free energy change (ΔΔG, kcal/mol) for variants seen in cohort of patients (positive values indicate that the substitution stabilizes the protein; negative values indicate destabilization).

Web Server	ΔΔG Phe368Cys (kcal/mol)	ΔΔG Ser531Leu (kcal/mol)	ΔΔG Ala505Ser (kcal/mol)
mCSM[13]	−1.83	+0.3	−0.7
SDM[14]	+0.72	+1.4	−2.25
DUET[15]	−1.71	+1.0	−0.48
SAAFEC[16]	−3.43	−4.6	−0.77
I-Mutant[17]	−1.89	N/A	−0.32

N/A, not applicable.

The previously reported pathogenic substitution p.Phe368Cys, occurs in a hydrophobic core of the SAM domain of ZAK, in a side chain that is buried and surrounded by several other hydrophobic residues. Substitution to cysteine is predicted to be unfavorable for such a hydrophobic environment and would affect the stability of the SAM domain. The folding free energy calculations (Table 2) support our analysis and indicate a large destabilizing effect caused by substitution p.Phe368Cys. Furthermore, the Phe residue at position 368 is highly conserved in all species tested in MSA (Figure 3e), indicating its functional importance in the ZAK protein.

In contrast, the Ser 531 side chain was totally exposed in the model, having no contact with neighboring groups. The predicted effect of p.Ser531Leu on protein stability was mixed, with some servers predicting destabilization and others stabilization. Thus, we presumed that residue 531 could tolerate these substitutions. Indeed, MSA (data not shown) indicated that residue 531 could be a charged residue (Arg) or a hydrophobic residue (Ala) in different species. Furthermore, the high allele frequency of this variant renders it unlikely to be implicated in a rare disease.

Finally, p.Ala505Ser is located in an alpha helix and is predicted to be exposed to water in the liquid matrix of cytoplasm. The folding free energy calculations (Table 2) consistently indicated that p.Ala505Ser is destabilizing, which is logical given the substitution of a hydrophobic amino acid with a hydrophilic amino acid. However, the magnitude of the change for p.Ala505Ser (median =−0.7 kcal/mol) was much smaller than that observed for p.Phe368Cys (median =−1.83 kcal/mol). Thus, missense variant p.Ala505Ser induces a mild destabilization on the total protein tertiary structure of ZAK.

Patient has deletion at chromosome 7q21.3-q22.1

To extend this analysis, the patient with the p.Ala505Ser variant underwent additional laboratory testing by the GGC cytogenetics laboratory. Microarray analysis revealed a copy number loss of approximately 2.01 Mb on chromosome 7q21.3-q22.1 between positions 96,070,437 and 98,080,551. This deletion includes the SHFM1 region of chromosome 7q and involves a total of 27 genes, 7 of which have Online Mendelian Inheritance in Man (OMIM) entries (SHFM1, DLX6, DLX5, TAC1, ASNS, LMTK2, BHLHB8). Within this region, the genes SHFM1, DLX5, and DLX6 have been previously implicated as SHFM-causing genes.

Discussion

The gene of interest in this study, ZAK, has been shown, in two consanguineous Pakistani pedigrees and in confirmatory animal models, to be relevant to SHFM with associated hearing loss.[7] The 38 patients selected for study in this cohort were screened for phenotypic match to patients with the reported pathogenic variant p.Phe368Cys in ZAK. All 38 patients in this cohort manifested the SHFM phenotype and had variable degrees of sensorineural hearing loss. Next, because the reported pedigrees[7] exhibit autosomal recessive inheritance of the ZAK mutation, patients were included in this study only if their pedigree suggested possible autosomal recessive inheritance. In summary, the cohort of unrelated patients selected for this study was an appropriate group of patients in which to sequence ZAK based on their similarities to previously reported cases.

The patients in the two Pakistani pedigrees that prompted this study contained a pathogenic substitution p.Phe368Cys in the ZAK protein. This residue lies within a SAM domain, which was previously shown to be critical to the SHFM-related functions of ZAK.[7] Sanger automated DNA sequencing revealed a rare homozygous missense mutation, c.1515 G>T, a previously unreported variant in ZAK, for one patient in our cohort. This variant encodes substitution p.Ala505Ser, which is located within the ZAK protein’s hydrophobic region of unknown function. To assess the significance of this variant, we performed in silico protein modeling, which predicted that substitution p.Ala505Ser destabilizes the ZAK protein, as indicated by the negative predicted binding free energy change values shown in Table 2. Notably, the amplitude of the destabilization for this novel variant (median = −0.7 kcal/mol) was lower than that predicted for the known pathogenic variant p.Phe368Cys (median =−1.83 kcal/mol). In summary, in silico modeling suggested that p.Ala505Ser elicits some destabilization of the tertiary structure of the ZAK protein.

These results must be interpreted in the context of low evolutionary conservation at this residue in the ZAK protein (Figure 3e). Alanine at residue 505 is not conserved, and variation with polar (including arginine) or charged residues is observed in other species. More distant sequence homology analysis revealed that zebrafish naturally express serine at residue 505. Although zebrafish lack limbs, the presence of p.Ala505Ser in nature suggests that the protein destabilization caused by this variant may not abrogate protein function. However, this is countered by the fact that serine occurs at a rare allele frequency at this locus in humans.

To search for a more convincing etiology of this patient’s SHFM, additional diagnostic laboratory studies were performed. Microarray analysis detected a deletion of chromosome 7q21.3-q22.1, involving the SHFM1 region. The SHFM1 region is an approximately 1.2-Mb critical region at 7q21.3 that contains the SHFM1, DLX5, and DLX6 genes, all of which are relevant to limb development.[18,19] Deletion, translocation, or disruption of this locus is reported to cause ectrodactyly and, in 35% of these patients, associated sensorineural hearing loss.[20] Tackels-Horne et al. (2001) were the first to associate hearing loss with the 7q21.3-q22.1 locus in two pedigrees of patients containing a variety of deletions or translocations involving the SHFM1, DLX5, and DLX6 genes.[19] In addition to associating this region with sensorineural hearing loss, the authors hypothesized that DLX5 or DLX6 may be important to developmental processes. This was later demonstrated to be the case, as Robledo et al. (2002) showed that DLX5 and DLX6 were implicated in craniofacial, axial, and appendicular skeletal development in a mouse model.[21] The low penetrance of hearing loss (∼35%) in patients with deletion of these genes suggests that other gene products might have functions similar to those of DLX5 and DLX6. Thus, we question whether destabilization in ZAK might have contributed to the patient’s SHFM phenotype and associated hearing loss.

This question led us to study the reported pathways, which are summarized in Figure 1. The similarity of known biological functions for the ZAK and p63 proteins provide a potential mechanistic explanation for their similar effects on limb formation and hearing loss. As thoroughly reviewed by Kantaputra and Carlson (2018),[8] multiple etiologies of SHFM converge upon regulation of fibroblast growth factor 8 (Fgf8), which regulates genetic imprinting mechanisms that influence cell polarization and stratification within the AER.[22] No clear mechanistic connection between ZAK and p63 is present in the literature, but Spielmann et al. (2016) determined, using quantitative PCR, that RNA levels of p63 were 60% lower in ZAK knockout mice than in WT mice (Figure 1).[7] Ultimately, this leads to abnormal transcriptional regulation of DLX5 and DLX6 by the protein ΔNp63α, altering expression of DLX5 and DLX6.[2] The influence of these pathways upon hearing is less well characterized. In summary, these pathways suggest that deficiency in ZAK enzyme activity, alongside deletion of DLX5 and DLX6, may have an additive effect on an SHFM phenotype with hearing loss.

The observation that knockout of ZAK was embryonic lethal but deletion of the SAM domain recapitulated the SHFM phenotype implicates the SAM domain as the critical domain for function.[7] Because p.Ala505Ser is not predicted to significantly alter SAM domain structure or function, this missense variant is likely insufficient to independently cause the SHFM phenotype. Whether this variant contributed, by a polygenic mechanism, to worsen the phenotype for this patient or to cause penetrance of hearing loss—a finding only seen in 35% of patients with SHFM1[20]—remains to be further confirmed but is plausible given the overlap between ZAK signaling and DLX5/6 pathways. We report this analysis to potentially aid interpretation of variations in the ZAK gene in future patients.

Conclusions

The mild destabilization effect of p.Ala505Ser on ZAK protein stability suggests the potential for p.Ala505Ser to contribute an additive effect on the chromosome 7q21.3-q22.1 deletion, but that the primary pathogenic locus is the deleted 7q21.3-q22.1 region containing SHFM1, DLX5, and DLX6. These results must also be considered in context of limitations of our study, which include that the rare variant p.Ala505Ser was only observed in one patient; thus, these results must be considered as a single case report.

Click here for additional data file.

Supplemental material, IMR879293 Supplemental Material for Rare missense variant p.Ala505Ser in the ZAK protein observed in a patient with split-hand/foot malformation from a non-consanguineous pedigree by Christopher Ronald Funk, Elizabeth S. Huey, Melanie M. May, Yunhui Peng, Ekaterina Michonova, Robert G. Best, Charles E. Schwartz and Anna V. Blenda in Journal of International Medical Research