Accelerated diversification of nonhuman primate malarias in Southeast Asia: adaptive radiation or geographic speciation?

Although parasitic organisms are found worldwide, the relative importance of host specificity and geographic isolation for parasite speciation has been explored in only a few systems. Here, we study Plasmodium parasites known to infect Asian nonhuman primates, a monophyletic group that includes the lineage leading to the human parasite Plasmodium vivax and several species used as laboratory models in malaria research. We analyze the available data together with new samples from three sympatric primate species from Borneo: The Bornean orangutan and the long-tailed and the pig-tailed macaques. We find several species of malaria parasites, including three putatively new species in this biodiversity hotspot. Among those newly discovered lineages, we report two sympatric parasites in orangutans. We find no differences in the sets of malaria species infecting each macaque species indicating that these species show no host specificity. Finally, phylogenetic analysis of these data suggests that the malaria parasites infecting Southeast Asian macaques and their relatives are speciating three to four times more rapidly than those with other mammalian hosts such as lemurs and African apes. We estimate that these events took place in approximately a 3-4-Ma period. Based on the genetic and phenotypic diversity of the macaque malarias, we hypothesize that the diversification of this group of parasites has been facilitated by the diversity, geographic distributions, and demographic histories of their primate hosts.

Introduction

Parasitism is one of the most common multispecies interactions observed in nature (Poulin 2005). Most research programs in evolutionary parasitology focus on long-term coevolutionary processes where the adaptation of parasites to their hosts is assumed to have an essential role in parasite speciation (Page 2003; Poulin et al. 2011). A complementary approach considers that, in addition to adaptation to a host or hosts, geographic factors come into play when explaining the diversity of parasite species (Huyse et al. 2005; Brooks et al. 2006; Hoberg et al. 2008; Ricklefs 2010). This second perspective is particularly valuable wherever there are not obvious phylogenetic concordances between host species and their parasites.

Plasmodium (phylum Apicomplexa, family Plasmodiidae) is a group of vector-borne parasitic protozoa that are better known by the handful of species that cause malaria in humans. The global health importance of human malarias, however, overshadows the extraordinary diversity in this genus represented by more than 200 described species parasitizing many vertebrate hosts including rodents and nonhuman primates (Garnham 1966; Coatney et al. 1971; Valkiūnas 2005; Telford 2009). Of particular interest are those Plasmodium found in nonhuman primates and their relationships to the extant human malarias (Ayala et al. 1998; Escalante et al. 1998; Singh et al. 2004; Escalante et al. 2005; Krief et al. 2010; Liu et al. 2010; Prugnolle et al. 2013; Pacheco et al. 2013). Less attention, however, has been given to the evolutionary processes leading to the staggering diversity observed in primate malarial parasites (Hayakawa et al. 2008; Pacheco et al. 2011; Pacheco et al. 2012b).

Based on a rate of speciation that was considered rapid for Plasmodium, it was proposed that such a “big-bang” pattern of parasite diversity was evidence for an adaptive radiation of the malarias infecting mammalian hosts facilitated by multiple host switches (Hayakawa et al. 2008). Although the hypothesis of adaptive radiation is appealing for parasites, such processes are not easy to study in any group of organisms, parasitic, or free-living. Specifically, showing rapid diversification as result of putative adaptations to a newly available ecological opportunity is complex, especially in a diverse and geographically widespread group of organisms (Hunter 1998; Glor 2010; Aguilée et al. 2013). Furthermore, apparent patterns of accelerated evolution could be consistent with many other processes, if they are observed at broad time and geographic scales (Glor 2010).

One group of malarial parasites that is especially well suited for exploring processes such as adaptive radiation and/or geographic speciation is composed of the parasites that infect nonhuman primates in South East Asia. These species are part of a monophyletic group that exhibits extraordinary phenotypic diversity, and its radiation is connected with the origin of the human parasite Plasmodium vivax (Escalante et al. 1998, 2005; Pacheco et al. 2011; Prugnolle et al. 2013). This group includes species that are zoonotic (Singh et al. 2004; Ta et al. 2014), and others that are considered model organisms in malaria research (Coatney et al. 1971; Cox-Singh and Singh 2008). In addition, this clade is found in a defined geographic region that shares a common geologic history with many parasite species exhibiting overlapping host ranges (the host species utilized by a parasite species) and spatial distributions (Coatney et al. 1971; Fooden 1994).

Here, we report new data on the diversity of Plasmodium species infecting macaques and orangutans at two study sites in Borneo and analyze them together with previously reported data (Hayakawa et al. 2008; Pacheco et al. 2011; Pacheco et al. 2012b). We find several lineages that suggest new parasite species in this biodiversity hotspot. Among those, we report divergent sympatric lineages in orangutans as well as evidence of parasite population genetic structures determined by geographic isolation. We also present evidence suggesting that the diversification rate of the Southeast Asian parasites is accelerated relative to other Plasmodium clades (e.g., the lemur parasites or the Laverania clade that infects African apes and humans).

Thus, taking into consideration the apparent accelerated rate of evolution in this monophyletic group and that many species share multiple hosts, we hypothesize that Plasmodium species diversity in Southeast Asia was driven, at least in part, by the history of the host species utilized by these parasites. However, determining whether such rapid diversification was adaptive in nature will require additional evidence beyond the observation of accelerated evolution.

Results

In this study, we characterized the diversity of nonhuman primate malarias found in blood sampled from one population of Pongo pygmaeus morio and two populations of Macaca nemestrina and Macaca fascicularis from Sabah, Malaysian Borneo.

Plasmodium Species Diversity and Phylogenetic Relationship

Table 1 describes the diversity of nonhuman primate malarias found in one population of orangutans (Po. p. morio) and two populations from different macaque species (M. nemestrina and M. fascicularis) present in Sabah, Malaysian Borneo. Blood smears were collected and some of the samples were positive by microscopy; however, the quality of the smears did not allow us to identify the parasites beyond the genus level. Thus, we relied on molecular diagnostics for species identification.

Table 1.

Diversity of Plasmodium Parasites in Primates from Sabah, Borneo (complete mtDNA genomes).

Host Species	Parasite Species	No. Positive Samples
Macaca nemestrina (Sundaland pig-tailed macaque)	Plasmodium knowlesi	2/15
	P. coatneyi	1/15
	Plasmodium sp.	1/15
	P. fieldi	3/15a
	P. inui	9/15a
	P. cynomolgi	1/15
Macaca fascicularis (long-tailed macaque)	P. knowlesi	4/26a
	P. coatneyi	1/26
	Plasmodium sp.	1/26
	P. fieldi	1/26
	P. inui	8/26a
	P. cynomolgi	3/26a
	Hepatocystis sp.	2/26a
Pongo pygmaeus (Bornean orangutan)	Plasmodium sp. (clade A)b	7/38a
	Plasmodium sp. (clade B)b	13/38a
	Plasmodium sp. (clade C)b	10/38a

aMixed infection.

bFor clade information see figure 1.

F

Phylogenetic tree of Plasmodium species based on complete mitochondrial genomes. Bayesian and ML methods yielded identical topologies and so only the Bayesian tree is shown. Both phylogenetic methods used a general time reversible model with gamma-distributed substitution rates and a proportion of invariant sites (GTR + Γ + I). The values above branches are posterior probabilities together with bootstrap values (in bold) as a percentage obtained for the ML tree (see Materials and Methods). Macaque parasites from Borneo are labeled in blue and Plasmodium species from orangutans are labeled in red. Numbers of individuals with each parasite are shown. Samples with * indicate mixed infection. The outgroup is indicated by the gray branches.

