Substitutions into amino acids that are pathogenic in human mitochondrial proteins are more frequent in lineages closely related to human than in distant lineages

Data

Here, we reanalyzed the phylogenetic data on a set of five mitochondrial protein-coding genes in 4,350 species of opisthokonts (Klink & Bazykin, 2017), as well as on each of the 12 mitochondrial protein-coding genes in several thousand species of metazoans (Breen et al., 2012; Klink & Bazykin, 2017). These two datasets strongly overlap (Table S1). We focused on the amino acids observed in the human mitochondrial genome as reference, polymorphic, or pathogenic variants. A joint alignment of five concatenated mitochondrial genes of opisthokonts and alignments of 12 mitochondrial proteins of metazoans (Breen et al., 2012) were obtained as described in Klink & Bazykin (2017). These alignments were used to reconstruct constrained phylogenetic trees, ancestral states and phylogenetic positions of substitutions as in Klink & Bazykin (2017). Briefly, we obtained taxonomy-based trees from the Interactive Tree of Life Project (Letunic & Bork, 2007), resolving multifurcations and estimating branch lengths using RAxML (Stamatakis, 2014). As the reference human allele, we used the revised Cambridge Reference Sequence (rCRS) of the human mitochondrial DNA (Andrews et al., 1999). As non-reference alleles, we used amino acid changing variants from “mtDNA Coding Region & RNA Sequence Variants” section of the MITOMAP database (Lott et al., 2013). As pathogenic alleles, we used amino acid changing variants with “reported” (i.e., considered by one or more publication as possibly pathogenic) or “confirmed” (i.e., supported by at least two independent laboratories as pathogenic) status from the “Reported Mitochondrial DNA Base Substitution Diseases: Coding and Control Region Point Mutations” section of MITOMAP.

Clustering of substitutions giving rise to the human amino acid around the human branch

For each amino acid site in a protein, we considered those amino acid variants that (i) constitute the reference allele in humans and had arisen in the human lineage at some point during its evolution, or (ii) had originated in humans as a derived polymorphic allele, or (iii) are annotated in humans as pathogenic alleles. For further consideration, we retained only such alleles from each class for which at least one homoplasic (i.e., giving rise to the same allele by way of parallelism, convergence, or reversal) and at least one divergent (i.e., giving rise to a different allele) substitution from the same ancestral variant was observed at this site elsewhere on the phylogeny outside of the human lineage. Substitutions, including reversals, that occurred anywhere on the path between the tree root and H. sapiens were excluded. While the homoplasic and divergent substitutions had to derive from the same ancestral variant, it could be either the same or a different variant than that ancestral to the variant observed in human.

For each such allele, we compared the phylogenetic distances between human and positions of homoplasic substitutions with the distances between human and positions of divergent substitutions, using a previously described procedure which controls for the differences in SPFLs between sites or in mutational probabilities of different substitutions (Klink & Bazykin, 2017). Briefly, for each ancestral amino acid at each site, we subsampled equal numbers of homoplasic substitutions to human (reference, non-reference or pathogenic) amino acids and of divergent substitutions to non-human amino acids. The subsample size was min(N_homo, N_diverg), where N_homo and N_diverg are the numbers of homoplasic and divergent substitutions originating from this amino acid at this site. We repeated this procedure for all ancestral amino acids at all sites, thus obtaining two equal-sized subsamples of homoplasic and divergent substitutions, and measured all distances in these resulting subsamples. We then pooled these values across all considered sites and groups of derived alleles. Next, we categorized them by the phylogenetic distance between the human and the position of the substitution, and calculated, for each bin of the phylogenetic distances, the ratio of the numbers of homoplasic (H) and divergent (D) substitutions (H/D). To obtain the mean values and 95% confidence intervals for the H/D statistic, we bootstrapped sites in 1,000 replicates, each time repeating the entire resampling procedure. As a control, we performed the same analyses using instead of the human variant a random non-human amino acid among those observed at this site, or using data obtained by simulating the evolution at each site along the same phylogeny and with gene-specific GTR+Gamma amino acid substitution matrices using the evolver program of the PAML package (Yang, 1997) as in Klink & Bazykin (2017).

Molecular dynamics simulation of single point mutations in position 91 of COX3

