Gene 552 (2014) 239–245

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

High prevalence of non-synonymous substitutions in mtDNA of cichlid fishes from Lake Victoria Kazumasa Shirai a, Nobuyuki Inomata b, Shinji Mizoiri c, Mitsuto Aibara d, Yohey Terai e, Norihiro Okada d,f, Hidenori Tachida g,⁎ a

Graduate School of Systems Life Sciences, Kyushu University, Fukuoka, Japan International College of Arts and Sciences, Fukuoka Women's University, Fukuoka, Japan Sekisui Chemical Tanzania Ltd., Dar es Salaam, Tanzania d Foundation for Advancement of International Science, Tsukuba, Japan e The Graduate University for Advanced Studies, Kanagawa, Japan f Department of Life Sciences, National Cheng Kung University, Tainan 701, Taiwan g Department of Biology, Faculty of Sciences, Kyushu University, Fukuoka, Japan b c

a r t i c l e

i n f o

Article history: Received 14 June 2014 Received in revised form 28 August 2014 Accepted 17 September 2014 Available online 18 September 2014 Keywords: Adaptive substitution Branch model Branch-site model Cichlids Population size Slightly deleterious mutation

a b s t r a c t When a population size is reduced, genetic drift may fix slightly deleterious mutations, and an increase in nonsynonymous substitution is expected. It has been suggested that past aridity has seriously affected and decreased the populations of cichlid fishes in Lake Victoria, while geographical studies have shown that the water levels in Lake Tanganyika and Lake Malawi have remained fairly constant. The comparably stable environments in the latter two lakes might have kept the populations of cichlid fishes large enough to remove slightly deleterious mutations. The difference in the stability of cichlid fish population sizes between Lake Victoria and the Lakes Tanganyika and Malawi is expected to have caused differences in the nonsynonymous/synonymous ratio, ω (=dN/dS), of the evolutionary rate. Here, we estimated ω and compared it between the cichlids of the three lakes for 13 mitochondrial protein-coding genes using maximum likelihood methods. We found that the lineages of the cichlids in Lake Victoria had a significantly higher ω for several mitochondrial loci. Moreover, positive selection was indicated for several codons in the mtDNA of the Lake Victoria cichlid lineage. Our results indicate that both adaptive and slightly deleterious molecular evolution has taken place in the Lake Victoria cichlids' mtDNA genes, whose nonsynonymous sites are generally conserved. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Most species have experienced reductions in population size for various reasons (e.g., climate changes, including droughts and floods; volcanic eruptions; competition and human activities). Such reductions may have large effects on molecular evolution. The nearly neutral theory asserts that fixations of slightly deleterious mutations occur ubiquitously in molecular evolution (Ohta, 2002; Akashi et al., 2012). When a population size is large, slightly deleterious mutations are effectively removed from the population by natural selection. Conversely, slightly deleterious mutations may be fixed in small populations by genetic drift (Ohta, 1973). Because nonsynonymous substitutions are more strongly affected by natural selection than are synonymous substitutions, fixations of slightly

Abbreviations: LM, Lake Malawi; LT, Lake Tanganyika; LV, Lake Victoria; mtDNA, mitochondrial DNA. ⁎ Corresponding author at: Department of Biology, Faculty of Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan. E-mail address: [email protected] (H. Tachida).

http://dx.doi.org/10.1016/j.gene.2014.09.039 0378-1119/© 2014 Elsevier B.V. All rights reserved.

deleterious mutations are expected to increase the nonsynonymous substitution rate. The observation that the nonsynonymous to synonymous substitution ratios, ω (dN/dS), are smaller in 49 nuclear genes in rodents, considered to have larger population sizes, than in those of primates and Artiodactyla, thought to have smaller population sizes (Ohta, 1995), supports this theory. The results of studies on mitochondrial DNA sequences in island birds (Johnson and Seger, 2001), mammals (Popadin et al., 2007) and dogs (Björnerfeldt et al., 2006) also agree with the predictions of this theory. Here, we examine this prediction of the nearly neutral theory, i.e., ω increases when population sizes become small, by studying fishes belonging to the Cichlidae family in African lakes. The Cichlidae is known as one of the most species-rich families of vertebrates. In particular, there are a large number of endemic species of cichlid fishes in three of the African Great Lakes: Lakes Victoria (LV), Tanganyika (LT) and Malawi (LM). The number of endemic cichlid species is approximately 700 in LV, 250 in LT and 700 in LM (Turner et al., 2001). The species flock in LT is the most diverse phenotypically and genetically of the three lakes (Nishida, 1991; Salzburger et al., 2002), and these species

240

K. Shirai et al. / Gene 552 (2014) 239–245