All blood samples from 15 pig-tailed macaques (M. nemestrina) and 26 long-tailed macaques (M. fascicularis) were positive for malarial parasites by polymerase chain reaction (PCR) using cytochrome b (cytb), and some of them had mixed infections. We were able to generate 34 parasite mitochondrial genomes (mtDNA) for different Plasmodium parasites, 16 from M. nemestrina and 18 from M. fascicularis. However, we were unable to obtain complete mtDNA sequences for several positive M. fascicularis samples, probably because of low parasitemia. Plasmodium inui was the most common parasite in these two macaque species (9/15 for M. nemestrina and 8/26 for M. fascicularis, table 1). Regardless of the relatively small sample size, we found the same six Plasmodium species (P. knowlesi, P. coatneyi, Plasmodium sp., P. fieldi, P. inui, and P. cynomolgi) in both species of Macaca (table 1).

Figure 1 shows a Bayesian phylogenetic analysis of a combined set of nearly complete mitochondrial genome sequences including those generated in this and several previous studies (supplementary table S1, Supplementary Material online). We used the African parasite P. gonderi and one found in drills in order to root the diversity of species observed in this study (the six Plasmodium species listed in table 1) together with other parasites already reported in nonhuman primates from southeast Asia (Escalante et al. 2005; Pacheco et al. 2013; Pacheco et al. 2012b). Overall, this phylogeny is similar to those obtained in previous studies (Pacheco et al. 2011, 2013; Pacheco et al. 2012b), but includes several lineages that may correspond to new Plasmodium species. Of particular interest is a lineage that was present in both macaque species which shares a recent common ancestor with P. coatneyi. This Plasmodium species is divergent from the other species in that clade and does not appear to belong to either P. knowlesi or P. coatneyi (table 2). Specifically, the divergence between Plasmodium sp. and P. coatneyi (0.023 ± 0.002, table 2) was comparable to that observed between other pairs of Plasmodium species, including the human parasite P. vivax and the macaque parasite P. cynomolgi (0.012 ± 0.002, table 2), and the chimpanzee parasite P. reichenowi and the human parasite P. falciparum (0.025 ± 0.002, table 2). This result may indicate that both species of macaques harbor a new Plasmodium species. In addition, we found Hepatocystis sp. in two individuals of M. fascicularis but not in M. nemestrina (table 1, sequences not reported).

Table 2.

Genetic Divergences (substitutions per site) among Different Plasmodium Species Using the Kimura 2-Parameter Model As Implemented in MEGA v5.2.2 (Tamura et al. 2011).

Species	Genetic distance (d ± SE)
	n	COX1	COX3	CYTB	COX1 + CYTB	Complete mtDNA
Plasmodium sp. (orangutan-clade A)	12	0.005 ± 0.001	0.006 ± 0.002	0.003 ± 0.001	0.004 ± 0.001	0.004 ± 0.0003
Plasmodium sp. (orangutan-clade B)	13	0.003 ± 0.001	0.004 ± 0.001	0.004 ± 0.001	0.003 ± 0.001	0.003 ± 0.0003
Plasmodium sp. (orangutan-clade C)	10	0.002 ± 0.001	0.003 ± 0.001	0.004 ± 0.001	0.003 ± 0.001	0.003 ± 0.0003
Orangutan clade A − orangutan clade B	12 vs. 13	0.011 ± 0.002	0.019 ± 0.005	0.008 ± 0.002	0.010 ± 0.001	0.009 ± 0.0010
Orangutan clade A − orangutan clade C	12 vs. 10	0.034 ± 0.004	0.043 ± 0.007	0.028 ± 0.004	0.031 ± 0.003	0.023 ± 0.0017
Orangutan clade B − orangutan clade C	13 vs. 10	0.034 ± 0.004	0.049 ± 0.008	0.028 ± 0.005	0.031 ± 0.003	0.024 ± 0.0018
Orangutan clade A + orangutan clade B	25	0.008 ± 0.001	0.012 ± 0.003	0.006 ± 0.001	0.007 ± 0.001	0.006 ± 0.0006
P. cynomolgi (all strains)	14	0.005 ± 0.001	0.006 ± 0.001	0.005 ± 0.001	0.004 ± 0.001	0.004 ± 0.0006
P. cynomolgi (clade A-Berok)	4	0.007 ± 0.002	0.007 ± 0.002	0.006 ± 0.002	0.006 ± 0.001	0.006 ± 0.0008
P. cynomolgi (clade B-B)	10	0.001 ± 0.001	0.001 ± 0.001	0.001 ± 0.001	0.001 ± 0.000	0.001 ± 0.0002
P. cynomolgi (clade A vs. clade B)	4 vs.10	0.010 ± 0.002	0.010 ± 0.003	0.009 ± 0.002	0.009 ± 0.002	0.009 ± 0.0011
P. vivax	109	0.001 ± 0.000	0.001 ± 0.001	0	0.001 ± 0.000	0.001 ± 0.0002
P. vivax–P. cynomolgi	109 vs. 14	0.014 ± 0.003	0.017 ± 0.004	0.013 ± 0.003	0.014 ± 0.002	0.012 ± 0.0015
P. inui (all)	31	0.013 ± 0.002	0.023 ± 0.004	0.015 ± 0.002	0.014 ± 0.001	0.013 ± 0.0010
P. inui (Macaca nemestrina)-(M. fascicularis)	9 vs. 8	0.004 ± 0.001	0.005 ± 0.001	0.003 ± 0.001	0.003 ± 0.000	0.004 ± 0.0000
P. inui (all)–P. sp. from orang clade A	31 vs. 12	0.023 ± 0.003	0.045 ± 0.007	0.018 ± 0.002	0.024 ± 0.003	0.023 ± 0.0010
P. inui (all)–P. sp. from orang clade B	31 vs. 13	0.025 ± 0.004	0.054 ± 0.008	0.027 ± 0.004	0.026 ± 0.002	0.025 ± 0.0020
P. inui (all)–P. sp. from orang clade A + B	31 vs. 25	0.024 ± 0.003	0.050 ± 0.007	0.026 ± 0.004	0.025 ± 0.002	0.024 ± 0.0020
P. knowlesi (all)	65	0.001 ± 0.000	0.002 ± 0.001	0.001 ± 0.000	0.001 ± 0.000	0.001 ± 0.0000
P. knowlesi (NCBI)–Sepilok	60 vs. 5	0.002 ± 0.000	0.002 ± 0.001	0.001 ± 0.000	0.001 ± 0.000	0.001 ± 0.0000
P. knowlesi–P. coatneyi	65 vs. 2	0.043 ± 0.005	0.034 ± 0.007	0.025 ± 0.005	0.035 ± 0.004	0.032 ± 0.0020
Plasmodium sp. (Mne–Mfa)	1 vs. 1	0	0	0.003 ± 0.001	0.001 ± 0.001	0.003 ± 0.0010
P. knowlesi–Plasmodium sp.	65 vs. 2	0.037 ± 0.005	0.037 ± 0.007	0.020 ± 0.004	0.030 ± 0.003	0.029 ± 0.0020
P. coatneyi–Plasmodium sp.	2 vs. 2	0.032 ± 0.005	0.037 ± 0.007	0.019 ± 0.004	0.026 ± 0.003	0.023 ± 0.0020
P. falciparum	101	0	0.000 ± 0.001	0.001 ± 0.001	0.001 ± 0.000	0.000 ± 0.0001
P. reichenowi–P. falciparum	2 vs. 101	0.032 ± 0.005	0.055 ± 0.009	0.025 ± 0.005	0.029 ± 0.003	0.025 ± 0.0018

