A Novel Missense Mutation in TNNI3K Causes Recessively Inherited Cardiac Conduction Disease in a Consanguineous Pakistani Family

Cardiac conduction disease (CCD), which causes altered electrical impulse propagation in the heart, is a life-threatening condition with high morbidity and mortality. It exhibits genetic and clinical heterogeneity with diverse pathomechanisms, but in most cases, it disrupts the synchronous activity of impulse-generating nodes and impulse-conduction underlying the normal heartbeat. In this study, we investigated a consanguineous Pakistani family comprised of four patients with CCD. We applied whole exome sequencing (WES) and co-segregation analysis, which identified a novel homozygous missense mutation (c.1531T>C;(p.Ser511Pro)) in the highly conserved kinase domain of the cardiac troponin I-interacting kinase (TNNI3K) encoding gene. The behaviors of mutant and native TNNI3K were compared by performing all-atom long-term molecular dynamics simulations, which revealed changes at the protein surface and in the hydrogen bond network. Furthermore, intra and intermolecular interaction analyses revealed that p.Ser511Pro causes structural variation in the ATP-binding pocket and the homodimer interface. These findings suggest p.Ser511Pro to be a pathogenic variant. Our study provides insights into how the variant perturbs the TNNI3K structure-function relationship, leading to a disease state. This is the first report of a recessive mutation in TNNI3K and the first mutation in this gene identified in the Pakistani population.