can be categorized into 12 tribes (Poll, 1986). Conversely, most of the species flocks in LV and LM belong to only one tribe: Haplochromini (Fryer and Iles, 1972; Greenwood, 1979; Meyer et al., 1990; Verheyen et al., 2003). Haplochromini forms a sister group to an endemic LT tribe Tropheini (Salzburger et al., 2005). This suggests that Haplochromini originated in LT. The species flocks in LM and LV are different in terms of the degree of divergence. The species flock in LM can be categorized genetically into three groups: pelagic species, rock-dwelling species and sand-dwelling species (Meyer et al., 1990; Shaw et al., 2000). In LV, the phylogenetic relationship of the species flock is more complicated than in the other two lakes (Samonte et al., 2007; Wagner et al., 2013; Takeda et al., 2013). Thus, the degree of the genetic and phenotypic divergence of the species flocks is different among the three lakes. In addition, these three lakes are very different in terms of their depths. Lake Victoria is shallow compared to the other two lakes; its maximum depth is only 79 m (Bootsma and Hecky, 2003). In comparison, the other two lakes are deep. The maximum depth is 1470 m in LT and 700 m in LM (Bootsma and Hecky, 2003). The different depths of the lakes have resulted in very different histories for the cichlids inhabiting them (Danley et al., 2012). Geological studies have revealed that these three lakes have experienced several declines in water level (Delvaux, 1995; Cohen et al., 1997; Lezzar et al., 1996). Lake Victoria dried up at least twice between 17,000 and 14,000 years ago (Johnson et al., 1996; Stager et al., 2002, 2011). Such dramatic changes are expected to have decreased the population size of cichlids in the lake. Indeed, this expectation is consistent with indications from microsatellite data that the population sizes of the cichlids in the LV region superflock started to decline approximately 18,000 years ago (Elmer et al., 2009). Although the other two lakes have lost much of their water in the past in several desiccations, they have kept their water levels high (N100 m, Scholz et al., 2007). Thus, the degree of reduction in cichlid population size in these deep lakes might have been less drastic than that of the cichlid fishes in LV. In summary, lake depth seems to have affected cichlid populations by influencing the magnitude of dry-ups in the lakes. More specifically, the less stable environment of LV, due to its low water level, is thought to have affected the sizes of its cichlid populations. These examples suggest that differences in habitat stability might have caused differences in molecular evolution among the cichlids in these three African Great Lakes. In LV, past desiccations have dried up the water, possibly reducing their cichlid population sizes drastically. Such a small population size might have allowed fixations of slightly deleterious mutations. Conversely, in LT and LM, the cichlid population sizes might have remained relatively stable and large because of the lakes' depths. Because natural selection operates effectively in large populations, slightly deleterious mutations must have been excluded from the populations in these two lakes. Consequently, we expect to observe a difference in ω between the fluctuating environment (the shallow lake (LV)) and the stable environments (the deep lakes (LT and LM)). In the present study, we sequenced 13 mitochondrial genes in cichlid fishes from LV, LT and LM. From these sequences, we estimated ω with the likelihood-based method using the branch model developed by Yang (2007). By comparing the estimates of ω for the cichlids in fluctuating environment with those in stable environments, we examined whether fixations of slightly deleterious mutations occurred more frequently in the cichlids inhabiting fluctuating environments. Although an increase of ω may be the result of fixations of slightly deleterious mutations, it may also occur due to positive natural selection. Indeed, previous studies have indicated action of positive selection in mitochondrial evolution (Ballard and Kreitman, 1994; Bazin et al., 2006). Although the proportion of adaptive substitutions was previously thought to be larger in larger populations, this predicted relationship has been shown to be model-dependent. A recent study on a model involving environmental changes showed only a moderate effect of

population size but large effect of environmental changes on adaptive substitutions (Lourenço et al., 2013). Considering expected environmental fluctuations in LV, an inflation of ω might occur via positive selection. Therefore, to evaluate the contribution of positive selection to increasing ω, we applied another likelihood-based method using the branch-site model developed by Yang (2007). In addition, we estimated the ratio of radical to conservative amino acid changes (the Kr/Kc ratio, see Hughes et al., 1990) to complement the analysis of ω. We found significant increases in ω in some of the mtDNA genes in the lineage of the LV cichlids. In addition, positive selection is shown to have operated for some codon positions in the mitochondrial genes with ω exceeding one in the LV lineage.

2. Materials and methods 2.1. Samples We chose samples so that most of the major cichlid lineages in the three African Great Lakes were included (Table S1). Haplochromis pyrrhocephalus was used as the representative for LV. We used only one species from LV because the species in LV are very closely related (Nagl et al., 1998, 2000; Seehausen et al., 2003; Terai et al., 2004). Indeed, a preliminary sequencing of five mtDNA genes in two species, H. pyrrhocephalus and Haplochromis (Paralabidochromis) sauvagei (Pfeffer, 1896) (this species was also known as Haplochromis sp. “rockkribensis” sensu Seehausen, 1996) showed that only one nonsynonymous substitution occurred between the genes in these two species (the sequences of the H. (Paralabidochromis) sauvagei are deposited with DDBJ/EMBL/GenBank (accession number AB968212–AB968216)). H. (Paralabidochromis) sauvagei is considered to have diverged earlier from the other cichlid species in LV (Samonte et al., 2007). Therefore, number of nonsynonymous substitutions between the mtDNA genes in the LV cichlids seems small. Three species, Labidochromis caeruleus, Dimidiochromis compressiceps and Rhamphochromis sp., were chosen to represent the three major LM lineages. We chose four species based on the classification of Salzburger et al. (2002); Boulengerochromis microlepis belonging to Tilapiini that is thought to be an old lineage in LT, Julidochromis transcriptus belonging to Lamprologini that is the largest tribe in LT, Cyprichromis leptosoma belonging to Cyprichromini in the H-lineage that includes Haplochromini and Tropheni, and Tropheus duboisi belonging to Tropheini that is most closely related to Haplochromini, as representatives of the LT lineages. In addition, a riverine species, Astatotilapia burtoni, was used. This species is distributed near the mouths of the rivers flowing into Lake Tanganyika. A. burtoni is more closely related to the LV flock than to the LM flock (Meyer et al., 1991). As an outgroup for all these species, Oreochromis niloticus was used.