Of the 38 orangutans sampled for this study, 23 were positive (60.5%) by PCR and 10 harbored mixed infections (table 1). Three samples had low parasitemia and the PCR did not yield bands that would allow us to clone them. Unfortunately, although positive blood smears were available, these were not of sufficient quality to determine whether these parasites coincide with any of the Plasmodium species previously described in orangutans by Garnham (1966) and Peters et al. (1976). We generated 30 parasite mitochondrial genomes from the remaining 20 samples; these haplotypes belonged to three different lineages (A–C, fig. 1). The phylogenetic analysis showed that two of these lineages (A and B) form a monophyletic group and share a common ancestor with P. inui, a species of macaque parasite with broad geographic distribution (Coatney et al. 1971; Fooden 1994).

Lineage C is at the base of a monophyletic group that includes P. hylobati from gibbons, P. inui from macaques, as well as lineages A and B from orangutans. This parasite phylogeny is not congruent with any primate phylogenetic hypotheses estimated using nuclear and/or mitochondrial genes (e.g., Perelman et al. 2011); such discrepancy indicates that host switches likely took place during the evolution of this Plasmodium clade. Specifically, P. inui seems to derive from parasites infecting Asian apes (Pacheco et al. 2012b). The average distance between the haplotypes belonging to lineages A and B in orangutan malarias (0.009 ± 0.001 for the complete mtDNA, table 2) is similar to the divergence between P. vivax and P. cynomolgi (0.012 ± 0.002, table 2). On the other hand, the average distances between clades A and C (0.023 ± 0.002) and between clades B and C (0.024 ± 0.002) are somewhat higher than those between clades A and B, between P. vivax and P. cynomolgi, and between P. reichenowi and P. falciparum. These results suggest that the haplotypes identified in figure 1 as A–C may belong to three different species. Clade A includes haplotypes described previously in a different population of orangutans (Pacheco et al. 2012b).

In the case of P. cynomolgi, we found two closely related lineages in both species of macaques. These lineages belong to one of the two clades (fig. 1) that have been described for P. cynomolgi using antigen-encoding genes (Pacheco et al. 2007, 2010; Pacheco et al. 2012a; Pacheco et al. 2012b), which contains the Berok, PT2 and Gombak strains. Indeed, the average genetic distance among these four P. cynomolgi strains (0.006 ± 0.0008) is relatively small if we consider its broad geographic range (see table 2). It is also worth noting that the divergence of these two clades of P. cynomolgi is comparable to that observed between P. vivax and P. cynomolgi (table 2), an observation that has been corroborated by genomic data (Tachibana et al. 2012).

Malaria Antigens: Apical Membrane Antigen-1 and Merozoite Surface Protein 1 42 kDa

In addition to mtDNA genomes, phylogenetic analyses were also performed separately on data from several macaque and orangutan parasites using nuclear genes encoding two major antigens: Apical membrane antigen-1 (AMA-1; fig. 2 and supplementary table S1, Supplementary Material online) and merozoite surface protein 1 42 kDa (MSP-142; fig. 3 and supplementary table S1, Supplementary Material online). Despite the smaller number of species for which data is available, the phylogenies estimated from these two loci are largely compatible with previous studies using genes encoding merozoite antigens (MSP-1: Pacheco et al. 2007; Pacheco et al. 2012b; RAP-1: Pacheco et al. 2010; MSP-8/MSP-10: Pacheco et al. 2012a). In particular, in both phylogenies we observe that several of the orangutan malaria lineages form a monophyletic group sharing a common ancestor with the macaque parasite, P. inui. One difference between the phylogenies generated using either mtDNA or these two antigen-encoding genes is in the relative position of P. hylobati. This is not surprising considering that the position of P. hylobati has modest support in the mtDNA (0.85 posterior probability and 63% bootstrap support, see fig 1). Although in the mitochondrial tree this lineage is located within the clade that includes the three orangutan lineages and P. inui, in the ama-1 and msp-142 phylogenies it is located at the base of the clade containing P. inui and the three orangutan parasite clades (A-C, fig. 1). This discrepancy could have three different explanations. First, because the mitochondrial and nuclear loci are unlinked, there could be random differences in the genealogies of these loci. Second, this could be due to diversifying- or frequency-dependent selection acting on the antigen loci. Finally, this could be an artifact of our inability to generate antigen sequences for all of the taxa included in the mitochondrial tree. Indeed, because we were unable to generate any sequences of ama-1 that can be unambiguously assigned to mtDNA clade C, the position of P. hylobati relative to the orangutan clades is unresolved (fig. 2).

F

Bayesian phylogenetic tree of Plasmodium species based on the gene encoding AMA-1. Bayesian and ML methods yielded identical topologies and so only the Bayesian tree is shown. We used a general time reversible model with gamma-distributed substitution rates (GTR + Γ). The values above branches are posterior probabilities together with bootstrap values (in bold) as a percentage obtained for the ML tree (see Materials and Methods). Macaque parasites from Borneo are labeled in blue and Plasmodium species from orangutans are labeled in red. Clone numbers are shown. Clade numbers are shown and not necessarily correspond to clades A–C (fig. 1). For more details see Result section. Samples with * indicate mixed infection. The outgroup is indicated by the gray branch.

F

Bayesian phylogenetic tree of Plasmodium species based on the gene encoding MSP-142 kDa. Bayesian and ML methods yielded identical topologies and so only the Bayesian tree is shown. We used a general time reversible model with gamma-distributed substitution rates (GTR + Γ). The values above branches are posterior probabilities together with bootstrap values (in bold) as a percentage obtained for the ML tree (see Materials and Methods). Macaque parasites from Borneo are labeled in blue and Plasmodium species from orangutans are labeled in red. Clone numbers are shown. Clade numbers are shown and not necessarily correspond to clades A–C (fig. 1). For more details see Result section. The outgroup is indicated by the gray branch.

In the case of ama-1, we were able to obtain sequences for P. knowlesi, Plasmodium sp., P. cynomolgi, P. inui from M. nemestrina, and Plasmodium sp. from orangutan. The genetic polymorphism of ama-1 per species is reported in the supplementary table S2, Supplementary Material online. The ama-1 allele found in this study for P. cynomolgi was most closely related to that contained in the Berok strain; this result is consistent with the mtDNA phylogeny. Again, the undescribed Plasmodium sp. isolated from macaques forms a clade along with P. coatneyi, and the average genetic distance estimated between the ama-1 allele for this parasite and P. coatneyi is very high (0.076 ± 0.073), supporting the hypothesis that the former organism is a new species. On the other hand, we only were able to generate msp-142 sequences for P. inui and the Plasmodium species detected in orangutans. The genetic polymorphism per species for msp-142 is reported in supplementary table S3, Supplementary Material online. Overall, using these additional msp-142 sequences, we recovered the same topology that was previously reported for this antigen (Pacheco et al. 2007; Pacheco et al. 2012b).

Phylogenetic-based analyses provide some insight into the role of selection at these two loci (Pond et al. 2011). In the case of ama-1, several lineages contain codons that appear to be under positive selection (ω >1; P < 0.05 corrected for multiple testing using Holm–Bonferroni method; Supplementary table S4 and fig. S1, Supplementary Material online). Furthermore, many of the nodes that may correspond to host switches also have significant uncorrected sequential likelihood ratio tests (LRTs). A similar pattern emerges from msp-142 where this method identifies four lineages that appear to be under episodic diversifying selection (supplementary table S4 and fig. S1, Supplementary Material online).

Haplotype Networks

A median joining network of 41 Plasmodium mtDNA haplotypes isolated from orangutans from Borneo is depicted in figure 4. This network includes both the 30 Plasmodium mtDNA haplotypes from Sabah and the 11 mtDNA haplotypes from the Kalimantan orangutans. As expected, the three clades identified in the mtDNA tree were clearly separated in the haplotype network. Furthermore, clade A itself (colored in green) contains two genetically differentiated populations (light green: Sabah; dark green: Kalimantan) and indeed the fixation index (FST) between these two populations is high (0.56). In another context, the Venn diagram depicted in figure 4 shows the number of orangutans with mixed infections of Plasmodium. At least 10 of the 20 orangutan samples harbored mixed infections with parasites from two distinct clades, but none was found to be simultaneously infected by all three species.