Introduction
Cardiac conduction disease (CCD) is a hypernym that encompasses numerous rare disorders in which the integrity of the heart conduction system is disrupted [1]. In CCD, impairment of cardiac regeneration [28]. To date, five single nucleotide variants have been reported in human TNNI3K. Among them, three missense, one splicing, and one nonsense mutation were reported in different patients from Germany, Canada, Salvador, and China [16,[25][26][27]. Here, we report a novel recessive missense variant in TNNI3K (NM_015978.2: c.1531T>C) in a Pakistani family, which was not present in the Human Gene Mutation Database (HGMD) or any other publicly accessible databases.

Human Subjects, DNA Extraction, and Clinical Evaluation
A four-generation consanguineous family with conduction heart disease was recruited from Swat, Khyber Pukhtunkhwa, Pakistan (Figure 1a). Written informed consent was obtained from all adult individuals and the parents of patients. Blood samples from two affected and four unaffected individuals were obtained in ethylenediaminetetraacetic acidcontaining BD vacutainers ® (BD-Plymouth, Plymouth, UK). DNA from blood samples was isolated using the Qiamp DNA FFPE Tissue kit (Qiagen, Hilden, Germany). Phenotypic assessment and diagnoses were conducted in two living symptomatic family members based on cardiological investigations and clinical data. Rhythm and conduction abnormalities were detected by a standard 24 h Holter electrocardiogram (ECG) and 12-lead surface ECG (speed 25 mm/s, 10 mm/mV). Echocardiography was performed to check for problems with the valves or chambers of the heart.

Sanger Sequencing Analysis
Primer pairs were designed using Primer 3 software. PCR products were directly sequenced using a commercial service (Microsynth AG, Sankt Gallen, Switzerland). Segregation analysis of DNA from all available family members was performed by bidirectional Sanger sequencing to validate mutations and sequences were analyzed using DNASTAR Lasergene12 (DNASTAR Inc., Madison, WI, USA). Primer sequences are listed in Table S1.

Molecular Dynamic Simulations
For the molecular dynamic (MD) simulation, the X-ray structure of wild-type (WT) TNNI3K was retrieved from the Protein Data Bank (PDB; http://www.rcsb.org/; download on 15 August 2020; PDB-ID: 4YFI) and converted into its biological oligomerization state [29,41]. The mutant structure of TNNI3K (containing p.Ser511Pro) was built based on the WT X-ray structure, with an introduction of the single-point mutation. The WT and mutant structures were preprocessed using the Protein Preparation Wizard (SCHRÖDINGER ® 2020-3) [42,43]. Bond orders were assigned, all hydrogen atoms were added to each structure, and termini were capped. Incomplete side chains and loops were calculated and water molecules deleted. Heterogenetic states were generated at pH 7 ± 2. Hydrogen bonding networks were optimized by minimization of sampled hydrogens [44]. Each structure was minimized to a gradient of 0.05 using an OPLS3e force field. Finally, the minimized structures were subjected to solvent-explicit all-atom MD simulations using the GPU-accelerated SCHRÖDINGER ® DESMOND software [45,46]. Each structure was integrated in a minimized periodic simple point charged water model with an orthorhombic box size of 10.0 Å × 10.0 Å × 10.0 Å and 0.15 M NaCl ions. Each solvated system was minimized and relaxed using the default protocol of SCHRÖDINGER ® Desmond. MD simulations were carried out with periodic boundary conditions in the NPT ensemble using OPLS3e force field parameters. Temperature and pressure were maintained at 300 K and 1 atmospheric pressure, respectively, using Nose-Hoover temperature coupling and isotropic scaling. The operation was followed by running a 1 µs NPT production simulation for each structure, saved at the 1 ns interval. Three independent MD simulations were conducted for each case.
Next, the WT and mutated structures were compared based on root mean square deviation (RMSD) and root mean square fluctuation (RMSF). Residue-residue contacts between WT and mutant structures were monitored over the course of the last 200 ns of the simulation time using the trajectory_asl_monitor.py script supported by SCHRÖDINGER ® . A contact was defined as a distance of < 5 Å between the centers of masses of defined contact residue pairs. For each residue contact pair, the fraction of the simulation time in which the residues were in contact was calculated as a mean of triplicate measurements. Interaction energy values between two monomers for the WT and mutant structures were computed using the Molecular Mechanics (MM)/Generalized Born Surface Area (GBSA) method implemented in SCHRÖDINGER ® . The average binding free energy (∆G_bind), based on MM/GBSA, was calculated using the SCHRÖDINGER ® Desmond commandline thermal_mmgbsa.py script. During the MM/GBSA calculation, the last 200 ns MD simulation trajectory (every 10th snapshot) was used as an input to compute the average binding free energy. All heatmaps and plots were generated using GraphPad Prism version 9.0.0 for Windows (GraphPad Software, San Diego, CA USA, www.graphpad.com). Figures were created using SCHRÖDINGER ® Meastro 2020-3 [47].

Clinical Features
Early onset cardiac conduction disease (CCD) was clinically diagnosed in four members of a Pakistani family (Figure 1a). The proband (IV-3) was a 22-year-old female with a cardiac conduction defect and no underlying structural heart anomalies. A 24 h Holter recording revealed mild sinus bradycardia with an average heart rate of 55 beats per minute (bpm) at rest (range: 24-110) and sinus tachycardia of 160 bpm during exercise. Standard 12-lead ECG recorded at rest showed changes in the right bundle branch block and left anterior fascicular block (RBBB/LAFB) with a bifascicular block. In individual IV-5, the 24 h Holter monitor recorded bradycardia with an average heart rate of 56 bpm and tachycardia of 160 bpm; however, the echocardiogram revealed a restrictive perimembranous ventricular septal defect with bilateral shunt. This patient died during this study due to sudden cardiac arrest (SCA). In addition, two other siblings (IV 1 and IV 2) who had a history of CCD died due to SCA at 26 and 23 years-of-age, respectively. Both parents and healthy siblings were assessed for CCD but no relevant clinical features were detected. Clinical data are presented in Table 1.

Exome Sequencing Reveals a Pathogenic TNNI3K Variant
To uncover the underlying genetic cause of CCD in the family, we performed exome sequencing in all available individuals, including two affected individuals, their parents, and two healthy siblings (Figure 1a). We assumed an autosomal-recessive mode of inheritance because of the consanguinity of the healthy parents. After applying filtering criteria as described in Section 2.2, a total of 81 rare genetic variants, including 79 nonsynonymous and two synonymous variants, were identified ( Figure S1). The identified variants spanned 69 protein coding genes and seven non-coding genes. Among the variants, eight co-segregated with the phenotype in the pedigree (Table S2). Candidate variants were further screened by filtering within regions containing runs of homozygosity (ROH) ( Table S3). Of the three resulting candidate variants, the TTLL7 (NM_001350215.1) missense variant is predicted to lead a serine to alanine change in a moderately conserved amino acid, which lies outside of protein domains with known functions; TTLL7 is expressed at high levels in the mammalian nervous system. The MET (NM_000245.3, NM_001324402.1) missense variant (rs545332056) changes an alanine to serine at position 319; this is a highly conserved amino acid but reported in ClinVar (CA4448049) as of uncertain significance. Interestingly, a homozygous missense mutation (Phe841Val) in MET was reported in a Pakistani family with hearing impairment; however, no such phenotype was detected in any affected individuals in the current study. The TNNI3K (NM_015978.2) serine residue at position 511 is highly conserved in TNNI3K orthologues; accordingly, p.Ser511Pro is classified as "damaging" based on SIFT and Polyphen2 scores ( Figure S2, Table S3). It also had very high dbNSFP and phred-scaled CADD (v1.3) scores of 0.79 and 26.9, respectively, consistent with pathogenicity (Table S3). Moreover, mutations in TNNI3K have also previously been linked to the conduction heart disease spectrum (MIM 616117) (Table S4). Therefore, the variant in TNNI3K is the most likely causative candidate for the phenotype under investigation. The identified variant is an extremely rare homozygous variant that has not been reported in a heterozygous/homozygous state in publicly available exome/genome databases. Furthermore, it was absent from 110 exomes from Pakistani families with various heart diseases and in exome sequence data from 50,000 UK Biobank participants.

p.Ser511Pro Influences TNNI3K Structural Fluctuation
To understand the effect of p.Ser511Pro on the structural and dynamic properties of the protein, we conducted three independent 1 µs MD stimulations on the wild type and mutant structures. The X-ray structure of TNNI3K with the lowest resolution (PDB-ID: 4YFI [29]), containing amino acids Glu441 to Ile726 for monomer I (chain A) and amino acids Gly451 to Met721 for monomer II (chain B), was used for the analysis of the wild type model (TNNI3K-WT; Figure 1d). The mutant model (TNNI3K-S511P) was obtained by introducing a point mutation (p.Ser511Pro) into TNNI3K-WT. The resulting RMSD Cα values of the MD simulations indicated that both models reached an equilibrium state within a simulation time of 1 µs ( Figure S3a). The distinct conformational states sampled for both models showed no major differences ( Figure S3b). Structural fluctuation of TNNI3K-P511S was comparable to that of TNNI3K-WT over the entire lengths of the models (average RMSF Cα : 1.0 Å (TNNI3K-WT) and 1.1 Å (TNNI3K-S511P); Figure 2a). The only pronounced differences were seen in the ATP-binding pocket of monomer I (Leu495-Leu496) and in the dimer interface of monomer II (Thr637, Met677, and His680-Ile682). All other regions exhibited changes smaller than the standard deviation or changes below 0.5 Å.

p.Ser511Pro Affects the Local Hydrogen Bonding Network and Protein Surface Hydrophobicity of TNNI3K
Next, we investigated the local effect of p.Ser511Pro, which is located in helix 2 (Lys500-Asn517). Ser511 formed one hydrogen bond from its backbone nitrogen to the backbone oxygen of Cys507 and a second hydrogen bond from its side chain hydroxyl oxygen to the backbone oxygen of Cys507 throughout the simulation time (Figure 2b). Interestingly, p.Ser511Pro cannot form a hydrogen bond to Cys507 because of its specific structure. Furthermore, the protein surface of TNNI3K-S511P was more hydrophobic than that of TNNI3K-WT (Figure 2c,d). Genes 2021, 12, x FOR PEER REVIEW 9 of 16

p.Ser511Pro Affects the Structure of Helix 2 and the ATP-Binding Pocket
To analyze the effects of p.Ser511Pro on each monomer structure, a comparative time-dependent intra-contact analysis was performed for TNNI3K-WT and TNNI3K-S511P between the amino acids of helix 2 and the amino acids of the corresponding monomer. As illustrated in Figure 3a,b, the largest contact difference of 97% (TNNI3K-WT: 98%; TNNI3K-S511P: 1%) in monomer I was found for the contact between Cys507 and Ile512, which are neighboring amino acids of S511 and S511P, respectively. For monomer II, the largest contact difference (91%) was also between Cys507 and Ile512 (TNNI3K-WT: 97%; TNNI3K-S511P: 6%). In monomer I, the distances between the Cα atoms of Ile512 and Cys507 were 6.02 ±0.18 Å for TNNI3K-WT and 7.05 ±0.40 Å for TNNI3K-S511P, respectively. Similar observations were made for monomer II (TNNI3K-WT: 6.03 ±0.19 Å; TNNI3K-S511P: 6.99 ±0.28 Å). These data imply that p.Ser511Pro causes minimal stretching of helix 2 in the mutant region (Cys507-Ile512) relative to TNNI3K-WT. Intriguingly, > ±50% differences were further observed only for contacts between helix 2 amino acids, namely Asp502, Val503, Phe506 and Cys514, in monomer I and amino acids of the ATPbinding pocket (Phe448, Leu452, Tyr497, and Phe535; Figure 3a,b). of helix 2 in monomer I and the residues of the corresponding monomer, and (b) residues of helix 2 in monomer II and the residues of the corresponding monomer. Each point represents the mean of the average residue-residue contact difference between TNNI3K-WT and TNNI3K-S511P over the last 200 ns from the three independent MD simulations (spectrum range from green to white to blue, from 100% to −100% contact; green means more contacts in TNNI3K-WT and blue means more contacts in TNNI3K-S511P). (c) Contact map between residue-residue pairs in the dimer interface for TNNI3K-WT vs. TNNI3K-S511P. The contact spectrum ranging from green to white to blue, representing 100% to −100% contact, where green indicates more contacts in TNNI3K-WT and blue indicates more contacts in TNNI3K-S511P. Each point represents the mean difference in the residue-residue contact over the last 200 ns from the three independent MD simulations. (d) MM/GBSA analysis applied for dimerization of TNNI3K. Every 10th snapshot of the last 200 ns of the MD simulation was used. Each violin plot represents the mean from the three independent simulations (95% confidence interval).

p.Ser511Pro Creates a Less Favorable Binding Site at the Dimer Interface
In addition, amino acids at the dimer interface (Ile584-Asn590 and Glu609-Leu613; Figure S4) were located within a 10 Å radius around helix 2. Therefore, differences of ≤±30% were measured in intra contact analysis. The largest difference (−30%) was for the contact between Ile512 and Phe612 in monomer I (TNNI3K-WT: 70%; TNNI3K-S511P: 100%; Figure 3a,b). The Cα atom distance between Ile512 and Phe612 in monomer I was 10.29 ± 1.48 Å for TNNI3K-WT and 7.67 ± 0.68 Å for TNNI3K-S511P. Interestingly, Ile512 and Phe612 in monomer II showed a difference of only 2% (TNNI3K-WT: 100%, TNNI3K-S511P: 98%) with a Cα distance of 10.01 ± 0.62 Å for TNNI3K-WT and 10.17 ± 0.75 Å for TNNI3K-S511P. These findings indicate that for monomer I, Phe612, located in the dimer interface loop (Glu609-Leu613), is brought closer to helix 2 by the p.Ser511Pro mutation. Therefore, we hypothesized that p.Ser511Pro might indirectly influence the dimerization of TNNI3K. To assess this, we conducted a comparative time-dependent inter-contact analysis between the dimer interface amino acids (monomer I/II: Ile584-Ans590, Glu609-Leu653, Leu668-His681, and Trp702-Arg709; Figure S4) of TNNI3K-WT and TNNI3K-S511P in monomer I, an amino acid close to Phe612, and Asn627 in monomer II (TNNI3K-WT: 89.50%; TNNI3K-S511P: 5.67%). By contrast, Glu609 in monomer II and Asn627 in monomer I had contact values of 39.00% in TNNI3K-WT and 31.33% in TNNI3K-S511P, representing a difference of 7.67%. This is not surprising as Phe612 and Glu512 only showed strong differences in monomer I and not in monomer II. Similar results were observed for the contact between Cys705 and Asp676 (Cys705 monomer I −Asp676 monomer II: 59.33%; Asp676 monomer I −Cys705 monomer II: 5%). By contrast, the contact between Arg587 monomer I and Thr622 monomer II showed similar effects to those of their respective counterparts (Arg587 monomer I −Thr622 monomer II: 60.66%; Thr622 monomer I −Arg587 monomer II: 42.33%). Conversely, for Arg641 monomer I and Gln639 monomer II, more contacts were measured in TNNI3K-S511P than in TNNI3K-WT (Arg641 monomer I. −Gln639 monomer II: −58.66%; Gln639 monomer I −Arg641 monomer II: −10.33%). For all other 17 contacts, the counterpart went in the opposite direction. Together, these results indicate that the mutation affects the dimer-binding interface. To better understand the observed effects of the mutations on the dimer interface, we also employed the MM/GBSA method to calculate the binding energy between the two monomers; this is a simple post-processing technique in which the free energy of a state is calculated from the internal energy (MM) of the molecule and its interaction with an implicit representation of the solvent (GBSA). Although this method does not estimate the unfavorable entropy contribution to the binding and hence is less accurate than, for example, free energy perturbation, thermodynamic integration, or Bennett acceptance ratio, it avoids the enormous costs associated with other techniques and is consequently used often for large systems [48][49][50]. We considered the MM/GBSA values determined here as scores indicating predictive power in the ranking rather than representing quantitative agreement with experimentally observed affinities. The calculated MM/GBSA values (Figure 3d) indicate that the observed structural change leads to a less favorable total free binding energy and thereby to a less favorable bound state.

Discussion
In the present study, we identified a recessive CCD-associated sequence alteration, c.1531T>C;(p.Ser511Pro), in TNNI3K in a consanguineous Pakistani family. This variant affects an amino acid that is completely evolutionarily conserved ( Figure S2) and which is predicted to be a causative mutation using in silico prediction tools. The missense mutation co-segregated with cardiac conduction disease in the family with complete penetrance (Figure 1a). It resides within the kinase domain of TNNI3K, which is highly conserved across species. This study aimed to gain insights into the effect of this mutation on the kinase domain and to assess its functional and structural consequences using long-term MD simulations.
The MD simulations provided detailed information about the structural fluctuation of the protein, the intra and intermolecular amino acid contacts, and the energetic properties of the mutant structure relative to the wild type. In the wild type model, Ser511 forms two hydrogen bonds with Cys507; one from its backbone nitrogen to the backbone oxygen of Cys507 and the other from its side-chain hydroxyl oxygen to the backbone nitrogen of Cys507. Pro511 cannot mimic the hydrogen bond network of Ser511 due to its unique structural property. This conformational change introduced by the mutation stretches helix 2 into the region containing the mutation (p.Cys507-Ile512) in both monomers.
Changes in both structural fluctuation ( Figure 2a) and amino acid contacts (Figure 3a,b) were also observed in the monomer I ATP-binding pocket adjacent to helix 2. As amino acid orientation and interactions in the ATP-binding pocket are critical for the activity of the complex [51], these findings suggest that the mutation alters the ATPase activity of TNNI3K and therefore causes a reduction in its catalytic phosphorylation activity. Interestingly, our hypothesis is in line with the findings of Podliesna et al. [26] who observed a reduction or loss of autophosphorylation activity in all previously reported human TNNI3K kinase domain mutants (p.Gly526Asp and p.Thr539Ala). They postulated that affected individuals would exhibit reduced TNNI3K protein levels and kinase activity, eventually leading to a cardiac excitation conduction defect [26].
Moreover, using MD, we observed that the mutation in monomer I brings Phe612, located within a loop (Ala605 to Glu634) in the dimer interface, closer to helix 2. Further intercontact analysis confirmed changes in the region of the dimer interface. Together with the results of MM/GBSA-based analysis, we hypothesize that the change in the amino acid interaction pattern at the dimer interface results in a more energetically unfavorable dimer arrangement in TNNI3K-S511P. Further analyses using a heterozygous dimer will be required to obtain a deeper understanding of the one-sided effects of the mutation in the monomer I.
The mutation-induced change in helix 2 also influenced the hydrophobicity of the TNNI3K protein surface. It is known that the surfaces of proteins in an aqueous environment are predominantly composed of hydrophilic amino acids [52]. Therefore, the substitution of the polar serine to a non-polar proline led to an increase in the hydrophobicity of the protein surface (Figure 2c,d). Changes in surface hydrophobicity are closely associated with protein misfolding and aggregation, which are pathological hallmarks of many human diseases such as Alzheimer's disease, Parkinson's disease, and heart diseases [16,53]. Interestingly, previous studies have reported that a heterozygous TNNI3K missense mutation leading to a Gly526Asp substitution causes the formation of insoluble protein aggregates in in vitro analysis and a hydrophobic spot in the kinase domain surface in in silico protein surface analysis. Additionally, a dominantly negative mutation in TNNI3K (Gly526Asp) has been shown to deplete TNNI3K in ventricular cardiomyocytes [16]. Here, we show that a bi-allelic p.Ser511Pro mutation resulted in a recessive conduction defect in a Pakistani family.
Hitherto, disease-causing TNNI3K variants have been reported to be due to heterozygous mutations (Table S4). Symptoms manifest with age, with disease onset in CCD patients ranging from 18 to 84 years [16,[25][26][27]. Interestingly, in the current study, the parents, who were obligate carriers, and heterozygous siblings were all asymptomatic; no subclinical evidence of cardiac conduction system dysfunction was observed. However, in homozygous patients, disease symptoms appeared at a very early age and patients had a shorter life expectancy than those in previous reports. Two patients (IV-I and IV-2) had already died at 26 and 23 years-of-age, respectively. A plausible explanation for these results may be autozygosity, which increases the homozygous occurrence of alleles associated with a dominant inheritance pattern. The genetic background and/or modifier loci could also have a role in increasing the severity of the disease [54][55][56].
To date, no studies have been conducted in the Pakistani population to investigate the genetics underlying the cardiac system, except for a report of a homozygous mutation in CACNA1D, which causes deafness with bradycardia syndrome [57]. Therefore, to the best of our knowledge, this is the first reported recessive TNNI3K mutation, as well as the first TNNI3K mutation of all detected in the Pakistani population.
Clinical implications: First, it is important to understand the underlying pathomechanisms of a disease; this is mainly facilitated by the identification of pathogenic mutations, especially in families, with a genetic background that is not frequently studied. Regardless, knowledge of a pathogenic mutation in a family helps in the early diagnosis of asymptomatic family members.

Conclusions
In summary, the p.Ser511Pro mutation identified in TNNI3K leads to a different hydrogen bond network in helix 2, which affects the contacts between amino acids of helix 2 and those of the ATP-binding pocket. It also affects the dimer interface, resulting in an unfavorable binding energy. As the p.Ser511Pro substitution results in greater hydrophobicity and loss of interaction with the ATP-binding pocket, Ser511 clearly plays an important role in orienting the TNNI3K dimer; however, these new findings require functional validation. Overall, our study contributes to the TNNI3K mutation spectrum and will inform genetic counseling for CCD in Pakistani families.