2.2. Sequence data We analyzed 13 mitochondrial protein-coding genes (Table S2). There are 13 protein-coding genes in the mtDNA of cichlid fishes. We used all 13 genes: ND1, ND2, ND3, ND4L, ND4, ND5, ND6, CO1, CO2, CO3, CYTB, ATP8 and ATP6. Complete genomes of mtDNA have been determined in T. duboisi and O. niloticus (Mabuchi et al., 2007; He et al., 2010; GenBank accession number AP006015 and GU370126). PCR and sequencing primers were designed from these sequences (Table S3). In addition, some published primers were used for the amplification of ND2, CO1 and CYTB (Kocher et al., 1989, 1995; Duftner et al., 2005; Chakrabarty, 2006). The details of the conditions for the PCR amplification and sequencing are listed in Data S1 and Table S4. We sequenced more than 90% of these genes, totaling 3324 codons. Indel variants were not used for the analysis. These sequences have been deposited with DDBJ/EMBL/GenBank (accession number AB915406–AB915509).

K. Shirai et al. / Gene 552 (2014) 239–245

2.3. Estimation of ω The obtained sequences were aligned using Muscle v3.5 (Edgar, 2004). A phylogenetic tree was reconstructed using the maximum likelihood method implemented in PhyML v3.0 (Guindon et al., 2010). The best-fit substitution model of the tree (TIM2 + I + G) was selected by jModelTest v2 (Darriba et al., 2012). Using the codeml program in PAML (Yang, 2007), ω (=dN/dS) was estimated from the aligned sequences with a maximum likelihood method. We used the F3 × 4 model, which uses the equilibrium codon frequencies calculated from the average nucleotide frequencies for the three codon positions (Yang and Nielsen, 1998). Three different values (0.5, 1.0 and 2.0) were used as initial ω in order to avoid local optima in the maximum likelihood estimations. The Bonferroni correction was conducted on multiple tests. Because the Bonferroni correction is very conservative, we present the results both with and without the correction. We used two branch models in codeml: the one ratio model, which assumes that all branches have the same ω, and the two ratios model, which assumes that a particular branch or group of branches has an ω different from those of the other branches. A likelihood ratio test between the two models can reveal whether a particular lineage or lineages have a significantly different ω from those of the other lineages.

241

because we wanted to know the effects of a reduction in population size on different classes of amino acid substitutions. Here, four classifications of amino acids (A: charge, aromaticity and volume; B: polarity and volume; C: charge and aromaticity; D: charge and polarity) were employed (Zhang, 2000; Hanada et al., 2006; Hanada et al., 2007), and we estimated the Kr/Kc ratio using the Zhang method (Zhang, 2000). The ancestral sequences were estimated for all ancestral nodes in the phylogenetic tree of the cichlids by a maximum likelihood method using PAML (Yang, 2007), and these sequences were used for the estimation of the Kr/Kc ratios.

3. Results 3.1. Phylogenetic analysis Fig. 1 shows the maximum likelihood tree reconstructed using all 13 mtDNA genes combined together. This tree agrees with trees from past studies that investigated the phylogenetic relationships among the African cichlid fishes (Meyer et al., 1990; Shaw et al., 2000; Salzburger et al., 2002). This phylogenetic tree was used to estimate ω and Kr/Kc. In the following, we call the terminal branches leading to H. pyrrhocephalus (LV) the LV lineage.

2.4. Detection of positive selection To detect positive selection, we used the model M2 (selection) on the branch-site model (Yang and Nielsen, 2002), implemented in PAML (Yang, 2007). The branch-site model estimates ω for branches and sites. The branches in which positive selection is suspected to occur are referred to as the foreground branches. The remaining branches are referred to as the background branches. In the background branches, each codon site has either ω0 or ω1 (conserved sites have ω0 b 1; neutral sites have ω1 = 1). In the foreground branches, each codon site has either ω0, ω1 or ω2 (sites affected by positive selection have ω2 N 1). The branch-site model assumes four site classes. In site classes 0 and 1, ω = ω0 and ω = ω1, respectively, in all branches. In site class 2a, ω = ω0 in the background branches, and ω = ω2 in the foreground branches. In site class 2b, ω = ω1 in the background branches, and ω = ω2 in the foreground branches. The model with ω2 = 1 was used as the null hypothesis to avoid false positive selection (Zhang et al., 2005). A likelihood ratio test between the two models can detect positive selection. If the likelihood ratio test rejects the null hypothesis and significantly shows that ω2 N 1, it indicates that the branches underwent positive selection. In addition, we used the Bayes empirical Bayes (BEB) approach (Yang et al., 2005) to identify the sites with positive selection in the genes for which the likelihood ratio test gave a significant result. We call the site positively selected if the posterior probability for the site is not less than 0.95.

3.2. Comparison of ω estimated by branch models (cichlids in LV vs. other cichlids) To compare ω between the LV lineage and the other lineages, an ω (ωLV) was assigned to the LV lineage (Fig. 1; branch A), and another independent ω (ωOther) was assigned to all the other branches in the two ratios model. Fig. 2 and Table S5 show the results of this estimation for the mtDNA genes. The estimate of ω for the LV lineage was about three times larger than that for the other lineages when the sequences of all mtDNA genes were concatenated in the analysis (all mtDNA genes: ωLV = 0.1371, ωOther = 0.0439, P = 5.6186 × 10−17; dN = 0.0163, dS = 0.1191, number of nonsynonymous substitutions = 114.4, number of synonymous substitutions = 352.8). When we analyzed individual mtDNA genes, the estimate of ω for the LV lineage was significantly higher than that for the other lineages for five genes after the Bonferroni correction (ND2; ωLV = 0.2537, ωOther = 0.0739, P = 0.0004: ND4L; ωLV = 0.7537, ωOther = 0.0351, P = 0.0015: CO2; ωLV = 0.2214, ωOther = 0.0217, P = 2.9861 × 10− 5: ATP6; ωLV = 0.3212, ωOther = 0.0341, P = 7.0655 × 10 − 8 : ATP8; ωLV = 1.8676, ωOther = 0.1433, P = 3.6111 × 10− 5 (if P b 0.0038, the difference is significant at the 5% level after Bonferroni correction)). Note that the estimate of ω for the LV lineage was more than one for ATP8, suggesting that positive selection occurred for this gene in this lineage.

