Next-generation sequencing of the whole mitochondrial genome identifies functionally deleterious mutations in patients with multiple sclerosis

Abstract

Multiple sclerosis (MS) is an immune-mediated disease of the central nervous system with genetics and environmental determinants. Studies focused on the neurogenetics of MS showed that mitochondrial DNA (mtDNA) mutations that can ultimately lead to mitochondrial dysfunction, alter brain energy metabolism and cause neurodegeneration. We analyzed the whole mitochondrial genome using next-generation sequencing (NGS) from 47 Saudi individuals, 23 patients with relapsing-remitting MS and 24 healthy controls to identify mtDNA disease-related mutations/variants. A large number of variants were detected in the D-loop and coding genes of mtDNA. While distinct unique variants were only present in patients or only occur in controls, a number of common variants were shared among the two groups. The prevalence of some common variants differed significantly between patients and controls, thus could be implicated in susceptibility to MS. Of the unique variants only present in the patients, 34 were missense mutations, located in different mtDNA-encoded genes. Seven of these mutations were not previously reported in MS, and predicted to be deleterious with considerable impacts on the functions and structures of encoded-proteins and may play a role in the pathogenesis of MS. These include two heteroplasmic mutations namely 10237T>C in MT-ND3 gene and 15884G>C in MT-CYB gene; and three homoplasmic mutations namely 9288A>G in MT-CO3 gene, 14484T>C in MT-ND6 gene, 15431G>A in MT-CYB gene, 8490T>C in MT-ATP8 gene and 5437C>T in MT-ND2 gene. Notably some patients harboured multiple mutations while other patients carried the same mutations. This study is the first to sequence the entire mitochondrial genome in MS patients in an Arab population. Our results expanded the mutational spectrum of mtDNA variants in MS and highlighted the efficiency of NGS in population-specific mtDNA variant discovery. Further investigations in a larger cohort are warranted to confirm the role of mtDNA MS.

About the Authors:

Ghada Al-Kafaji

Roles Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

ORCID logo https://orcid.org/0000-0002-8397-3689

Halla F. Bakheit

Roles Methodology

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

Faisal AlAli

Roles Methodology

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

Mina Fattah

Roles Methodology

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

Saad Alhajeri

Roles Methodology

Affiliation: Ministry of Health, Kuwait City, State of Kuwait

Maram A. Alharbi

Roles Data curation, Methodology

Affiliation: College of Forensic Sciences, Naif Arab University for Security Sciences, Riyadh, Kingdom of Saudi Arabia

Abdulqader Daif

Roles Data curation

Affiliation: King Saud University Medical City, Riyadh, Kingdom of Saudi Arabia

Manahel Mahmood Alsabbagh

Roles Methodology

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

ORCID logo https://orcid.org/0000-0002-8013-760X

Materah Salem Alwehaidah

Roles Formal analysis

Affiliation: Department of Medical Laboratory, Faculty of Allied Health, Kuwait University, Kuwait City, State of Kuwait

Moiz Bakhiet

Roles Investigation, Writing – original draft

Affiliation: Department of Molecular Medicine and Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders, College of Medicine and Medical Sciences, Arabian Gulf University, Manama, Kingdom of Bahrain

Introduction

Multiple sclerosis (MS) is a progressive inflammatory, immune-mediated disease of the central nervous system characterized by demyelination, axonal loss and neurodegeneration [1]. MS is the most disabling neurological condition of young adults [2] and more commonly occurs in females than in males [3]. There are over two million people affected with MS worldwide [4] and a significant increase in the number of MS cases in the Arab Gulf region [5, 6]. MS can be categorized based on the phase and severity of disease progression into: clinically isolated syndrome (CIS), relapsing-remitting MS (RRMS), secondary progressive MS (SPMS), and primary progressive MS (PPMS). Approximately 85% of MS patients present with RRMS, which is marked by alternate periods of acute demyelination (relapses) and periods of neurological recovery and stability (remissions). After a period of 15–20 years, most patients with RRMS progress into a SPMS with worsened neurological function and signified by few or no acute relapses [7].

Etiologically, combined contributions of genetics and environmental factors are considered as necessary influences for MS risk and development [8]. In recent years, a considerable evidence placed mitochondrial dysfunction at the centre of the pathogenicity of MS in which mitochondrial DNA (mtDNA) defects, and mitochondrial structural and functional changes contribute to the disease progression [9, 10]. Additionally, some patients with Leber Hereditary Optic Neuropathy, a mitochondrial inherited disorder caused by mtDNA mutations present with MS symptoms including inflammatory demyelination [11, 12], suggesting a major role of mtDNA mutations in the predisposition to MS.

The mitochondria are important organelles present within almost all eukaryotic cells. The primary function of mitochondria is generation of energy in the form of adenosine triphosphate (ATP) through the process of oxidative phosphorylation (OXPHOS). Reactive oxygen species (ROS) are generated by mitochondria as by-products of respiration and OXPHOS. The mitochondria are unique in that they have their own genome (mtDNA), which is solely maternally inherited [13]. Human mtDNA is a double-stranded circular molecule of about 16,500 nucleotides and contains 37 genes coding for 22 tRNAs and 2 rRNAs (12S, 16S) which are essential for protein synthesis within mitochondria, and 13 protein subunits of the electron transport chain (ETC) complexes, all of which are involved in the OXPHOS. The 13 mtDNA gene-encoded proteins include seven subunits of NADH dehydrogenase of complex I (ND1, ND2, ND3, ND4, ND4L, ND5 and ND6); cytochrome b-c1 (CYB) of ubiquinol-cytochrome c reductase of complex III; three subunits of cytochrome c oxidase of complex IV (COX I, COX II and COX III); and two subunits of mitochondrial ATP synthase of complex V (ATP6 and ATP8). In addition, the mitochondrial control region (or D-loop) is a non-coding region and contains essential transcription and replication elements responsible for the expression of mitochondrial genome. In human, the mtDNA is present in multiple copies per cell, which gives rise to an important feature of mitochondrial genetics; homoplasmy (when all copies of mtDNA are identical) and heteroplasmy (when there is a mixture of normal and mutated mtDNA copies) [14].

It is well documented that the mtDNA is more susceptible to oxidative damage and acquires higher rates of mutation than the nuclear DNA (nDNA), possibly due to absence of protective histone, inadequate DNA repair pathways and high ROS exposure [15]. Mutations in the mtDNA are known to cause diverse mitochondrial disorders, and are also linked to other complex traits such as ageing, cancer and neurodegenerative diseases, all are associated with defects in oxidative energy metabolism [16, 17]. Specifically, impaired mitochondrial energy metabolism is considered as a key pathogenic factor in MS and plays an important role in axonal degeneration [9, 10, 18]. From this perspective, neurogenetic studies in relation to mtDNA provided an evidence of the implication of mtDNA mutations/variants in the pathogenesis or risk of MS. Specifically, variants in mtDNA-encoded complex I have been linked to MS in different ethnic populations [19, 20]. In our previous study in Saudi patients with MS, we identified four novel mutations in the mtDNA-encoded ND4 gene, which reduced stability of complex I [21]. However, most of the studies on MS have used targeted sequencing for specific mtDNA genes/ gene regions of interest [19–21] or selected mitochondrial polymorphisms for disease association [22–24]. While analysis of whole mtDNA genome highlighted specific haplogroup-defining variants in Caucasian patients with MS, but did not examine variants across the mitochondrial genome [25].