F

Venn diagram of mixed infections and median joining network of Plasmodium sp. mtDNA haplotypes from orangutans from Borneo. Branch lengths are proportional to divergence and node sizes are proportional to total haplotype frequencies. Colors correspond to different clades: green = clade A; red = clade B; and blue = clade C.

Figure 5 shows a median joining network for the 16 P. inui mtDNA haplotypes from both species of Macaca (eight haplotypes from M. nemestrina and eight from M. fascicularis). No evidence of population structure between hosts was found in this analysis; in fact, the fixation index (FST) is effectively 0 (−0.024). In addition, we also identified a divergent P. inui mtDNA haplotype that was carried by a single individual of M. nemestrina (see also fig. 1).

F

Median joining network of Plasmodium inui mtDNA haplotypes from macaques (Macaca nemestrina and M. fascicularis) from Sabah, Borneo. Branch lengths are proportional to divergence and node sizes are proportional to total haplotype frequencies. Colors correspond to different species of Macaca: blue = M. nemestrina and maroon = M. fascicularis.

The Ages and Branching Rates of the Asian Primate Malarias

The mitochondrial genome phylogeny was used to estimate the ages of the major clades of South Asian primate malarias under three different calibration scenarios (fig. 6 and supplementary fig. S2, Supplementary Material online). These times, along with the 95% credibility intervals (Crls), were calculated using BEAST v1.7.5 (Drummond et al. 2012) and are shown in table 3. (See supplementary fig. S2, Supplementary Material online, for the node numbers.) As expected, there is extensive overlap among the Crls obtained under the three different scenarios used for calibration. When the phylogeny is calibrated by assuming that the Asian and African primate parasites split in the same time frame that their hosts became geographically isolated (6–14.2 Ma) and that lemur parasites originated as a monophyletic group at least 20 Ma (see Materials and Methods), the Asian primate malarias are estimated to have originated between 5.38 and 8.52 Ma. The split between P. inui and the orangutan Plasmodium lineages in clades A and B was around 2.92 and 4.80 Ma, which is in agreement with an earlier study of clade A (2.46–4.38 Ma, Pacheco et al. 2012b). Furthermore, clades A and B diverged approximately 1.0–2.3 Ma (table 3); this is clearly older than time estimates for the split of the Bornean and Sumatran orangutans using nuclear single nucleotide polymorphisms (SNPs) but earlier than estimates for the same two species using the mitochondrial genome (Ma et al. 2013). On the other hand, the common ancestor of the orangutan Plasmodium lineages in clade C appears to have existed between 3.71 and 5.98 Ma, well after the time of the split between orangutan and gibbons, or the most recent common ancestor of Asian apes with macaques.

F

Timetree of the divergence of malarial parasites from Borneo, Sabah. Divergence times were estimated with BEAST using the most inclusive scenario based on the minimum divergence of Macaca and Papio using fossils (6–14.2 Ma) and a minimum of 20 Ma for the origin of lemur parasites. Times are shown in Ma. Macaque parasites from Borneo are labeled in blue and Plasmodium species from orangutans are labeled in red. The numbers of the nodes described in table 3 are provided and depicted with an orange dot.

Table 3.

Divergence Times of Major Splits in the Malaria Phylogeny As Estimated by Beast.

Calibrations (Ma)		Node 13: Min = 6, Max = 14.2; Node 19: Min 20		Node 13: Min = 6, max = 14.2; Node 14: Min = 23.5; Node 19: Min 20		Node 13: Min = 6, Max = 8; Node 19: Min 20
Divergence	Node	Node Age	95%CrI	Node Age	95%CrI	Node Age	95%CrI
Split P. vivax–P. cynomolgi	3	4.2	2.76–5.75	4.43	2.97–5.96	3.68	2.50–4.89
Origin of P. inui	6	2.92	2.20–3.77	3.08	2.33–3.90	2.55	2.01–3.21
Origin of Plasmodium sp. from orangutans (clade A and Ba)	8	1.55	1.00–2.24	1.64	1.05–2.32	1.37	0.91–1.94
Split P. inui–Plasmodium spp. from orangutans	9	3.78	2.92–4.80	3.98	3.10–4.97	3.31	2.66–4.06
Split Plasmodium spp. (clade Ca)–(P. inui–P. spp.–P. hylobati)	11	4.73	3.71–5.98	4.97	3.90–6.12	4.13	3.36–5.03
Split (P. vivax-P. cynomolgi)–(P. inui–P. spp.–P. hylobati)	12	5.65	4.47–7.10	5.94	4.73–7.28	4.91	4.03–5.92
Split P. knowlesi–P. coatneyi	17	5.54	4.17–7.10	5.80	4.50–7.37	4.77	3.70–5.90
Origin of Southern Asia primates malaria	19	6.84	5.38–8.52	7.18	5.78–8.78	5.88	4.82–7.09
Split of Plasmodium spp. from Papio–Macaca divergence	21	11.57	9.33–14.28	12.18	9.99–14.30	8.97	7.24–11.03
Radiation of rodent malarias	33	13.83	10.37–18.21	14.62	10.96–18.84	12.38	9.57–15.69
Radiation of African ape malarias	38	16.81	12.55–22.17	17.70	13.22–22.80	14.94	11.18–18.97
Origin of Plasmodium in mammals	39	38.14	31.44–47.17	40.38	33.80–48.54	34.71	30.02–40.54

Note.—Calibrations, point time estimates and their associated 95% Crls are shown in millions of years. Node numbers are listed in supplementary figure S2, Supplementary Material online. See Materials and Methods for more details.

aClades are shown in figure 1.

Previous phylogenetic studies have indicated that P. coatneyi and P. knowlesi are sister taxa (e.g., Escalante et al. 2005), and, indeed, these two species are usually referred as “closely related” in the malaria literature. In this study, the age of the most recent common ancestor of P. knowlesi and P. coatneyi was estimated to be 5.54 Ma (CrI 4.17–7.10, table 3), which lies within the time interval proposed for the radiation of macaques and other nonhuman primates in the Malaysian archipelago (Delson 1980; Perelman et al. 2011; Springer et al. 2012).

Visual inspection of the Plasmodium phylogeny (fig. 6) suggests that the speciation rate may have been elevated within the clade of malaria parasites infecting Asian macaques. In fact, this hypothesis is consistent with the results obtained when we used the program MEDUSA to compare different birth–death models on a sample of 50 trees describing the relationships between the Plasmodium species known to infect mammalian hosts (table 4 and supplementary fig. S3, Supplementary Material online). In every case, the corrected Akaike information criterion (AICc) calculated by MEDUSA was minimized by a two component birth–death model with accelerated branching within the clade of Asian macaque-infecting parasites. Although in five cases the reduction in AICc was too small to reject a simpler model with homogeneous birth and death rates across the entire tree, in the remaining 45 cases, the two-component model was identified as the best-fitting model by the stepwise selection procedure implemented in MEDUSA (table 4). Maximum likelihood (ML) estimates of the branching and extinction rates within the two-component models suggest that the parasites infecting Asian macaques have speciated approximately three to four times more rapidly than the other lineages contained in the tree. Across the 50 trees analyzed, the diversification rates (r = b–d; events/million years) range from 0.1245 to 0.3257 (mean 0.2055) within the Asian clade and 0.0394 to 0.0728 (mean 0.0580) within the rest of the tree. Furthermore, these differences appear to be almost entirely due to changes in the branching rate rather than the extinction rate. Indeed, in every case the extinction rate was estimated to be less than 0.2% of the branching rate.

Table 4.

Result of the Test of Different Species Divergence Rates.