2.5. Estimation of Kr/Kc The ratio of the radical replacement rate (Kr) to the conservative replacement rate (Kc) can be used to measure selection pressure in addition to the ω ratio (Hughes et al., 1990). This is because the effects of radical replacements on the protein structure are likely stronger than those of conservative replacements, and thus, radical replacements should be more strongly selected against (Smith, 2003). Therefore, if a population size decreases, some radical replacements may become slightly deleterious, and therefore, the Kr/Kc ratio would increase (Smith, 2003; Wernegreen, 2011). One of the reasons why the Kr/Kc ratio is used instead of ω to detect slightly deleterious substitutions is that the saturation of synonymous substitutions can cause the overestimation of ω if distantly related species are used to estimate substitutions (Smith and Smith, 1996). In the present study, only closely related species were used; nevertheless, we estimated the Kr/Kc ratio

Fig. 1. A phylogenetic tree reconstructed from 13 mitochondrial loci by maximum likelihood method. Only its topology is shown. Branch A is the lineage of the cichlids in Lake Victoria (H. pyrrhocephalus).

242

K. Shirai et al. / Gene 552 (2014) 239–245 Table 2 The sites affected by positive selection in H. pyrrhocephalus. Locus

Number of amino acid site

Ancestral amino acid

Amino acid of H. pyrrhocephalus

Posterior probability (ω N 1)

ND2 ND5 CO2 ATP8

123 269 117 25 30 38 42

D V P P A T T

S T F I Q L L

0.968 0.971 0.971 0.979 0.984 0.973 0.979

Ancestral amino acids were estimated by maximum likelihood method.

Fig. 2. Estimates of ω in the cichlids in Lake Victoria and the other areas. ω was estimated with the two ratios model (LV cichlids vs. other cichlids) of the branch model. Gray bars show the ω of the cichlids in Lake Victoria (H. pyrrhocephalus). White bars show the ω of the other cichlids. All mtDNA is the sequence of the 13 mitochondrial genes combined together. The other 13 loci shown are mitochondrial loci. The asterisks show the significance (P b 0.05) of the likelihood ratio test results from the two ratios model (LV cichlids vs. other cichlids) and the one ratio model.

for the LV lineage were higher than those for the other lineages, as in the case of ω. Thus, radical amino acid substitutions occurred more frequently in the LV than in the other lineages, indicating weaker selection in the LV lineage, as was shown by the change in ω.

3.3. Detection of positive selection by branch-site model 4. Discussion To detect positive selection, the branch-site model was used on the 13 mtDNA protein-coding genes. First, we assigned the branch representing the LV lineage (Fig. 1; branch A) as the foreground branch. To increase the power of detection and to avoid false positives (Zhang et al., 2005; Nozawa et al., 2009), we concatenated all 13 genes and analyzed the concatenated sequence. Table 1 shows the results. There were codon sites with an ω significantly larger than one, suggesting positive selection had occurred for some of the mitochondrial genes in the LV lineage. The BEB method was used to identify putative selected codons in the mtDNA genes of the LV lineage for which the branch-site model suggested positive selection. Table 2 shows the results. Among 112 codons that have changed in the LV lineage, seven codons in four mitochondrial genes (ND2, ND5, CO2 and ATP8) were identified as positively selected by this method (posterior probability ≥ 0.950). Notably, ATP8 contained four of the seven positively selected codons.

The sizes of the cichlid populations in LV are thought to fluctuate because of unstable environment, while those in LT and LM are relatively stable because LT and LM are deep. These considerations led us to expect ω would be higher in the LV cichlid lineage than in the LT and LM lineages because nonsynonymous mutations, considered to be mostly deleterious, accumulate when population size is small (Ohta, 1973; Akashi et al., 2012). Indeed, we found ω was significantly higher in the LV than in the LT and LM lineages for several mtDNA genes. Surprisingly, ω was significantly larger than one for seven codon sites of mtDNA in the LV lineage, indicating positive selection in the mtDNA of this lineage.

3.4. Estimation of ω excluding positively selected codons In three of the five genes that showed a significantly higher ω in the LV lineage, positively selected codons were found. To estimate ω excluding the effects of inflation caused by positive selection, we estimated the ω of the LV lineage for these three genes and for all mtDNA genes combined using the branch model excluding the positively selected codons. Fig. 3 and Table S6 show the results. Although the positively selected codons were excluded, the LV lineage still showed a significantly higher ω for the three genes (P b 0.0167: significant after Bonferroni correction) and in all mtDNA genes combined. 3.5. Estimation of Kr/Kc The Kr/Kc ratio was estimated for all mitochondrial genes combined using four classifications (A: charge, aromaticity and volume; B: polarity and volume; C: charge and aromaticity, D: charge and polarity). Fig. 4 and Table S7 show the estimates of Kr/Kc for the LV (branch A in Fig. 1) and other lineages. In three of the classifications, the Kr/Kc ratios Table 1 Results of analyses using branch site model. Species

ω2 (ω ≥ 1)

Proportion of ω2 sites

2ΔlnL

P-value

H. pyrrhocephalus

3.1025

0.0341

6.2670

0.0123⁎

⁎ Significant without Bonferroni correction P b 0.05.