The amino acid at position 91 of COX3 is adjacent to the pore in the complex IV of the respiratory chain which is thought to be required for the transport of oxygen. Therefore, estimation of the effect of the mutation in this position requires a complete description of the cytochrome c oxidase complex in membrane environment. We applied coarse grain description with martini forcefield (Marrink et al., 2007). Human cytochrome C oxidase was modelled with homology modelling approach using Modeller 9.18 (Eswar et al., 2006). The membrane environment was rebuild from a random position of lipids and restrained protein structure as in MemProt MD database (Stansfeld et al., 2015). Each mutant was subjected to 0.5 mks simulation with two replicas in GROMACS 2016 (Abraham et al., 2015) with time step 10fs with PME electrostatics (Wennberg et al., 2015). Simulation results were analyzed with MDAnalysis Python module (Michaud-Agrawal et al., 2011). Water molecules within the pore were counted in a cylindrical selection with radius 2 nm, height 5 nm and center position determined dynamically as the center of mass of the two helices of the dimer harboring mutations. Counts of water molecules and their standard deviations were estimated from the last 100 ns of the trajectory.

Substitutions giving rise to reference human amino acids are more frequent in species closely related to human

The phylogenetic distribution of homoplasic (parallel, convergent or reversing) substitutions is relevant for understanding which amino acids are permitted at a particular species, and which are selected against. In particular, an excess of such substitutions giving rise to a particular allele in closely related lineages implies that the fitness conferred by this allele in these species is higher than in other species. Here, using phylogenies reconstructed from two datasets of mitochondrial protein coding genes obtained earlier (Klink & Bazykin, 2017), we asked how homoplasic substitutions giving rise to the human variant are positioned phylogenetically relative to the human branch, compared with other (divergent) substitutions.

The reference human allele is also observed, on average, in over half of other considered species (71.2% of species of opisthokonts for the 5-genes dataset, and 60% of all species of metazoans for the 12-genes dataset). In 7.1% (10%) of these species, this allele did not share common ancestry with the human, but instead originated independently in an average of 30 (21) independent homoplasic substitutions per site (Table S2 and Figs. S1–S2).

We asked whether the phylogenetic positions of homoplasic substitutions giving rise to the human allele are biased towards the human branch, compared to the positions of divergent substitutions giving rise to other alleles. This analysis controls for the biases associated with pooling sites and amino acid variants (see ‘Materials and Methods’). Although it does not allow us to analyze very conserved sites, most sites were sufficiently variable for this and all subsequent analyses (Tables S1–S3).

In most proteins, the mean phylogenetic distances from human to the branches at which substitutions to the human reference amino acid occurred independently due to a homoplasy were ∼10% shorter than to the branches at which substitutions to another amino acid occurred (Fig. 1A; Fig. S3). No such decrease was observed for a random amino acid among those that were observed at this site in non-human species, or in simulated data (Fig. 1A; Fig. S3). The number of homoplasic substitutions giving rise to the human reference amino acid relative to the number of divergent substitutions towards non-human amino acids (H/D ratio) uniformly decreases with the evolutionary distance from the human (Fig. 1B; Fig. S4). At most genes, the relative number of homoplasic substitutions giving rise to the human allele drops 1.5–3-fold with phylogenetic distance from the human branch (Fig. 1B; Fig. S4), implying that such substitutions preferentially occur in lineages more closely related to humans.

The excess of homoplasic changes to the human reference amino acid at small phylogenetic distances from human is not an artefact of differences in mutation rates between amino acids in distinct clades, since these rates are similar in all considered species and cannot lead to such clustering (Klink & Bazykin, 2017). It is also not an artefact of differences in codon usage bias between species, as it was still observed when we considered only “accessible” amino acid pairs where the derived amino acid could be reached through a single nucleotide substitution from any codon of the ancestral amino acid (Fig. S3).

Substitutions to amino acids corresponding to variant alleles at human polymorphic sites more frequently occur in species closely related to humans

Next, we considered human SNPs in mitochondrial proteins in the MITOMAP database. These SNPs are representative of human polymorphism across multiple populations, although as any such database, it preferentially represents variants that are more frequent. We analyzed the phylogenetic distribution of substitutions in non-human species giving rise to the amino acid that is also observed as the non-reference (usually minor) allele in humans (Table S3).

Similarly to the human reference amino acid variant, the homoplasic substitutions giving rise to non-reference alleles were clustered on the phylogeny near humans, compared to the divergent substitutions giving rise to a variant never observed in humans (Fig. 2; Figs. S5–S6). Again, the mean phylogenetic distance from human to a substitution giving rise to the human non-reference amino acid was ∼10% lower than to other substitutions (Fig. 2A; Fig. S5), and the density of homoplasic substitutions giving rise to such an allele dropped significantly with phylogenetic distance from the human branch (p < 0.001, Fig. 2B; Fig. S6). As before, this clustering was also observed if only accessible pairs of amino acids were considered, while no systematic differences were observed for random amino acids (Fig. S5).

Substitutions to amino acids corresponding to human pathogenic variants more frequently arise in species closely related to humans