nd1	nd2	# Taxa	mod1	mod2	P value	Evidence	Contents
42	43	18	exp	v. r.	0.00923	33.89	+P. malariae
43	44	17	exp.	exp	0.00313	76.24	+P. gonderi, Plasmodium sp. (AY800112)
44	45	15	exp.	v. r.	0.00197	152.56	Asian clade
45	46	14	exp.	exp	0.00304	109.34	−P. fragile
46	47	11	exp.	exp	0.00886	36.77	−P. coatneyi, P. knowlesi, Plasmodium sp. (pig-tailed 27C)

Note.—The results were obtained in the PRC test in software MEDUSA (see Materials and Methods). Table shows splits with significant support (P < 0.05) under an LRT. The change in branch length distribution is inferred to occur between nodes nd1 and nd2 (see supplementary fig. S3, Supplementary Material online, for ML with labeled nodes). # taxa denotes the number of taxa within the subtree, whereas mod1 and mod2 specify the branch length distributions for the subtree and its complement, respectively. The evidence ratio is shown under the column labeled evidence. The exponential (exp) and variable rates (v.r.) models were tested.

Discussion

Although there has been increasing awareness of the importance of parasitic organisms for biodiversity (Poulin 2005), most studies continue to focus on the harm caused by parasites to their hosts (Nunn and Altizer 2006) and pay relatively little attention to the role of parasites as intrinsic elements of an ecosystem. Nonetheless (Brooks and McLennan 1993; Page 2003; Poulin 2005; Brooks et al. 2006; Poulin et al. 2011), it is notable that research programs in evolutionary parasitology addressing host specialization (host range) and parasite speciation have followed intertwined paths. Although host evolutionary history is considered a major driver in parasite speciation (Poulin et al. 2011), the discussion has been recently enriched by highlighting the importance of geographic patterns in the host or hosts populations (Nieberding et al 2008; Ricklefs 2010). However, a common framework that encompasses host-range and different forms of geographic speciation is still a work in progress (Huyse et al. 2005). Malarial parasites and other Apicomplexa allow such basic issues to be addressed since there is an overall knowledge of their basic biology and their species diversity.

Here, we start to explore the relative importance of geography and host diversity and evolutionary history in the observed pattern of speciation of a diverse monophyletic group of parasites with overlapping geographical distributions. Indeed, processes such as adaptive radiation can be formally investigated under such circumstances (Glor 2010; Aguilée et al. 2013); otherwise, phylogenetic patterns could be consistent with other processes making any assertion of adaptive radiation untestable (Hunter 1998; Hoberg et al. 2008; Glor 2010; Aguilée et al. 2013). It follows that demonstrating adaptive radiation in a genus such as Plasmodium with a complex evolutionary history and wide geographic distribution based solely on its apparent “rapid speciation” is particularly difficult (Pacheco et al. 2011; Pacheco et al. 2012b). Although we are far from understanding the complex evolutionary history of the nonhuman primate malarias in Southeast Asia, this study unveils several interesting patterns.

First, we found that the two sympatric macaque species found in Borneo share the same set of Plasmodium parasites. Thus, even with limited data, we do not observe evidence of host specificity. This result is perhaps surprising when we consider that the two Macaca species probably diverged 4–5 Ma (Tosi et al. 2003; Fabre et al. 2009) and now have broad geographic distributions. Second, we did not observe that populations of P. inui infecting these two host species were genetically differentiated. This observation suggests that transmission of this multihost parasite occurs frequently between host individuals belonging to different species, at least on the time scale resolved by mtDNA. Third, our results indicate that some of the parasite populations are geographically structured within the island of Borneo. This is evidenced by the orangutan parasites from clade A (fig. 4). This result is expected and simply reflects the spatial distribution of the host population.

The most important finding concerning the orangutan malarias is that there are multiple sympatric clades of Plasmodium. The phylogenetic relationships between lineages C and A and B in a clade that also includes the gibbon parasite P. hylobati (found in Hylobates moloch from Western Java, Indonesia) and the macaque parasite P. inui indicates that host switches are relatively common. However, because we have only limited data on parasites from gibbons, the direction of such host switches cannot currently be established. Nevertheless, P. inui seems to have originated through a host switch from Asian apes into macaques as previously suggested using a more limited data set (Pacheco et al. 2012b). Asian ape parasites, thus far, have only been found in single host species, despite multiple sampling efforts (Coatney et al. 1971). In contrast, P. inui infects a broad range of macaque species and Presbytis and also has a broad geographic distribution. Thus, it seems interesting that a parasite clade with many species with apparently high host specificity in Asian apes has given rise to a generalist parasite with a broad host range.

In the case of orangutan malarias, the occurrence of such sympatric lineages within a host could be the result of sympatric speciation, secondary contacts between originally isolated host populations with divergent parasite lineages, or secondary acquisition via host switches as has been proposed in avian parasites (Ricklefs 2010). Although we cannot rule out the possibility of sympatric speciation, the life cycle of malarial parasites involving sexual reproduction in the mosquito vector makes it likely that different parasite lineages infecting the same host species in the same area will frequently mate. Under these conditions, sympatric speciation is expected to occur only if there is very strong selection promoting divergence and reproductive isolation (Gavrilets 2004; Giraud et al. 2006). For this reason, we believe that our data are best explained by scenarios involving secondary contacts. The finding of sympatric divergent lineages in orangutans (A and B) is the third case of “canonically” sympatric parasite species or subspecies (same host and sympatric in a geographic sense) described in primate malarias where the host demographic history seems to provide a suitable explanation for this pattern (Sutherland et al. 2010; Pacheco et al. 2013). Indeed, the two lineages, A and B, shared their most recent common ancestor 1–2 Ma, a time frame that lies within the demographic history of the extant orangutans, a genus that evolved under complex climatic and geographic events that changed the connectivity of its populations (Nieberding et al. 2008; Locke et al. 2011).

Finally, our data suggest that speciation has occurred at an elevated rate within the group of malaria parasites that infect Southeast Asian nonhuman primates (table 4). This rate seems to be 3–4 time faster than in the other groups included in this study (e.g., lemur malarias or the Leverania lineages observed in African Apes) in a time period of approximately 3–4 Ma from the origin of the southeast Asia parasite clade to the split of P. inui and the orangutan parasites; those events also include the divergence of the lineage leading to the human parasite P. vivax from its closest relative in macaques, P. cynomolgi (table 4). Furthermore, we found no evidence of a significant extinction rate in the sampled species of mammalian malarias. Although this inference is contingent on the assumption that we have sampled all extant Plasmodium species infecting mammalian hosts, we obtained qualitatively similar results when we repeated the MEDUSA analyses but specified that the lemur parasite clade contained an additional three species that have not yet been sequenced. Indeed, there are partial mtDNA sequences that suggest the existence of three such species (unpublished data). It is also important to highlight that additional sampling of nonhuman primates in Southeast Asia may increase the number of species contained within this clade (e.g., parasites in gibbons, see Coatney et al. 1971), in which case the evidence for accelerated branching would be even stronger.