In this study, we aimed to analyze the whole mitochondrial genome using next-generation sequencing (NGS) from blood samples of Saudi subjects with RRMS and healthy controls to identify mtDNA disease-related mutations/variants. We hypothesized that mtDNA mutations/variants would occur at higher rate in cases compared to controls and implicated in susceptibility to MS. We also hypothesized that mutations in mtDNA-encoded genes occur in cases to a deleterious level and play a role in the pathogenesis of MS. NGS data can yield a large degree of population-specific mtDNA variants and permit extensive interpretation of potential disease-causing mutations.

Methods

Subjects

A cohort of 47 Saudi individuals consisting 23 patients with RRMS and 24 healthy controls were included in this study. Written consent forms were obtained from each participant. The contents of the consent forms and the study were reviewed and approved by the Scientific and Ethics Committee in King Saud University, College of Medicine, Kingdom of Saudi Arabia (KSA), and the Medical Research and Ethics Committee in the College of Medicine and Medical Sciences, Arabian Gulf University, Kingdom of Bahrain. The RRMS patients were recruited from the Neurology Outpatient Clinic at King Khalid Hospital, King Saud University, KSA. All patients were assessed independently by a neurologist and the diagnosis of RRMS was done according to the McDonald criteria [26] with at least two previous relapses in CNS regions, confirmed by a neurological examination, a magnetic resonance imaging (MRI), and electrophysiological studies. The healthy control individuals who had no neurological conditions or history of autoimmune and inflammatory disease were enrolled in this study from King Khalid hospital Blood Bank, KSA. Demographic and clinical data including age, gender, body mass index (BMI) and blood pressure were recorded for both patients and controls. Disease duration, disability status, and medications were recorded for all patients. The disability status was evaluated using the Kurtzke Expanded Disability Status Scale (EDSS).

Extraction of genomic DNA

Peripheral blood samples (5 ml) from patients and healthy controls were collected in ethylenediaminetetraacetic acid (EDTA) tubes. Genomic DNA was extracted from the peripheral blood using the QIAMP DSP DNA kit (Qiagen, Hilden, Germany) as previously described [21]. Briefly, 200 μl blood was mixed with 20 μl protease and 200 μl of lysis buffer. The mixture was incubated at 56°C for 10 min, and centrifuged at 20,000 x g for 1 min at 4°C. Then 200 μl of absolute ethanol was added and centrifuged at 6000 ×g for 1 min. This was followed by washing steps using 500 μl of washing buffer and centrifugation at 6,000 x g for 1 min and then at 20,000 x g for 3 min at room temperature. For genomic DNA elution, 200 μl of elution buffer was added, incubated at room temperature for 1 min and centrifuged at 6,000 x g for 1 min at room temperature. For quality control, all DNA samples were quantified and checked for purity using the NanoDrop ND-1000 ultraviolet-visible light spectrophotometer (Thermo Fisher Scientic, Inc.).

Amplicon generation