Finally, we considered human alleles annotated as disease-causing in the MITOMAP database (Lott et al., 2013). Since only a handful of mutations is thus annotated in each gene (Table S4), the variance in the estimates of the H/D ratio is, as expected, large. Still, in the opisthokont dataset (Fig. 3), as well as in five of the twelve genes of the metazoan dataset (Figs. S7–S8), the human pathogenic variant also arose independently more frequently in the phylogenetic vicinity of humans. The opposite pattern, i.e., biased occurrence of the human pathogenic variant in phylogenetically remote species, has not been observed in any of the genes. This trend was even stronger for the six pathogenic mutations confirmed by two or more independent studies (“confirmed” status in MITOMAP; Fig. 3A). It was observed for the mutations that had been observed as homoplasmic, but not for the mutations that had always only been seen as heteroplasmic (i.e., that were never fixed within an individual) and therefore could confer complete loss of function (Fig. 3A). As before, this result is not due to preferential usage of codons more likely to mutate into the human variant in species closely related to human (Fig. S7).

Human pathogenic variants are biochemically similar to normal human variants

To understand what drives the substitutions to the human variant, either normal or pathogenic, in species closely related to humans, we analyzed the identity of these variants. Both normal (Fig. 4A) and pathogenic (Fig. 4B) human amino acids were more similar in their biochemical properties according to the Miyata matrix (Miyata, Miyazawa & Yasunaga, 1979) to the human reference variant than amino acids observed in non-human species. In turn, amino acids observed in non-human species were more similar to the human reference variant than amino acids never observed at this site in any species.

Individual mutations

To illustrate the observed phylogenetic clustering, we schematically plotted the distribution over the opisthokont phylogeny of substitutions at the six amino acid sites that carry pathogenic mutations with “confirmed” status. Visual inspection of these plots confirms that the substitutions giving rise to pathogenic alleles tend to be clustered in the vicinity of the human, compared to other substitutions of the same ancestral amino acids (Fig. S9). For the metazoan dataset, we also plot three select individual amino acid sites with known pathogenic mutations which are considered below. In all these sites, a higher density of substitutions to a particular amino acid in vicinity of the human branch may imply that its relative fitness is higher in species closely related to human.

The V →A mutation at ND1 site 113 has been reported to cause bipolar disorder, and decreases the mitochondrial membrane potential and reduces ND1 activity in experiments (Munakata et al., 2004). According to the mtDB database (Ingman & Gyllensten, 2006), the A allele persists in human population at 0.5% frequency. However, we observed that the same allele originated independently in three clades of vertebrates: old World monkeys (Cercopithecidae), flying lemurs (Cynocephalidae) and turtles (Geoemydidae), while most of the other substitutions of V at this site occurred in invertebrates (Fig. 5). As a result, the mean phylogenetic distance between human and the parallel V →A substitutions is 2.35 (median 0.75), while it is 4.12 (median 4.6) for substitutions of V to other amino acids.

The V →A mutation at COX3 site 91 has been reported to cause Leigh disease (Mkaouar-Rebai et al., 2011). In metazoans, A allele at this site has originated independently seven times, including six times from V and once from I. All but one of these substitutions occurred in mammals, while tens of substitutions of V and I giving rise to other amino acids occurred throughout metazoans (Fig. 6). As a result, the mean phylogenetic distance from human to V →A substitutions was 0.6 (median 0.7), while the distance to other mutations from V was 1.8 (median 2.1); the corresponding numbers for I were 0.6 (median 0.6) and 2.7 (median 2.2).

To better understand the possible reasons for the unexpected pattern of clustering of the deleterious variant in the phylogenetic vicinity of human, we used molecular dynamics simulations to predict the effect of each mutation on the structure and function of the human protein. As site 91 is positioned within the wall of a pore that is thought to be a channel for oxygen transport (Shinzawa-Itoh et al., 2007), we estimated the pore size that would correspond to each amino acid that occured at this site elswhere on the phylogeny if it arose in the mammalian context. All amino acids led to an increase of the pore size, compared to the normal V allele. Such an increase is expected to permit water molecules to enter the pore, and to impede or prevent oxygen transport. Among the eight observed amino acids, the human pathogenic variant A alters the pore size to the smallest extent, while variants observed in other species substantially increase it, potentially interfering with function (Fig. 7).

Finally, the I →V mutation at ND6 site 33 has been reported to cause type two diabetes (Tawata et al., 2000) and has population frequency of 0.1% according to the mtDB. However, this substitution has occurred repeatedly in parallel in mammals and amphibians, while other substitutions of I were frequent in invertebrates (Fig. 8). The mean phylogenetic distance from human to parallel substitutions to V was 2.4 (median 1.5), while it was as high as 12.9 (median 14.8) for other amino acids.