Mitochondrial complex 1 gene analysis in keratoconus

Ethics statement

All patients provided informed consent before participation in this study. The study was approved by Institutional Review Board (IRB00006862; All India Institute of Medical Sciences, New Delhi, India).

Clinical examination and selection of cases

A total of twenty keratoconus patients (Table 1) presented (during April 2009 to April 2010) at the Dr. R. P. Centre for Ophthalmic Sciences (AIIMS, New Delhi, India) were enrolled in this study. Clinical evaluation involved Ultrasonic Pachymetry, videokeratography (VKG), Orbscan, visual testing, fundoscopy, slitlamp-biomicroscopy, and retinoscopy. Of these patients, 14 were males and 6 were females. The mean age of presentation was 17.2 years. Diagnosis of keratoconus involved the presence of characteristic topographic features, such as inferior or central corneal steepening, or an asymmetric bowtie pattern with skewing of the radial axes, and the presence of one or more of the following characteristic, clinical features in one or both eyes: conical corneal deformation, munsen sign, corneal stromal thinning, a Fleischer ring or Vogt striae. Family history up to three generations was collected and pedigrees were drawn. All 20 cases were sporadic without any family history. All keratoconus cases secondary to causes like trauma, surgery, Ehlers Danlos syndrome, osteogenesis imperfecta, and pellucid marginal degeneration were excluded from the study. Twenty ethnically matched normal individuals without any ocular disorder were enrolled as controls. Health information was obtained from controls through the questionnaire; all underwent ophthalmological examination. Five milliliters of blood was collected by venipuncture in EDTA (EDTA) vaccutainers (Greiner Bio-One GmbH, Frickenhausen, Germany) from both patients and controls. DNA was extracted from whole blood samples by the inorganic method. For the population study, controls were taken from published data defining the lineage of the Indian population [14,15].

Polymerase chain reaction (PCR) amplification and sequence analysis of the mitochondrial DNA coding region