Fig. 3. Estimates of ω in the cichlids in Lake Victoria and the other areas, excluding the sites affected by positive selection. ω was estimated with the two ratios model (LV cichlids vs. other cichlids) of the branch model. Gray bars show the ω of the LV cichlids (H. pyrrhocephalus). White bars show the ω of the other cichlids. All mtDNA is the sequence of the 13 mitochondrial genes combined together. The other three genes are mitochondrial. The asterisks show the significance (P b 0.05) of the likelihood ratio from the two ratios model (LV cichlids vs. other cichlids) and the one ratio model.

K. Shirai et al. / Gene 552 (2014) 239–245

Fig. 4. Kr/Kc calculated by the Zhang method (Zhang, 2000) using four classifications. White and gray bars show the Kr/Kc ratios in the other and LV lineages, respectively. A is the classification based on charge, aromaticity and volume. B is classified by polarity and volume. C is classified by charge and aromaticity. D is classified by charge and polarity (see Zhang, 2000; Hanada et al., 2006; Hanada et al., 2007).

4.1. Accumulation of slightly deleterious mutations in Lake Victoria The estimates of ω for ND2, ND4L, CO2, ATP6, ATP8 and in all mtDNA genes combined in the LV lineage were significantly higher than those in the other lineages. Although positive selection was indicated for seven codon sites, the analysis excluding these codon sites also showed a significant elevation of ω in the LV lineage. One possible explanation for this observation is the accumulation of slightly deleterious mutations, as we originally hypothesized. Because Lake Victoria is shallow, its cichlid populations were likely severely affected by past desiccations. These desiccations must have drastically reduced the population size of the cichlids in the lake. Under such conditions, slightly deleterious mutations would be fixed. This explains why the ω for five mitochondrial genes in the LV lineage was significantly higher than those in the other lineages. The higher Kr/Kc ratios in the LV lineage for the mtDNA genes also support this hypothesis of fixations of slightly deleterious mutations. Fixations of slightly deleterious mutations in populations with reduced sizes have also been observed in other organisms (Akashi et al., 2012; Smith, 2003; Wernegreen, 2011). In mtDNA, because there is no recombination, slightly deleterious mutations are not removed by recombination, and therefore, they accumulate in mtDNA. To restore the fitness, compensatory mutations may play a key role (Rand et al., 2004). For example, Osada and Akashi (2012) have suggested that mutations in the nuclear genes that code for subunits of mitochondrial complexes compensate for slightly deleterious mutations in some of the mtDNA genes in mammals. Although these compensatory substitutions occurred in nuclear genes, Oliveira et al. (2008) suggested that substitutions compensating for slightly deleterious mutations fixed by drift or Wolbachia sweeps occurred in mtDNA genes in the parasitic wasp Nasonia. Such compensatory processes might have occurred in the mtDNA in the LV lineage and also may have contributed to its high ω values. 4.2. Adaptive evolution of mtDNA genes In this study, positive selection was indicated for seven codon sites of the mtDNA in the LV lineage. In fact, these adaptive substitutions occurred in very conservative sites. To know the extent of the conservation of amino acids in these sites, we added the sequence data from two South American cichlids (Astronotus ocellatus and Hypselecara temporalis; Mabuchi et al., 2007; Azuma et al., 2008; GenBank accession number AP009127 and AP009506) to ours and aligned their amino acid

243

sequences (Table 3). The divergence time between the cichlids in East Africa and those in South America has been estimated as approximately 100 to 120 MYA (Azuma et al., 2008). Six of the seven positively selected amino acids are conserved in all cichlids except for the LV species H. pyrrhocephalus. The changes found in H. pyrrhocephalus were radical in classification A of Hanada et al. (2007). The radical amino acid substitutions in sites highly conserved among cichlid species must have large fitness effects. The 13 mitochondrial genes encode subunits of the complexes involved in oxidative phosphorylation. Although mtDNA protein-coding genes are generally conserved, several studies have reported adaptive evolution of mtDNA genes. An example is the positive selection detected for CO1 in billfish. It has been suggested that this substitution might be related to the unique endothermy of billfish (Dalziel et al., 2006; Little et al., 2012). Also Shen et al. (2010) found accelerated evolution of respiratory enzymes including those encoded by mtDNA in the ancestral lineage leading to bats, which they considered to be related to evolution of flights. These studies indicate that mtDNA genes, although usually conserved, might be changed in adaptation to new environments. Furthermore, a survey of genetic variation in a broad-range of species indicated action of positive selection on mtDNA (Bazin et al., 2006) though other complicating factors need to be considered to fully understand mtDNA evolution (Meiklejohn et al., 2007). The relationship between the rate of adaptive evolution and effective size (Ne) was originally thought to be positive but there are many complications in empirically showing this relationship (see Lanfear et al., 2014; Gossmann et al., 2012). On the other hand, Lourenço et al. (2013) have theoretically shown that the rate of adaptive evolution is positively related to the rate of environmental changes and organismal complexity. This suggests that the fluctuating environment of LV might have promoted adaptive molecular evolution at some of the mtDNA genes of the cichlids in the lake. When positively selected nucleotides were fixed in the population, slightly deleterious mutations in mtDNA genes would be also fixed by hitchhiking. The hitchhiking, if it occurs in mtDNA, affects all nucleotides of mtDNA because mtDNA lacks recombination. This may have contributed to the increase of ω found in the mtDNA of the LV lineage. 5. Conclusion The present study revealed that the LV lineage experienced a significant increase in ω and had positively selected codon sites with ω N 1 for some mtDNA genes. This result is consistent with our original hypothesis that cichlids in a fluctuating environment have a higher ω caused by the accumulation of slightly deleterious mutations. However, the increased ω could be explained by other factors (adaptive evolution, compensatory substitutions and hitchhiking). Determining to what extent respective factors affected the evolution in mtDNA of cichlids in LV is Table 3 Evolution of putative positively selected amino acid sites among cichlids in East Africa and South America. Locus