In addition to accelerated cladogenesis, there are other lines of evidence needed to determine whether this pattern is the result of adaptations to new hosts and/or environments. This aspect is far more complex to address. Here, we have found some evidence for episodic selection in two proteins involved in the invasion of the red blood cell. Considering the number of proteins that participate in the invasion of the red blood cell, these two antigens only provide a first glimpse of the many molecular adaptations that have presumably emerged in this group of parasites. This is already evident when we consider that some parasites like P. knowlesi and P. coatneyi exhibit antigenic variation while other species in this clade do not (Lapp et al. 2013). In addition, there is an extraordinary diversity of life history traits in this group (Coatney et al. 1971), including the preferred type of red blood cell (e.g., reticulocytes for P. cynomolgi whereas P. knowlesi invades all, including mature erythrocytes) and their periodicity (parasite replication cycle) including a 24 h cycle for P. knowlesi (the only one in mammals with such a short cycle) and a 72 h for P. inui (a convergent trait with an unrelated parasite, P. malariae) among other traits. At this stage, we can only speculate that such phenotypic diversity may allow the coexistence of multiple species by reducing competition. Indeed, we observed a relatively high number of mixed infections in this study (table 1). Nevertheless, all this phenotypic diversity and variation at antigen-encoding genes reported here (supplementary table S2 and S3, Supplementary Material online) and elsewhere emerged in a relatively short period of 3–4 Ma (Pacheco et al. 2011; Pacheco et al. 2012b, Rice et al. 2014, and this investigation); these facts are consistent with, but of course do not prove, that this relatively speciose clade of parasites is the product of an adaptive radiation.

In summary, contrary to previous studies, the evidence does not suggest that host switches have played a central role in the diversification of this group of parasites. On the one hand, we find no evidence of host specialization in the macaque parasites. In particular, we did not observe differentiation of P. inui lineages between the two host species nor did we find parasite species restricted to a single macaque species. On the other hand, although the orangutan parasites do appear to be host specific, phylogenetic analyses suggest that at least two of these lineages have diversified within orangutans. The two lineages that do appear to have originated via host switches are P. inui which emerged from lineages found in Asian apes (Pacheco et al. 2012b) and the lineage leading to P. vivax which now infects humans and African apes (Prugnolle et al. 2013). We hypothesize that the host ranges of these parasite species underwent several historical changes as result of each host species demographic history that affected their population densities and spatial distributions. Importantly, this notion of a host range evolutionary dynamic does not assume that parasites affected their hosts' speciation or distribution (e.g., via differential virulence or other coevolutionary processes). However, host–parasite antagonisms could be part of the parasite evolutionary dynamics as suggested by the observation of episodic selection in the antigens included in this investigation. Although hypothesizing adaptive radiation is tempting in light of the rapid diversification of host species, we are unable to draw such a conclusion in this case due to the confounding effects of biogeographical processes. Overall, elucidating the relative importance of adaptive radiation and biogeographical processes is a matter of great importance. Understanding the evolutionary dynamic of the host ranges in parasite speciation will provide valuable information about the origin of the observed biodiversity of parasites.

Materials and Methods

Study Sites

The island of Borneo is considered a biodiversity hotspot (Meijaard and Nijman 2003). Just in the state of Sabah (73,618 km2 in the northern portion of Borneo), there are more than ten species of nonhuman primates. Approximately 11,000 Eastern Bornean orangutans (Po. p. morio) are found throughout the major lowland forests and floodplains (Ancrenaz et al. 2004). An unknown number of Sundaland pig-tailed macaque (M. nemestrina) and common long-tailed macaque (M. fasicularis) are found almost ubiquitously throughout the state. The samples were collected in two locations separated by approximately 300 km in Sabah: The Sepilok Orangutan Rehabilitation Centre (SORC) and a satellite facility operated by the Sabah Wildlife Department.

The SORC is located 22 km outside of the city of Sandakan (the second largest city in Sabah with a human population of ∼500,000) on the northern edge of the 5,529 hectare Kabili-Sepilok Virgin Jungle Forest Reserve, bordered by oil palm plantations, poultry farms, and fruit orchards, with at least six hotels located within 1 km of the Centre (supplementary fig. S4, Supplementary Material online). Established in 1964 as a sanctuary for the rehabilitation of orphaned, injured, and/or confiscated orangutans (Po. p. morio) and other endangered species, SORC is operated by the Sabah Wildlife Department. Following a 6-month quarantine period, orangutans are taught how to transverse the forest and forage for food. Following extensive health inspections, these animals are eventually relocated or released into the surrounding forest. To facilitate public education and generate operational funds, the public is allowed to view two daily feedings of the free-ranging animals (Ambu 2007).

The satellite facility operated by the Sabah Wildlife Department included in this study is located at the Shangri-La Rasa Ria Resort on Pantai Dalit Beach, in Tuaran. It is approximately 30 km north of Kota Kinabalu, surrounded by 400 hectares of tropical forest, and contains an Orangutan Education Centre, part of the Rasa Ria Nature Reserve and home to four orangutans at the time of sampling. These animals were originally housed at SORC, but brought to Rasa Ria for further rehabilitation. Like SORC, the public is allowed to view two daily feedings of the animals.

Sample Collection

As part of routine husbandry practices by the Sabah Wildlife Department, veterinary staff at the SORC, and Rasa Ria’s medical clinics inspected and sampled animals for health status. This included 34 orangutans (Po. p. morio, 17 female, 17 male), ages 1.6–18.6 years, at the SORC. Eleven of these animals were collected from the clinic’s indoor/outdoor nursery, whereas the remainders were free-ranging animals in the surrounding forest. An additional four orangutans (two female and two male), ages 2–4 years, were collected and sampled at the Rasa Ria indoor/outdoor nursery. For sampling purposes, all but one orangutan (the largest male) were manually restrained for clinical inspection and sample collection. None of these animals was injured or febrile at the time.

At SORC, samples were collected from 15 pig-tailed macaques (M. nemestrina, seven females and eight males, all adults) and 26 long-tailed macaques (M. fascicularis, 10 females, and 16 males, all adults). These animals were inspected and sampled after approximately 2 mg/kg of Zoletil 100 (zolazepam-tiletamine; Virbac, South Africa) was delivered intramuscularly via syringe dart through a Vario CO2 rifle or blowpipe from Telinject (Saugus, California). Respiration and heart rate were monitored while anesthetized. For the macaques and orangutans alike, blood was collected from the femoral vein using standard venipuncture technique with sterile, single-use Vacutainer productions (Beckton-Dickson). Multiple aliquots of whole blood were stored in cryovials, frozen at 0 °C for 2 weeks before shipment to the United States where they were stored at −80 °C. Aliquots of whole blood were also stored in cryovials using RNAlater (500 µl of blood mixed with 1.3 ml RNAlater; Life Technologies Corporation), frozen at 0 ° for 2 weeks before shipment to the United States where they were stored at −80 °C. All the samples were collected in July 2010 and November 2011. No animal was sampled more than once.

Samples were imported from Sabah to Indiana University using CDC Public Health Service import permit numbers 2011-03-048 and 2013-06-111 and CITES permit numbers 12US77006A/9 and 9062. Permission to conduct all research in Sabah was granted by the Sabah Wildlife Department, and all animal handling was done by trained staff of the Sabah Wildlife Department. Ethical permission was granted by the Bloomington Institutional Animal Care and Use Committee at Indiana University.

DNA Extraction and Malaria Molecular Diagnostics

Genomic DNA was extracted from whole blood using QIAamp DNA Blood Mini kit (Qiagen, GmbH, Hilden, Germany) and each sample was screened for Plasmodium parasites by nested PCR, using primers for a 1,200 bp fragment of the Cytochrome b (cytb) gene that have been used in previous studies (Pacheco et al. 2011, 2013; Pacheco et al. 2012b). The cytb external primers were Forward 5′-TGT AAT GCC TAG ACG TAT TCC/Reverse 5′-GT CAA WCA AAC ATG AAT ATA GAC and the internal primers were Forward 5′-T CTA TTA ATT TAG YWA AAG CAC/Reverse 5′-G CTT GGG AGC TGT AAT CAT AAT. The primary PCR amplifications were carried out in a 50 µl volume reaction using 20 ng of total genomic DNA, 3 mM MgCl2, 1 × PCR buffer, 1.25 mM of each deoxynucleoside triphosphate, 0.4 mM of each primer, and 0.03 U/µl AmpliTaq polymerase (Applied Biosystems, Roche-USA). The primary PCR conditions were: A partial denaturation at 94 °C for 4 min and 36 cycles with 1 min at 94 °C, 1 min at 53 °C and 2 min extension at 72 °C, and a final extension of 10 min at 72 °C was added in the last cycle.

