Mitochondrial genome variation in male LHON patients with the m.11778G > A mutation

Leber hereditary optic neuropathy (LHON) is a mitochondrial disorder with symptoms limited to a single tissue, optic nerve, resulting in vision loss. In the majority of cases it is caused by one of three point mutations in mitochondrial DNA (mtDNA) but their presence is not sufficient for disease development, since ~50% of men and ~10% women who carry them are affected. Thus additional modifying factors must exist. In this study, we use next generation sequencing to investigate the role of whole mtDNA variation in male Polish patients with LHON and m.11778G > A, the most frequent LHON mutation. We present a possible association between mtDNA haplogroup K and variants in its background, a combination of m.3480A > G, m.9055G > A, m.11299 T > C and m.14167C > T, and LHON mutation. These variants may have a negative effect on m.11778G > A increasing its penetrance and the risk of LHON in the Polish population. Surprisingly, we did not observe associations previously reported for m.11778G > A and LHON in European populations, particularly for haplogroup J as a risk factor, implying that mtDNA variation is much more complex. Our results indicate possible contribution of novel combination of mtDNA genetic factors to the LHON phenotype. Electronic supplementary material The online version of this article (10.1007/s11011-020-00605-3) contains supplementary material, which is available to authorized users.


Introduction
Maternally transmitted pathogenic variants of mitochondrial DNA (mtDNA) may lead to dysfunction of the oxidative phosphorylation system and the development of mitochondrial diseases, a heterogeneous group of usually multi-organ disorders that can occur at any age. Symptoms usually affect post-mitotic tissues with high energy demands, such as muscles and nerves (Greaves et al. 2012;Schon et al. 2012). The m.11778G > A point mutation was one of the first mtDNA pathogenic variants to be associated with human disease (Wallace et al. 1988). It was described in patients with hereditary predisposition to blindness, a disease known as Leber hereditary optic neuropathy (LHON) (Wallace et al. 1988).
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s11011-020-00605-3) contains supplementary material, which is available to authorized users. LHON is one of the most peculiar mtDNA diseases. Although mitochondrial disorders are generally rare, LHON is the most common in this group (Chinnery et al. 2000;Kim et al. 2018). Sudden, painless and progressive loss of vision is the only symptom in the majority of cases. This clinical presentation exclusively limited to a single tissue results from retinal ganglion cell damage and loss, and characteristic degeneration of the optic nerve, however, the cause(s) for this phenomenon are still under study (Koilkonda and Guy 2011;Yu-Wai-Man et al. 2011). Over 90% of LHON cases are caused by one of three mtDNA point mutations, m.3460G > A, m.11778G > A or m.14484 T > C, in genes encoding subunits 1, 4 and 6, respectively, of respiratory chain complex I (Mackey et al. 1996;Tońska et al. 2010;Piotrowska et al. 2015;Caporali et al. 2017) and the m.11778G > A variant is the most frequent LHON mutation (accounting for 70-90% of all cases), associated with disease worldwide (Lott et al. 2013;Meyerson et al. 2015;Kim et al. 2018). Although necessary, their presence is not sufficient for disease development, implying that additional modifying factors, genetic, epigenetic and/or environmental, must exist (Yu-Wai- Man et al. 2011;Giordano et al. 2014;Piotrowska et al. 2015;Meyerson et al. 2015;Bianco et al. 2016;Caporali et al. 2017). Besides incomplete penetration, LHON also shows gender bias as symptoms affect~50% of men and only~10% of women who carry one of the three most common LHON mutations (Yu-Wai- Man et al. 2009;Tońska et al. 2010; Koilkonda and Guy 2011;Piotrowska et al. 2015;Meyerson et al. 2015).
One of the most frequently studied and well known modifying genetic factors of LHON are mtDNA haplogroups that are sets of specific common mitochondrial genetic variants segregating together and reflecting mtDNA evolutionary history. Increased risk of vision loss was observed when the m.14484 T > C mutation was present on the background of haplogroup J1, m.11778G > A on J1 and J2 and m.3460G > A on K, whereas haplogroup H had a protective effect in the case of m.11778G > A (Brown et al. 1997;Man 2004;Carelli et al. 2006;Hudson et al. 2007;Piotrowska et al. 2015;Meyerson et al. 2015;Caporali et al. 2017). These observations compellingly suggested that one or more common mtDNA variants defining these haplogroups act synergistically with primary LHON mutations to modify the risk of developing the disease.
The above-mentioned associations were found in European populations, however, studies on mtDNA variation show that there are significant differences between distinct populations in the frequency of individual mtDNA variants and haplogroups, highlighting the necessity for carefully designed, ethnically matched case-control studies. Thus extrapolation of results from one population to another might cause serious errors. With the advent of next-generation sequencing (NGS) methods, high-throughput analysis of whole genomes, including mitochondrial ones, was enabled, allowing for powerful and sophisticated screening for genetic modifiers of hereditary traits. In this study, we use NGS to investigate the mtDNA variation in male Polish LHON patients with the m.11778G > A mutation to confirm previously reported or to find novel mitochondrial genetic modifying factors. We found that haplogroup K and mtDNA variants defining it, m.3480A > G , m.9055G > A , m.11299 T > C and m.14167C > T, might increase the risk of LHON in male patients with the m.11778G > A mutation in the Polish population.