The mitochondrial complex 1 (ND1, ND2, ND3, ND4, ND4L, ND5, and ND6 [ND stands for NADH dehydrogenase]) was amplified in all patients and controls using 12 pairs of primers using cycling conditions as described by Kumar and associates [16] and presented in Table 2. Briefly, PCR amplifications for all primer sets were performed in a 40-μl volume containing 1.0 μl of 20 μM stock solution for each primer (Eurofins Genomics India pvt Ltd, Bangalore, India), 100 ng of genomic DNA, 1 unit of Taq polymerase (Banglore Genei, Bengaluru, Karnataka, India), 0.1 mM of each deoxynucleotide triphosphate (dNTP), and 4 μl of 10× PCR buffer (with 15 mM MgCl₂) by means of 30 cycles of amplification, each consisting of 30 s denaturation at 94 °C, 30 s annealing at 55 °C, and 1 min extension at 72 °C. Finally, an extension for 5 min at 72 °C was performed. Amplified PCR products were purified using a gel/PCR DNA fragments extraction kit (catalog number DF100; Geneaid Biotech Ltd., Sijhih City, Taiwan). Purified PCR products of both primers (forward and reverse) were sent for sequencing to MCLAB (Molecular Cloning Laboratories, South San Francisco, CA). All sequence variants from both KC patients and controls were compared to the Human Mitochondrial reference sequence NC_012920 provided by the National Center for Biotechnology Information (NCBI) using ClustalW2 (multiple sequence alignment program for DNA; European Molecular Biology Laboratory (EMBL)-European Bioinformatics Institute (EBI).

Computational assessment of missense mutations

For prediction of pathogenic characteristics of all non-synonymous mtDNA changes two homology based programs PolyPhen-2 (Polymorphism Phenotyping) and SIFT (Sorting Intolerant From Tolerant) were used.

PolyPhen structurally analyzes an amino acid polymorphism and predicts whether that amino acid change is likely to be deleterious to protein function [17-19]. Polyphen-2 is more advanced version of the earlier version PolyPhen [20]. The prediction is based on the position-specific independent counts (PSIC) score derived from multiple sequence alignments of observations in case of functional domain of protein and predicted hydrophobic and transmembrane (PHAT) matrix element difference in case of transmembrane region of protein. PolyPhen scores of above 0.85 indicate the polymorphism is probably damaging to protein function. Scores of above 0.15 are possibly damaging, and scores of less than 0.15 are classified as benign.

SIFT is a sequence homology-based tool that sorts intolerant from tolerant amino acid substitutions and predicts whether an amino acid substitution in a protein will have a phenotypic effect [21-23]. SIFT is based on the premise that protein evolution is correlated with protein function. Positions important for function should be conserved in an alignment of the protein family, whereas unimportant positions should appear diverse in an alignment. Positions with normalized probabilities less than 0.05 are predicted to be deleterious and, those greater than or equal to 0.05 are predicted to be tolerated.

Haplogroup and phylogenetic analysis

To check the fidelity of our conclusion, the evolutionary information and the significance of mutations should be known. For haplogroups (Hg) analysis we have carefully chosen two hundred healthy individual samples from same area for comparison analysis and these were also treated as controls. For all control samples, sequences of the control region were determined from position 16024 to 00300, using the ABI Prism Dye Terminator cycle-sequencing protocols developed by Applied Biosystems (Perkin-Elmer, Foster City, CA), to provide an initial presumed Hg assignment and cases were haplogrouped by complete coding region sequences. The C-track length variation at positions 16182 and 16183 in HVS-I and the indels at positions 00309 and 00315 in HVS-II were excluded from further analyses. Hg assignment was then confirmed, based on control and coding region Hg defining polymorphisms determined by means of direct sequencing.

The NETWORK 4.5.1.6 program was run for placing all the mutations of control samples in their respective phylogenetic tree using the protocol as described at the Fluxus Engineering Website. The matrilineal lineages of the case were drawn manually in the reduced median network of control samples, to create the topology map we have applied the reduced median algorithm (r=1), followed by the median-joining algorithm (epsilon=2).

Principal component analysis

To minimize errors both strands were double-sequenced. Principal component analysis (PCA) of mtDNA was performed. For this experiment we compiled our data of mitochondrial haplogroups. Controls samples of five different populations were taken from published data defining Indian lineage for comparative analysis representing each sector of India i.e., Northern India, North-Western India, Western India, Eastern India, and Southern India.

The MVSP software package (Kovach WL, Services KC. MVSP - A multi-variate statistical package for Windows ver 3.13m. 2004) was used to identify the principal components (PCs) of mitochondrial variations that lead to form a haplogroup for every individual. To express the relative importance of top two eigenvectors in the resulting PCA plot, two axes were scaled by a factor equal to the square root of the corresponding eigen value. This experiment was repeated to confirm the outcomes.

Sequence variation in Complex I genes

DNA sequencing of Complex I genes revealed a total 84 nucleotide variation in patients (Table 3) and 29 variations in controls (Table 4). Of the 84 nucleotide variations in patients, 18 (21.42%) were non synonymous (Table 5), 52 (61.90%) were synonymous, 9 (10.71%) variations were in RNA genes and 3 (3.57%) were in non-coding region. Of 84 nucleotide variations in KC patients; five variations (3918G>A, 5348C>T, 12007G>A, 12372G>A and 12561G>A) were also present in controls. Out of 29 nucleotide variations found in controls 5 were non-synonymous. Maximum nucleotide variations were in ND5 (n=28) followed by ND4 (n=15); ND2 (n=13); ND3 (n=7); ND1 (n=3); ND4L (n=3), tRNAs (n=9), rRNAs (n=2) and 3 in non-coding regions.

Of 84 variations, 2 variations were frame-shift (11273G>A, 5300T>T). In one patient (KC 16) a single base deletion of guanine was observed at mtDNA position 11273. This caused a frame shift mutation after codon 172 (Gly>Ala) and introduced a stop codon at position 174 which resulted in a 173 amino acids truncated protein.This variation was homoplasmic (Figure 1).

In patient KC 2, we found a 2 base pair (CA) deletion at genomic position 5300 and 5301. This frame shift mutation altered the amino acid reading frame in ND2 protein at position 277. This CA deletion produced a truncated protein of 287 amino acids (Figure 2).

In silico analysis

SIFT analysis revealed two pathogenic changes (p.L71I and p.T9A) and PolyPhen revealed two pathogenic changes (p.N30L and p.I172V). The polyphen score of p.T9A was not available (no result for this mutation was available through PolyPhen; Table 5).

Principle component analysis

The tight cluster in Principal component Analysis (PCA) plot comprises the north-western, western, and north Indian population whereas the southern Indian and eastern Indian population is caught in a loose cluster (Figure 3). The controls were taken from published data defining Indian lineage for PCA and Haplogroup Network. We have treated the patients as a sub group of individuals having genetic structure different from normal Indians e.g., population. The inferences from PCA plot strongly supports our motive behind the planning of experiment, interestingly the patient population has not shown any relevant genetic affinity with other macropopulations of India.

Evolutionary insight of Mt complex I sequence variations

The genetic diversity in India is very complex. Several mutations from even control regions have been classified into the associative agent for various diseases [34]. The degree of haplotype sharing between populations is to investigate the combined frequency of the shared haplotypes in two population groups. Thus, among the northern and the southern population groups the combined frequency of the haplotypes present also in the other group is significantly lower than that which we observed in the case of random groups. This is not surprising because West Eurasian-specific mtDNA haplogroups are rather frequent in northwest India [35]. Because the Indo-European and the Dravidic speakers of India are largely concentrated to the northern and southern parts of the subcontinent, respectively, the differences arising from geographic division of the Indian populations also correspond to different linguistic groupings [36]. In this study, we found that all the mutations were apparently North-Indian specific with some novel mutations. The sequencing of Complex 1 revealed 84 mutations, of which 14, including 2 frame shift mutations and 4 non-synonymous mutations, were novel and exclusively observed in KC patients. Interestingly, most of the patients and their maternal relatives were clustered under the haplogroups (T, C4a2a, R2’TJ, M21’Q1a, M12’G2a2a, M8’CZ, M7a2a, U5b1, U1a3) which are present as negligible frequency in normal Indian population, whereas only few patients were found to be a part of the haplogroups whose origin is contentious i.e., U7 (Indo-European), R2 and R31. We have found three patients who fall under Indian haplogroups (M4, M4’63, R31a1) but they also carry the same sets of novel synonymous and non-synonymous mutations i.e., 4769, 4985, 5580, and 12850 (Figure 4). We have found some novel mutations in addition to each individual’s lineages and they are different from each other. This finding suggests the positive/causative role of different combinations of the mitochondrial coding mutations in this disease, as the normal population, completely lack these mutations.

The patients harbored some novel mutations at the different sites in mitochondria i.e., 4769, 4985, 5580, and 12850. These variants have never been reported in any of the population studies whereas they were present in every patient. Nevertheless, it is impossible from the evolutionary point of view that these sets of mutations in the individuals from different haplogroup. By keeping in mind about mutation rates in the coding region and its natural selection [37], we propose that these variants could theoretically influence the patient’s phenotype. However, the variants present in coding regions of mitochondrial gene are not conserved in course of evolution. The patients were apparently homoplasmic (only one type of mtDNA was present). To determine whether the maternal inherited gene pool of keratoconus patients is truly closer to any major populations in India, we have constructed the haplogroup frequency based PCA plot for mtDNA (Figure 3). Indeed, this analysis shows ambiguously that the three Indian populations clusters tightly among themselves viz. North, North West and West populations and two populations are to be found in a loose cluster viz South Indian and East Indian, whereas the keratoconus population matches with none of them in the mtDNA PCA plots. However the genetic data indicates that the keratoconus patients comprise several different haplotypes, if they are compared to normal populations around them. Most of the patients are in the clades which are nonspecific to Indian lineages. This information suggests that keratoconus patients are among those who are recent migrants into India and some genes in mitochondria have acquired mutations which are not filtered by purifying selection. Our results explain that the patients are genetically unrelated to each other due to the present maternal lineages which were diversified in the history of evolution. This fact suggests that the polymorphisms which are playing pivotal roles in causing the disease are recently accumulated in the mitochondrial coding regions of an individual patient. We have found that mutations specifically found in KC patients can affect transcription, translation or have synergistic effect with other variants in causing the disease. It has been reported many times about synergistic effect of different mutations in mitochondria that can cause many severe diseases [38]. Nevertheless, it is impossible from the evolutionary point of view these sets of mutations to occur in the normal individual from different haplogroups. By keeping in mind about mutation rates in the coding region and its natural selection, we propose that these variants could theoretically influence the patient’s disease. Non-synonymous mutations and frame shift mutations adversely affect C1 synergetics resulting in increase ROS production and mitochondrial dysfunction. KC corneas are unable to process ROS due to depleted or low ATP levels and increased ROS production and thereby undergo oxidative damage. These corneas have increased levels of malondialdehyde (MDA), which can results in altered protein function leading to cascade of events, including apoptosis that can damage the corneal tissues.

Thus this pilot study highlights the role of sequence variation in mitochondrial complex I gene in keratoconus patients. Such cases with elevated free radicals levels and oxidative damage to cornea may benefit immensely by antioxidant therapy.

Conclusions

Keratoconus corneas are known to suffer from oxidative damage. It is important to analyze the cause of raised ROS levels. Hence in this study we analyzed mitochondrial complex I genes, as sequence variations in this gene complex are associated with depleted ATP and increased ROS levels. We found a higher number of non-synonymous sequence variations in KC patients in comparison to controls.