The nested PCRs were also made in a 50 µl volume reaction using only 1 µl of the primary PCRs, 1.5 mM MgCl2, 1 × PCR buffer, 1.25 mM of each deoxynucleoside triphosphate, 0.4 mM of each primer, and 0.03 U/µl AmpliTaq polymerase. The nested PCR conditions were: A partial denaturation at 94 °C for 4 min and 25 cycles with 1 min at 94 °C, 1 min at 56 °C and 2 min extension at 72 °C, and a final extension of 10 min at 72 °C was added in the last cycle. Both strands for all the cytb fragments were direct sequenced, using an Applied Biosystems 3730 capillary sequencer, and identified as Plasmodium using BLAST (Altschul et al. 1997).

Molecular Data

Three loci were amplified from malaria positive samples. First, we amplified approximately 5,800 bp of the parasites mitochondrial genomes (mtDNA). The mtDNA is a linear genome with three genes encoding Cytochrome B, Cox I, and Cox III (Pacheco et al. 2011). We use the mtDNA genome because no pattern consistent with positive selection has been identified, it has a consistent A+T content across Plasmodium species (in contrast with the nuclear genome), and there is a relatively rich database of Plasmodium mitochondrial data that facilitates phylogenetic studies (Escalante et al. 1998; Pacheco et al. 2011; Pacheco et al. 2012b). In addition, we sequenced two nuclear genes encoding major antigens: AMA-1 and MSP-142. These two malarial antigens are essential in the invasion of the host red blood cell and have been widely studied (Chesne-Seck et al. 2005; Pacheco et al. 2007). As in the case of the mtDNA, there is a rich data set on these antigens that facilitates species identification.

The mtDNA was amplified using the oligos Forward 5′-GA GGA TTC TCT CCA CAC TTC AAT TCG TAC TTC/Reverse 5′-CAG GAA AAT WAT AGA CCG AAC CTT GGA CTC with TaKaRa LA TaqTM Polymerase (TaKaRa Mirus Bio Inc). The PCR conditions were: A partial denaturation at 94 °C for 1 min and 30 cycles with 30 s at 94 °C and 7 min at 68 °C, a final extension of 10 min at 72 °C was added. To detect mixed infections, in all cases, at least two independent PCR products were purified using QIAquick Gel extraction kit (Qiagen, GmbH), cloned in the pGEM-T Easy Vector systems (Promega, USA), and a minimum of four clones and both strands were sequenced from each individual.

The AMA-1 was amplified by using degenerated primers Forward 5′-AT GAA TAA AAT ATA CTR CAT AMT/Reverse 5′-TC AGT AGT AWG GCT TCT CCA. PCR was carried out in a 50 µl volume reaction using 20 ng of total genomic DNA, 2.5 mM MgCl2, 1 × PCR buffer, 1.25 mM of each deoxynucleoside triphosphate, 0.4 mM of each primer, and 0.03 U/µl AmpliTaq Gold polymerase, and the PCR conditions were: A partial denaturation at 94 °C for 4 min and 36 cycles with 1 min at 94 °C, 1 min at 52 °C and 2 min extension at 72 °C, and a final extension of 10 min at 72 °C was added in the last cycle. For each sample, at least two independent PCR products were also purified from the agarose gel using QIAquick gel extraction kit (Qiagen, GmbH), cloned, and both strands for a minimum of five clones were sequenced from each individual.

Finally, the MSP-142 was amplified by using the primers Forward 5′-GAC CAA GTA ACA ACG GGA G/Reverse 5′-C AAA GAG TGG CTC AGA ACC in 50 µl volume with 20 ng/µl of total genomic DNA, 2.5 mM MgCl2, 1 × PCR buffer, 1.25 mM of each deoxynucleoside triphosphate, 0.4 mM of each primer, and 0.03 U/µM of AmpliTaq Gold DNA polymerase (Applied Biosystems, Roche-USA). The PCR conditions were a partial denaturation at 94 °C for 4 min and 35 cycles of 1 min at 94 °C, 1 min at 55 °C, and 2 min extension at 72 °C, and a final extension of 10 min was added in the last cycle. Like in the case of the mtDNA, at least two independent PCR products per sample were purified from the agarose gel using QIAquick gel extraction kit (Qiagen, GmbH), cloned, and both strands were sequenced for a minimum of five clones for each sample.

Genetic diversity within different Plasmodium species/lineages was quantified by calculating the parameter π (nucleotide diversity) for the mtDNA and the two antigens (AMA-1 and MSP-142). This statistic estimates the average number of substitutions between any two sequences as implemented in MEGA v5.2.2 (Tamura et al. 2011). Standard errors of π estimates were calculated with 1,000 bootstrap replicates.

The sequences reported in this investigation were submitted to GenBank under the accession numbers KJ569800 to KJ569917 (see supplementary table S1, Supplementary Material online, for details).

Phylogenetic Analyses

Independent alignments for nucleotide sequences of the mtDNA (5,717 bp excluding gaps), ama-1 (1,663 bp excluding gaps), and msp-142 (1,057 bp excluding gaps) genes were made using ClustalX v2.0.12 and Muscle as implemented in SeaView v4.3.5 with manual editing. First, phylogenetic relationships were inferred for the mtDNA using both the ML method implemented in PhyML v3.0 (Guindon and Gascuel 2003) and Bayesian methods using MrBayes v3.1.2 with the default priors (Ronquist and Huelsenbeck 2003). The reliability of the nodes in the ML tree was assessed by the bootstrap method with 200 pseudoreplications. In the Bayesian analysis, each of the three mitochondrial genes plus the noncoding regions were used as separate partitions (Pacheco et al. 2011). Bayesian support for the nodes was inferred in MrBayes by sampling every 100 generations from two independent chains lasting 8 × 106 Markov Chain Monte Carlo (MCMC) steps. The chains were assumed to have converged once the average standard deviation of the posterior probability was below 0.01, and the value of the potential scale reduction factor was between 1.00 and 1.02 (Ronquist and Huelsenbeck 2003). Fifty percent of the sample was then discarded as a burn-in once convergence was reached. Both phylogenetic methods used a general time reversible model with gamma-distributed substitution rates and a proportion of invariant sites (GTR + Γ + I) and indeed this model had the lowest Bayesian Information Criterion scores as estimated by MEGA v5.2.2 (Nei and Kumar 2000; Tamura et al. 2011). In order to compare the different Plasmodium species found in this study, we estimated the genetic divergences of the mtDNA genomes among and within different species using the Kimura 2-parameter model as implemented in MEGA v5.2.2 (Tamura et al. 2011).

Phylogenetic relationships were also inferred for the ama-1 and msp-142 kDa genes using both the ML method implemented in PhyML v3.0 and Bayesian methods with the default priors (Ronquist and 2003). In the case of the Bayesian inference, two independent chains were sampled every 100 generations in runs lasting 6 × 106 MCMC steps for ama-1 gene and 3 × 106 for msp-142. As before, 50% of the sample was discarded as a burn-in period once the chain had achieved convergence. We used a general time reversible model with gamma-distributed substitution rates (GTR + Γ) for both genes because it best fit the data as estimated by MEGA v5.2.2 (Nei and Kumar 2000; Tamura et al. 2011).

Supplementary table S1, Supplementary Material online, provides a complete list of the species and sequences obtained in this study; it also includes those previously published that were used as part of our phylogenetic analyses. Because the same Plasmodium species (almost identical mtDNA with the exception of few SNPs) were found in both species of Macaca, we only amplified the genes encoding the two major antigens (ama-1 and msp-142) from the 15 positive samples of pig-tailed macaques (M. nemestrina). Information about basic biology, geographic distribution and host-range of these species can be found elsewhere (Coatney et al. 1971).