The mtDNA of all genomic DNA samples were amplified using long-range PCR (Qiagen Long range PCR kit, Cat# 206403) with two sets of primers: Primer Set 1, Forward L644–GACGGGCTCACATCACCCCATAA, and Reverse H8982-GCGTACGGCCAGGGCTATTGGT. Primer Set 2, Forward L8789-GCCACAACTAACCTCCTCGGACTCCT and Reverse H877- GGTGGCTGGCACGAAATTGACC. A 50 ng DNA was taken as input for the PCR amplifications. The amplified products were checked using 1% Agarose gel electrophoresis.

Next generation sequencing of whole mitochondrial genome

Library preparation and purity were performed according to the manufacturer’s protocol. Briefly, PCR products amplified from both primer sets were pooled together. The KAPA HTP Library Preparation Kit (Roche NimbleGen) was used to prepare libraries for whole genome sequencing. First, the DNA was fragmented to ~250 bp using Covaris M220 Focused-Ultrasonicator and the sheared DNA was checked on TapeStation using the HSD1000 DNA ScreenTapes (Agilent) for fragment distribution. Sheared DNA was then end repaired and adapter ligated in a series of enzymatic steps. The Adapter ligated products were then purified by Agencourt AmpureXP Beads (BeckmanCoulter) followed by size selection and then enriched using the following thermal conditions: Initial denaturation 98°C for 45 sec; 6 cycles of 98°C for 15 sec, 60°C for 30 sec,72°C 30 sec; final extension of 72°C for 1 min. PCR products were then purified by Agencourt AmpureXP Beads (BeckmanCoulter) and the final libraries were checked for fragment size distribution on TapeStation using the D1000 DNA ScreenTapes (Agilent). Prepared libraries were quantified using the Qubit HS Assay (Invitrogen). The obtained libraries were pooled and diluted to final optimal loading concentration before cluster amplification on Illumina flow cell. Once the cluster generation was completed, the cluster flow cell was loaded on Illumina HiSeq X instrument following the manufacturer’s instructions (Illumina) to generate 150 bp paired end reads.

Sequence data analysis

Paired end sequencing data were exported to FASTQ file. Sequence reads were trimmed using custom script to remove adapters and bases where the quality value was less than 20. The trimmed reads were aligned using Burrows–Wheeler Aligner (BWA) aligner. Alignment was performed to hg19 version of the genome available from UCSC genome browser. The aligned reads were filtered for low-quality mapping score, unusual insert-size and cross-chromosome mapping. Reads were re-aligned around variants and Base Quality Score Recalibration (BQSR) was done using GATK-lite program. Revised Cambridge Reference Sequence (rCRS) of the Human Mitochondrial DNA (NC_012920.1) was applied as reference mitochondrial sequence. Variants were annotated using MITOMAP and Human Mitochondrial Genome Database (mtDB). A total sequencing coverage of 10,000 X of the mitochondrial genome for this work was achieved. At such coverage heteroplasmy at 5% levels would be detectable. We determined whether the mutations were novel or known by investigating their presence in the MITOMAP, mtDB, ClinVar databases as well as Google search for individual studies reporting mtDNA mutations.

Functional and structural prediction analysis

The effect of mutations-caused amino acid changes on protein function was predicted by a combination of tools that use sequence homology, evolutionary conservation, and protein structural information [27]. These tools include:

1. Polymorphism Phenotyping (PolyPhen): Prediction based on protein sequence alignment and structure-based method, and uses a machine learning classification to categorize variants as benign, possibly damaging, and probably damaging with scores of < 0.446, > 0.446, and > 0.908 respectively.

2. Sorting Intolerant From Tolerant (SIFT): Uses sequence homology based on multiple sequence alignment conservation approach and how an amino acid substitution affects protein function. SIFT categories are benign and tolerated with scores of > 0.05 and < 0.05 respectively.

3. Combined Annotation Dependent Depletion (CADD): Integrates multiple information sources including conservation, structure-based features and functional information, and uses a machine learning approach to categorize variants as benign or deleterious. CADD scores range from 1–99 (variants with higher scores are more likely to be deleterious).

4. Mutation Assessor: Prediction based on evolutionary conservation of the affected amino acid in protein homologs. Mutation Assessor categories are neutral, low, medium and high with score ranging from 0–1 (variants with higher scores are more likely to be deleterious).

Moreover, the impact of single nucleotide substitutions-caused amino acid changes on three-dimensional (3-D) protein structure was determined using PyMOL (Molecular Graphics System Version 2.0, Schrodinger LLC). The 3-D protein structure was retrieved from https://www.ebi.ac.uk/pdbe.

Statistical analysis

The SPSS software (version 23; IBM Corp., Armonk, NY, USA) was used for data analysis. For differences in the demographic and clinical parameters between the patients and healthy controls, the Kolmogorov-Smirnov test was first used to evaluate the normal distribution of the data. Then the comparisons of variables between the two groups were done using a Chi-square test for categorical variables or equivalent non-parametric Wilcoxon signed-rank test and Mann-Whitney test for normally distributed variables. Fischer exact test was used to assess the differences in variant frequencies between patients and controls. Chi-Square test was used to determine the association of variants in patients and controls and the odds ratios (OD), relative risk (RR) and 95% confidence interval (95% CI) were reported. A two sided P value < 0.05 was considered to be statistically significant.

Results

Demographic and clinical data

Table 1 shows the demographic and clinical data of RRMS patients (n = 23) and healthy control individuals (n = 24). The mean age of patients was 28 years (standard deviation [SD], 7.5) and ranged from 18–44 years. The mean age of healthy controls was 31 years (SD, 7.5) and ranged from 22–52 years. The gender ratio (Male:Female; M:F) of patient was 5M:18F and the gender ratio of healthy controls was 5M:19F. There was no significant difference in the mean age (P = 0.14) or sex distribution between the two groups (P = 0.66). Also no significant differences were found between the two groups for mean BMI (P = 0.23) or mean systolic blood pressure (P = 0.16) and mean diastolic blood pressure (P = 0.37). The mean disease duration for the patient group was 5.3 years (SD, 4) and ranged from 1–15 years. The mean EDSS for the patient group was 3.9 (SD, 1.4) and ranged from 2–6.5. The patients were under the following treatment: of Avonex, Betaferon, Glienya, Rebif and Tysabri (n = 3, 8, 3, 5 and 4 respectively).

[Figure omitted. See PDF.]

Table 1. Demographic and clinical data of RRMS patients and healthy controls.

https://doi.org/10.1371/journal.pone.0263606.t001

Next generation sequencing of whole mitochondrial genome

We analyzed the whole mitochondrial genome from 47 subjects, including 23 RRMS patients and 24 healthy controls using NGS. All 47 mtDNA sequences passed quality controls. The overall alignment and the passed alignment percentage (alignment to hg19) in all samples were around 94.29% and 85.20% respectively. The average coverage of the samples on the panel was around 99.54% and the on-target percentage for the samples was around 86.72%. The data sets of mtDNA sequences generated in the present study were deposited at NCBI (http://www.ncbi.nlm.nih.gov/), Sequence Read Archive Repository (accession SRA data ref. no. PRJNA781092). The reference BioSample accession numbers are SAMN23235107 and SAMN23235153.

Distribution of variants in mitochondrial genome among patients and controls

Overall 3248 mtDNA variants were found across all samples. On average, 69 variants were detected per sample. Total number of variants detected in patients was higher as compared to controls. Of total, 1686 variants were found in patients and 1562 variants were found in controls. The variants were compared at the genomic position level for patients and controls, and were found to be distributed around all regions of mitochondrial genome in the two groups (Fig 1A). The difference in variant distribution across mitochondrial genome indicates that a large number of variants were located in the D-loop region, and the total number of coding and non-coding variants were higher in patients than in controls (Fig 1A). The majority of mtDNA variants across all samples were homoplasmic (85% in patients and 86% in controls) as compared with the heteroplasmic variants (15% in patients and 14% in controls) (Fig 1B).

[Figure omitted. See PDF.]

Fig 1.

A. Distribution of variants across mitochondrial genome in RRMS patients and healthy controls. A large number of variants were located in the D-loop region, and the total number of coding and non-coding variants were higher in patients than in controls. B. Percentage (%) of homoplasmic and heteroplasmic variants in RRMS patients and healthy controls. The majority of mtDNA variants across all samples were homoplasmic (85% in patients and 86% in controls) as compared with the heteroplasmic variants (15% in patients and 14% in controls).

https://doi.org/10.1371/journal.pone.0263606.g001

Unique and shared mtDNA variants in patients and controls

Comparison of mtDNA variants among patients and controls revealed a number of unique variants only observed in patients or only observed in controls, which occurred in one or multiple samples. Of the unique variants in patients, 177 were synonymous and 141 were nonsynonymous including 34 missense and 107 silent variants (Fig 2). Of the unique variants observed in controls, 152 were synonymous and 118 were nonsynonymous including 29 missense and 89 silent variants (Fig 2). Details of these variants are displayed in Table 2.

[Figure omitted. See PDF.]

Fig 2. Unique variants in RRMS patients and healthy controls.

The unique variant only found in patients include 177 synonymous and 141 nonsynonymous (34 missense and 107 silent). The unique variant only found in controls include 152 synonymous and 118 nonsynonymous (29 missense and 89 silent) variants.

https://doi.org/10.1371/journal.pone.0263606.g002

[Figure omitted. See PDF.]

Table 2. mtDNA variants observed only in RRMS patients or only in healthy controls.

https://doi.org/10.1371/journal.pone.0263606.t002

Moreover, a large number of variants were found to be shared among patients and controls (Table 3). The majority of the observed variants showed no significant differences in their prevalence between patients and controls (P>0.05). Whereas the prevalence of 8 variants differed significantly between the two groups (P<0.05). Of these, two synonymous variants namely 152T>C and 263A>G (referred as rs285351 polymorphism) were located in the D-loop region of the mtDNA. The variant 152T>C was detected in 23% of patients and 13% of controls (P = 0.024). Whereas the variant 263A>G was found in 100% of patients and 50% of controls (P<0.01). Other 6 variants were found in protein-coding genes: The variant 4216T>C represented a non-synonymous polymorphism referred as rs1599988 in MT-ND1 gene and was observed in 39% of patients compared to 13% of controls (P = 0.048). The variant 7028C>T represented a synonymous polymorphism referred as rs2015062 in MT-CO1 gene and was observed in 96% of patients compared to 42% of controls (P<0.01). The variant 10398A>G represented a non-synonymous polymorphism referred as rs2853826 in MT-ND3 gene and was found in 30% of patients and 63% of controls (P = 0.041). The variant 11719G>A represented a synonymous polymorphism in MT-ND4 gene and was found in 83% of patients and 42% of controls (P = 0.006). The variant 13708G>A represented a non-synonymous polymorphism referred as rs28359178 in MT-ND5 gene was detected in 39% of patients compared to 8% of controls (P = 0.017). Finally, the variant 14766C>T represented a non-synonymous polymorphisms referred as rs193302980 in MT-CYB gene and was observed in 96% of patients and 50% of controls (P = 0.003).

[Figure omitted. See PDF.]

Table 3. Common mtDNA variants and their prevalence among RRMS patients and healthy controls.

https://doi.org/10.1371/journal.pone.0263606.t003

Nonsynonymous mtDNA variants identified only in patients

We next focused on nonsynonymous mtDNA variants that present only in MS patients as a possible genetic predisposition to disease. A total of 34 nonsynonymous variants categorized as missense mutations were found only present in patients but not in controls, some were detected in one patient and others were found in more than one patient (Table 4). They were located within protein-coding genes, therefore potentially functional. Of these, the 12031G>A mutation in MT-ND4 gene was never reported in the MTOMAP or other databases and thus classified as a novel mutation. While 2 mutations, namely 3865A>G in MT-ND1 gene and 3944T>C in MT-ND1 gene were previously reported in MS. The remaining 31 mutations were previously described to be associated with diseases including mitochondrial disorders, neurological disorders and other diseases (Table 4).

[Figure omitted. See PDF.]

Table 4. Nonsynonymous mtDNA variants identified only in RRMS patients.

https://doi.org/10.1371/journal.pone.0263606.t004

Prediction of impacts of missense mutations on protein function

To analyse the impacts of amino acid substitutions upon missense mutations on protein function, a combination of tools were used including PolyPhen, SIFT, CADD, and Mutation Assessor (Dong et al., 2015). Accordingly, the missense mutations were defined as deleterious when predicted to be damaging, probably/possibly damaging by these methods. The analysis predicted 7 mutations to be deleterious, all of which were not previously reported in MS (Table 5). Of these, 4 mutations namely 10237T>C in MT-ND3 gene, 9288A>G in MT-CO3 gene, 14484T>C in MT-ND6 gene and 15431G>A in MT-CYB gene were all found in one patient (P2). Two mutations, namely 8490T>C in MT-ATP8 gene and 15884G>C in MT-CYB gene were found in another patient (P4). The same mutation 15884G>C in MT-CYB gene was also found in two other patients (P12 and P18). Whereas, the mutation 5437C>T in MT-ND2 gene was in two additional patients (P14 and P23).

[Figure omitted. See PDF.]

Table 5. Prediction of impacts of missense mutations on protein function.

https://doi.org/10.1371/journal.pone.0263606.t005

Prediction of impacts of missense mutations on protein structure

To analyze the impacts of deleterious missense mutations on protein structure, we have further studied the 3-D structure of encoded proteins for investigating the structural changes caused by specific mutation. The PyMOL (Molecular Graphics System Version 2.0, Schrodinger LLC) was used and the 3-D protein structure was retrieved from https://www.ebi.ac.uk/pdbe. The PyMOL was used to induce the deleterious mutations observed in mtDNA-encoded genes of the OXPHOS system including MT-ND2 and MT-ND6 of Complex I, MT-CO3 of Complex IV, and MT-CYB of Complex III of OXPHOS system on 3-D retrieved models: 5xtd, 5xtd, 5z62 and 5xte, respectively. Each mutation was examined at the level of polar bonds and clashes. However, the effect of mutations in MT-ATP8 gene of Complex V could not be illustrated at the protein structural level as the human 3-D structure are not available for Complex V proteins.

The mutation 10237T>C of MT-ND3 gene has change in amino acid at Ile60Thr (S1 Fig). The wildtype Ile60 (red) (A) forms three polar bonds (yellow dashed lines) with Phe56, Leu-63 and Leu-64 (orange). Whereas the mutated Thr60 (red) (B) forms five polar bonds (yellow dashed lines) with Phe56, Leu57, Thr61, Leu63 and Leu64 (orange), and causes physical clashes with the surrounding amino acids (C). The mutation 9288A>G of MT-CO3 gene has change in amino acid at Thr28Ala (S2 Fig). The wild type Thr28 (red) (A) forms four polar bonds (yellow dashed lines) with Ala24, Leu25 and Ser29 (all in orange). Whereas the mutated Ala28 (red) (B) shows a total loss of polar bonds, and causes physical clashes with the surrounding amino acids (C). The mutation 14484T>C of MT-ND6 gene has change in amino acid at Met64Val (S3 Fig). The wildtype Met64 (red) (A) does not form any polar bonds with the surrounding amino acids. Whereas the mutated Val64 (red) (B) forms a single polar bond with Gly-68 (orange) and causes physical clashes with the surrounding amino acids (red) (C). The mutation 15431G>A of MT-CYB gene has change in amino acid at Ala229Thr (S4 Fig). Both the wildtype Ala229 (red) (A) and the mutated Thr229 (red) (B) form two polar bonds (yellow dashed lines) with the surrounding amino acids (orange). However, the mutated Thr229 causes physical clashes with the surrounding amino acids (red) (C). For the mutation 15884G>C of MT-CYB that has change in amino acid at Ala380Pro, the 3-D structure could not be done. Although alanine is the last amino acid within the protein, it was not shown in the 3-D structure, whereas the last shown amino acid was tryptophan 379. The mutation 5437C>T of MT-ND2 gene has change in amino acid at Thr323Ile (S5 Fig). The wildtype Thr323 (red) (A) forms a single polar bond (yellow dashed line) with Pro322 (orange), while the mutated Ile323 (red) (B) shows a total loss of polar bonds and causes physical clashes with the surrounding amino acids (C).

Discussion

Mutations in the mtDNA and disruption of mitochondrial function alter brain energy metabolism and cause neurodegeneration [9, 10, 18]. Previous studies on MS have used specific primers for targeted regions of mtDNA [19–21] or investigated specific mitochondrial polymorphisms and haplogroups for disease association [22–25]. In this study, we analyzed the whole mitochondrial genome to examine mtDNA mutations/variants in blood samples from Saudi patients with RRMS and healthy controls using NGS. Overall we identified 3248 variants, 1686 in patients and 1562 in controls that were distributed around all regions of mitochondrial genome (Fig 1A). We observed a large number of variants located in the D-loop region with a higher number of variants in patients than in controls (Fig 1A). There was an excess of homoplasmic variants (85% in patients and 86% in controls) compared to heteroplasmic variants (15% in patients and 14% in controls) (Fig 1B). When we analyzed the variants only observed in patients or only observed in controls, we found distinct unique variants in patients or controls (Fig 2 and Table 2). There were 177 synonymous variants only present in patients versus 152 synonymous variants only present in the controls. While there were 141 nonsynonymous (34 missense and 107 silent) variants among patients versus 118 nonsynonymous (29 missense and 89 silent) variants only present in controls. We also found a large number of common synonymous and nonsynonymous variants that were shared among patients and controls (Table 3). Common genetic variants which frequently occur in the human mtDNA are expected to be under selective constraints [28]. Such variants have functional qualities and can disturb mitochondrial function since natural selection acts on phenotype [29]. Indeed, mtDNA variants can influence mitochondrial activities [30], mtDNA transcription [31] and generation of ROS [32]. The impact of mtDNA variants has been also implicated in aging and age-related diseases [33]. Because of the degeneracy of the genetic code in which more than one codon may specify the same amino acid, many mutations in coding sequences, especially at the third base of codons, do not affect protein sequence or structure and function of proteins, and are therefore considered silent. However, it has been shown that silent mutations created by synonymous DNA changes can have a significant impact on mRNA stability and folding, and consequently gene expression [34–37]. They also influence the rate of translation as well as posttranslational modification of proteins [31, 38].

Some of the common mtDNA variants in our study showed differences in their prevalence among patients and controls (P<0.05) (Table 3). In the D-loop region, the synonymous variants 152T>C and 263A>G occurred at significantly higher frequencies in patients (23% and 100% respectively) compared to controls (13% and 50% respectively). The mitochondrial non-coding D-loop is known to accumulate mutations at a higher rate than other regions of mtDNA. Variants in the D-loop region have important contribution to the proper functioning of mitochondria through an effect on mtDNA replication, transcription or translation [31]. While variants in the mtDNA D-loop have been described in diseases such as cancer [39], cardiovascular disease [40] and Huntington’s disease [41], both 152T>C and 263A>G variants found in this study were not previously reported in MS.

In protein-coding genes also differed significantly in their frequencies among patients and controls. Namely, 4216T>C in MT-ND1 gene (referred as rs1599988), 13708G>A in MT-ND5 gene (referred as rs28359178) and 14766C>T in MT-CYB gene (referred as rs193302980) are all missense mutations, observed more frequent in patients (39%, 39% and 96% respectively) than in controls (13%, 8% and 50% respectively). 4216T>C and 13708G>A were previously reported in MS patients and considered as diseases predisposing markers [42, 43]. Whereas 14766C>T variant was linked to 13708G>A variant, a potential marker associated with the increased risk of MS in European cohorts [43]. Other variants namely 7028C>T in MT-CO1 gene (referred as rs2015062) and 11719G>A in MT-ND4 gene are silent polymorphisms, found more frequently in patients (96 and 83% respectively) than in controls (40 and 42% respectively). Both 7028C>T and 11719G>A were reported in MS patients from different populations [21, 43]. Contrary to these results, 7028C>T was not found to be associated with MS in a Russian cohort [44]. The variability in the occurrence and prevalence of these variants in MS patients in different studies may be due to genetic and racial factors, which differ among populations. Nevertheless, we could not exclude the possibility that these variants may be possible risk factors for MS. In addition, 10398A>G (referred as rs2853826) in MT-ND3 gene represented a missense mutation and was found less frequently in patients (30%) than in controls (63%). This mutation was reported in previous studies to reduce the risk to Parkinson’s disease [45, 46], and could have a similar effect in MS.

Analysis of nonsynonymous variants detected only in patients but not in controls (Table 4) revealed 34 missense mutations. According to MITOMAP, mtDB and ClinVar as well as Google search, one mutation namely 12031G>A in MT-ND4 gene was assessed as a novel mutation since it was not found in the above mentioned data sets and was non haplogroup associated with Phylotree. Whereas 2 mutations namely 3865A>G in MT-ND1 gene and 3944T>C in MT-ND1 gene were previously reported in MS, and the 31 mutations were previously described in several diseases such as mitochondrial disorders and neurological diseases. Pathogenic variants in the mtDNA due to changes in the amino acid sequences of protein coding genes have been associated with susceptibility to complex diseases [47]. The substitution of amino acid may cause change in the hydrophobicity/ hydrophilicity and influence the charges or binding properties of a given protein, thus affects the protein structure and function [48]. The PolyPhen, SIFT, CADD, and Mutation Assessor which are widely used to predict the effects of nonsynonymous mutations on protein function [27], confirmed 7 deleterious missense mutations (Table 5). The mutations were identified in different mtDNA-encoded genes of the ETC complexes. Three of these mutations were found in NADH dehydrogenase subunits: The transition mutation 5437C>T in MT-ND2 is homoplasmic and causes changes in amino acid from Thr323Ile. The transition mutation 10237T>C in MT-ND3 gene is heteroplasmic and causes changes in amino acid from Ile60Thr. The transition mutation 14484T>C in MT-ND6 gene is homoplasmic and causes changes in amino acid from Met64Val. Threonine (Thr) is a neutral uncharged polar amino acids, whereas isoleucine (Ile), methionine (Met) and valine (Val) are hydrophobic uncharged nonpolar. NADH dehydrogenase (complex I) is the first enzyme of the respiratory chain, which catalyzes the transfer of electrons from NADH to ubiquinone. It a multi-subunit complex in which 7 subunits (ND1, ND2, ND3, ND4, ND4L, ND5 and ND6) are encoded by the mtDNA genes. There is a causative connection between defects in complex I and neurodegenerative diseases such as Parkinson’s disease [49] and MS [50], and variants in mtDNA-encoded complex I genes are implicated in the pathogenicity or risk of MS [19–21]. The identified deleterious mutations in NADH dehydrogenase subunits are likely to impair the function of the whole enzyme as deduced from their effects on complex I.

Another deleterious mutation was found in Cytochrome c oxidase, namely 9288A>G in MT-CO3 gene. This transition mutation is homoplasmic and causes changes in amino acid from Thr28Ala. Threonine (Thr) is neutral uncharged polar amino acids, whereas alanine (Ala) is hydrophobic uncharged nonpolar. Cytochrome c oxidase (complex IV) catalyzes the final step in mitochondrial ETC, functioning in the transfer of electrons from cytochrome c to oxygen that is then converted to water. It consists of multi-subunit essential for assembly and respiratory function of the enzyme complex [51]. Three of the biggest subunits (COX I, COX II and COX III) are encoded by the mtDNA and form the functional core of the enzyme complex. Mutations in mtDNA-encoded COX subunits can lead to complex IV deficiency by affecting complex structure and assembly [51]. Functional defects of mitochondrial COX have been reported in MS lesions [52] and our study is the first to report the harmful effect of 9288A>G mutation of MT-CO3 gene in MS patients.

Three other mutations were found in MT-CYB gene: The transition mutation 15431G>A is homoplasmic and causes changes in amino acid from Ala229Thr. The transition mutation 15674T>C is also homoplasmic and causes changes in amino acid from Ser310Pro. Whereas the transversion mutation 15884G>C is heteroplasmic and causes changes in amino acid from Ala380Pro. Alanine (Ala) is a hydrophobic uncharged nonpolar amino acid, threonine (Thr) and serine (Ser) are neutral uncharged polar amino acids, and proline (Pro) is a neutral uncharged nonpolar amino acid. The MT-CYB gene is a component of the ubiquinol-cytochrome c reductase (cytochrome b-c1 complex or complex III), the third complex in the ETC which plays a critical role in ATP generation. It contributes to the generation of a proton gradient across the mitochondrial membrane that is then used for ATP synthesis. The deleterious mutations identified in MT-CYB gene may have an important impact on complex III, particularly the heteroplasmic transversion mutation 15884G>C. The guanine base (G) in genomic DNA is highly susceptible to oxidative stress due to having the lowest oxidation potential. Transversion mutations such as G>C frequently occur under oxidative stress [53], an important factor in the progression of MS [54] and other neurodegenerative diseases [55].

Additionally the mutation 8490T>C was found in MT-ATP8 gene of mitochondrial ATP synthase of complex V. This transition mutation is homoplasmic and causes changes in amino acid from Met42Thr. Methionine (Met) is a hydrophobic uncharged nonpolar amino acid, whereas threonine (Thr) is a neutral uncharged amino acids. Mitochondrial ATP synthase composed of 17 subunits, of which 2 (ATP6 and ATP8) are encoded by the mtDNA. Complex V is responsible for last step of OXPHOS and synthesizes ATP from ADP in the mitochondrial matrix using the energy provided by the proton electrochemical gradient [56]. While the 8490T>C mutation in MT-ATP8 gene was not previously reported in MS, different pathogenic mutations in MT-ATP8 gene were described in human disorders such as epilepsy, MIDD, brain pseudo-atrophy, episodic weakness and neurological disorders [57]. Mutations in MT-ATP8 gene not only disturb the ATP synthase function but can also affect the assembly process and stability of complex V [57].

The structure of a protein is closely linked to its stability, function, and interaction. Mutations that predominantly impact the function or structure of proteins show their effects through several factors including impairment of the catalytic machinery, disruption of a functional site, modification of inter-subunit interfaces, disruption of protein folding and destabilization of the protein [58, 59]. In addition to the harmful effects of mtDNA missense mutations on protein function, these mutations disturbed the structure of encoded-proteins and cause physical clashes with the surrounding amino acids as indicated by PyMOL analysis (S1–S5 Figs). When the substituting residue does not fit in the protein structure as a result of clashes with the surrounding amino acids, a drastic conformation change occurs as the consequence, which leads to local or global structural alteration in the protein [60]. Moreover, clashes of amino acids can drastically alter the folding and stability of proteins [59]. Therefore, these mutations appear functional and structural impact on encoded-proteins.

In our study, the effect of the mutation 8490T>C (Met42Thr) in MT-ATP8 gene could not be illustrated at the protein structural level. However as mentioned above, replacement of amino acid from methionine to threonine due to this mutation can change the polarity of amino acids and hence affect the function and structure of Complex V.

In light of all this information, it can be postulated that the identified deleterious missense mutations in MS patients may be disease-related since they are located in structurally/functionally important sites and could impair the encoded-proteins involved in the OXPHOS system in mitochondria.

The mtDNA is particularly susceptible to acquiring somatic mutations due to its close proximity to the ROS production site, lack of protective histone and insufficient repair systems [15]. The accumulation of mtDNA mutations impairs the ETC, triggers oxidative damage, and subsequently a decline in mitochondrial function with more ROS production. Mitochondrial dysfunction as a result of mtDNA mutations can lead to to energy failure, a hallmark of neurodegenerative disease and contributes to axonal degeneration in MS [10, 18, 55]. Given the important role of the ETC complexes, particularly complexes I and complex III in ROS generation [61, 62], the detected deleterious mutations in these complexes may be implicated in MS by affecting the mitochondrial ROS production. However, it is difficult to identify ROS-causing or -pathogenic mutations because of the co-occurrence of multiple mtDNA mutations in several mitochondrial complexes.

When we further analyzed the patient-specific mutations (Table 5), we observed that some patients harboured multiple mutations, while other patients carried the same mutations. In this respect, patient P2 carried four mutations; the heteroplasmic mutation 10237T>C in MT-ND3 gene, and the homoplasmic mutations 9288A>G in MT-CO3 gene, 1484T>C in MT-ND6 gene and 15431G>A in MT-CYB gene. Patient P4 carried the homoplasmic mutation 8490T>C in MT-ATP8 gene and the heteroplasmic mutation 15884G>C in MT-CYB gene. The same heteroplasmic mutation 15884G>C in MT-CYB gene was also carried by two other patients (P12 and P18). Whereas two patients (P14 and P23) carried the homoplasmic mutation 5437C>T in MT-ND2 gene. Studying other family members of those patients, especially patient P2 who carried four harmful mutations which can lead to respiratory defects, is significant to further understand the role of mitochondrial neurogenetics in MS. Nevertheless, when mtDNA mutations are studied in a complex disease as MS, other modifying factors including environmental factors and the possible effects of the nuclear genes should be considered for the phenotypic manifestation of mutations.

Conclusions

Taken together, out of a large number of variants, we found distinct unique variants that were only present in the patients or only occur in the controls. We also found a number of common variants and the prevalence of some of them differed significantly between patients and controls. Additionally we identified 34 missense mutations only present in the patients. Of them, seven mutations (including homoplasmic and heteroplasmic) were not previously reported in MS, and predicted to be deleterious with considerable impacts on the functions and structures of encoded-proteins. We also observed that some patients harbored multiple deleterious mutations while other patients carried the same deleterious mutations. While the identified mtDNA sequence variants are not sufficient to cause a classical disease, they may be associated with a cascade involving altered energy output, and thus could contribute to the susceptibility or pathogenicity of MS. To the best of our knowledge, this is the first study to perform sequence analysis of the whole mtDNA of MS patients in an Arab population. Our results expanded the mutational spectrum of mtDNA variants in MS and highlighted the efficiency of NGS in population-specific mtDNA variant discovery. A more extensive analysis using a larger sample size is required to confirm the link between mitochondrial neurogenetics and MS.

Supporting information

S1 Fig. 10237T>C (Ile60Thr) mutation of MT-ND3 gene.

The wildtype Ile60 (red) (A) forms three polar bonds (yellow dashed lines) with Phe56, Leu-63 and Leu-64 (all in orange). Whereas the mutated Thr60 (red) (B) forms five polar bonds (yellow dashed lines) with Phe56, Leu57, Thr61, Leu63 and Leu64 (orange), and causes physical clashes with the surrounding amino acids (C).

10.1371/journal.pone.0263606.s001

(PDF)

S2 Fig. 9288A>G (Thr28Ala) mutation of MT-CO3 gene.

The wild type Thr28 (red) (A) forms four polar bonds (yellow dashed lines) with Ala24, Leu25 and Ser29 (all in orange). Whereas the mutated Ala28 (red) (B) shows a total loss of polar bonds, and causes physical clashes with the surrounding amino acids (C).

10.1371/journal.pone.0263606.s002

(PDF)

S3 Fig. 14484T>C (Met64Val) mutation of MT-ND6 gene.

The wildtype Met64 (red) (A) does not form any polar bonds with the surrounding amino acids. Whereas the mutated Val64 (red) (B) forms a single polar bond with Gly-68 (orange) and causes physical clashes with the surrounding amino acids (red) (C).

10.1371/journal.pone.0263606.s003

(PDF)

S4 Fig. 15431G>A (Ala229Thr) mutation of MT-CYB gene.

Both the wildtype Ala229 (red) (A) and the mutated Thr229 (red) (B) form two polar bonds (yellow dashed lines) with the surrounding amino acids (all in orange). However, the mutated Thr229 causes physical clashes with the surrounding amino acids (red) (C).

10.1371/journal.pone.0263606.s004

(PDF)

S5 Fig. 5437C>T (Thr323Ile) mutation of MT-ND2 gene.

The wildtype Thr323 (red) (A) forms a single polar bond (yellow dashed line) with Pro322 (orange), while the mutated Ile323 (red) (B) shows a total loss of polar bonds and causes physical clashes with the surrounding amino acids (C).

10.1371/journal.pone.0263606.s005

(PDF)

Acknowledgments

We would like to thank the technical support with sample transportation in the Kingdom of Saudi Arabia and the technical staff of Al-Jawhara Centre for Molecular Medicine, Genetics, and Inherited Disorders at the College of Medicine and Medical Sciences, Arabian Gulf University, Kingdom of Bahrain.

Citation: Al-Kafaji G, Bakheit HF, AlAli F, Fattah M, Alhajeri S, Alharbi MA, et al. (2022) Next-generation sequencing of the whole mitochondrial genome identifies functionally deleterious mutations in patients with multiple sclerosis. PLoS ONE 17(2): e0263606. https://doi.org/10.1371/journal.pone.0263606

References

1. Frohman EM, Racke MK, Raine CS. Multiple sclerosis—the plaque and its pathogenesis. N Engl J Med. 2006;354: 942–955. pmid:16510748

2. Reich DS, Lucchinetti CF, Calabresi PA. Multiple sclerosis. N Engl J Med. 2018;11; 378(2): 169–180. pmid:29320652

3. Harbo HF, Gold R, Tintore M. Sex and gender issues in multiple sclerosis. Ther Adv Neurol Disord. 2013;6(4): 237–248. pmid:23858327

4. Browne P, Chandraratna D, Angood C, Tremlett H, Baker C, Taylor BV, et al. Atlas of multiple sclerosis 2013: A growing global problem with widespread inequity. Neurology. 2014;83(11): 1022–4. pmid:25200713

5. Bohlega S, Inshasi J, Al Tahan AR, Madani AB, Qahtani H and Rieckmann P. Multiple sclerosis in the Arabian Gulf countries: A consensus statement. J Neurol. 2013;260: 2959–2963. pmid:23504049

6. AlJumah M, Bunyan R, Al Otaibi H, Al Towaijri G, Karim A, Al Malik Y, et al. Rising prevalence of multiple sclerosis in Saudi Arabia, a descriptive study. BMC Neurol. 2020;8;20(1): 49. pmid:32035478

7. Loma I, Heyman R. Multiple sclerosis: Pathogenesis and treatment. Curr Neuropharmacol. 2011;9: 409–416. pmid:22379455

8. Waubant E, Lucas R, Mowry E, Graves J, Olsson T, Alfredsson L, et al. Environmental and genetic risk factors for MS: an integrated review. Ann Clin Transl Neurol. 2019;6(9): 1905–1922. pmid:31392849

9. Mao P, Reddy PH. Is multiple sclerosis a mitochondrial disease? Biochim Biophys Acta 2010;1802: 66–79. pmid:19607913

10. Mahad D, Lassmann H, Turnbull D. Review: Mitochondria and disease progression in multiple sclerosis. Neuropathol Appl Neurobiol. 2008;34: 577–589. pmid:19076696

11. Harding AE, Riordan-Eva P and Govan GG: Mitochondrial DNA diseases: Genotype and phenotype in Leber’s hereditary optic neuropathy. Muscle Nerve Suppl. 1995;3: S82–S84. pmid:7603533

12. Bargiela D, Chinnery PF. Mitochondria in neuroinflammation—Multiple sclerosis (MS), leber hereditary optic neuropathy (LHON) and LHON-MS. Neurosci Lett. 2019;710: 132932. pmid:28668384

13. Garcia I, Jones E, Ramos M, Innis-Whitehouse W, Gilkerson R. The little big genome: The organization of mitochondrial DNA. Front Biosci. 2017;22: 710–721. pmid:27814641

14. Taylor RW, Turnbull DM. Mitochondrial DNA mutations in human disease. Nat Rev Genet. 2005;6: 389–402. pmid:15861210

15. Santos JH, Hunakova L, Chen Y, Bortner C, Van Houten B. Cell sorting experiments link persistent mitochondrial DNA damage with loss of mitochondrial membrane potential and apoptotic cell death. J Biol Chem. 2003;278: 1728–1734. pmid:12424245

16. Schon EA, DiMauro S, Hirano M. Human mitochondrial DNA: roles of inherited and somatic mutations. Nat Rev Genet. 2012;13(12): 878–90. pmid:23154810

17. Cha MY, Kim DK, Mook-Jung I. The role of mitochondrial DNA mutation on neurodegenerative diseases. Exp Mol Med. 2015;47: e150. pmid:25766619

18. Campbell GR, Worrall JT, Mahad DJ. The central role of mitochondria in axonal degeneration in multiple sclerosis. Mult Scler. 2014;20(14): 1806–13. pmid:25122475

19. Vyshkina T, Sylvester A, Sadiq S, Bonilla E, Canter JA, Perl A, et al. Association of common mitochondrial DNA variants with multiple sclerosis and systemic lupus erythematosus. Clin Immunol. 2008;129: 31–35. pmid:18708297

20. Zonouzi PA, Ghorbian S, Abkar M, Zonouzi PA and Azadi A: Mitochondrial complex I gene variations as a potential risk factor in the pathogenesis of multiple sclerosis. J Neurol Sci. 2014;345: 220–223. pmid:25172194

21. Alharbi MA, Al-Kafaji G, Ben Khalaf N, Messaoudi SA, Taha S, Daif A, et al. Four novel mutations in mitochondrial ND4 gene of complex I in Saudi patients with multiple sclerosis. Biomed Rep. 2019;11(6): 257–268. pmid:31798871

22. Andalib S, Talebi M, Sakhinia E, Farhoudi M, Sadeghi-Bazargani H, Gjedde A. Mitochondrial DNA T4216C and A4917G variations in multiple sclerosis. J Neurol Sci. 2015;15;356(1–2): 55–60. pmid:26201854

23. Andalib S, Emamhadi M, Yousefzadeh-Chabok S, Salari A, Sigaroudi AE, Vafaee MS. MtDNA T4216C variation in multiple sclerosis: a systematic review and meta-analysis. Acta Neurol Belg. 2016;116(4): 439–443. pmid:27457485

24. Tranah GJ, Santaniello A, Caillier SJ, D’Alfonso S, Boneschi FM, Hauser SL, et al. Mitochondrial DNA sequence variation in multiple sclerosis. Neurology. 2015;28;85(4): 325–30. pmid:26136518

25. Ban M, Elson J, Walton A, Turnbull D, Compston A, Chinnery P, et al. Investigation of the role of mitochondrial DNA in multiple sclerosis susceptibility. PLoS One 2008;3: e2891. pmid:18682780

26. Polman CH, Reingold SC, Banwell B, Clanet M, Cohen JA, Filippi M, et al. Diagnostic criteria for multiple sclerosis: 2010 revisions to the McDonald criteria. Ann Neurol. 2011;69: 292–302. pmid:21387374

27. Dong C, Wei P, Jian X, Gibbs R, Boerwinkle E, Wang K, et al. Comparison and integration of deleteriousness prediction methods for nonsynonymous SNVs in whole exome sequencing studies. Hum Mol Genet. 2015;15;24(8): 2125–37. pmid:25552646

28. Wallace DC. The mitochondrial genome in human adaptive radiation and disease: on the road to therapeutics and performance enhancement. Gene. 2005;354: 169–180. pmid:16024186

29. Levin L, Zhidkov I, Gurman Y, Hawlena H, Mishmar D. Functional recurrent mutations in the human mitochondrial Phylogeny: Dual roles in evolution and disease. Genome Biol. Evol. 2013;5(5): 876–890. pmid:23563965

30. Kazuno A, Munakata K, Nagai T, Shimozono S, Tanaka M, Yoneda M., et al. Identification of mitochondrial DNA polymorphisms that alter mitochondrial matrix pH and intracellular calcium dynamics. PLoS Genet. 2006;2: e128. pmid:16895436

31. Suissa S, Wang Z, Poole J, Wittkopp S, Feder J, Shutt TE, et al., Ancient mtDNA genetic variants modulate mtDNA transcription and replication. PLoS Genet. 2009;5: e1000474. pmid:19424428

32. Moreno-Loshuertos R, Acín-Pérez R, Fernández-Silva P, Movilla N, Pérez-Martos A, de Cordoba SR, et al. Differences in reactive oxygen species production explain the phenotypes associated with common mouse mitochondrial DNA variants. Nat Genet. 2006;38: 1261–1268. pmid:17013393

33. Mishmar D, Zhidkov I. Evolution and disease converge in the mitochondrion. Biochim Biophys Acta. 2010;1797: 1099–1104. pmid:20074547

34. Chamary JV, Parmley JL, Hurst LD. Hearing silence: non-neutral evolution at synonymous sites in mammals. Nature Reviews Genetics. 2006;7(2): 98–108. pmid:16418745

35. Nackley G, Shabalina SA, Tchivileva IE, Satterfield K, Korchynskyi O. Makarov SS, et al. Human catechol-O-methyltransferase haplotypes modulate protein expression by altering mRNA secondary structure. Science. 2006; 314(5807): 1930–1933. pmid:17185601

36. Komar AA. SNPs, silent but not invisible. Science. 2007; 315(5811): 466–467. pmid:17185559

37. Shabalina SA, Spiridonov NA, Kashina A. Sounds of silence: synonymous nucleotides as a key to biological regulation and complexity. Nucleic Acids Res. 2013;41(4): 2073–2094. pmid:23293005

38. Kimchi-Sarfaty C, Oh JM, Kim I-W, Sauna ZE, Calcagno A M, Ambudkar SV. et al. A silent” polymorphism in the MDR1 gene changes substrate specificity. Science. 2007; 315(5811): 525–528. pmid:17185560

39. Luo Z, Fan P, Zhang S, Li Z. Single nucleotide polymorphisms in the D-loop region predicts earlyage-at-onset of malignant fibrous histiocytoma. Mitochondrial DNA. 2019; 4(1)1: 2078–2083.

40. Umbria M, Ramos A, PilarAluja M, Santos C. The role of control region mitochondrial DNA mutations in cardiovascular disease: stroke and myocardial infarction. Sci Rep. 2020;10: 2766. pmid:32066781

41. Mousavizadeh K, Rajab P, Alaee M, Dadgar S, Houshman M. Usage of mitochondrial D-loop variation to predict risk for Huntington disease. Mitochondrial DNA. 2015;26(4). pmid:24471944

42. Mihailova SM, Ivanova MI, Quin LM, Naumova EJ. Mitochondrial DNA variants in Bulgarian patients affected by multiple sclerosis. Eur J Neurol. 2007;14(1): 44–47. pmid:17222112

43. Yu X, Koczan D, Sulonen A, Akkad DA, Kroner A, Comabella M, et al. mtDNA nt13708A variant increases the risk of multiple sclerosis. PLoS One. 2008; 3(2): e1530. pmid:18270557

44. Kozin MS, Kulakova OG, Kiselev IS, Balanovsky OP, Boyko AN, Favorova OO. Variants of mitochondrial genome and risk of multiple sclerosis development in Russians. Acta Naturae. 2018;10(4): 79–86. pmid:30713765

45. Liou CW, Chuang JH, Chen JB, Tiao MM, Wang PW, Huang ST, et al. Mitochondrial DNA variants as genetic risk factors for Parkinson disease. 2016;23(8):1289–1300.

46. Huerta C, Castro MG., Coto E, Bla´zquez M, Ribacoba R, Guisasola LM, et al. Mitochondrial DNA polymorphisms and risk of Parkinson’s disease in Spanish population. J Neurol Sci. 2005;236: 49–54. pmid:15975594

47. Wallace DC, Brown MD, Lott MT. Mitochondrial DNA variation in human evolution and disease. Gene 1999;238: 211–230. pmid:10570998

48. Moilanen JS, Kari M. Phylogenetic network and physicochemical properties of nonsynonymous mutations in the protein coding genes of human mitochondrial DNA. Mol Biol Evol. 2003;20:1195–1210. pmid:12777521

49. Winklhofer KF, Haass C. Mitochondrial dysfunction in Parkinson’s disease. Biochim Biophys Acta. 2010;1802: 29–44. pmid:19733240

50. Kumleh HH, Riazi GH, Houshmand M, Sanati MH, Gharagozli K, Shafa M. Complex I deficiency in Persian multiple sclerosis patients. J Neurol Sci. 2006;243: 65–69. pmid:16413582

51. Rak M, Benit P, Chretien D, Bouchereau J, Schiff M, El-Khoury R, et al. Mitochondrial cytochrome c oxidase deficiency. Clin Sci (Lond). 2016;130(6): 393–470. pmid:26846578

52. Mahad D, Ziabreva I, Lassmann H, Turnbull D. Mitochondrial defects in acute multiple sclerosis lesions. Brain. 2008;131(7):1722–1735. pmid:18515320

53. Kino K, Sugiyama H. Possible cause of GWCCCWG transversion mutation by guanineoxidation product, imidazolone. Chem Biolo. 2001; 8(4):369–378. pmid:11325592

54. Ohl K, Tenbrock K, Kipp M. Oxidative stress in multiple sclerosis: Central and peripheral mode of action. Exp Neurol. 2016;277: 58–67. pmid:26626971

55. Guo C, Sun L, Chen X, Zhang D. Oxidative stress, mitochondrial damage and neurodegenerative diseases. Neural Regen Res. 2013;25;8(21): 2003–2014. pmid:25206509

56. Jonckheere A, Smeitink J, Rodenbur R. Mitochondrial ATP synthase: architecture, function and pathology. J Inherit Metab Dis. 2012;35: 211–225. pmid:21874297

57. Dautant A, Meier T, Hahn A, Tribouillard-Tanvier D, Rago J-P, Kucharczyk R. ATP synthase diseases of mitochondrial genetic origin. 2018;9:239.

58. Ittisoponpisan S, Islam SA, Khanna T, Alhuzimi E, David A, Sternberg MJE. Can predicted protein 3D structures provide reliable insights into whether missense variants are disease associated? J Mol Biol. 2019; 17;431(11): 2197–2212. pmid:30995449

59. Rudenko G, Bonten E, Hol W, d’Azzo A. The atomic model of the human protective protein/cathepsin A suggests a structural basis for galactosialidosis. PNAS. 1998;95(2): 621–625. pmid:9435242

60. Calyseva J, Vihinen M. PON-SC–program for identifying steric clashes caused by amino acid substitutions. BMC Bioinformatics. 2017;18: 531. pmid:29187139

61. Bleier L, Dröse S. Superoxide generation by complex III: From mechanistic rationales to functional consequences. Biochimica et Biophysica Acta. 2013;1827: 1320–1331. pmid:23269318

62. Hernansanz-Agustin P, Enriquez JA. Generation of reactive oxygen species by mitochondria. Antioxidants. 2021;10(3): 415. pmid:33803273