ND2

ND5

CO2

ATP8

Distribution

No. of amino acid site

123

269

117

25

30

38

42

New world

H. temporalis A. ocellatus O. niloticus B. microlepis J. transcriptus C. leptosoma T. duboisi Rhamphochromis sp. L. caeruleus D. compressiceps A. burtoni H. pyrrhocephalus

D . . . . . . . . . . S

V . . . . . . . . . . T

P . . . . . . . . . . F

P . . . . . . . . . . I

A . . . . . . . . . . Q

T A . . . . . . . . . L

T . . . . . . . . . . L

East Africa

A dot indicates that the amino acid is the same as that of H. temporalis at the site.

244

K. Shirai et al. / Gene 552 (2014) 239–245

difficult with only the DNA sequence information obtained here. This question may be answered through analyses of the structures of mitochondrial proteins or through estimating the fitness of these cichlids. Although the cause of the large ω is currently unknown, this study revealed some unique features of the molecular evolution of the mtDNA of cichlids in LV. In addition to African cichlids, comparisons of molecular evolutionary rates between species in stable and fluctuating environments may reveal other similar cases. Acknowledgments We thank two referees for giving us useful comments on an earlier draft of this manusript. This work was partially supported by a Grantin-Aid for Scientific Research on Priority Areas no. 16057201 to H. T and Grants-in-Aid for Scientific Research on Priority Areas nos. 14087101 and 14087202 by the Ministry of Education, Culture, Sports, Science, and Technology of Japan to N. O. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.09.039. References Akashi, H., Osada, N., Ohta, T., 2012. Weak selection and protein evolution. Genetics 192, 15–31. Azuma, Y., Kumazawa, Y., Miya, M., Mabuchi, K., Nishida, M., 2008. Mitogenomic evaluation of the historical biogeography of cichlids toward reliable dating of teleostean divergences. BMC Evol. Biol. 8, 215. Ballard, J.W.O., Kreitman, M., 1994. Unraveling selection in the mitochondrial genome of Drosophila. Genetics 138, 757–772. Bazin, E., Glémin, S., Galtier, N., 2006. Population size does not influence mitochondrial genetic diversity in animals. Science 312, 570–572. Björnerfeldt, S., Webster, M.T., Vilà, C., 2006. Relaxation of selective constraint on dog mitochondrial DNA following domestication. Genome Res. 16, 990–994. Bootsma, H.A., Hecky, R.E., 2003. A comparative introduction to the biology and limnology of the African Great Lakes. J. Great Lakes Res. 29, 3–18. Chakrabarty, P., 2006. Systematics and historical biogeography of Greater Antillean Cichlidae. Mol. Phylogenet. Evol. 39, 619–627. Cohen, A.S., Lezzar, K.E., Tiercelin, J.J., Sorgehan, M., 1997. New palaeogeographic and lake level reconstructions of Lake Tanganyika: implications for tectonic, climatic and biological evolution in a rift lake. Basin Res. 9, 107–132. Dalziel, A.C., Moyes, C.D., Fredriksson, E., Lougheed, S.C., 2006. Molecular evolution of cytochrome c oxidase in high-performance fish (Teleostei: Scombroidei). J. Mol. Evol. 62, 319–331. Danley, P.D., Husemann, M., Ding, B., DiPietro, L.M., Beverly, E.J., Peppe, D.J., 2012. The impact of the geologic history and paleoclimate on the diversification of East African cichlids. Int. J. Evol. Biol. 2012. http://dx.doi.org/10.1155/2012/574851. Darriba, D., Taboada, G.L., Doallo, R., Posada, D., 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat. Methods 9, 772. Delvaux, D., 1995. Age of Lake Malawi (Nyasa) and water level fluctuations. Mus. R. Afr. Centr. Tervuren. (Belg.) Dept. Geol. Min. Rapp. Ann. 1995–1996, 99–108. Duftner, N., Koblmüller, S., Sturmbauer, C., 2005. Evolutionary relationships of the limnochromini, a tribe of benthic deepwater cichlid fish endemic to Lake Tanganyika, East Africa. J. Mol. Evol. 60, 277–289. Edgar, R.C., 2004. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinforma. 5, 113. Elmer, K.R., Reggio, C., Wirth, T., Verheyen, E., Salzburger, W., Meyer, A., 2009. Pleistocene desiccation in East Africa bottlenecked but did not extirpate the adaptive radiation of Lake Victoria haplochromine cichlid fishes. Proc. Natl. Acad. Sci. U. S. A. 106, 13404–13409. Fryer, G., Iles, T.D., 1972. The Cichlid Fishes of the Great Lakes of Africa: Their Biology and Evolution. Oliver & Boyd, Edinburgh. Gossmann, T.I., Keightley, P.D., Eyre-Walker, A., 2012. The effect of variation in the effective population size on the rate of adaptive molecular evolution in eukaryotes. Genome Biol. Evol. 4, 658–667. Greenwood, P.H., 1979. Towards a phyletic classification of the ‘genus’ Haplochromis (Pisces, Cichlidae) and related taxa. Bull. Br. Mus. Nat. Hist. (Zool.) 35, 265–322. Guindon, S., Dufayard, J.F., Lefort, V., Anisimova, M., Hordijk, W., Gascuel, O., 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59, 307–321. Hanada, K., Gojobori, T., Li, W.H., 2006. Radical amino acid change versus positive selection in the evolution of viral envelope proteins. Gene 385, 83–88. Hanada, K., Shiu, S.H., Li, W.H., 2007. The nonsynonymous/synonymous substitution rate ratio versus the radical/conservative replacement rate ratio in the evolution of mammalian genes. Mol. Biol. Evol. 24, 2235–2241.