The median joining network for a complete set of 41 Plasmodium sp. mtDNA haplotypes from orangutans was estimated using Network v4.6.1.0 (Fluxus Technologies 2011), with transversions weighted twice as large as transitions and the epsilon parameter set equal to 0. In this analysis, we included all 30 Plasmodium sp. mtDNA haplotypes obtained from orangutan in this study (Sabah) and also 11 mtDNA haplotypes from Kalimantan orangutans that have been reported in somewhere else (Pacheco et al. 2012b). Also, we estimated a median joining network for a complete set of 16 P. inui mtDNA haplotypes from both species of Macaca (eight haplotypes from M. nemestrina and eight from M. fascicularis), as was described before. In addition, we also used DnaSP v.5 (Librado and Rozas 2009) to estimate the fixation index (FST) between the two populations (Sabah and Kalimantan, Borneo) of Plasmodium sp. from orangutan-clade A (for more details about the clades see the Results section) and between the two populations of P. inui from M. nemestrina and M. fascicularis (Sabah).

Signatures Consistent with Natural Selection

Although we cannot attest that msp-1 or ama-1 genes are particularly important in parasite speciation, their role in the invasion of the red blood cell could indicate adaptation to a new host. We use a branch-site method or REL (Pond et al. 2011) as implemented in HyPhy (Pond et al. 2005). This approach aims to identify lineages where a proportion of codons evolve with ω >1 by using an LRT. This method does not make assumptions about which lineages are under selection so it is ideal for this type of exploratory analyses. Because these methods assume the true phylogeny, we used the consensus phylogeny derived from a Bayesian analysis that only included lineages that we consider “good species” or clearly supported lineages (1.0 posterior probability on both antigens). Although many of these lineages are accepted species, others have just been proposed in this investigation.

Estimation of Divergence Times

BEAST v1.7.5 (Drummond et al. 2012) was used to estimate time trees for the mtDNA sequences (supplementary table S1, Supplementary Material online) using four partitions of the data: each gene (cox1, cox3, cytb) plus the noncoding regions (Pacheco et al. 2011, 2013; Pacheco et al. 2012b). Relaxed clock methods were applied with a general time reversible model with gamma-distributed substitution rates and a proportion of invariant sites (GTR + Γ + I), with heterogeneity both among sites and across the four partitions. Two independent chains were run for each analysis until convergence appeared to have been achieved. Three scenarios were explored, each of them incorporating different calibration points, two based on host fossils and the third based on a biogeographical landmark for the host.

The first calibration point assumes that African parasites found in Mandrillus spp. and Cercocebus spp. diverged from those Plasmodium spp. found in Southeast Asia macaques when Macaca branched from Papio (Mu et al. 2005; Hayakawa et al. 2008; Pacheco et al. 2011). Fossils identified as Macaca spp. indicate that such an event took place 6–8 Ma as minimum and maximum boundaries, respectively (Delson 1980). These time points are used as calibration with an exponential prior assuming that the divergence of the African and Asian parasites took place with 95 % of probability between 6 and 8 Ma (Mu et al. 2005). However, these fossils are used as minimum times when studying primates (Perelman et al. 2011). A more inclusive alternative calibration was constructed around this event assuming that the divergence at the same node took place at any time between 6 and 14.2 Ma using a uniform prior; this incorporates molecular estimates for the same Papio–Macaca divergence event (Pacheco et al. 2011). The use of 14.2 Ma as a maximum is consistent with older fossils reported for Macaca spp. (see Paleobiology Database at http://www.paleodb.org/, last accessed November 17, 2014) and considers also the fact that P. gonderi is a parasite of Chlorocebus (a Cercopithecini). In either case, these calibrations assume that the parasite lineages did not diverge before their host taxa had diverged, giving an upper bound of either 14.2 or 8 Ma. They also assume that the parasite lineages did not diverge after the host taxa had become isolated on different continents, that is, parasite diversification did not involve intercontinental dispersal separate from that of the host taxa, giving a lower bound of 6 Ma.

A second fossil-based calibration is the minimum of 23.5 Ma for the human/Macaca split (Benton and Donoghue 2007), which is assumed to be the time when the monophyletic group that includes P. malariae (a parasite found in humans) and all the Asian parasites originated. This event is used as calibration with an exponential prior set conservatively with 97.5% of the prior probability below to 65 Ma, so it allows for the clade to be as old as the origin of primates, or even older, if the data support that (Perelman et al. 2011; Springer et al. 2012). In addition to these calibrations that incorporate information from fossils, an additional calibration point considers the origin of lemur malarias (a biogeographical event), a monophyletic group that radiated in Madagascar (Pacheco et al. 2011). The colonization of Madagascar by terrestrial mammals, including lemurs, happened during a relatively short time frame in the Cenozoic (∼65−20 Ma) during which ocean currents allowed for the eastward transport of vegetation rafts from Africa to Madagascar (Poux et al. 2005; Ali and Huber 2010). Thus, times younger than approximately 20 Ma are unlikely to be biologically realistic for the origin of lemur malarias (Pacheco et al 2011). This event was used as calibration with an exponential prior set conservatively with 97.5% of the prior probability below to 65 Ma (Perelman et al. 2011; Springer et al. 2012).

The three scenarios considered in this investigation utilize the calibrations described above as follow: 1) a combination of the relaxed 6–14.2 Ma calibration for Papio–Macaca divergence (Delson 1980; Mu et al. 2005; Hayakawa et al. 2008) with a minimum of 20 Ma for the origin of lemur parasites (Poux et al. 2005; Ali and Huber 2010); 2) a combination of the 6–14.2 Ma calibration with a minimum of 23.5 Ma for the human/Macaca split and 20 Ma for the origin of the lemur parasites; and 3) the 6–8 Ma narrowly defined around the fossils of the Papio–Macaca divergence (Delson 1980) with a minimum of 20 Ma for the origin of the lemur lineage (see Pacheco et al. 2011, 2013; Pacheco et al. 2012b). Although the time estimates under these scenarios are expected to overlap, their comparison is important since the calibrations rely on different assumptions (Pacheco et al. 2011; Ramiro et al. 2012). An extensive discussion on different timing methods and calibrations was reported elsewhere (Pacheco et al. 2011).

Species Diversification Rates

Two methods were used to test whether the branching rate (cladogenesis) has changed across the Plasmodium mtDNA phylogeny. In both cases, we retained only those 38 lineages thought to correspond to good species, while excluding closely related haplotypes that appear to have been sampled from within species. We first carried out the nonparametric cladogenesis test (Nater et al. 1992) which looks for evidence of accelerated branching within a phylogeny. This was applied to a ML tree inferred for these 38 taxa. We then used the program MEDUSA (Alfaro et al. 2009) to analyze a sample of 50 trees from a Bayesian phylogenetic analysis using an inhomogeneous birth–death model. Given a rooted tree, MEDUSA fits a sequence of increasingly complicated birth–death models to that tree and uses a stepwise procedure using a corrected AICc to select the best-fitting model. The models differ in the number of components into which the tree is partitioned, with each component being assigned its own birth and death rates. In our analyses, we considered models with up to ten components. The trees used in this analysis were generated by two independent runs of BEAST v1.7.5 using the same prior distributions and calibration points as previously described, with one tree being sampled every million generations from generations 20–50 million in the first chain and from generations 32–50 million in the second chain. Monophyly of the parasites infecting mammalian hosts was enforced in these analyses and the clade containing four species with avian or reptilian hosts was pruned from each tree prior to analysis with MEDUSA.

Supplementary Material

Supplementary tables S1–S4 and figures S1–S4 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).