Material
In this work only male cases were studied as men are primarily affected by this disorder. Since LHON shows incomplete penetration and gender bias, women could not be easily enrolled as they very rarely have the disease. The whole mtDNA analysis was performed in a total of 89 Polish male individuals. The patient group consisted of 47 unrelated men (mean age 28.3 ± 9.70) with diagnosed Leber Hereditary Optic Neuropathy and confirmed presence of the m.11778G > A mtDNA mutation. LHON samples included in this study were collected in the Medical University of Warsaw and the Centers for Medical Genetics GENESIS. The control group consisted of 42 unrelated men (mean age 71.0 ± 9.25) treated at the Department of Diagnostics and Microsurgery of Glaucoma, Medical University of Lublin, because of simple senile cataract. LHON (with the presence of any of three most common causal mutations in mtDNA), glaucoma and other intraocular pathologies were excluded. In all cases total DNA collected from peripheral blood was used. Mitochondrial sequence data of control individuals included in this study were extracted from datasets obtained previously in different study (Piotrowska-Nowak et al. 2018) and re-used for the analysis here.

mtDNA sequencing and analysis
High-throughput sequencing of mitochondrial genomes was performed as described previously (Piotrowska-Nowak et al. 2018. In brief, whole mtDNA was amplified using long-range PCR. PCR products were subsequently purified and quantified before proceeding to DNA library preparation for next generation sequencing using Illumina's Nextera XT protocol and sequencing on MiSeq instrument. The NGS data were processed using CLC Genomics Workbench software (CLC bio, Qiagen). Bioinformatic analysis workflow consisted of quality control of sequencing reads, mapping to the human mtDNA reference sequence (rCRS, GenBank sequence NC_012920), variant detection, annotation and evaluation based on the strategy described in detail previously (Piotrowska-Nowak et al. 2018). mtDNA haplogroup assignment was performed for each subject using all variants identified in their mtDNA and the Mitomaster tool (Lott et al. 2013). The variant population frequency was obtained using a current human mtDNA sequence dataset deposited in GenBank and available from MITOMAP database (Lott et al. 2013). The population frequency cutoff value for rare mtDNA variants was set to ≤0.5%. Prediction of the impact of amino acid substitutions and analysis of patient's sequence variant load was performed using a MutPred algorithm as described previously (Pienaar et al. 2017;Venter et al. 2017;Piotrowska-Nowak et al. 2019). Bioinformatic analysis of NGS data in search of large mtDNA deletions was performed using CLC Genomics Workbench InDels and Structural Variants tool and eKLIPse tool (Goudenège et al. 2018).

Statistical analysis
The statistical analysis of mtDNA variation between LHON patients and control subjects was performed using IBM SPSS Statistics and GraphPad Prism software. Mitochondrial haplogroup distribution and frequency of individual mtDNA SNVs were compared using Fisher's exact test together with odds ratio (OR) and 95% confidence intervals (95% CI) calculations of the strength of association. To assess the differences between rare (≤0.5%) and common (>0.5%) variant frequency, in the ratio of transitions to transversions and non-synonymous to synonymous variants in particular mtDNA regions Fisher's exact test and OR with 95% CI were also applied. Independent t-tests were used to compare mean MutPred derived variant loads. The number of subjects having rare variants or variants with MutPred score above 0.5 was compared using Fisher's exact test. The differences were considered statistically significant when the p value was <0.05.

Proven pathogenic variants
In this study we investigated the mtDNA sequence variation in male LHON patients with the m.11778G > A mutation and control male subjects. The mitochondrial genome was covered by an average sequencing depth of 5141 ± 1661X. In all patients the m.11778G > A mutation was found in the homoplasmic state (i.e. in all mtDNA molecules), except for four men in whom high levels of heteroplasmy (i.e. mixture of mtDNA molecules with different sequence) was detected, 72%, 75%, 84% and 92% respectively (subject IDs: 9594, L7-2018, L12-2007 and 10911, respectively). No other proven pathogenic variants were found either in LHON patients or control individuals. However, three distinct mtDNA variants reported previously to associate with LHON (Fauser et al. 2002;Abu-Amero and Bosley 2006;Dai et al. 2018) were identified in three different patients and one control subject, in all in the homoplasmic state, and are shown in more detail in Table 1. All three variants define haplogroups on the background of which they were detected in patients (Oven and Kayser 2009). The presence of large deletions of mtDNA, that is above 1 kb, was also investigated. We did not observe any additional products in LR-PCR indicating the presence of mtDNA molecules with large deletions. Reliable rearrangements in mtDNA of studied participants were also not detected with bioinformatic analysis of NGS data.

mtDNA haplogroup distribution
Based on a full set of sequence variants detected with NGS, each subject's mtDNA was assigned to the appropriate haplogroup. A total of 12 different mtDNA haplogroups were identified in the studied groups. Their observed frequencies in both cohorts are shown in Table 2. Comparison of mtDNA haplogroup distribution showed no statistically significant differences between men with LHON and control subjects. Only a borderline significant difference was noted for haplogroup K and its marginally higher than expected prevalence in LHON patients compared to controls (11% vs 0%, p = 0.057). For the individual mtDNA haplogroup details see Supplementary  Table S1 in Online Resource 1.

SNV association analysis
We also compared the frequency of individual variants in whole mtDNA between LHON patients and healthy individuals. A statistically significant difference was observed for five variants depicted in Table 3. All of them were found in the homoplasmic state and were more prevalent in the patient group and thus associated with increased risk of LHON in the studied men. Four out of five detected variants are markers of haplogroup U8b or haplogroup K which is phylogenetically derived from U8b (based on the mtDNA phylogenetic tree (Oven and Kayser 2009)). Moreover, those variants were identified together, exclusively on the background of the haplogroups they define, largely K, in the same six LHON subjects (except for one case of m.11299 T > C in a patient with haplogroup H). A m.73A > G change is a common variant present in many lineages of the human mtDNA tree (Oven and Kayser 2009;Lott et al. 2013), yet it was overrepresented in the LHON patient group in this study.

Regional variant distribution
We have divided the mitochondrial genome into seven regions based on the function of the encoded products, as described previously (Piotrowska-Nowak et al. 2018): respiratory chain complex I, III and IV coding regions, ATP synthase coding region, rRNA and tRNA coding regions and non-coding region. For each of the studied groups in each region of mtDNA, as well as in the entire mtDNA, the number of transitions (T S ) and transversions (T V ) together with their ratio (T S /T V ) was determined. The obtained results are presented in Table 4. No statistically significant differences in the ratio of transitions to transversions were noted for all categories tested.
Subsequently we determined the number of rare (≤0.5% population frequency) and common (>0.5% population frequency) variants in each region of mtDNA. Their contribution in the sequence variation of particular mtDNA regions in both groups is shown in Fig. 1. Comparison of variant distribution showed a statistically significant difference in the frequency of rare variants in the region covering genes encoding subunits of respiratory chain complex IV (p = 0.048, Table 5). Specifically, rare variants in this region were more frequent in control individuals than in LHON patients (34% vs 19%). Healthy men were more likely to have at least one rare variant in this region compared to LHON patients (45% vs 30%), however, the difference was not large enough to be statistically significant (p = 0.187, OR = 0.514, 95% CI: 0.215-1.228). A large but not significant difference observed in the frequency of rare tRNA variants between both groups was associated with very few variants detected.
We also determined the number of non-synonymous and synonymous variants in each protein coding region of mtDNA and compared their ratios between both groups. Statistical analysis did not show any significant differences between both groups (Table 6).
To avoid false positive results caused by the presence of the m.11778G > A mutation, we excluded it from calculations of transition, rare and non-synonymous variants in complex I region and both 'overall' categories.

MutPred variant load
Variant load calculations were based on MutPred pathogenicity scores as described in detail elsewhere (Pienaar et al. 2017;Venter et al. 2017;Piotrowska-Nowak et al. 2019). In brief, MutPred scores for all non-synonymous variants detected in each subject's mtDNA sequence were summed and the obtained value, total variant load,  served as indicator of predicted mildly deleterious variants in participants. Additionally, as variants with high (>0.5) MutPred score are more likely to be deleterious to protein function, they were used to calculate second variant load, >0.5 threshold, as a clearer indication of a pathogenic load. m.11778G > A mutation was excluded from the calculations. Total and MutPred score > 0.5 variant load distributions in LHON patients and healthy subjects are shown in Fig. 2. No significant differences in the distribution of both variant load scores was observed. The number of subjects carrying variants with MutPred score > 0.5 did not differ significantly between the studied groups, however LHON patients were somewhat more likely to have at least one scored high pathogenicity variant (patients 62%, controls 43%, p = 0.091, OR = 2.148, 95% CI: 0.920-5.017). For the exact variant load values with the number of non-synonymous substitutions see Supplementary Table S1 in Online Resource 1.

Discussion
LHON disease was first described almost 150 years ago (Leber 1871). Although our understanding of its pathogenesis has increased remarkably since then, there are still some aspects to be resolved, such as sex bias and incomplete penetration. In search of genetic modifying factors, we investigated whole mtDNA variation in male Polish LHON patients with m.11778G > A primary mutation.
Screening for the presence of other proven pathogenic variants in mtDNA, both point mutations and large deletions, brought negative results, implying that such deleterious changes do not contribute to disease development. Nevertheless, we identified three non-synonymous variants reported previously to associate with LHON (Fauser et al. 2002;Abu-Amero and Bosley 2006;Dai et al. 2018). However, their role in disease development in a synergistic mechanism, as an additional, secondary mutation, in three Table 4 Number of transitions (T S ), transversions (T V ), and their ratio in particular mtDNA regions in male LHON patients and control subjects. Statistical analysis was performed using the Fischer's exact test. Statistical significance was assumed at the level of p < 0.05. 'Overall' category represents calculations on the entire mtDNA variation and 'overall unique'on the entire and unique (not iterated) mtDNA variation mtDNA region Patient group Control group p OR 95% CI  patients carrying them is debatable since their status was not yet confirmed with functional studies and thus their pathogenicity is unclear. Such variants accompanying primary LHON mutations are reported to be found also in unaffected control subjects, similarly as in one case in this study, and thus are suggested to be rather simple mtDNA variations with neutral effect, but it cannot be completely ruled out that their reported co-occurrence may confer phenotypic variability (Brown et al. 1997;Koilkonda and Guy 2011;Dai et al. 2018).

Number of T S Number of T V T S /T V Number of T S Number of T V T S /T V
Haplogroup and SNP association analysis are commonly used methods applied in the search for mitochondrial genetic factors predisposing to or protecting from disease development. The modifying effect of mitochondrial genetic background on the penetrance of LHON mutations is well documented (Brown et al. 1997Carelli et al. 2006;Hudson et al. 2007;Yu-Wai-Man et al. 2011;Piotrowska et al. 2015;Caporali et al. 2017). In this study mtDNA haplogroup K was overrepresented in LHON men compared to control individuals. Although the observed frequency difference was on the border of statistical significance, probably due to generally low frequency of haplogroup K in Polish population (about 4%, as reported in the latest large study by Jarczak et al. (2019)) and relatively small sample size in our study groups, subsequent SNP analysis revealed that four variants, markers of haplogroup U8b or K, which derives from U8b (Oven and Kayser 2009), associate with increased risk of LHON and have moderate to high effect size (OR > 2.5). These observations lead us to propose that haplogroup K and its characterizing variants, particularly specific combination of m.3480A > G, m.9055G > A, m.11299 T > C and m.14167C > T variants, may have a   (Brown et al. 1997;Man 2004;Carelli et al. 2006;Hudson et al. 2007;Yu-Wai-Man et al. 2011;Meyerson et al. 2015). Moreover, so far haplogroup K was associated with increased risk of LHON for the m.3460G > A primary mutation in European patients (Hudson et al. 2007;Yu-Wai-Man et al. 2011;Meyerson et al. 2015), however, in previous studies this mutation was distributed randomly among mtDNA haplogroups (Torroni et al. 1997;Brown et al. 1997;Carelli et al. 2006). Discrepancies between those and our study might result from marked geographic variation of mtDNA haplogroups. However, it is conceivable that haplogroup K may have similar effects when associated with previously reported m.3460G > A and m.11778G > A described here in different populations. Further studies are needed to unravel this question.
Synergistic and deleterious effect of haplogroup K-related variants on the pathogenic potential of mtDNA mutation may come from subtle conformational changes shifting the assembly kinetics and stability of respiratory chain complexes, as was shown in LHON mutant cybrid cells belonging to different mtDNA haplogroups (Dudkina et al. 2005;Hudson et al. 2007;Pello et al. 2008;Yu-Wai-Man et al. 2011). Particularly interesting is the m.9055G > A variant which leads to an amino acid change in subunit 6 of ATP synthase and possibly may affect the efficiency of energy production in mitochondria. Although the remaining identified SNPs are synonymous variants, it is worth noting that all localize exclusively in genes encoding subunits of the respiratory chain complex I that is fundamental for proper functioning of the respiratory chain and that is affected by primary LHON mutations, including m.11778G > A. It has been recognized that silent sequence changes can impact the secondary structure or stability of mRNA and thus affect protein expression Kimchi-Sarfaty 2011, 2013) what can be true also in this case. Table 6 Non-synonymous (NS) and synonymous (S) variant distribution and their ratios were determined in different mtDNA regions in LHON patients and control individuals. Statistical analysis was performed using the Fischer's exact test. Statistical significance was assumed at the level of p < 0.05. 'Overall' category represents calculations on the entire mtDNA variation and 'overall unique'on the entire and unique (not iterated) mtDNA variation  (Brown et al. 2002). We also identified one common non-coding variant, m.73A > G, to be associated with increased risk of LHON in men with m.11778G > A in this study. It is an ancestral polymorphism present in multiple mtDNA lineages. It localizes in the main non-coding region of mtDNA referred to as the control region as it covers the main regulatory sequences associated with the replication and transcription of the mtDNA molecule (Anderson et al. 1981;Falkenberg et al. 2007). Therefore a non-coding control region variant, such as m.73A > G, will not directly affect the oxidative phosphorylation system altering energy or ROS production but, by changing a regulatory motif or being adjacent to one, it may impact the replication and transcription of the mitochondrial genome e.g. slightly influencing mtDNA copy number or gene expression (Suissa et al. 2009;Lott et al. 2013;Umbria et al. 2018). This together with other risk factors and m.11778G > A mutation may possibly modulate the LHON phenotype.
DNA sequence variants of a harmful nature are removed from the population through purifying selection. It is therefore assumed that pathogenic variants of mtDNA, even those mildly deleterious, will be rare and more frequently observed on the younger branches of the mtDNA phylogenetic tree (Elson et al. 2004;Pereira et al. 2011;Soares et al. 2013;Wei et al. 2017;Venter et al. 2017;Piotrowska-Nowak et al. 2019). According to the hypothesis "common disease -rare variant" (CDRV), a single rare variant may have mild deleterious effect, but such multiple changes may contribute to the development of the disease in a synergistic way by cumulative effect (Elson et al. 2006;Schork et al. 2009;Pienaar et al. 2017;Venter et al. 2017;Piotrowska-Nowak et al. 2019). Regarding mitochondria this seems particularly attractive because all genes in the mitochondrial genome encode products that are involved in the same basic biochemical process, oxidative phosphorylation. In order to verify this hypothesis, we compared the frequency of rare mtDNA variants between patients and controls in individual, arbitrarily determined, mtDNA regions. Analysis of sequence variation showed increased frequency of rare variants in the region covering genes encoding subunits of respiratory chain complex IV in healthy subjects when compared to LHON patients. This suggests that rare variants may also play a protective role in the development of disease, what was already reported previously (Singh et al. 2018).
Not all variants of low population frequency will actually be harmful, thus we also predicted the pathogenicity of variants in silico using the MutPred algorithm and compared variant cumulative impact between two groups. The usefulness and strength of the MutPred tool in estimating the harmfulness of nonsynonymous mtDNA variants has been previously demonstrated (Pereira et al. 2011(Pereira et al. , 2012Pienaar et al. 2017;Venter et al. 2017;Piotrowska-Nowak et al. 2019). In this study we did not observe any significant differences in MutPred variant load between patients and controls. However, this can be attributed to the small number of individuals in study groups and thus should be verified in larger cohorts.
The burden of possibly harmful variants was also investigated by analysis of transversions and non-synonymous variants, known to bear deleterious potential. We did not find any significant differences between the studied groups in any of the mtDNA regions tested what implies lack of association between transversion and non-synonymous variation load and risk of LHON associated with m.11778G > A in men in this study.
In summary, in this study we investigated the whole mtDNA variation in male patients with LHON and its association with the m.11778G > A mutation in the Polish population. Most interestingly, we present a possible association between mtDNA haplogroup K and variants in its background and m.11778G > A mutation. Although the presented data are preliminary with a limited sample size not allowing to make firm conclusions, our results indicate possible contribution of novel combination of mtDNA genetic factors to the LHON phenotype, e.g. by increasing the penetrance of the mutation in this mtDNA background. Surprisingly, we did not observe associations previously reported for m.11778G > A and LHON in European populations, particularly for haplogroup J as a risk factor, what implies that mtDNA variation is much more complex. Further investigation in larger cohorts are required to verify these important findings.
Code availability Not applicable.
Authors' contributions AP-N thought of the work design, carried out the experimental work, provided sequence data, conducted analysis, and wrote the manuscript. MRK referred patients, collected, screened and selected samples, and revised the manuscript. EK-J obtained ethics committee approval, referred patients, collected samples, and revised the manuscript. AMA, MK and MO referred patients, collected samples, and revised the manuscript. KT screened and selected patient samples, managed collaboration and revised the manuscript. EB thought of the work concept and design, obtained ethics committee approval, managed project, collaboration and funding, and revised the manuscript providing critical comments. All authors read and approved the final manuscript. The funding source had no involvement in study design; in the collection, analysis and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethics approval All procedures performed in studies involving human participants were in accordance with the ethical standards of the local research committees (Ethics Committee of the Medical University of Warsaw, approval number KB/187/2015, and Medical University of Lublin, approval number 127/12) and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.
Consent to participate and for publication Informed consent was obtained from all individual participants included in the study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.