He, A., Luo, Y., Yang, H., Liu, L., Li, S., Wang, C., 2010. Complete mitochondrial DNA sequences of the Nile tilapia (Oreochromis niloticus) and Blue tilapia (Oreochromis aureus): genome characterization and phylogeny applications. Mol. Biol. Rep. http://dx. doi.org/10.1007/s11033-010-0324-7. Hughes, A.L., Ota, T., Nei, M., 1990. Positive Darwinian selection promotes charge profile diversity in the antigen-binding cleft of class I major-histocompatibility-complex molecules. Mol. Biol. Evol. 7, 515–524. Johnson, K.P., Seger, J., 2001. Elevated rates of nonsynonymous substitution in island birds. Mol. Biol. Evol. 18, 874–881. Johnson, T.C., Scholz, C.A., Talbot, M.R., Kelts, K., Ricketts, R.D., Ngobi, G., et al., 1996. Late Pleistocene desication of Lake Victoria and rapid evolution of cichlid fishes. Science 273, 1091–1093. Kocher, T.D., Thomas, W.K., Meyer, A., Edwards, S.V., Paabo, S., Villablanca, F.X., et al., 1989. Dynamics of mitochondrial DNA evolution in animals: amplification and sequencing with conserved primers. Proc. Natl. Acad. Sci. U. S. A. 86, 6196–6200. Kocher, T.D., Conroy, J.A., Mckaye, K.R., Stauffer, J.R., Lockwood, S.F., 1995. Evolution of NADH dehydrogenase subunit 2 in East African cichlid fishes. Mol. Phylogenet. Evol. 4, 420–432. Lanfear, R., Kokko, H., Eyre-Walker, A., 2014. Population size and the rate of evolution. Trends Genet. 29, 33–41. Lezzar, K.E., Tiercelin, J.J., Batist, M., Cohen, A.S., Bandora, T., Rensbergen, P., et al., 1996. New seismic stratigraphy and Late Tertiary history of the North Tanganyika Basin, East African Rift system, deduced from multichannel and high-resolution reflection seismic data and piston core evidence. Basin Res. 8, 1–28. Little, A.G., Lougheed, S.C., Moyes, C.D., 2012. Evolution of mitochondrial-encoded cytochrome oxidase subunits in endothermic fish: the importance of taxon-sampling in codon-based models. Mol. Phylogenet. Evol. 63, 679–684. Lourenço, J.M., Glémin, S., Galtier, N., 2013. The rate of molecular adaptation in a changing environment. Mol. Biol. Evol. 30, 1292–1301. Mabuchi, K., Miya, M., Azuma, Y., Nishida, M., 2007. Independent evolution of the specialized pharyngeal jaw apparatus in cichlid and labrid fishes. BMC Evol. Biol. 7, 10. Meiklejohn, C.D., Montooth, K.L., Rand, D.M., 2007. Positive and negative selection on the mitochondrial genome. Trends Genet. 23, 259–263. Meyer, A., Kocher, T.D., Basasibwaki, P., Wilson, A.C., 1990. Monophyletic origin of Lake Victoria cichlid fishes suggested by mitochondrial DNA sequences. Nature 347, 550–553. Meyer, A., Kocher, T.D., Wilson, A.C., 1991. African fishes. Nature 350, 467–468. Nagl, S., Tichy, H., Mayer, W.E., Takahata, N., Klein, J., 1998. Persistence of neutral polymorphisms in Lake Victoria cichlid fish. Proc. Natl. Acad. Sci. U. S. A. 95, 14238–14243. Nagl, S., Tichy, H., Mayer, W.E., Takezaki, N., Takahata, N., Klein, J., 2000. The origin and age of the haplochromine species flock in Lake Victoria, East Africa. Proc. R. Soc. Lond. Ser. B 267, 1049–1061. Nishida, M., 1991. Lake Tanganyika as an evolutionary reservoir of old lineages of East African fishes: inferences from allozyme data. Experientia 47, 974–979. Nozawa, M., Suzuki, Y., Nei, M., 2009. Reliabilities of identifying positive selection by the branch-site and the site-prediction methods. Proc. Natl. Acad. Sci. U. S. A. 106, 6700–6705. Ohta, T., 1973. Slightly deleterious mutant substitutions in evolution. Nature 246, 96–98. Ohta, T., 1995. Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral theory. J. Mol. Evol. 40, 56–63. Ohta, T., 2002. Near-neutrality in evolution of genes and gene regulation. Proc. Natl. Acad. Sci. U. S. A. 99, 16134–16137. Oliveira, D.C.S.G., Raychoudhury, R., Lavrov, D.V., Werren, J.H., 2008. Rapidly evolving mitochondrial genome and directional selection in mitochondrial genes in the parasitic wasp Nasonia (Hymenoptera: Pteromalidae). Mol. Biol. Evol. 25, 2167–2180. Osada, N., Akashi, H., 2012. Mitochondrial-nuclear interactions and accelerated compensatory evolution: evidence from the primate cytochrome c oxidase complex. Mol. Biol. Evol. 29, 337–346. Pfeffer, G.J., 1896. Die Fische Ost-Afrikas. In: Mæbius, F. (Ed.), Deutsch-Ost-Afrika, Band III. Die Thierwelt Ost-Afrikas und der Nachbargebiete. Wirbelthiere 4, pp. i–xviii (+1-72, Separatum: 1-72). Poll, M., 1986. Classification des Cichlidae du lac Tanganyika, Tribus, genres et espéces. Brussels (Belgium). Acad R. Belg. Mem. Cl. Sci. 45, 1–163. Popadin, K., Polishchuk, L.V., Mamirova, L., Knorre, D., Gunbin, K., 2007. Accumulation of slightly deleterious mutations in mitochondrial protein-coding genes of large versus small mammals. Proc. Natl. Acad. Sci. U. S. A. 104, 13390–13395. Rand, D.M., Haney, R.A., Fry, A.J., 2004. Cytonuclear coevolution: the genomics of cooperation. Trends Ecol. Evol. 19, 645–653. Salzburger, W., Meyer, A., Baric, S., Verheyen, E., Sturmbauer, C., 2002. Phylogeny of the Lake Tanganyika cichlid species flock and its relationship to the Central and East African haplochromine cichlid fish faunas. Syst. Biol. 51, 113–135. Salzburger, W., Mack, T., Verheyen, E., Meyer, A., 2005. Out of Tanganyika: genesis, explosive speciation, key-innovations and phylogeography of the haplochromine cichlid fishes. BMC Evol. Biol. 5, 17. Samonte, I.E., Satta, Y., Sato, A., Tichy, H., Takahata, N., Klein, J., 2007. Gene flow between species of Lake Victoria haplochromine fishes. Mol. Biol. Evol. 24, 2069–2080. Scholz, C.A., Johnson, T.C., Cohen, A.S., King, J.W., Peck, J.A., Overpeck, J.T., et al., 2007. East African megadroughts between 135 and 75 thousand years ago and bearing on earlymodern human origins. Proc. Natl. Acad. Sci. U. S. A. 104, 16416–16421. Seehausen, O., 1996. Lake Victoria Rock Cichlids: Taxonomy, Ecology, and Distribution. Verduijn Cichlids, Zevenhuizen. Seehausen, O., Koetsier, E., Schneider, M.V., Chapman, L.J., Chapman, C.A., Knight, M.E., et al., 2003. Nuclear markers reveal unexpected genetic variation and a CongoleseNilotic origin of the Lake Victoria cichlid species flock. Proc. R. Soc. Lond. Ser. B 270, 129–137.

K. Shirai et al. / Gene 552 (2014) 239–245 Shaw, P.W., Turner, G.F., Idid, M.R., Robinson, R.L., Carvalho, G.R., 2000. Genetic population structure indicates sympatric speciation of Lake Malawi pelagic cichlids. Proc. R. Soc. Lond. Ser. B 267, 2273–2280. Shen, Y.Y., Liang, L., Zhu, Z.H., Zhou, W.P., Irwin, D.M., Zhang, Y.P., 2010. Adaptive evolution of energy metabolism genes and the origin of flight in bats. Proc. Natl. Acad. Sci. U. S. A. 107, 8666–8671. Smith, N.G.C., 2003. Are radical and conservative substitution rates useful statistics in molecular evolution? J. Mol. Evol. 57, 467–478. Smith, J.M., Smith, N.H., 1996. Synonymous nucleotide divergence: what is “saturation”? Genetics 142, 1033–1036. Stager, J.C., Mayewski, P.A., Meeker, L.D., 2002. Cooling cycles, Heinrich event 1, and the desiccation of Lake Victoria. Palaeogeogr. Palaeoclimatol. Palaeoecol. 183, 169–178. Stager, J.C., Ryves, D.B., Chase, B.M., Pausata, F.S., 2011. Catastrophic drought in the Afro-Asian monsoon region during Heinrich event 1. Science 331, 1299–1302. Takeda, M., Kusumi, J., Mizoiri, S., Aibara, M., Mzighani, S.I., Sato, T., et al., 2013. Genetic structure of pelagic and littoral cichlid fishes from Lake Victoria. PLoS One 8, e74088. Terai, Y., Takezaki, N., Mayer, W.E., Tichy, H., Takahata, N., Klein, J., et al., 2004. Phylogenetic relationships among East African haplochromine fishes as revealed by short inter spersed elements (SINEs). J. Mol. Evol. 58, 64–78. Turner, G.F., Seehausen, O., Knight, K.E., Allender, C.J., Robinson, R.L., 2001. How many species of cichlid fishes are there in African lakes? Mol. Ecol. 10, 793–806.

245

Verheyen, E., Salzburger, W., Snoeks, J., Meyer, A., 2003. Origin of the super-flock of cichlid fishes from Lake Victoria, East Africa. Science 300, 325–329. Wagner, C.E., Keller, I., Wittwer, S., Selz, O.M., Mwaiko, S., Greuter, L., et al., 2013. Genomewide RAD sequence data provide unprecedented resolution of species boundaries and relationships in the Lake Victoria cichlid adaptive radiation. Mol. Ecol. 22, 787–798. Wernegreen, J.J., 2011. Reduced selective constraint in endosymbionts: elevation in radical amino acid replacements occurs genome-wide. PLoS One 6, e28905. Yang, Z., 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591. Yang, Z., Nielsen, R., 1998. Synonymous and nonsynonymous rate variation in nuclear genes of mammals. J. Mol. Evol. 46, 409–418. Yang, Z., Nielsen, R., 2002. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol. Biol. Evol. 19, 908–917. Yang, Z., Wong, W.S.W., Nielsen, R., 2005. Bayes empirical Bayes inference of amino acid sites under positive selection. Mol. Biol. Evol. 22, 1107–1118. Zhang, J., 2000. Rates of conservative and radical nonsynonymous nucleotide substitutions in mammalian nuclear genes. J. Mol. Evol. 50, 56–68. Zhang, J., Nielsen, R., Yang, Z., 2005. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol. Biol. Evol. 22, 2472–2479.